-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwwm_st4.F90
322 lines (291 loc) · 13.5 KB
/
wwm_st4.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
#include "wwm_functions.h"
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE ST4_PRE (IP, WALOC, PHI, DPHIDN, SSINE, DSSINE, SSDS, DSSDS, SSNL4, DSSNL4, SSINL)
USE DATAPOOL
USE W3SRC4MD
IMPLICIT NONE
INTEGER, INTENT(IN) :: IP
REAL(rkind), INTENT(IN) :: WALOC(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: PHI(NUMSIG,NUMDIR), DPHIDN(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: SSINE(NUMSIG,NUMDIR), DSSINE(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: SSDS(NUMSIG,NUMDIR),DSSDS(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: SSNL4(NUMSIG,NUMDIR),DSSNL4(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: SSINL(NUMSIG,NUMDIR)
INTEGER :: IS, ID, ITH, IK, IS0
REAL(rkind) :: AWW3(NSPEC)
REAL(rkind) :: VDDS(NSPEC), VSDS(NSPEC), BRLAMBDA(NSPEC)
REAL(rkind) :: WN2(NUMSIG*NUMDIR), WHITECAP(1:4)
REAL(rkind) :: VSIN(NSPEC), VDIN(NSPEC)
REAL(rkind) :: ETOT, FAVG, FMEAN1, WNMEAN, AS, SUMWALOC, FAVGWS
REAL(rkind) :: TAUWAX, TAUWAY, AMAX, FPM, WIND10, WINDTH
REAL(rkind) :: HS,SME01,SME10,KME01,KMWAM,KMWAM2
DO IS = 1, NUMSIG
DO ID = 1, NUMDIR
AWW3(ID + (IS-1) * NUMDIR) = WALOC(IS,ID) * CG(IS,IP)
END DO
END DO
DO IK=1, NK
WN2(1+(IK-1)*NTH) = WK(IK,IP)
END DO
DO IK=1, NK
IS0 = (IK-1)*NTH
DO ITH=2, NTH
WN2(ITH+IS0) = WN2(1+IS0)
END DO
END DO
!
! wind input
!
TAUWX(IP) = ZERO
TAUWY(IP) = ZERO
SSINL = ZERO
NUMSIG_HF(IP) = NUMSIG
AS = 0.
BRLAMBDA = ZERO
IF (MESIN .GT. 0) THEN
CALL SET_WIND( IP, WIND10, WINDTH )
CALL SET_FRICTION( IP, WALOC, WIND10, WINDTH, FPM )
LLWS=.TRUE.
#ifdef DEBUGSRC
WRITE(740+myrank,*) '1: input value USTAR=', UFRIC(IP), ' USTDIR=', USTDIR(IP)
#endif
CALL W3SPR4 ( AWW3, CG(:,IP), WK(:,IP), EMEAN(IP), FMEAN(IP), FMEAN1, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS(IP))
#ifdef DEBUGSRC
WRITE(740+myrank,*) '1: out value USTAR=', UFRIC(IP), ' USTDIR=', USTDIR(IP)
WRITE(740+myrank,*) '1: out value EMEAN=', EMEAN(IP), ' FMEAN=', FMEAN(IP)
WRITE(740+myrank,*) '1: out value FMEAN1=', FMEAN1, ' WNMEAN=', WNMEAN
WRITE(740+myrank,*) '1: out value CD=', CD(IP), ' Z0=', Z0(IP)
WRITE(740+myrank,*) '1: out value ALPHA=', ALPHA_CH(IP), ' FMEANWS=', FMEANWS(IP)
#endif
IF (EMEAN(IP) .LT. THR .AND. WIND10 .GT. THR) CALL SIN_LIN_CAV(IP,WINDTH,FPM,SSINL)
CALL W3SIN4 ( IP, AWW3, CG(:,IP), WN2, WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, VSIN, VDIN, LLWS, BRLAMBDA)
#ifdef DEBUGSRC
WRITE(740+myrank,*) '1: WINDTH=', WINDTH, ' Z0=', Z0(IP), ' CD=', CD(IP)
WRITE(740+myrank,*) '1: UFRIC=', UFRIC(IP), 'WIND10=', WIND10, ' RHOAW=', RHOAW
WRITE(740+myrank,*) '1: TAUWX=', TAUWX(IP), ' TAUWY=', TAUWY(IP)
WRITE(740+myrank,*) '1: TAUWAX=', TAUWAX, ' TAUWAY=', TAUWAY
WRITE(740+myrank,*) '1: W3SIN4min/max/sum(VSIN)=', minval(VSIN), maxval(VSIN), sum(VSIN)
WRITE(740+myrank,*) '1: W3SIN4min/max/sum(VDIN)=', minval(VDIN), maxval(VDIN), sum(VDIN)
#endif
CALL W3SPR4 ( AWW3, CG(:,IP), WK(:,IP), EMEAN(IP), FMEAN(IP), FMEAN1, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS(IP))
CALL W3SIN4 ( IP, AWW3, CG(:,IP), WN2, WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, VSIN, VDIN, LLWS, BRLAMBDA)
#ifdef DEBUGSRC
WRITE(740+myrank,*) '2: W3SIN4min/max/sum(VSIN)=', minval(VSIN), maxval(VSIN), sum(VSIN)
WRITE(740+myrank,*) '2: W3SIN4min/max/sum(VDIN)=', minval(VDIN), maxval(VDIN), sum(VDIN)
#endif
CALL CONVERT_VS_VD_WWM(IP, VSIN, VDIN, SSINE, DSSINE)
ENDIF
IF (MESNL .GT. 0) THEN
CALL MEAN_WAVE_PARAMETER(IP,WALOC,HS,ETOT,SME01,SME10,KME01,KMWAM,KMWAM2)
CALL DIASNL4WW3(IP, KMWAM, WALOC, SSNL4, DSSNL4)
END IF
IF (MESDS .GT. 0) THEN
CALL W3SDS4(AWW3,WK(:,IP),CG(:,IP),UFRIC(IP),USTDIR(IP),DEP(IP),VSDS,VDDS,BRLAMBDA,WHITECAP)
#ifdef DEBUGSRC
WRITE(740+myrank,*) '2: W3SDS4min/max/sum(VSDS)=', minval(VSDS), maxval(VSDS), sum(VSDS)
WRITE(740+myrank,*) '2: W3SDS4min/max/sum(VDDS)=', minval(VDDS), maxval(VDDS), sum(VDDS)
#endif
CALL CONVERT_VS_VD_WWM(IP, VSDS, VDDS, SSDS, DSSDS)
ENDIF
!
PHI = SSINL + SSINE + SSNL4 + SSDS
DPHIDN = DSSINE + DSSNL4 + DSSDS
!
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE ST4_POST (IP, WALOCOLD, WALOC)
USE DATAPOOL
USE W3SRC4MD
IMPLICIT NONE
INTEGER, INTENT(IN) :: IP
REAL(rkind), INTENT(IN) :: WALOCOLD(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: WALOC(NUMSIG,NUMDIR)
REAL(rkind) :: SSINE(NUMSIG,NUMDIR),DSSINE(NUMSIG,NUMDIR)
REAL(rkind) :: SSDS(NUMSIG,NUMDIR),DSSDS(NUMSIG,NUMDIR)
REAL(rkind) :: SSINL(NUMSIG,NUMDIR)
INTEGER :: IS, ID, IK, ITH, ITH2, IS0, NKH, NKH1
REAL(rkind) :: PHI(NUMSIG,NUMDIR), DPHIDN(NUMSIG,NUMDIR)
REAL(rkind) :: AWW3(NSPEC), AWW3OLD(NSPEC), WN2(NUMSIG*NUMDIR), BRLAMBDA(NSPEC)
REAL(rkind) :: SSINE_WW3(NSPEC), DSSINE_WW3(NSPEC), TMP_DS(NUMSIG)
REAL(rkind) :: SSDS_WW3(NSPEC), DSSDS_WW3(NSPEC), SSNL4_WW3(NSPEC), DSSNL4_WW3(NSPEC)
REAL(rkind) :: SSBF_WW3(NSPEC), DSSBF_WW3(NSPEC), SSNL4(NUMSIG,NUMDIR), DSSNL4(NUMSIG,NUMDIR)
REAL(rkind) :: SSBF(NUMSIG,NUMDIR), DSSBF(NUMSIG,NUMDIR)
REAL(rkind) :: ETOT, FAVG, FMEAN1, WNMEAN, AS, FAVGWS
REAL(rkind) :: TAUWAX, TAUWAY, AMAX, WIND10, WINDTH, EBAND, B1BAND
REAL(rkind) :: WHITECAP(1:4), SUMWALOC, FPM, FH1, FH2, TAUBBL(2)
REAL(rkind) :: PHIAW, CHARN, PHINL, PHIBBL, TAUWIX, TAUWIY, TAUWNX, TAUWNY, TAUOX, TAUOY
REAL(rkind) :: FACTOR, FACTOR2, MWXFINISH, MWYFINISH, EFINISH, DIFF, CONST2, TEMP2
REAL(rkind) :: A1BAND, PHIOC, HSTOT, PIBBL, FAGE, FHIGH, FACHFA
! IF (MINVAL(WALOC) .LT. 0) STOP 'POST 1'
! IF (MINVAL(WALOCOLD) .LT. 0) STOP 'POST 2'
DO IS = 1, NUMSIG
DO ID = 1, NUMDIR
AWW3(ID + (IS-1) * NUMDIR) = WALOC(IS,ID) * CG(IS,IP)
AWW3OLD(ID + (IS-1) * NUMDIR) = WALOCOLD(IS,ID) * CG(IS,IP)
END DO
END DO
DO IK=1, NK
WN2(1+(IK-1)*NTH) = WK(IK,IP)
END DO
DO IK=1, NK
IS0 = (IK-1)*NTH
DO ITH=2, NTH
WN2(ITH+IS0) = WN2(1+IS0)
END DO
END DO
!
! wind input
!
AS = ZERO
NUMSIG_HF(IP) = NUMSIG
CALL SET_WIND( IP, WIND10, WINDTH )
CALL SET_FRICTION( IP, WALOC, WIND10, WINDTH, FPM )
CALL W3SPR4 ( AWW3, CG(:,IP), WK(:,IP), EMEAN(IP), FMEAN(IP), FMEAN1, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS(IP))
IF (EMEAN(IP) .LT. THR .AND. WIND10 .GT. THR) CALL SIN_LIN_CAV(IP,WINDTH,FPM,SSINL)
CALL W3SIN4 ( IP, AWW3, CG(:,IP), WN2, WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, SSINE_WW3, DSSINE_WW3, LLWS, BRLAMBDA)
CALL W3SPR4 ( AWW3, CG(:,IP), WK(:,IP), EMEAN(IP), FMEAN(IP), FMEAN1, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS(IP))
!
! dissipation
!
CALL W3SDS4(AWW3,WK(:,IP),CG(:,IP),UFRIC(IP),USTDIR(IP),DEP(IP),SSDS_WW3,DSSDS_WW3,BRLAMBDA,WHITECAP)
!
! snl4
!
IF (MESNL .GT. 0) CALL DIASNL4WW3(IP, WNMEAN, WALOC, SSNL4, DSSNL4)
!
! and last but not least the bottom friction dissipation ...
!
IF (MESBF .GT. 0) CALL SDS_BOTF(IP,WALOC,SSBF,DSSBF)
!
! Resio & Perrie ... scaling ...
!
FAGE = FFXFA*TANH(0.3*WIND10*FMEANWS(IP)*PI2/G9)
FH1 = (FFXFM+FAGE) * FMEAN1
FH2 = FFXPM / UFRIC(IP)
FHIGH = MIN ( SIG(NK) , MAX ( FH1 , FH2 ) )
NKH = MIN ( NK , INT(FACTI2+FACTI1*LOG(MAX(1.E-7,FHIGH))) )
NKH1 = MIN ( NK , NKH+1 )
NKH = MAX ( 2 , MIN ( NKH1 , INT ( FACTI2 + FACTI1*LOG(MAX(1.E-7_rkind,FHIGH)) ) ) )
!
! Add tail
!
FACHF = 5.
FACHFA = XFR**(-FACHF-2)
DO IK=NKH+1,NK
DO ITH=1, NTH
AWW3(ITH+(IK-1)*NTH) = AWW3(ITH+(IK-2)*NTH) * FACHFA + 0. ! Adding a magic zero without a comment ...
END DO
END DO
!
! convert to wwm convetion ...
!
DO IS = 1, NUMSIG
DO ID = 1, NUMDIR
WALOC(IS,ID) = AWW3(ID + (IS-1) * NUMDIR) / CG(IS,IP)
END DO
END DO
CALL CONVERT_VS_VD_WWM(IP, SSINE_WW3, DSSINE_WW3, SSINE, DSSINE)
CALL CONVERT_VS_VD_WWM(IP, SSDS_WW3, DSSDS_WW3, SSDS, DSSDS)
IF (MESNL .GT. 0) CALL CONVERT_WWM_VS_VD(IP, SSNL4_WW3, DSSNL4_WW3, SSNL4, DSSNL4)
IF (MESBF .GT. 0) CALL CONVERT_WWM_VS_VD(IP, SSBF_WW3, DSSBF_WW3, SSBF, DSSBF)
!
! compute momentum to BBL
!
DO IK=1, NK
CONST2=DDEN(IK)/CG(IK,IP)*G9/(SIG(IK)/WK(IK,IP))
DO ITH=1,NTH
IS=ITH+(IK-1)*NTH
TEMP2=CONST2*DSSBF_WW3(IS)*AWW3(IS)
TAUBBL(1) = TAUBBL(1) - TEMP2*ECOS(IS)
TAUBBL(2) = TAUBBL(2) - TEMP2*ESIN(IS)
END DO
END DO
!
! adding the fluxes from waves to ocean ...
!
WHITECAP(3) = ZERO
PHIAW = ZERO
CHARN = ZERO
PHINL = ZERO
PHIBBL = ZERO
TAUWIX = ZERO
TAUWIY = ZERO
TAUWNX = ZERO
TAUWNY = ZERO
HSTOT = ZERO
DO IK = 1, NK
FACTOR = DDEN(IK)/CG(IK,IP)
FACTOR2 = FACTOR*G9*WK(IK,IP)/SIG(IK)
DO ITH = 1, NTH
IS = (IK-1) * NTH + ITH
PHIAW = PHIAW + SSINE_WW3(IS) * DT4S * FACTOR / MAX ( ONE , (ONE-DT4S*DSSINE_WW3(IS)))
PIBBL = PHIBBL - SSBF_WW3(IS) * DT4S * FACTOR / MAX ( ONE , (ONE-DT4S*DSSBF_WW3(IS)))
PHINL = PHINL + SSNL4_WW3(IS) * DT4S * FACTOR / MAX ( ONE , (ONE-DT4S*DSSNL4_WW3(IS)))
IF (SSINE_WW3(IS) .GT. ZERO) THEN
WHITECAP(3) = WHITECAP(3) + AWW3(IS) * FACTOR
ELSE
! computes the upward energy flux (counted > 0 upward)
CHARN = CHARN - (SSINE_WW3(IS))* DT4S * FACTOR / MAX ( ONE , (ONE-DT4S*DSSINE_WW3(IS)))
END IF
HSTOT = HSTOT + AWW3(IS) * FACTOR
END DO
END DO
WHITECAP(3) = 4.*SQRT(WHITECAP(3))
HSTOT = 4.*SQRT(HSTOT)
TAUWIX = TAUWIX + TAUWX(IP) * RHOAW * DT4S
TAUWIY = TAUWIY + TAUWY(IP) * RHOAW * DT4S
TAUWNX = TAUWNX + TAUWAX * RHOAW * DT4S
TAUWNY = TAUWNY + TAUWAY * RHOAW * DT4S
!
! The wave to ocean flux is the difference between initial energy
! and final energy, plus wind input plus the SNL flux to high freq.,
! minus the energy lost to the bottom boundary layer (BBL)
!
EFINISH = 0.
MWXFINISH = 0.
MWYFINISH = 0.
DO IK = 1, NK
EBAND = 0.
A1BAND = 0.
B1BAND = 0.
DO ITH = 1, NTH
DIFF = AWW3OLD(ITH+(IK-1)*NTH)-AWW3(ITH+(IK-1)*NTH) ! Check that the difference is taken in the right direction!
EBAND = EBAND + DIFF
A1BAND = A1BAND + DIFF*ECOS(ITH)
B1BAND = B1BAND + DIFF*ESIN(ITH)
END DO
EFINISH = EFINISH + EBAND * DDEN(IK) / CG(IK,IP)
MWXFINISH = MWXFINISH + A1BAND * DDEN(IK) / CG(IK,IP) * WK(IK,IP)/SIG(IK)
MWYFINISH = MWYFINISH + B1BAND * DDEN(IK) / CG(IK,IP) * WK(IK,IP)/SIG(IK)
END DO
!
! Transfoe rmation in momentum flux in m^2 / s^2
!
TAUOX = (G9*MWXFINISH+TAUWIX-TAUBBL(1))/DT4S
TAUOY = (G9*MWYFINISH+TAUWIY-TAUBBL(2))/DT4S
TAUWIX = TAUWIX/DT4S
TAUWIY = TAUWIY/DT4S
TAUWNX = TAUWNX/DT4S
TAUWNY = TAUWNY/DT4S
!
! Transformation in wave energy flux in W/m^2=kg / s^3
!
PHIOC = RHOW*G9*(EFINISH+PHIAW-PHIBBL)/DT4S
PHIAW = RHOW*G9*PHIAW /DT4S
PHINL = RHOW*G9*PHINL /DT4S
PHIBBL= RHOW*G9*PHIBBL/DT4S
#ifdef SCHISM
OUTT_INTPAR(IP,32) = TAUOX ! x-component of the wave-ocean momentum flux (tauox in m2.s-2)
OUTT_INTPAR(IP,33) = TAUOX ! y-component of the wave-ocean momentum flux (tauoy in m2.s-2)
OUTT_INTPAR(IP,34) = PHIOC ! Wave-to-ocean TKE flux (phioc in W.m-2)
OUTT_INTPAR(IP,35) = PHIBBL ! phibbl / wave_phibbl : Energy flux due to bottom friction (phioc in W.m-2)
#endif
! IF (MINVAL(WALOC) .LT. 0) STOP 'POST 3'
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************