-
Notifications
You must be signed in to change notification settings - Fork 1
/
wwm_sdiss_ardh_vec_local.F90
461 lines (410 loc) · 14.4 KB
/
wwm_sdiss_ardh_vec_local.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
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
SUBROUTINE SDISS_ARDH_VEC_LOCAL (IPP, F, FL, SL, F1MEAN, XKMEAN, &
& PHIEPS, TAUWD, M, SSDS, DSSDS)
! ----------------------------------------------------------------------
!**** *SDISSIP_ARDH_VEC* - COMPUTATION OF DISSIPATION SOURCE FUNCTION.
! LOTFI AOUF METEO FRANCE 2013
! FABRICE ARDHUIN IFREMER 2013
!* PURPOSE.
! --------
! COMPUTE DISSIPATION SOURCE FUNCTION AND STORE ADDITIVELY INTO
! NET SOURCE FUNCTION ARRAY. ALSO COMPUTE FUNCTIONAL DERIVATIVE
! OF DISSIPATION SOURCE FUNCTION.
!** INTERFACE.
! ----------
! *CALL* *SDISSIP (F, FL, IJS, IJL, SL, F1MEAN, XKMEAN,*
! PHIEPS, TAUWD, M LCFLX)
! *F* - SPECTRUM.
! *FL* - DIAGONAL MATRIX OF FUNCTIONAL DERIVATIVE
! *IJS* - INDEX OF FIRST GRIDPOINT
! *IJL* - INDEX OF LAST GRIDPOINT
! *SL* - TOTAL SOURCE FUNCTION ARRAY
! *F1MEAN* - MEAN FREQUENCY BASED ON 1st MOMENT.
! *XKMEAN* - MEAN WAVE NUMBER BASED ON 1st MOMENT.
! *PHIEPS* - ENERGY FLUX FROM WAVES TO OCEAN INTEGRATED OVER
! THE PROGNOSTIC RANGE.
! *TAUWD* - DISSIPATION STRESS INTEGRATED OVER
! THE PROGNOSTIC RANGE.
! *MIJ* - LAST FREQUENCY INDEX OF THE PROGNOSTIC RANGE.
! *LCFLX* - TRUE IF THE CALCULATION FOR THE FLUXES ARE
! PERFORMED.
! METHOD.
! -------
! SEE REFERENCES.
! EXTERNALS.
! ----------
! NONE.
! REFERENCE.
! ----------
! ARDHUIN et AL. JPO DOI:10.1175/20110JPO4324.1
! ----------------------------------------------------------------------
!USE YOWFRED , ONLY : FR, TH, FRATIO, DELTH, GOM, COSTH, SINTH, DFIM
!USE YOWPCONS , ONLY : RAD ,G ,ZPI ,ROWATER ,YEPS
!USE YOWMEAN , ONLY : EMEAN, FMEAN
!USE YOWMPP , ONLY : NINF ,NSUP
!USE YOWPARAM , ONLY : NANG ,NFRE
!USE YOWSHAL , ONLY : TFAK ,INDEP, TCGOND
!USE YOWSPEC , ONLY : U10NEW, THWNEW, USNEW
!USE YOWSTAT , ONLY : ISHALLO
!USE YOWTEST , ONLY : IU06 ,ITEST
!USE PARKIND1 ,ONLY : JPIM ,JPRB
!USE YOMHOOK ,ONLY : LHOOK ,DR_HOOK
USE DATAPOOL, ONLY : FR, WETAIL, FRTAIL, WP1TAIL, ISHALLO, RKIND, &
& DFIM, DFIMOFR, DFFR, DFFR2, WK, ZALP, IAB, SWELLFT, &
& IUSTAR, IALPHA, USTARM, TAUT, ONETHIRD, RKIND, ONE, &
& DELUST, DELALP, BETAMAX, XKAPPA, IDAMPING, TAUWSHELTER, &
& FRATIO, EMEAN, USNEW, THWNEW, DEGRAD, LCFLX, MSC, MDC, &
& SINTH, COSTH, &
& ROWATER => RHOW, &
& ROAIR => RHOA, &
& TH => SPDIR, &
& DELTH => DDIR, &
& G => G9, &
& ZPI => PI2, &
& EPSMIN => SMALL, &
& NANG => MDC, &
& NFRE => MSC, &
& INDEP => DEP, &
& CG0 => CG, &
& WK0 => WK
! ----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: K , M, IPP
INTEGER :: I, J, I1, J1
INTEGER :: IS, SDSNTH, DIKCUMUL
INTEGER :: I_INT, J_INT, M2, K2
INTEGER :: NSMOOTH(NFRE)
INTEGER :: ISDSDTH
INTEGER :: ISSDSBRFDF
INTEGER :: MIJ
INTEGER, ALLOCATABLE :: SATINDICES(:,:)
REAL(rkind) :: XK(NFRE), CG(NFRE)
REAL(rkind) :: ALFAMEAN
REAL(rkind) :: TPIINV, TMP00, TMP01, TMP02, TMP03, TMP04
REAL(rkind) :: COSWIND
REAL(rkind) :: DTURB
REAL(rkind) :: EPSR
REAL(rkind) :: W, P0, MICHE, DELTA1, DELTA2
REAL(rkind) :: SCDFM
REAL(rkind) :: SDISS, YEPS
REAL(rkind) :: SDSBR, SDSBR2
REAL(rkind) :: SATURATION2,FACSAT
REAL(rkind) :: SSDSC3, SSDSC4, SSDSC5, SSDSC6, SDSCOS
REAL(rkind) :: SSDSHF, SSDSLF, X, DTEMP, TEMP
REAL(rkind) :: SSDSBRF1, XFR, SXFR,SSDSBR,SDSC1,SSDSP,SSDSC2
REAL(rkind) :: DSIP_
REAL(rkind) :: ROG
REAL(rkind) :: SSDSC1, SSDSC , SSDSISO
!REAL(KIND=JPRB) :: ZHOOK_HANDLE
REAL(rkind) :: SIG(NFRE)
REAL(rkind) :: BSIGBAJ
REAL(rkind) :: F1MEAN, XKMEAN, WNMEAN2
REAL(rkind) :: PHIEPS, TAUWD, CM
REAL(rkind) :: FACTOR, FACTURB
REAL(rkind) :: RENEWALFREQ(NANG)
REAL(rkind), DIMENSION(NFRE) :: CONSTFM
REAL(rkind) :: BTH0(NFRE) !saturation spectrum
REAL(rkind) :: BTH0S(NFRE) !smoothed saturation spectrum
REAL(rkind) :: BTH(NANG,NFRE) !saturation spectrum
REAL(rkind) :: BTHS(NANG,NFRE) !smoothed saturation spectrum
REAL(rkind),DIMENSION(NANG,NFRE) :: F,FL,SL,A, D
REAL(rkind),DIMENSION(NANG,NFRE) :: SSDS, DSSDS
REAL(rkind), ALLOCATABLE, DIMENSION(:) :: C_,C_C,C2_,C2_C2,DSIP_05_C2
REAL(rkind) , ALLOCATABLE :: SATWEIGHTS(:,:)
REAL(rkind) , ALLOCATABLE :: CUMULW(:,:,:,:)
LOGICAL :: LLTEST
!IF (LHOOK) CALL DR_HOOK('SDISSIP_ARDH_VEC',0,ZHOOK_HANDLE)
! ----------------------------------------------------------------------
W=26.
TPIINV = 1./ZPI
MIJ = MSC
ROG = ROWATER*G
IF (LCFLX) THEN
PHIEPS = 0.
TAUWD = 0.
! !!!! CONSTFM is only defined up to M=MIJ
SCDFM = 0.5*DELTH*(1.-1./FRATIO)
DO M=1,MIJ-1
CONSTFM(M) = ROG*DFIM(M)
ENDDO
CONSTFM = ROG*SCDFM*FR(MIJ)
ENDIF
XFR = 1.1 !! ??? JEAN BIDLOT: what is XFR is it FRATIO ????
SSDSBRF1 = 0.5
SXFR = 0.5*(FRATIO-1/FRATIO)
! SDSBR = 1.40E-3
! SDSBR = 1.20E-3 ! from Babanin (personal communication)
SDSBR = 9.0E-4
SDSBR2 = 0.8
SSDSBR = SDSBR
! SDSC1 = -2.1 !! This is Bidlot et al. 2005, Otherwise WAM4 uses -4.5
SDSC1 = 0.
SSDSC1 = SDSC1
SSDSC5 = 0.
SSDSISO = 1.
SSDSP = 2.
! SSDSC2 = -3.0E-5 ! Retuned by Fabrice from VdW
SSDSC2 = -2.8E-5
SSDSC3 = 0.8
SSDSC4 = 1.
SSDSC6 = 0.3
ISDSDTH = 80
SDSCOS = 2.
SSDSHF = 1.
SSDSLF = 1.
EPSR=SQRT(SSDSBR)
TMP00 = SSDSC1*ZPI
YEPS = ROAIR/ROWATER
TMP01 = SSDSC5/G*YEPS
DO M=1, NFRE
SIG(M) = ZPI*FR(M)
END DO
ALFAMEAN = (XKMEAN**2)*EMEAN(IPP)
FACTOR = TMP00 *F1MEAN*(ALFAMEAN*ALFAMEAN)
FACTURB = TMP01*USNEW(IPP)*USNEW(IPP)
WNMEAN2 = MAX( 1.E-10 , XKMEAN )
IF (ISHALLO.EQ.0) THEN
DO M=1, NFRE
XK(M) = WK0(M,IPP)!TFAK(INDEP,M)
CG(M) = CG0(M,IPP)!TCGOND(INDEP,M)
END DO
ELSE
DO M=1, NFRE
XK(M) = (SIG(M)**2)/G
CG(M) = G/(2*ZPI*FR(M))
END DO
ENDIF
DO M=1, NFRE
!cdir outerunroll=4
DO K=1, NANG
A(K,M) = TPIINV*CG(M)*F(K,M)
ENDDO
END DO
IF (ISDSDTH.NE.180) THEN
SDSNTH = MIN(NINT(ISDSDTH*DEGRAD/(DELTH)),NANG/2-1)
ALLOCATE(SATINDICES(NANG,SDSNTH*2+1))
ALLOCATE(SATWEIGHTS(NANG,SDSNTH*2+1))
DO K=1,NANG
DO I_INT=K-SDSNTH, K+SDSNTH
J_INT=I_INT
IF (I_INT.LT.1) J_INT=I_INT+NANG
IF (I_INT.GT.NANG) J_INT=I_INT-NANG
SATINDICES(K,I_INT-(K-SDSNTH)+1)=J_INT
SATWEIGHTS(K,I_INT-(K-SDSNTH)+1)=COS(TH(K)-TH(J_INT))**SDSCOS
END DO
END DO
END IF
! calcul de cumulw
SSDSBRF1 = 0.5
SXFR = 0.5*(FRATIO-1/FRATIO)
! initialise CUMULW
ALLOCATE(CUMULW(NFRE,NANG,NFRE,NANG))
DO I=1,NFRE
DO J=1,NANG
DO I1=1,NFRE
DO J1=1,NANG
CUMULW(I,J,I1,J1)=0.
ENDDO
ENDDO
ENDDO
END DO
IF (SSDSC3.NE.0.) THEN
! DIKCUMUL is the integer difference in frequency bands
! between the "large breakers" and short "wiped-out waves"
DIKCUMUL = NINT(SSDSBRF1/(XFR-1.))
ALLOCATE(C_(NFRE))
ALLOCATE(C_C(NFRE))
ALLOCATE(C2_(NFRE-DIKCUMUL))
ALLOCATE(C2_C2(NFRE-DIKCUMUL))
ALLOCATE(DSIP_05_C2(NFRE-DIKCUMUL))
DO M=1,NFRE
C_(M)=G/SIG(M) ! Valid in deep water only
C_C(M)=C_(M)*C_(M)
END DO
DO M=1,NFRE
DO K=1,NANG
DO M2=1,M-DIKCUMUL
C2_(M2)=G/SIG(M2)
C2_C2(M2)=C2_(M2)*C2_(M2)
DSIP_ = SIG(M2)*SXFR
DSIP_05_C2(M2)=DSIP_/(0.5*C2_(M2))
DO K2=1,NANG
CUMULW(M,K,M2,K2)=SQRT(C_C(M)+C2_C2(M2) &
& -2*C_(M)*C2_(M2)*COSTH(1+ABS(K2-K)))*DSIP_05_C2(M2)
END DO
END DO
END DO
END DO
DEALLOCATE(C_)
DEALLOCATE(C_C)
DEALLOCATE(C2_)
DEALLOCATE(C2_C2)
DEALLOCATE(DSIP_05_C2)
! Multiplies by lambda(k,theta)=1/(2*pi**2) and
! and the coefficient that transforms SQRT(B) to Banner et al. (2000)'s epsilon
! 2.26 is equal to 5.55 (Banner & al. 2000) times 1.6**2 / 2pi where
! 1.6 is the ratio between Banner's epsilon and SQRT(B)
TMP02 = 2*TPIINV*2.26
DO I=1,NFRE
DO J=1,NANG
DO I1=1,NFRE
DO J1=1,NANG
CUMULW(I,J,I1,J1)=CUMULW(I,J,I1,J1)*TMP02
END DO
ENDDO
END DO
END DO
END IF
DO M=1, NFRE
FACSAT=(XK(M)**3)*DELTH
BTH0(M)=SUM(A(1:NANG,M))*FACSAT
END DO
! DO K=1,NANG
DO M=1, NFRE
DO K=1,NANG
FACSAT=(XK(M)**3)*DELTH
! integrates around full circle
BTH(K,M) = SUM(SATWEIGHTS(K,1:SDSNTH*2+1)*A(SATINDICES(K,1:SDSNTH*2+1),M))*FACSAT
END DO
BTH0(M) = MAXVAL(BTH(1:NANG,M))
END DO
!/ST3 SDSBR = 1.20E-3 ! Babanin (personnal communication)
ISSDSBRFDF = 22 ! test pour DC
ISSDSBRFDF = 0
!/ST3 SDSBRF1 = 0.5
!/ST3 SDSBRF2 = 0.
IF (ISSDSBRFDF.GT.0.AND.ISSDSBRFDF.LT.NFRE/2) THEN
!cdir collapse
BTH0S=BTH0
!cdir collapse
NSMOOTH=ONE
!cdir collapse
BTHS=BTH
!cdir outerunroll=4
DO M=1, ISSDSBRFDF
BTH0S (1+ISSDSBRFDF)=BTH0S (1+ISSDSBRFDF)+BTH0(M)
NSMOOTH(1+ISSDSBRFDF)=NSMOOTH(1+ISSDSBRFDF)+1
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)+BTH(K,M)
END DO
ENDDO
DO M=2+ISSDSBRFDF,1+2*ISSDSBRFDF
!cdir nodep
BTH0S (1+ISSDSBRFDF)=BTH0S (1+ISSDSBRFDF)+BTH0(M)
NSMOOTH(1+ISSDSBRFDF)=NSMOOTH(1+ISSDSBRFDF)+1
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)+BTH(K,M)
END DO
ENDDO
DO M=ISSDSBRFDF,1,-1
!cdir nodep
BTH0S (M)=BTH0S (M+1)-BTH0(M+ISSDSBRFDF+1)
NSMOOTH(M)=NSMOOTH(M+1)-1
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)-BTH(K,M)
END DO
ENDDO
DO M=2+ISSDSBRFDF,NFRE-ISSDSBRFDF
!cdir nodep
BTH0S (M)=BTH0S (M-1)-BTH0(M-ISSDSBRFDF-1)+BTH0(M+ISSDSBRFDF)
NSMOOTH(M)=NSMOOTH(M-1)
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)-BTH(K,M)+BTH(K,M)
END DO
ENDDO
!cdir novector
DO M=NFRE-ISSDSBRFDF+1,NFRE
!cdir nodep
BTH0S (M)=BTH0S (M-1)-BTH0(M-ISSDSBRFDF)
NSMOOTH(M)=NSMOOTH(M-1)-1
!cdir collapse
DO K=1,NANG
BTHS(K,M)=BTHS(K,M)-BTH(K,M)
END DO
END DO
! final division by NSMOOTH
!cdir collapse
DO M=1,NFRE
BTH0(M)=MAX(0.,BTH0S(M)/NSMOOTH(M))
END DO
DO M=1,NFRE
!cdir outerunroll=4
DO K=1,NANG
BTH(K,M)=MAX(0.,BTHS(K,M)/NSMOOTH(M))
END DO
END DO
END IF
! DELTA1 = 0.4
! DELTA2 = 0.6
DELTA1 = 0.
DELTA2 = 0.
MICHE = 1.0
TMP03 = 1.0/(SSDSBR*MICHE)
DO M=1, NFRE
LLTEST = (SSDSC3.NE.0.AND.M.GT.DIKCUMUL)
IF (XKMEAN.NE.0) THEN
X = WK0(M,IPP)/XKMEAN!TFAK(INDEP(IJ),M)/XKMEAN(IJ)
BSIGBAJ = FACTOR*( (1.-DELTA2)*X + DELTA2*X**2)
ELSE
BSIGBAJ = 0
ENDIF
IF (ISHALLO.EQ.0) THEN
CM=WK0(M,IPP)/SIG(M)!TFAK(INDEP(IJ),M)/SIG(M)
ELSE
!AR: below is a nice bug always the last index ... in the original, should be IJ
CM=SIG(M)/G
ENDIF
DO K=1,NANG
RENEWALFREQ(K)=0.
! Correction of saturation level for shallow-water kinematics
! Cumulative effect based on lambda (breaking probability is
! the expected rate of sweeping by larger breaking waves)
IF (LLTEST) THEN
DO M2=1,M-DIKCUMUL
!AR: below the index is wrong ...
DO K2=1,NANG
IF (BTH0(M2).GT.SSDSBR) THEN
! Integrates over frequencies M2 and directions K2 to
! Integration is performed from M2=1 to a frequency lower than M: IK-DIKCUMUL
RENEWALFREQ(K)=RENEWALFREQ(K)+ CUMULW(M,K,M2,K2) &
& *(MAX(SQRT(BTH(K2,M2))-EPSR,0.))**2
ENDIF
END DO
END DO
ENDIF
SATURATION2=TANH(10*(((BTH(K,M)/SSDSBR)**0.5)-SDSBR2))
COSWIND=(COSTH(K)*COS(THWNEW(IPP))+SINTH(K)*SIN(THWNEW(IPP))) ! Ãvérifier K ?
DTURB=-2.*SIG(M)*XK(M)*FACTURB*COSWIND ! Theory -> stress direction
P0=SSDSP ! -0.5*SSDSC3*(1-TANH(W*USTAR*XK(M)/SIG(M)-0.1)) ! for SDSC3=1 this is vdW et al.
TMP04 = SSDSC3*RENEWALFREQ(K)
! DTEMP=SSDSC2 * SIG(M) &
DTEMP=SSDSC2 * SIG(M) &
& * ( SSDSC6 *(MAX(0.,BTH0( M)*TMP03-SSDSC4))**P0 &
& + (1-SSDSC6)*(MAX(0.,BTH (K,M)*TMP03-SSDSC4))**P0)&
& - (TMP04+DTURB) !terme cumulatif
D(K,M) = DTEMP + BSIGBAJ*SSDSLF *0.5*(1-SATURATION2) &
& + BSIGBAJ*SSDSHF *0.5*(SATURATION2+1)
WRITE(111116,'(10F15.8)') D(K,M),DTEMP,SIG(M),RENEWALFREQ(K),BTH0(M),BTH(K,M)
END DO
!cdir outerunroll=4
DO K=1, NANG
SL(K,M) = SL(K,M)+D(K,M)*F(K,M)
FL(K,M) = FL(K,M)+D(K,M)
SSDS(K,M) = D(K,M)*F(K,M)
DSSDS(K,M) = D(K,M)
IF (LCFLX.AND.M.LE.MIJ) THEN
SDISS = D(K,M)*F(K,M)
PHIEPS = PHIEPS+SDISS*CONSTFM(M)
TAUWD = TAUWD+CM*SDISS*CONSTFM(M)
ENDIF
END DO
END DO
IF (ALLOCATED(CUMULW)) DEALLOCATE (CUMULW)
IF (ALLOCATED(SATWEIGHTS)) DEALLOCATE (SATWEIGHTS)
IF (ALLOCATED(SATINDICES)) DEALLOCATE (SATINDICES)
!IF (LHOOK) CALL DR_HOOK('SDISSIP_ARDH_VEC',1,ZHOOK_HANDLE)
RETURN
END SUBROUTINE SDISS_ARDH_VEC_LOCAL