Page 1 Source Listing W3SPR4 2014-09-16 16:52 w3src4md.f90 1 !/ ------------------------------------------------------------------- / 2 MODULE W3SRC4MD 3 !/ 4 !/ +-----------------------------------+ 5 !/ | WAVEWATCH III SHOM | 6 !/ ! F. Ardhuin ! 7 !/ | FORTRAN 90 | 8 !/ | Last update : 08-Sep-2011 | 9 !/ +-----------------------------------+ 10 !/ 11 !/ 30-Aug-2010 : Origination. ( version 3.14-Ifremer ) 12 !/ 02-Nov-2010 : Addding fudge factor for low freq. ( version 4.03 ) 13 !/ 02-Sep-2011 : Clean up and time optimization ( version 4.04 ) 14 !/ 04-Sep-2011 : Estimation of whitecap stats. ( version 4.04 ) 15 !/ 16 ! 1. Purpose : 17 ! 18 ! The 'SHOM/Ifremer' source terms based on P.A.E.M. Janssen's wind input 19 ! and dissipation functions by Ardhuin et al. (2009,2010) 20 ! and Filipot & Ardhuin (2010) 21 ! The wind input is converted from the original 22 ! WAM codes, courtesy of P.A.E.M. Janssen and J. Bidlot 23 ! 24 ! 2. Variables and types : 25 ! 26 ! Name Type Scope Description 27 ! ---------------------------------------------------------------- 28 ! ---------------------------------------------------------------- 29 ! 30 ! 3. Subroutines and functions : 31 ! 32 ! Name Type Scope Description 33 ! ---------------------------------------------------------------- 34 ! W3SPR4 Subr. Public Mean parameters from spectrum. 35 ! W3SIN4 Subr. Public WAM4+ input source term. 36 ! INSIN4 Subr. Public Corresponding initialization routine. 37 ! TABU_STRESS, TABU_TAUHF, TABU_TAUHF2 38 ! Subr. Public Populate various tables. 39 ! CALC_USTAR 40 ! Subr. Public Compute stresses. 41 ! W3SDS4 Subr. Public Dissipation (Ardhuin & al. / Filipot & Ardhuin) 42 ! ---------------------------------------------------------------- 43 ! 44 ! 4. Subroutines and functions used : 45 ! 46 ! Name Type Module Description 47 ! ---------------------------------------------------------------- 48 ! ---------------------------------------------------------------- 49 ! 50 ! 5. Remarks : 51 ! 52 ! 6. Switches : 53 ! 54 ! 7. Source code : 55 !/ 56 !/ ------------------------------------------------------------------- / 57 !/ Page 2 Source Listing W3SPR4 2014-09-16 16:52 w3src4md.f90 58 PUBLIC 59 !/ 60 !/ Public variables 61 !/ 62 INTEGER, PARAMETER :: NHMAX = 25 63 INTEGER(KIND=4) :: NH(3), THO(2,3,NHMAX) 64 REAL :: HA(NHMAX,3), HD(NHMAX,3), HA2(NHMAX,3) 65 !air kinematic viscosity (used in WAM) 66 INTEGER, PARAMETER :: ITAUMAX=200,JUMAX=200 67 INTEGER, PARAMETER :: IUSTAR=100,IALPHA=200, ILEVTAIL=50 68 REAL :: TAUT(0:ITAUMAX,0:JUMAX), DELTAUW, DELU 69 ! Table for H.F. stress as a function of 2 variables 70 REAL :: TAUHFT(0:IUSTAR,0:IALPHA), DELUST, DELALP 71 ! Table for H.F. stress as a function of 3 variables 72 REAL :: TAUHFT2(0:IUSTAR,0:IALPHA,0:ILEVTAIL) 73 ! Table for swell damping 74 REAL :: DELTAIL 75 REAL, PARAMETER :: UMAX = 50. 76 REAL, PARAMETER :: TAUWMAX = 2.2361 !SQRT(5.) 77 REAL, DIMENSION(:) , ALLOCATABLE :: XSTRESS,YSTRESS 78 INTEGER :: DIKCUMUL 79 ! Size of wave height table for integrating the PDF of wave heights 80 INTEGER, PARAMETER :: NKHI=100 81 REAL, PARAMETER :: PI=3.14157, G=9.81 82 REAL, PARAMETER :: FAC_KD1=1.01, FAC_KD2=1000., KHSMAX=2., KHMAX=2. 83 REAL, PARAMETER ::KDMAX=200000. 84 ! variables for negative wind input (beta from ST2) 85 ! 86 INTEGER, PARAMETER, PRIVATE :: NRSIGA = 400 87 INTEGER, PARAMETER, PRIVATE :: NRDRAG = 20 88 REAL, PARAMETER, PRIVATE :: SIGAMX = 40. 89 REAL, PARAMETER, PRIVATE :: DRAGMX = 1.E-2 90 ! 91 REAL, PRIVATE :: DSIGA, DDRAG, & 92 BETATB(-NRSIGA:NRSIGA+1,NRDRAG+1) 93 !/ 94 LOGICAL, SAVE , PRIVATE :: FIRST = .TRUE. 95 !/ 96 CONTAINS 97 !/ ------------------------------------------------------------------- / 98 SUBROUTINE W3SPR4 (A, CG, WN, EMEAN, FMEAN, FMEAN1, WNMEAN, & 99 AMAX, U, UDIR, USTAR, USDIR, TAUWX, TAUWY, CD, Z0,& 100 CHARN, LLWS, FMEANWS) 101 !/ 102 !/ +-----------------------------------+ 103 !/ | WAVEWATCH III SHOM | 104 !/ ! F. Ardhuin ! 105 !/ | H. L. Tolman | 106 !/ | FORTRAN 90 | 107 !/ | Last update : 13-Jun-2011 | 108 !/ +-----------------------------------+ 109 !/ 110 !/ 03-Oct-2007 : Origination. ( version 3.13 ) 111 !/ 13-Jun-2011 : Adds f_m0,-1 as FMEAN in the outout ( version 4.04 ) 112 !/ 113 ! 1. Purpose : 114 ! Page 3 Source Listing W3SPR4 2014-09-16 16:52 w3src4md.f90 115 ! Calculate mean wave parameters for the use in the source term 116 ! routines. 117 ! 118 ! 2. Method : 119 ! 120 ! See source term routines. 121 ! 122 ! 3. Parameters : 123 ! 124 ! Parameter list 125 ! ---------------------------------------------------------------- 126 ! A R.A. I Action density spectrum. 127 ! CG R.A. I Group velocities. 128 ! WN R.A. I Wavenumbers. 129 ! EMEAN Real O Energy 130 ! FMEAN1 Real O Mean frequency (fm0,-1) used for reflection 131 ! FMEAN Real O Mean frequency for determination of tail 132 ! WNMEAN Real O Mean wavenumber. 133 ! AMAX Real O Maximum of action spectrum. 134 ! U Real I Wind speed. 135 ! UDIR Real I Wind direction. 136 ! USTAR Real I/O Friction velocity. 137 ! USDIR Real I/O wind stress direction. 138 ! TAUWX-Y Real I Components of wave-supported stress. 139 ! CD Real O Drag coefficient at wind level ZWND. 140 ! Z0 Real O Corresponding z0. 141 ! CHARN Real O Corresponding Charnock coefficient 142 ! LLWS L.A. I Wind sea true/false array for each component 143 ! FMEANWS Real O Mean frequency of wind sea, used for tail 144 ! ---------------------------------------------------------------- 145 ! 146 ! 4. Subroutines used : 147 ! 148 ! STRACE Service routine. 149 ! 150 ! 5. Called by : 151 ! 152 ! W3SRCE Source term integration routine. 153 ! W3OUTP Point output program. 154 ! GXEXPO GrADS point output program. 155 ! 156 ! 6. Error messages : 157 ! 158 ! 7. Remarks : 159 ! 160 ! 8. Structure : 161 ! 162 ! See source code. 163 ! 164 ! 9. Switches : 165 ! 166 ! !/S Enable subroutine tracing. 167 ! !/T Enable test output. 168 ! 169 ! 10. Source code : 170 ! 171 !/ ------------------------------------------------------------------- / Page 4 Source Listing W3SPR4 2014-09-16 16:52 w3src4md.f90 172 USE CONSTANTS 173 USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, DTH, DDEN, WWNMEANP, & 174 WWNMEANPTAIL, FTE, FTF, SSTXFTF, SSTXFTWN,& 175 SSTXFTFTAIL, SSWELLF 176 ! 177 IMPLICIT NONE 178 !/ 179 !/ ------------------------------------------------------------------- / 180 !/ Parameter list 181 !/ 182 REAL, INTENT(IN) :: A(NTH,NK), CG(NK), WN(NK), U, UDIR 183 REAL, INTENT(IN) :: TAUWX, TAUWY 184 LOGICAL, INTENT(IN) :: LLWS(NSPEC) 185 REAL, INTENT(INOUT) :: USTAR ,USDIR 186 REAL, INTENT(OUT) :: EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, & 187 CD, Z0, CHARN, FMEANWS 188 !/ 189 !/ ------------------------------------------------------------------- / 190 !/ Local parameters 191 !/ 192 INTEGER :: IKP1 ! wind sea peak index 193 INTEGER :: IS, IK, ITH, I1, ITT 194 195 REAL :: TAUW, EBAND, EMEANWS, RDCH, FXPMC, & 196 WNP, UNZ, FP, & 197 R1, CP, EB(NK),EB2(NK),ALFA(NK) 198 !/ 199 !/ ------------------------------------------------------------------- / 200 !/ 201 ! 202 UNZ = MAX ( 0.01 , U ) 203 USTAR = MAX ( 0.0001 , USTAR ) 204 ! 205 EMEAN = 0. 206 EMEANWS= 0. 207 FMEANWS= 0. 208 FMEAN = 0. 209 FMEAN1 = 0. 210 WNMEAN = 0. 211 AMAX = 0. 212 ! 213 ! 1. Integral over directions and maximum --------------------------- * 214 ! 215 DO IK=1, NK 216 EB(IK) = 0. 217 EB2(IK) = 0. 218 DO ITH=1, NTH 219 IS=ITH+(IK-1)*NTH 220 EB(IK) = EB(IK) + A(ITH,IK) 221 IF (LLWS(IS)) EB2(IK) = EB2(IK) + A(ITH,IK) 222 AMAX = MAX ( AMAX , A(ITH,IK) ) 223 END DO 224 END DO 225 ! 226 ! 2. Integrate over directions -------------------------------------- * 227 ! 228 DO IK=1, NK Page 5 Source Listing W3SPR4 2014-09-16 16:52 w3src4md.f90 229 ALFA(IK) = 2. * DTH * SIG(IK) * EB(IK) * WN(IK)**3 230 EB(IK) = EB(IK) * DDEN(IK) / CG(IK) 231 EB2(IK) = EB2(IK) * DDEN(IK) / CG(IK) 232 EMEAN = EMEAN + EB(IK) 233 FMEAN = FMEAN + EB(IK) /SIG(IK) 234 FMEAN1 = FMEAN1 + EB(IK) *(SIG(IK)**(2.*WWNMEANPTAIL)) 235 WNMEAN = WNMEAN + EB(IK) *(WN(IK)**WWNMEANP) 236 EMEANWS = EMEANWS+ EB2(IK) 237 FMEANWS = FMEANWS+ EB2(IK)*(SIG(IK)**(2.*WWNMEANPTAIL)) 238 END DO 239 ! 240 ! 3. Add tail beyond discrete spectrum and get mean pars ------------ * 241 ! ( DTH * SIG absorbed in FTxx ) 242 ! 243 EBAND = EB(NK) / DDEN(NK) 244 EMEAN = EMEAN + EBAND * FTE 245 FMEAN = FMEAN + EBAND * FTF 246 FMEAN1 = FMEAN1 + EBAND * SSTXFTFTAIL 247 WNMEAN = WNMEAN + EBAND * SSTXFTWN 248 EBAND = EB2(NK) / DDEN(NK) 249 EMEANWS = EMEANWS + EBAND * FTE 250 FMEANWS = FMEANWS + EBAND * SSTXFTFTAIL 251 ! 252 ! 4. Final processing 253 ! 254 FMEAN = TPIINV * EMEAN / MAX ( 1.E-7 , FMEAN ) 255 IF (FMEAN1.LT.1.E-7) THEN 256 FMEAN1=TPIINV * SIG(NK) 257 ELSE 258 FMEAN1 = TPIINV *( MAX ( 1.E-7 , FMEAN1 ) & 259 / MAX ( 1.E-7 , EMEAN ))**(1/(2.*WWNMEANPTAIL)) 260 ENDIF 261 WNMEAN = ( MAX ( 1.E-7 , WNMEAN ) & 262 / MAX ( 1.E-7 , EMEAN ) )**(1/WWNMEANP) 263 IF (FMEANWS.LT.1.E-7.OR.EMEANWS.LT.1.E-7) THEN 264 FMEANWS=TPIINV * SIG(NK) 265 ELSE 266 FMEANWS = TPIINV *( MAX ( 1.E-7 , FMEANWS ) & 267 / MAX ( 1.E-7 , EMEANWS ))**(1/(2.*WWNMEANPTAIL)) 268 END IF 269 270 ! 271 ! 5. Cd and z0 ----------------------------------------------- * 272 ! 273 TAUW = SQRT(TAUWX**2+TAUWY**2) 274 275 Z0=0. 276 CALL CALC_USTAR(U,TAUW,USTAR,Z0,CHARN) 277 UNZ = MAX ( 0.01 , U ) 278 CD = (USTAR/UNZ)**2 279 USDIR = UDIR 280 ! 281 ! 6. Final test output ---------------------------------------------- * 282 ! 283 RETURN 284 ! 285 ! Formats Page 6 Source Listing W3SPR4 2014-09-16 16:52 w3src4md.f90 286 ! 287 !/ 288 !/ End of W3SPR3 ----------------------------------------------------- / 289 !/ 290 END SUBROUTINE ENTRY POINTS Name w3src4md_mp_w3spr4_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References A Dummy 98 R(4) 4 2 0 ARG,IN 220,221,222 ALFA Local 197 R(4) 4 1 0 229 AMAX Dummy 99 R(4) 4 scalar ARG,OUT 211,222 CD Dummy 99 R(4) 4 scalar ARG,OUT 278 CG Dummy 98 R(4) 4 1 0 ARG,IN 230,231 CHARN Dummy 100 R(4) 4 scalar ARG,OUT 276 CONSTANTS Module 172 172 CP Local 197 R(4) 4 scalar DDEN Local 173 R(4) 4 1 1 PTR 173,230,231,243,248 DTH Local 173 R(4) 4 scalar PTR 173,229 EB Local 197 R(4) 4 1 0 216,220,229,230,232,233,234,235,24 3 EB2 Local 197 R(4) 4 1 0 217,221,231,236,237,248 EBAND Local 195 R(4) 4 scalar 243,244,245,246,247,248,249,250 EMEAN Dummy 98 R(4) 4 scalar ARG,OUT 205,232,244,254,259,262 EMEANWS Local 195 R(4) 4 scalar 206,236,249,263,267 FMEAN Dummy 98 R(4) 4 scalar ARG,OUT 208,233,245,254 FMEAN1 Dummy 98 R(4) 4 scalar ARG,OUT 209,234,246,255,256,258 FMEANWS Dummy 100 R(4) 4 scalar ARG,OUT 207,237,250,263,264,266 FP Local 196 R(4) 4 scalar FTE Local 174 R(4) 4 scalar PTR 174,244,249 FTF Local 174 R(4) 4 scalar PTR 174,245 FXPMC Local 195 R(4) 4 scalar I1 Local 193 I(4) 4 scalar IK Local 193 I(4) 4 scalar 215,216,217,219,220,221,222,228,22 9,230,231,232,233,234,235,236,237 IKP1 Local 192 I(4) 4 scalar IS Local 193 I(4) 4 scalar 219,221 ITH Local 193 I(4) 4 scalar 218,219,220,221,222 ITT Local 193 I(4) 4 scalar LLWS Dummy 100 L(4) 4 1 0 ARG,IN 221 MAX Func 202 scalar 202,203,222,254,258,259,261,262,26 6,267,277 NK Local 173 I(4) 4 scalar PTR 173,182,197,215,228,243,248,256,26 4 NSPEC Local 173 I(4) 4 scalar PTR 173,184 NTH Local 173 I(4) 4 scalar PTR 173,182,218,219 R1 Local 197 R(4) 4 scalar RDCH Local 195 R(4) 4 scalar Page 7 Source Listing W3SPR4 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References SIG Local 173 R(4) 4 1 1 PTR 173,229,233,234,237,256,264 SQRT Func 273 scalar 273 SSTXFTF Local 174 R(4) 4 scalar PTR 174 SSTXFTFTAIL Local 175 R(4) 4 scalar PTR 175,246,250 SSTXFTWN Local 174 R(4) 4 scalar PTR 174,247 SSWELLF Local 175 R(4) 4 1 1 PTR 175 TAUW Local 195 R(4) 4 scalar 273,276 TAUWX Dummy 99 R(4) 4 scalar ARG,IN 273 TAUWY Dummy 99 R(4) 4 scalar ARG,IN 273 TPIINV Param 254 R(4) 4 scalar 254,256,258,264,266 U Dummy 99 R(4) 4 scalar ARG,IN 202,276,277 UDIR Dummy 99 R(4) 4 scalar ARG,IN 279 UNZ Local 196 R(4) 4 scalar 202,277,278 USDIR Dummy 99 R(4) 4 scalar ARG,INOUT 279 USTAR Dummy 99 R(4) 4 scalar ARG,INOUT 203,276,278 W3GDATMD Module 173 173 W3SPR4 Subr 98 WN Dummy 98 R(4) 4 1 0 ARG,IN 229,235 WNMEAN Dummy 98 R(4) 4 scalar ARG,OUT 210,235,247,261 WNP Local 196 R(4) 4 scalar WWNMEANP Local 173 R(4) 4 scalar PTR 173,235,262 WWNMEANPTAIL Local 174 R(4) 4 scalar PTR 174,234,237,259,267 Z0 Dummy 99 R(4) 4 scalar ARG,OUT 275,276 Page 8 Source Listing W3SPR4 2014-09-16 16:52 w3src4md.f90 291 !/ ------------------------------------------------------------------- / 292 SUBROUTINE W3SIN4 (A, CG, K, U, USTAR, DRAT, AS, USDIR, Z0, CD, & 293 TAUWX, TAUWY, TAUWNX, TAUWNY, ICE, S, D, LLWS, IX, IY) 294 !/ 295 !/ +-----------------------------------+ 296 !/ | WAVEWATCH III SHOM | 297 !/ ! F. Ardhuin ! 298 !/ | H. L. Tolman | 299 !/ | FORTRAN 90 | 300 !/ | Last update : 16-May-2010 | 301 !/ +-----------------------------------+ 302 !/ 303 !/ 09-Oct-2007 : Origination. ( version 3.13 ) 304 !/ 16-May-2010 : Adding sea ice ( version 3.14_Ifremer ) 305 !/ 306 ! 1. Purpose : 307 ! 308 ! Calculate diagonal and input source term for WAM4+ approach. 309 ! 310 ! 2. Method : 311 ! 312 ! WAM-4 : Janssen et al. 313 ! WAM-"4.5" : gustiness effect (Cavaleri et al. ) 314 ! SAT : high-frequency input reduction for balance with 315 ! saturation dissipation (Ardhuin et al., 2008) 316 ! SWELL : negative wind input (Ardhuin et al. 2008) 317 ! 318 ! 3. Parameters : 319 ! 320 ! Parameter list 321 ! ---------------------------------------------------------------- 322 ! A R.A. I Action density spectrum (1-D). 323 ! CG R.A. I Group speed *) 324 ! K R.A. I Wavenumber for entire spectrum. *) 325 ! U Real I WIND SPEED 326 ! USTAR Real I Friction velocity. 327 ! DRAT Real I Air/water density ratio. 328 ! AS Real I Air-sea temperature difference 329 ! USDIR Real I wind stress direction 330 ! Z0 Real I Air-side roughness lengh. 331 ! CD Real I Wind drag coefficient. 332 ! USDIR Real I Direction of friction velocity 333 ! TAUWX-Y Real I Components of the wave-supported stress. 334 ! TAUWNX Real I Component of the negative wave-supported stress. 335 ! TAUWNY Real I Component of the negative wave-supported stress. 336 ! ICE Real I Sea ice fraction. !/Stefan: ICE is DUMMY argument; remove later. 337 ! S R.A. O Source term (1-D version). 338 ! D R.A. O Diagonal term of derivative. *) 339 ! ---------------------------------------------------------------- 340 ! *) Stored as 1-D array with dimension NTH*NK 341 ! 342 ! 4. Subroutines used : 343 ! 344 ! STRACE Subroutine tracing. ( !/S switch ) 345 ! PRT2DS Print plot of spectrum. ( !/T0 switch ) 346 ! OUTMAT Print out matrix. ( !/T1 switch ) 347 ! Page 9 Source Listing W3SIN4 2014-09-16 16:52 w3src4md.f90 348 ! 5. Called by : 349 ! 350 ! W3SRCE Source term integration. 351 ! W3EXPO Point output program. 352 ! GXEXPO GrADS point output program. 353 ! 354 ! 6. Error messages : 355 ! 356 ! 7. Remarks : 357 ! 358 ! 8. Structure : 359 ! 360 ! See source code. 361 ! 362 ! 9. Switches : 363 ! 364 ! !/S Enable subroutine tracing. 365 ! !/T Enable general test output. 366 ! !/T0 2-D print plot of source term. 367 ! !/T1 Print arrays. 368 ! 369 ! 10. Source code : 370 ! 371 !/ ------------------------------------------------------------------- / 372 USE CONSTANTS 373 USE W3GDATMD, ONLY: NK, NTH, NSPEC, XFR, DDEN, SIG, SIG2, TH, & 374 ESIN, ECOS, EC2, ZZWND, AALPHA, BBETA, ZZALP,& 375 TTAUWSHELTER, SSWELLF, DDEN2, DTH, SSINTHP, & 376 ZZ0RAT, FLICES 377 ! 378 IMPLICIT NONE 379 !/ 380 !/ ------------------------------------------------------------------- / 381 !/ Parameter list 382 !/ 383 REAL, INTENT(IN) :: A(NSPEC) 384 REAL, INTENT(IN) :: CG(NK), K(NSPEC),Z0,U, CD 385 REAL, INTENT(IN) :: USTAR, USDIR, AS, DRAT, ICE !/Stefan: ICE is DUMMY argument; remove later. 386 REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC), TAUWX, TAUWY, TAUWNX, TAUWNY 387 LOGICAL, INTENT(OUT) :: LLWS(NSPEC) 388 INTEGER, INTENT(IN) :: IX, IY 389 !/ 390 !/ ------------------------------------------------------------------- / 391 !/ Local parameters 392 !/ 393 INTEGER :: IS,IK,ITH, IOMA, ICL 394 REAL :: FACLN1, FACLN2, ULAM, CLAM, OMA, & 395 RD1, RD2, LAMBDA, COSFAC 396 REAL :: COSU, SINU, TAUX, TAUY, USDIRP, USTP 397 REAL :: TAUPX, TAUPY, UST2, TAUW, TAUWB 398 REAL , PARAMETER :: EPS1 = 0.00001, EPS2 = 0.000001 399 REAL :: Usigma !standard deviation of U due to gustiness 400 REAL :: USTARsigma !standard deviation of USTAR due to gustiness 401 REAL :: BETA, mu_janssen, omega_janssen, & 402 CM,ZCO,UCO,UCN,ZCN, & 403 Z0VISC, Z0NOZ, EB, & 404 EBX, EBY, AORB, AORB1, FW, UORB, M2, TH2, & Page 10 Source Listing W3SIN4 2014-09-16 16:52 w3src4md.f90 405 RE, FU, FUD, SWELLCOEFV, SWELLCOEFT 406 REAL :: HSBLOW, ABJSEA, FACTOR 407 REAL :: PTURB, PVISC, SMOOTH 408 REAL XI,DELI1,DELI2 409 REAL XJ,DELJ1,DELJ2 410 REAL XK,DELK1,DELK2 411 REAL :: CONST, CONST0, CONST2, TAU1 412 REAL X,ZARG,ZLOG,UST 413 REAL COSWIND,XSTRESS,YSTRESS,TAUHF 414 REAL TEMP, TEMP2 415 INTEGER IND,J,I,ISTAB 416 REAL DSTAB(3,NSPEC), DVISC, DTURB 417 REAL STRESSSTAB(3,2),STRESSSTABN(3,2) 418 !/ 419 !/ ------------------------------------------------------------------- / 420 !/ 421 ! 422 ! 1. Preparations 423 ! 424 ! 1.a estimation of surface roughness parameters 425 ! 426 Z0VISC = 0.1*nu_air/MAX(USTAR,0.0001) 427 Z0NOZ = MAX(Z0VISC,ZZ0RAT*Z0) 428 FACLN1 = U / LOG(ZZWND/Z0NOZ) 429 FACLN2 = LOG(Z0NOZ) 430 ! 431 ! 1.b estimation of surface orbital velocity and displacement 432 ! 433 UORB=0. 434 AORB=0. 435 436 DO IK=1, NK 437 EB = 0. 438 EBX = 0. 439 EBY = 0. 440 DO ITH=1, NTH 441 IS=ITH+(IK-1)*NTH 442 EB = EB + A(IS) 443 END DO 444 ! 445 ! At this point UORB and AORB are the variances of the orbital velocity and surface elevation 446 ! 447 UORB = UORB + EB *SIG(IK)**2 * DDEN(IK) / CG(IK) 448 AORB = AORB + EB * DDEN(IK) / CG(IK) !deep water only 449 END DO 450 451 UORB = 2*SQRT(UORB) ! significant orbital amplitude 452 AORB1 = 2*AORB**(1-0.5*SSWELLF(6)) ! half the significant wave height ... if SWELLF(6)=1 453 RE = 4*UORB*AORB1 / NU_AIR ! Reynolds number 454 ! 455 ! Defines the swell dissipation based on the "Reynolds number" 456 ! 457 IF (SSWELLF(4).GT.0) THEN 458 IF (SSWELLF(7).GT.0.) THEN 459 SMOOTH = 0.5*TANH((RE-SSWELLF(4))/SSWELLF(7)) 460 PTURB=(0.5+SMOOTH) 461 PVISC=(0.5-SMOOTH) Page 11 Source Listing W3SIN4 2014-09-16 16:52 w3src4md.f90 462 ELSE 463 IF (RE.LE.SSWELLF(4)) THEN 464 PTURB = 0. 465 PVISC = 1. 466 ELSE 467 PTURB = 1. 468 PVISC = 0. 469 END IF 470 END IF 471 ELSE 472 PTURB=1. 473 PVISC=1. 474 END IF 475 476 ! 477 IF (SSWELLF(2).EQ.0) THEN 478 FW=MAX(ABS(SSWELLF(3)),0.) 479 FU=0. 480 FUD=0. 481 ELSE 482 FU=ABS(SSWELLF(3)) 483 FUD=SSWELLF(2) 484 AORB=2*SQRT(AORB) 485 XI=(ALOG10(MAX(AORB/Z0NOZ,3.))-ABMIN)/DELAB 486 IND = MIN (SIZEFWTABLE-1, INT(XI)) 487 DELI1= MIN (1. ,XI-FLOAT(IND)) 488 DELI2= 1. - DELI1 489 FW =FWTABLE(IND)*DELI2+FWTABLE(IND+1)*DELI1 490 END IF 491 ! 492 ! 2. Diagonal 493 ! 494 ! Here AS is the air-sea temperature difference in degrees. Expression given by 495 ! Abdalla & Cavaleri, JGR 2002 for Usigma. For USTARsigma ... I do not see where 496 ! I got it from, maybe just made up from drag law ... 497 ! 498 UST=USTAR 499 ISTAB=3 500 TAUX = UST**2* COS(USDIR) 501 TAUY = UST**2* SIN(USDIR) 502 ! 503 ! Loop over the resolved part of the spectrum 504 ! 505 STRESSSTAB(ISTAB,:)=0. 506 STRESSSTABN(ISTAB,:)=0. 507 ! 508 ! Coupling coefficient times density ratio DRAT and fraction of free surface (1-ICE) 509 ! Stefan Zieger: Coupling is permanently disabled and MIN(0.,MAX(1.,1.-ICE)) has 510 ! FIXME no effect because it is always 0.0; remove later. 511 ! Implemented in W3SRCE for all source terms. 512 ! 513 IF (FLICES) THEN 514 CONST0=MIN(0.,MAX(1.,1.-ICE))*BBETA*DRAT/(kappa**2) 515 ELSE 516 CONST0=BBETA*DRAT/(kappa**2) 517 END IF 518 DO IK=1, NK Page 12 Source Listing W3SIN4 2014-09-16 16:52 w3src4md.f90 519 TAUPX=TAUX-ABS(TTAUWSHELTER)*STRESSSTAB(ISTAB,1) 520 TAUPY=TAUY-ABS(TTAUWSHELTER)*STRESSSTAB(ISTAB,2) 521 ! With MIN and MAX the bug should disappear.... but where did it come from? 522 USTP=MIN((TAUPX**2+TAUPY**2)**0.25,MAX(UST,0.3)) 523 USDIRP=ATAN2(TAUPY,TAUPX) 524 COSU = COS(USDIRP) 525 SINU = SIN(USDIRP) 526 IS=1+(IK-1)*NTH 527 CM=K(IS)/SIG2(IS) !inverse of phase speed 528 UCN=USTP*CM+ZZALP !this is the inverse wave age 529 ! the stress is the real stress (N/m^2) divided by 530 ! rho_a, and thus comparable to USTAR**2 531 ! it is the integral of rho_w g Sin/C /rho_a 532 ! (air-> waves momentum flux) 533 CONST2=DDEN2(IS)/CG(IK) & !Jacobian to get energy in band 534 *GRAV/(SIG(IK)/K(IS)*DRAT) ! coefficient to get momentum 535 CONST=SIG2(IS)*CONST0 536 ! CM parameter is 1 / C_phi 537 ! Z0 corresponds to Z0+Z1 of the Janssen eq. 14 538 ZCN=ALOG(K(IS)*Z0) 539 ! below is the original WAM version (OK for deep water) g*z0/C^2 540 ! ZCN=ALOG(G*Z0b(I)*CM(I)**2) 541 ! 542 ! precomputes swell factors 543 ! 544 SWELLCOEFV=-SSWELLF(5)*DRAT*2*K(IS)*SQRT(2*NU_AIR*SIG2(IS)) 545 SWELLCOEFT=-DRAT*SSWELLF(1)*16*SIG2(IS)**2/GRAV 546 ! 547 DO ITH=1,NTH 548 IS=ITH+(IK-1)*NTH 549 COSWIND=(ECOS(IS)*COSU+ESIN(IS)*SINU) 550 IF (COSWIND.GT.0.01) THEN 551 X=COSWIND*UCN 552 ! this ZARG term is the argument of the exponential 553 ! in Janssen 1991 eq. 16. 554 ZARG=KAPPA/X 555 ! ZLOG is ALOG(MU) where MU is defined by Janssen 1991 eq. 15 556 ! MU= 557 ZLOG=ZCN+ZARG 558 559 IF (ZLOG.LT.0.) THEN 560 ! The source term Sp is beta * omega * X**2 561 ! as given by Janssen 1991 eq. 19 562 ! for a faster performance EXP(X)*X**4 should be tabulated 563 DSTAB(ISTAB,IS) = CONST*EXP(ZLOG)*ZLOG**4*UCN*UCN*COSWIND**SSINTHP 564 LLWS(IS)=.TRUE. 565 ELSE 566 DSTAB(ISTAB,IS) = 0. 567 LLWS(IS)=.FALSE. 568 END IF 569 ! 570 ! Added for consistency with ECWAM implsch.F 571 ! 572 IF (28.*CM*USTAR*COSWIND.GE.1) THEN 573 LLWS(IS)=.TRUE. 574 END IF 575 ELSE ! (COSWIND.LE.0.01) Page 13 Source Listing W3SIN4 2014-09-16 16:52 w3src4md.f90 576 DSTAB(ISTAB,IS) = 0. 577 LLWS(IS)=.FALSE. 578 END IF 579 ! 580 IF ((SSWELLF(1).NE.0.AND.DSTAB(ISTAB,IS).LT.1E-7*SIG2(IS)) & 581 .OR.SSWELLF(3).GT.0) THEN 582 ! 583 DVISC=SWELLCOEFV 584 DTURB=SWELLCOEFT*(FW*UORB+(FU+FUD*COSWIND)*USTP) 585 ! 586 DSTAB(ISTAB,IS) = DSTAB(ISTAB,IS) + PTURB*DTURB + PVISC*DVISC 587 END IF 588 ! 589 ! Sums up the wave-supported stress 590 ! 591 ! Wave direction is "direction to" 592 ! therefore there is a PLUS sign for the stress 593 TEMP2=CONST2*DSTAB(ISTAB,IS)*A(IS) 594 IF (DSTAB(ISTAB,IS).LT.0) THEN 595 STRESSSTABN(ISTAB,1)=STRESSSTABN(ISTAB,1)+TEMP2*ECOS(IS) 596 STRESSSTABN(ISTAB,2)=STRESSSTABN(ISTAB,2)+TEMP2*ESIN(IS) 597 ELSE 598 STRESSSTAB(ISTAB,1)=STRESSSTAB(ISTAB,1)+TEMP2*ECOS(IS) 599 STRESSSTAB(ISTAB,2)=STRESSSTAB(ISTAB,2)+TEMP2*ESIN(IS) 600 END IF 601 END DO 602 END DO 603 ! 604 D(:)=DSTAB(3,:) 605 XSTRESS=STRESSSTAB (3,1) 606 YSTRESS=STRESSSTAB (3,2) 607 TAUWNX =STRESSSTABN(3,1) 608 TAUWNY =STRESSSTABN(3,2) 609 S = D * A 610 ! 611 ! ... Test output of arrays 612 ! 613 ! Computes the high-frequency contribution 614 ! the difference in spectal density (kx,ky) to (f,theta) 615 ! is integrated in this modified CONST0 616 CONST0=DTH*SIG(NK)**5/((GRAV**2)*tpi) & 617 *TPI*SIG(NK) / CG(NK) !conversion WAM (E(f,theta) to WW3 A(k,theta) 618 TEMP=0. 619 DO ITH=1,NTH 620 IS=ITH+(NK-1)*NTH 621 COSWIND=(ECOS(IS)*COSU+ESIN(IS)*SINU) 622 TEMP=TEMP+A(IS)*(MAX(COSWIND,0.))**3 623 END DO 624 625 TAUPX=TAUX-ABS(TTAUWSHELTER)*XSTRESS 626 TAUPY=TAUY-ABS(TTAUWSHELTER)*YSTRESS 627 USTP=(TAUPX**2+TAUPY**2)**0.25 628 USDIRP=ATAN2(TAUPY,TAUPX) 629 630 UST=USTP 631 ! finds the values in the tabulated stress TAUHFT 632 XI=UST/DELUST Page 14 Source Listing W3SIN4 2014-09-16 16:52 w3src4md.f90 633 IND = MAX(1,MIN (IUSTAR-1, INT(XI))) 634 DELI1= MAX(MIN (1. ,XI-FLOAT(IND)),0.) 635 DELI2= 1. - DELI1 636 XJ=MAX(0.,(GRAV*Z0/MAX(UST,0.00001)**2-AALPHA) / DELALP) 637 J = MAX(1 ,MIN (IALPHA-1, INT(XJ))) 638 DELJ1= MAX(0.,MIN (1. , XJ-FLOAT(J))) 639 DELJ2=1. - DELJ1 640 IF (TTAUWSHELTER.GT.0) THEN 641 XK = CONST0*TEMP / DELTAIL 642 I = MIN (ILEVTAIL-1, INT(XK)) 643 DELK1= MIN (1. ,XK-FLOAT(I)) 644 DELK2=1. - DELK1 645 TAU1 =((TAUHFT2(IND,J,I)*DELI2+TAUHFT2(IND+1,J,I)*DELI1 )*DELJ2 & 646 +(TAUHFT2(IND,J+1,I)*DELI2+TAUHFT2(IND+1,J+1,I)*DELI1)*DELJ1)*DELK2 & 647 +((TAUHFT2(IND,J,I+1)*DELI2+TAUHFT2(IND+1,J,I+1)*DELI1 )*DELJ2 & 648 +(TAUHFT2(IND,J+1,I+1)*DELI2+TAUHFT2(IND+1,J+1,I+1)*DELI1)*DELJ1)*DELK1 649 ELSE 650 TAU1 =(TAUHFT(IND,J)*DELI2+TAUHFT(IND+1,J)*DELI1 )*DELJ2 & 651 +(TAUHFT(IND,J+1)*DELI2+TAUHFT(IND+1,J+1)*DELI1)*DELJ1 652 END IF 653 TAUHF = CONST0*TEMP*UST**2*TAU1 654 TAUWX = XSTRESS+TAUHF*COS(USDIRP) 655 TAUWY = YSTRESS+TAUHF*SIN(USDIRP) 656 ! 657 ! Reduces tail effect to make sure that wave-supported stress 658 ! is less than total stress, this is borrowed from ECWAM Stresso.F 659 ! 660 TAUW = SQRT(TAUWX**2+TAUWY**2) 661 UST2 = MAX(USTAR,EPS2)**2 662 TAUWB = MIN(TAUW,MAX(UST2-EPS1,EPS2**2)) 663 IF (TAUWB.LT.TAUW) THEN 664 TAUWX=TAUWX*TAUWB/TAUW 665 TAUWY=TAUWY*TAUWB/TAUW 666 END IF 667 ! 668 RETURN 669 ! 670 ! Formats 671 ! 672 !/ 673 !/ End of W3SIN3 ----------------------------------------------------- / 674 !/ 675 END SUBROUTINE Page 15 Source Listing W3SIN4 2014-09-16 16:52 Entry Points w3src4md.f90 ENTRY POINTS Name w3src4md_mp_w3sin4_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References A Dummy 292 R(4) 4 1 0 ARG,IN 442,593,609,622 AALPHA Local 374 R(4) 4 scalar PTR 374,636 ABJSEA Local 406 R(4) 4 scalar ABMIN Param 485 R(4) 4 scalar 485 ABS Func 478 scalar 478,482,519,520,625,626 ALOG Func 538 scalar 538 ALOG10 Func 485 scalar 485 AORB Local 404 R(4) 4 scalar 434,448,452,484,485 AORB1 Local 404 R(4) 4 scalar 452,453 AS Dummy 292 R(4) 4 scalar ARG,IN ATAN2 Func 523 scalar 523,628 BBETA Local 374 R(4) 4 scalar PTR 374,514,516 BETA Local 401 R(4) 4 scalar CD Dummy 292 R(4) 4 scalar ARG,IN CG Dummy 292 R(4) 4 1 0 ARG,IN 447,448,533,617 CLAM Local 394 R(4) 4 scalar CM Local 402 R(4) 4 scalar 527,528,572 CONST Local 411 R(4) 4 scalar 535,563 CONST0 Local 411 R(4) 4 scalar 514,516,535,616,641,653 CONST2 Local 411 R(4) 4 scalar 533,593 CONSTANTS Module 372 372 COS Func 500 scalar 500,524,654 COSFAC Local 395 R(4) 4 scalar COSU Local 396 R(4) 4 scalar 524,549,621 COSWIND Local 413 R(4) 4 scalar 549,550,551,563,572,584,621,622 D Dummy 293 R(4) 4 1 0 ARG,OUT 604,609 DDEN Local 373 R(4) 4 1 1 PTR 373,447,448 DDEN2 Local 375 R(4) 4 1 1 PTR 375,533 DELAB Local 485 R(4) 4 scalar 485 DELALP Local 636 R(4) 4 scalar 636,1144,1161,1283,1300 DELI1 Local 408 R(4) 4 scalar 487,488,489,634,635,645,646,647,64 8,650,651 DELI2 Local 408 R(4) 4 scalar 488,489,635,645,646,647,648,650,65 1 DELJ1 Local 409 R(4) 4 scalar 638,639,646,648,651 DELJ2 Local 409 R(4) 4 scalar 639,645,647,650 DELK1 Local 410 R(4) 4 scalar 643,644,648 DELK2 Local 410 R(4) 4 scalar 644,646 DELTAIL Local 641 R(4) 4 scalar 641,1284,1308 DELUST Local 632 R(4) 4 scalar 632,1143,1158,1282,1297 DRAT Dummy 292 R(4) 4 scalar ARG,IN 514,516,534,544,545 DSTAB Local 416 R(4) 4 2 0 563,566,576,580,586,593,594,604 DTH Local 375 R(4) 4 scalar PTR 375,616 DTURB Local 416 R(4) 4 scalar 584,586 Page 16 Source Listing W3SIN4 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References DVISC Local 416 R(4) 4 scalar 583,586 EB Local 403 R(4) 4 scalar 437,442,447,448 EBX Local 404 R(4) 4 scalar 438 EBY Local 404 R(4) 4 scalar 439 EC2 Local 374 R(4) 4 1 1 PTR 374 ECOS Local 374 R(4) 4 1 1 PTR 374,549,595,598,621 EPS1 Param 398 R(4) 4 scalar 662 EPS2 Param 398 R(4) 4 scalar 661,662 ESIN Local 374 R(4) 4 1 1 PTR 374,549,596,599,621 EXP Func 563 scalar 563 FACLN1 Local 394 R(4) 4 scalar 428 FACLN2 Local 394 R(4) 4 scalar 429 FACTOR Local 406 R(4) 4 scalar FLICES Local 376 L(4) 4 scalar PTR 376,513 FLOAT Func 487 scalar 487,634,638,643 FU Local 405 R(4) 4 scalar 479,482,584 FUD Local 405 R(4) 4 scalar 480,483,584 FW Local 404 R(4) 4 scalar 478,489,584 FWTABLE Local 489 R(4) 4 1 301 489 GRAV Param 534 R(4) 4 scalar 534,545,616,636 HSBLOW Local 406 R(4) 4 scalar I Local 415 I(4) 4 scalar 642,643,645,646,647,648 IALPHA Param 637 I(4) 4 scalar 70,72,637,1144,1148,1156,1283,1288 ,1298 ICE Dummy 293 R(4) 4 scalar ARG,IN 514 ICL Local 393 I(4) 4 scalar IK Local 393 I(4) 4 scalar 436,441,447,448,518,526,533,534,54 8 ILEVTAIL Param 642 I(4) 4 scalar 72,642,1284,1307 IND Local 415 I(4) 4 scalar 486,487,489,633,634,645,646,647,64 8,650,651 INT Func 486 scalar 486,633,637,642 IOMA Local 393 I(4) 4 scalar IS Local 393 I(4) 4 scalar 441,442,526,527,533,534,535,538,54 4,545,548,549,563,564,566,567,573, 576,577,580,586,593,594,595,596,59 8,599,620,621,622 ISTAB Local 415 I(4) 4 scalar 499,505,506,519,520,563,566,576,58 0,586,593,594,595,596,598,599 ITH Local 393 I(4) 4 scalar 440,441,547,548,619,620 IUSTAR Param 633 I(4) 4 scalar 70,72,633,1143,1148,1157,1282,1288 ,1296 IX Dummy 293 I(4) 4 scalar ARG,IN IY Dummy 293 I(4) 4 scalar ARG,IN J Local 415 I(4) 4 scalar 637,638,645,646,647,648,650,651 K Dummy 292 R(4) 4 1 0 ARG,IN 527,534,538,544 KAPPA Param 514 R(4) 4 scalar 514,516,554 LAMBDA Local 395 R(4) 4 scalar LLWS Dummy 293 L(4) 4 1 0 ARG,OUT 564,567,573,577 LOG Func 428 scalar 428,429 M2 Local 404 R(4) 4 scalar MAX Func 426 scalar 426,427,478,485,514,522,622,633,63 4,636,637,638,661,662 MIN Func 486 scalar 486,487,514,522,633,634,637,638,64 2,643,662 Page 17 Source Listing W3SIN4 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References MU_JANSSEN Local 401 R(4) 4 scalar NK Local 373 I(4) 4 scalar PTR 373,384,436,518,616,617,620 NSPEC Local 373 I(4) 4 scalar PTR 373,383,384,386,387,416 NTH Local 373 I(4) 4 scalar PTR 373,440,441,526,547,548,619,620 NU_AIR Param 426 R(4) 4 scalar 426,453,544 OMA Local 394 R(4) 4 scalar OMEGA_JANSSEN Local 401 R(4) 4 scalar PTURB Local 407 R(4) 4 scalar 460,464,467,472,586 PVISC Local 407 R(4) 4 scalar 461,465,468,473,586 RD1 Local 395 R(4) 4 scalar RD2 Local 395 R(4) 4 scalar RE Local 405 R(4) 4 scalar 453,459,463 S Dummy 293 R(4) 4 1 0 ARG,OUT 609 SIG Local 373 R(4) 4 1 1 PTR 373,447,534,616,617 SIG2 Local 373 R(4) 4 1 1 PTR 373,527,535,544,545,580 SIN Func 501 scalar 501,525,655 SINU Local 396 R(4) 4 scalar 525,549,621 SIZEFWTABLE Param 486 I(4) 4 scalar 486 SMOOTH Local 407 R(4) 4 scalar 459,460,461 SQRT Func 451 scalar 451,484,544,660 SSINTHP Local 375 R(4) 4 scalar PTR 375,563 SSWELLF Local 375 R(4) 4 1 1 PTR 375,452,457,458,459,463,477,478,48 2,483,544,545,580,581 STRESSSTAB Local 417 R(4) 4 2 6 505,519,520,598,599,605,606 STRESSSTABN Local 417 R(4) 4 2 6 506,595,596,607,608 SWELLCOEFT Local 405 R(4) 4 scalar 545,584 SWELLCOEFV Local 405 R(4) 4 scalar 544,583 TANH Func 459 scalar 459 TAU1 Local 411 R(4) 4 scalar 645,650,653 TAUHF Local 413 R(4) 4 scalar 653,654,655 TAUHFT Local 650 R(4) 4 2 20301 650,651,1148,1180,1288,1309,1325 TAUHFT2 Local 645 R(4) 4 3 1035351 645,646,647,648,1310,1332 TAUPX Local 397 R(4) 4 scalar 519,522,523,625,627,628 TAUPY Local 397 R(4) 4 scalar 520,522,523,626,627,628 TAUW Local 397 R(4) 4 scalar 660,662,663,664,665 TAUWB Local 397 R(4) 4 scalar 662,663,664,665 TAUWNX Dummy 293 R(4) 4 scalar ARG,OUT 607 TAUWNY Dummy 293 R(4) 4 scalar ARG,OUT 608 TAUWX Dummy 293 R(4) 4 scalar ARG,OUT 654,660,664 TAUWY Dummy 293 R(4) 4 scalar ARG,OUT 655,660,665 TAUX Local 396 R(4) 4 scalar 500,519,625 TAUY Local 396 R(4) 4 scalar 501,520,626 TEMP Local 414 R(4) 4 scalar 618,622,641,653 TEMP2 Local 414 R(4) 4 scalar 593,595,596,598,599 TH Local 373 R(4) 4 1 1 PTR 373 TH2 Local 404 R(4) 4 scalar TPI Param 616 R(4) 4 scalar 616,617 TTAUWSHELTER Local 375 R(4) 4 scalar PTR 375,519,520,625,626,640 U Dummy 292 R(4) 4 scalar ARG,IN 428 UCN Local 402 R(4) 4 scalar 528,551,563 UCO Local 402 R(4) 4 scalar ULAM Local 394 R(4) 4 scalar UORB Local 404 R(4) 4 scalar 433,447,451,453,584 USDIR Dummy 292 R(4) 4 scalar ARG,IN 500,501 USDIRP Local 396 R(4) 4 scalar 523,524,525,628,654,655 Page 18 Source Listing W3SIN4 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References USIGMA Local 399 R(4) 4 scalar UST Local 412 R(4) 4 scalar 498,500,501,522,630,632,636,653 UST2 Local 397 R(4) 4 scalar 661,662 USTAR Dummy 292 R(4) 4 scalar ARG,IN 426,498,572,661 USTARSIGMA Local 400 R(4) 4 scalar USTP Local 396 R(4) 4 scalar 522,528,584,627,630 W3GDATMD Module 373 373 W3SIN4 Subr 292 X Local 412 R(4) 4 scalar 551,554 XFR Local 373 R(4) 4 scalar PTR 373 XI Local 408 R(4) 4 scalar 485,486,487,632,633,634 XJ Local 409 R(4) 4 scalar 636,637,638 XK Local 410 R(4) 4 scalar 641,642,643 XSTRESS Local 413 R(4) 4 scalar 605,625,654 YSTRESS Local 413 R(4) 4 scalar 606,626,655 Z0 Dummy 292 R(4) 4 scalar ARG,IN 427,538,636 Z0NOZ Local 403 R(4) 4 scalar 427,428,429,485 Z0VISC Local 403 R(4) 4 scalar 426,427 ZARG Local 412 R(4) 4 scalar 554,557 ZCN Local 402 R(4) 4 scalar 538,557 ZCO Local 402 R(4) 4 scalar ZLOG Local 412 R(4) 4 scalar 557,559,563 ZZ0RAT Local 376 R(4) 4 scalar PTR 376,427 ZZALP Local 374 R(4) 4 scalar PTR 374,528 ZZWND Local 374 R(4) 4 scalar PTR 374,428 Page 19 Source Listing W3SIN4 2014-09-16 16:52 w3src4md.f90 676 !/ ------------------------------------------------------------------- / 677 SUBROUTINE INSIN4(FLTABS) 678 !/ 679 !/ +-----------------------------------+ 680 !/ | WAVEWATCH III NOAA/NCEP | 681 !/ | SHOM | 682 !/ | F. Ardhuin | 683 !/ | FORTRAN 90 | 684 !/ | Last update : 30-Aug-2010 | 685 !/ +-----------------------------------+ 686 !/ 687 !/ 30-Aug-2010 : Origination. ( version 3.14-Ifremer ) 688 ! 689 ! 1. Purpose : 690 ! 691 ! Initialization for source term routine. 692 ! 693 ! 2. Method : 694 ! 695 ! 3. Parameters : 696 ! 697 ! ---------------------------------------------------------------- 698 ! FLTABS Logical 699 ! ---------------------------------------------------------------- 700 ! 701 ! 4. Subroutines used : 702 ! 703 ! Name Type Module Description 704 ! ---------------------------------------------------------------- 705 ! STRACE Subr. W3SERVMD Subroutine tracing. 706 ! ---------------------------------------------------------------- 707 ! 708 ! 5. Called by : 709 ! 710 ! Name Type Module Description 711 ! ---------------------------------------------------------------- 712 ! W3SIN4 Subr. W3SRC3MD Corresponding source term. 713 ! ---------------------------------------------------------------- 714 ! 715 ! 6. Error messages : 716 ! 717 ! None. 718 ! 719 ! 7. Remarks : 720 ! 721 ! 8. Structure : 722 ! 723 ! See source code. 724 ! 725 ! 9. Switches : 726 ! 727 ! !/S Enable subroutine tracing. 728 ! 729 ! 10. Source code : 730 ! 731 !/ ------------------------------------------------------------------- / 732 USE CONSTANTS, ONLY: TPIINV, RADE, GRAV Page 20 Source Listing INSIN4 2014-09-16 16:52 w3src4md.f90 733 USE W3ODATMD, ONLY: NDSE 734 USE W3SERVMD, ONLY: EXTCDE 735 USE W3DISPMD, ONLY: WAVNU2 736 USE W3GDATMD, ONLY: SIG, DSIP, NK, NTH, TTAUWSHELTER, & 737 SSDSDTH, SSDSCOS, TH, DTH, XFR, ECOS, ESIN, & 738 SSDSC, SSDSBRF1, SSDSBCK, SSDSBINT, SSDSPBK, & 739 SSDSABK, SSDSHCK, IKTAB, DCKI, SATINDICES, & 740 SATWEIGHTS, CUMULW, NKHS, NKD, NDTAB, QBI 741 !/ 742 IMPLICIT NONE 743 !/ 744 !/ ------------------------------------------------------------------- / 745 !/ Parameter list 746 !/ 747 LOGICAL, INTENT(IN) :: FLTABS 748 !/ 749 !/ ------------------------------------------------------------------- / 750 !/ 751 INTEGER SDSNTH, ITH, I_INT, J_INT, IK, IK2, ITH2 , IS, IS2 752 INTEGER IKL, ID, ICON, IKD, IKHS, IKH, TOTO 753 REAL C, C2 754 REAL DIFF1, DIFF2, K_SUP(NK), BINF, BSUP, K(NK), CGG, PROF 755 REAL KIK, DHS, KD, KHS, KH, B, XT, GAM, DKH, PR, W, EPS 756 REAL DKD, DELTAFIT, NHI, H, IH, DH 757 REAL, DIMENSION(:,:) , ALLOCATABLE :: SIGTAB 758 REAL, DIMENSION(:,:) , ALLOCATABLE :: K1, K2 759 !/ 760 !/ ------------------------------------------------------------------- / 761 !/ Local parameters 762 !/ 763 !/ 764 !/ ------------------------------------------------------------------- / 765 !/ 766 ! 767 ! 1. Initializations ------------------------------------------------ * 768 ! 769 ! These precomputed tables are written in mod_def.ww3 770 ! 771 IF (FLTABS) THEN 772 CALL TABU_STRESS 773 CALL TABU_TAUHF(SIG(NK) * TPIINV) !tabulate high-frequency stress 774 IF (TTAUWSHELTER.GT.0) THEN 775 WRITE(NDSE,*) 'Computing 3D lookup table... please wait ...' 776 CALL TABU_TAUHF2(SIG(NK) * TPIINV) !tabulate high-frequency stress 777 END IF 778 END IF 779 ! 780 ! 2. SPONTANEOUS BREAKING 781 ! 2.a Precomputes the indices for integrating the spectrum to get saturation (TEST 4xx ) 782 ! 783 IF (SSDSDTH.LT.180) THEN 784 SDSNTH = MIN(NINT(SSDSDTH/(DTH*RADE)),NTH/2-1) 785 SATINDICES(:,:)=1 786 SATWEIGHTS(:,:)=0. 787 DO ITH=1,NTH 788 DO I_INT=ITH-SDSNTH, ITH+SDSNTH 789 J_INT=I_INT Page 21 Source Listing INSIN4 2014-09-16 16:52 w3src4md.f90 790 IF (I_INT.LT.1) J_INT=I_INT+NTH 791 IF (I_INT.GT.NTH) J_INT=I_INT-NTH 792 SATINDICES(I_INT-(ITH-SDSNTH)+1,ITH)=J_INT 793 SATWEIGHTS(I_INT-(ITH-SDSNTH)+1,ITH)= & 794 COS(TH(ITH)-TH(J_INT))**SSDSCOS 795 END DO 796 END DO 797 ELSE 798 SATINDICES(:,:)=1 799 SATWEIGHTS(:,:)=1. 800 END IF 801 !/ ------------------------------------------------------------------- / 802 ! 803 ! Precomputes QBI and DCKI (TEST 500) 804 ! 805 IF (SSDSBCK.GT.0) THEN 806 ! 807 ! Precomputes the indices for integrating the spectrum over frquency bandwidth 808 ! 809 BINF=(1-SSDSBINT) ! Banner et al 2002: Hp=4*sqrt(int_0.7^1.3fp E df), SSDSBINT=0.3 810 BSUP=(1+SSDSBINT) 811 KIK=0. 812 ! 813 ! High frequency tail for convolution calculation 814 ! 815 ALLOCATE(K1(NK,NDTAB)) 816 ALLOCATE(K2(NK,NDTAB)) 817 ALLOCATE(SIGTAB(NK,NDTAB)) 818 819 SIGTAB=0. !contains frequency for upper windows boundaries 820 IKTAB=0 ! contains indices for upper windows boundaries 821 822 DO ID=1,NDTAB 823 TOTO=0 824 PROF=REAL(ID) 825 DO IKL=1,NK ! last window starts at IK=NK 826 CALL WAVNU2(SIG(IKL), PROF, KIK, CGG, 1E-7, 15, ICON) 827 K1(IKL,ID)=KIK ! wavenumber lower boundary (is directly related to the frequency indices, IK) 828 K2(IKL,ID)=((BSUP/BINF)**2.)*K1(IKL,ID)! wavenumber upper boundary 829 SIGTAB(IKL,ID)=SQRT(GRAV*K2(IKL,ID)*TANH(K2(IKL,ID)*ID)) ! corresponding frequency upper boundary 830 IF(SIGTAB(IKL,ID) .LE. SIG(1)) THEN 831 IKTAB(IKL,ID)=1 832 END IF 833 IF(SIGTAB(IKL,ID) .GT. SIG(NK)) THEN 834 IKTAB(IKL,ID)=NK+TOTO ! in w3sds4 only windows with IKSUP<=NK will be kept 835 TOTO=1 836 END IF 837 DO IK=1,NK-1 838 DIFF1=0. 839 DIFF2=0. 840 IF(SIG(IK)=SIGTAB(IKL,ID)) THEN 841 DIFF1=SIGTAB(IKL,ID)-SIG(IK) ! seeks the indices of the upper boundary 842 DIFF2=SIG(IK+1)-SIGTAB(IKL,ID)! the indices of lower boudary = IK 843 IF (DIFF1TAUW. 1009 ! ---------------------------------------------------------------------- 1010 INTEGER I,J,ITER 1011 REAL ZTAUW,UTOP,CDRAG,WCD,USTOLD,TAUOLD 1012 REAL X,UST,ZZ0,ZNU,F,DELF,ZZ00 1013 ! 1014 DELU = UMAX/FLOAT(JUMAX) 1015 DELTAUW = TAUWMAX/FLOAT(ITAUMAX) 1016 DO I=0,ITAUMAX 1017 ZTAUW = (REAL(I)*DELTAUW)**2 1018 DO J=0,JUMAX 1019 UTOP = FLOAT(J)*DELU 1020 CDRAG = 0.0012875 1021 WCD = SQRT(CDRAG) 1022 USTOLD = UTOP*WCD 1023 TAUOLD = MAX(USTOLD**2, ZTAUW+EPS1) 1024 DO ITER=1,NITER 1025 X = ZTAUW/TAUOLD 1026 UST = SQRT(TAUOLD) 1027 ZZ00=AALPHA*TAUOLD/GRAV 1028 IF (ZZ0MAX.NE.0) ZZ00=MIN(ZZ00,ZZ0MAX) 1029 ! Corrects roughness ZZ00 for quasi-linear effect 1030 ZZ0 = ZZ00/(1.-X)**XM 1031 !ZNU = 0.1*nu_air/UST ! This was removed by Bidlot in 1996 1032 !ZZ0 = MAX(ZNU,ZZ0) 1033 F = UST-KAPPA*UTOP/(ALOG(ZZWND/ZZ0)) 1034 DELF= 1.-KAPPA*UTOP/(ALOG(ZZWND/ZZ0))**2*2./UST & 1035 *(1.-(XM+1)*X)/(1.-X) 1036 UST = UST-F/DELF 1037 TAUOLD= MAX(UST**2., ZTAUW+EPS1) 1038 END DO 1039 TAUT(I,J) = SQRT(TAUOLD) 1040 END DO 1041 END DO 1042 I=ITAUMAX 1043 J=JUMAX Page 29 Source Listing TABU_STRESS 2014-09-16 16:52 w3src4md.f90 1044 ! 1045 ! Force zero wind to have zero stress (Bidlot 1996) 1046 ! 1047 DO I=0,ITAUMAX 1048 TAUT(I,0)=0.0 1049 END DO 1050 RETURN 1051 END SUBROUTINE TABU_STRESS ENTRY POINTS Name w3src4md_mp_tabu_stress_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References AALPHA Local 999 R(4) 4 scalar PTR 999,1027 ALOG Func 1033 scalar 1033,1034 CDRAG Local 1011 R(4) 4 scalar 1020,1021 CONSTANTS Module 998 998 DELF Local 1012 R(4) 4 scalar 1034,1036 DELTAUW Local 1015 R(4) 4 scalar 1015,1017,1416 DELU Local 1014 R(4) 4 scalar 1014,1019,1420 EPS1 Param 1002 R(4) 4 scalar 1023,1037 F Local 1012 R(4) 4 scalar 1033,1036 FLOAT Func 1014 scalar 1014,1015,1019 GRAV Param 1027 R(4) 4 scalar 1027 I Local 1010 I(4) 4 scalar 1016,1017,1039,1042,1047,1048 ITAUMAX Param 1015 I(4) 4 scalar 68,1015,1016,1042,1047,1417 ITER Local 1010 I(4) 4 scalar 1024 J Local 1010 I(4) 4 scalar 1018,1019,1039,1043 JUMAX Param 1014 I(4) 4 scalar 68,1014,1018,1043,1421 KAPPA Param 1033 R(4) 4 scalar 1033,1034 MAX Func 1023 scalar 1023,1037 MIN Func 1028 scalar 1028 NITER Param 1001 I(4) 4 scalar 1024 REAL Func 1017 scalar 1017 SQRT Func 1021 scalar 1021,1026,1039 TABU_STRESS Subr 931 772 TAUOLD Local 1011 R(4) 4 scalar 1023,1025,1026,1027,1037,1039 TAUT Local 1039 R(4) 4 2 40401 1039,1048,1424,1425 TAUWMAX Param 1015 R(4) 4 scalar 1015,1415 UMAX Param 1014 R(4) 4 scalar 1014 UST Local 1012 R(4) 4 scalar 1026,1033,1034,1036,1037 USTOLD Local 1011 R(4) 4 scalar 1022,1023 UTOP Local 1011 R(4) 4 scalar 1019,1022,1033,1034 W3GDATMD Module 999 999 WCD Local 1011 R(4) 4 scalar 1021,1022 X Local 1012 R(4) 4 scalar 1025,1030,1035 XM Param 1002 R(4) 4 scalar 1030,1035 ZNU Local 1012 R(4) 4 scalar ZTAUW Local 1011 R(4) 4 scalar 1017,1023,1025,1037 Page 30 Source Listing TABU_STRESS 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References ZZ0 Local 1012 R(4) 4 scalar 1030,1033,1034 ZZ00 Local 1012 R(4) 4 scalar 1027,1028,1030 ZZ0MAX Local 999 R(4) 4 scalar PTR 999,1028 ZZWND Local 999 R(4) 4 scalar PTR 999,1033,1034 Page 31 Source Listing TABU_STRESS 2014-09-16 16:52 w3src4md.f90 1052 !/ ------------------------------------------------------------------- / 1053 SUBROUTINE TABU_TAUHF(FRMAX) 1054 !/ 1055 !/ +-----------------------------------+ 1056 !/ | WAVEWATCH III NOAA/NCEP | 1057 !/ | F. Ardhuin | 1058 !/ | FORTRAN 90 | 1059 !/ | Last update 2006/08/14 | 1060 !/ +-----------------------------------+ 1061 !/ 1062 !/ 27-Feb-2004 : Origination in WW3 ( version 2.22.SHOM ) 1063 !/ the resulting table was checked to be identical to the original f77 result 1064 !/ 14-Aug-2006 : Modified following Bidlot ( version 2.22.SHOM ) 1065 !/ 18-Aug-2006 : Ported to version 3.09 1066 ! 1067 ! 1. Purpose : 1068 ! 1069 ! Tabulation of the high-frequency wave-supported stress 1070 ! 1071 ! 2. Method : 1072 ! 1073 ! SEE REFERENCE FOR WAVE STRESS CALCULATION. 1074 ! FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990. 1075 ! See tech. Memo ECMWF 03 december 2003 by Bidlot & Janssen 1076 ! 1077 ! 3. Parameters : 1078 ! 1079 ! Parameter list 1080 ! ---------------------------------------------------------------- 1081 ! FRMAX Real I maximum frequency. 1082 ! ---------------------------------------------------------------- 1083 ! 1084 ! 4. Subroutines used : 1085 ! 1086 ! STRACE Service routine. 1087 ! 1088 ! 5. Called by : 1089 ! 1090 ! W3SIN3 Wind input Source term routine. 1091 ! 1092 ! 6. Error messages : 1093 ! 1094 ! 7. Remarks : 1095 ! 1096 ! 8. Structure : 1097 ! 1098 ! See source code. 1099 ! 1100 ! 9. Switches : 1101 ! 1102 ! !/S Enable subroutine tracing. 1103 ! !/T Enable test output. 1104 ! 1105 ! 10. Source code : 1106 ! 1107 !/ ------------------------------------------------------------------- / 1108 USE CONSTANTS Page 32 Source Listing TABU_TAUHF 2014-09-16 16:52 w3src4md.f90 1109 USE W3GDATMD, ONLY: AALPHA, BBETA, ZZALP, XFR, FACHFE, ZZ0MAX 1110 ! 1111 IMPLICIT NONE 1112 !/ 1113 !/ ------------------------------------------------------------------- / 1114 !/ Parameter list 1115 !/ 1116 REAL, intent(in) :: FRMAX ! maximum frequency 1117 !/ 1118 !/ ------------------------------------------------------------------- / 1119 !/ Local parameters 1120 !/ 1121 ! USTARM R.A. Maximum friction velocity 1122 ! ALPHAM R.A. Maximum Charnock Coefficient 1123 ! WLV R.A. Water levels. 1124 ! UA R.A. Absolute wind speeds. 1125 ! UD R.A. Absolute wind direction. 1126 ! U10 R.A. Wind speed used. 1127 ! U10D R.A. Wind direction used. 1128 ! 10. Source code : 1129 ! 1130 !/ ------------------------------------------------------------------- / 1131 REAL :: USTARM, ALPHAM 1132 REAL :: CONST1, OMEGA, OMEGAC 1133 REAL :: UST, ZZ0,OMEGACC, CM 1134 INTEGER, PARAMETER :: JTOT=250 1135 REAL, ALLOCATABLE :: W(:) 1136 REAL :: ZX,ZARG,ZMU,ZLOG,ZZ00,ZBETA 1137 REAL :: Y,YC,DELY 1138 INTEGER :: I,J,K,L 1139 REAL :: X0 1140 ! 1141 USTARM = 5. 1142 ALPHAM = 20.*AALPHA 1143 DELUST = USTARM/REAL(IUSTAR) 1144 DELALP = ALPHAM/REAL(IALPHA) 1145 CONST1 = BBETA/KAPPA**2 1146 OMEGAC = TPI*FRMAX 1147 ! 1148 TAUHFT(0:IUSTAR,0:IALPHA)=0. !table initialization 1149 ! 1150 ALLOCATE(W(JTOT)) 1151 W(2:JTOT-1)=1. 1152 W(1)=0.5 1153 W(JTOT)=0.5 1154 X0 = 0.05 1155 ! 1156 DO L=0,IALPHA 1157 DO K=0,IUSTAR 1158 UST = MAX(REAL(K)*DELUST,0.000001) 1159 ZZ00 = UST**2*AALPHA/GRAV 1160 IF (ZZ0MAX.NE.0) ZZ00=MIN(ZZ00,ZZ0MAX) 1161 ZZ0 = ZZ00*(1+FLOAT(L)*DELALP/AALPHA) 1162 OMEGACC = MAX(OMEGAC,X0*GRAV/UST) 1163 YC = OMEGACC*SQRT(ZZ0/GRAV) 1164 DELY = MAX((1.-YC)/REAL(JTOT),0.) 1165 ! For a given value of UST and ALPHA, Page 33 Source Listing TABU_TAUHF 2014-09-16 16:52 w3src4md.f90 1166 ! the wave-supported stress is integrated all the way 1167 ! to 0.05*g/UST 1168 DO J=1,JTOT 1169 Y = YC+REAL(J-1)*DELY 1170 OMEGA = Y*SQRT(GRAV/ZZ0) 1171 ! This is the deep water phase speed 1172 CM = GRAV/OMEGA 1173 !this is the inverse wave age, shifted by ZZALP (tuning) 1174 ZX = UST/CM +ZZALP 1175 ZARG = MIN(KAPPA/ZX,20.) 1176 ZMU = MIN(GRAV*ZZ0/CM**2*EXP(ZARG),1.) 1177 ZLOG = MIN(ALOG(ZMU),0.) 1178 ZBETA = CONST1*ZMU*ZLOG**4 1179 ! Power of Y in denominator should be FACHFE-4 1180 TAUHFT(K,L) = TAUHFT(K,L)+W(J)*ZBETA/Y*DELY 1181 END DO 1182 END DO 1183 END DO 1184 DEALLOCATE(W) 1185 RETURN 1186 END SUBROUTINE TABU_TAUHF ENTRY POINTS Name w3src4md_mp_tabu_tauhf_ Page 34 Source Listing TABU_TAUHF 2014-09-16 16:52 Symbol Table w3src4md.f90 SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References AALPHA Local 1109 R(4) 4 scalar PTR 1109,1142,1159,1161 ALOG Func 1177 scalar 1177 ALPHAM Local 1131 R(4) 4 scalar 1142,1144 BBETA Local 1109 R(4) 4 scalar PTR 1109,1145 CM Local 1133 R(4) 4 scalar 1172,1174,1176 CONST1 Local 1132 R(4) 4 scalar 1145,1178 CONSTANTS Module 1108 1108 DELY Local 1137 R(4) 4 scalar 1164,1169,1180 EXP Func 1176 scalar 1176 FACHFE Local 1109 R(4) 4 scalar PTR 1109 FLOAT Func 1161 scalar 1161 FRMAX Dummy 1053 R(4) 4 scalar ARG,IN 1146 GRAV Param 1159 R(4) 4 scalar 1159,1162,1163,1170,1172,1176 I Local 1138 I(4) 4 scalar J Local 1138 I(4) 4 scalar 1168,1169,1180 JTOT Param 1134 I(4) 4 scalar 1150,1151,1153,1164,1168 K Local 1138 I(4) 4 scalar 1157,1158,1180 KAPPA Param 1145 R(4) 4 scalar 1145,1175 L Local 1138 I(4) 4 scalar 1156,1161,1180 MAX Func 1158 scalar 1158,1162,1164 MIN Func 1160 scalar 1160,1175,1176,1177 OMEGA Local 1132 R(4) 4 scalar 1170,1172 OMEGAC Local 1132 R(4) 4 scalar 1146,1162 OMEGACC Local 1133 R(4) 4 scalar 1162,1163 REAL Func 1143 scalar 1143,1144,1158,1164,1169 SQRT Func 1163 scalar 1163,1170 TABU_TAUHF Subr 1053 773 TPI Param 1146 R(4) 4 scalar 1146 UST Local 1133 R(4) 4 scalar 1158,1159,1162,1174 USTARM Local 1131 R(4) 4 scalar 1141,1143 W Local 1135 R(4) 4 1 1 ALC 1150,1151,1152,1153,1180,1184 W3GDATMD Module 1109 1109 X0 Local 1139 R(4) 4 scalar 1154,1162 XFR Local 1109 R(4) 4 scalar PTR 1109 Y Local 1137 R(4) 4 scalar 1169,1170,1180 YC Local 1137 R(4) 4 scalar 1163,1164,1169 ZARG Local 1136 R(4) 4 scalar 1175,1176 ZBETA Local 1136 R(4) 4 scalar 1178,1180 ZLOG Local 1136 R(4) 4 scalar 1177,1178 ZMU Local 1136 R(4) 4 scalar 1176,1177,1178 ZX Local 1136 R(4) 4 scalar 1174,1175 ZZ0 Local 1133 R(4) 4 scalar 1161,1163,1170,1176 ZZ00 Local 1136 R(4) 4 scalar 1159,1160,1161 ZZ0MAX Local 1109 R(4) 4 scalar PTR 1109,1160 ZZALP Local 1109 R(4) 4 scalar PTR 1109,1174 Page 35 Source Listing TABU_TAUHF 2014-09-16 16:52 w3src4md.f90 1187 1188 !/ ------------------------------------------------------------------- / 1189 SUBROUTINE TABU_TAUHF2(FRMAX) 1190 !/ 1191 !/ +-----------------------------------+ 1192 !/ | WAVEWATCH III NOAA/NCEP | 1193 !/ | F. Ardhuin | 1194 !/ | FORTRAN 90 | 1195 !/ | Last update 2006/08/14 | 1196 !/ | Last update 2013/01/24 | 1197 !/ +-----------------------------------+ 1198 !/ 1199 !/ 15-May-2007 : Origination in WW3 ( version 3.10.SHOM ) 1200 !/ 24-Jan-2013 : Allows to read in table ( version 4.08 ) 1201 ! 1202 ! 1. Purpose : 1203 ! 1204 ! Tabulation of the high-frequency wave-supported stress as a function of 1205 ! ustar, alpha (modified Charnock), and tail energy level 1206 ! 1207 ! 2. Method : 1208 ! 1209 ! SEE REFERENCE FOR WAVE STRESS CALCULATION. 1210 ! FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990. 1211 ! See tech. Memo ECMWF 03 december 2003 by Bidlot & Janssen 1212 ! 1213 ! 3. Parameters : 1214 ! 1215 ! Parameter list 1216 ! ---------------------------------------------------------------- 1217 ! FRMAX Real I maximum frequency. 1218 ! ---------------------------------------------------------------- 1219 ! 1220 ! 4. Subroutines used : 1221 ! 1222 ! STRACE Service routine. 1223 ! 1224 ! 5. Called by : 1225 ! 1226 ! W3SIN3 Wind input Source term routine. 1227 ! 1228 ! 6. Error messages : 1229 ! 1230 ! 7. Remarks : 1231 ! 1232 ! 8. Structure : 1233 ! 1234 ! See source code. 1235 ! 1236 ! 9. Switches : 1237 ! 1238 ! !/S Enable subroutine tracing. 1239 ! !/T Enable test output. 1240 ! 1241 ! 10. Source code : 1242 ! 1243 !/ ------------------------------------------------------------------- / Page 36 Source Listing TABU_TAUHF2 2014-09-16 16:52 w3src4md.f90 1244 USE CONSTANTS 1245 USE W3GDATMD, ONLY: AALPHA, BBETA, ZZALP, XFR, FACHFE, & 1246 TTAUWSHELTER, ZZ0MAX 1247 ! 1248 IMPLICIT NONE 1249 !/ 1250 !/ ------------------------------------------------------------------- / 1251 !/ Parameter list 1252 !/ 1253 REAL, intent(in) :: FRMAX ! maximum frequency 1254 !/ 1255 !/ ------------------------------------------------------------------- / 1256 !/ Local parameters 1257 !/ 1258 ! USTARM R.A. Maximum friction velocity 1259 ! ALPHAM R.A. Maximum Charnock Coefficient 1260 ! WLV R.A. Water levels. 1261 ! UA R.A. Absolute wind speeds. 1262 ! UD R.A. Absolute wind direction. 1263 ! U10 R.A. Wind speed used. 1264 ! U10D R.A. Wind direction used. 1265 ! 10. Source code : 1266 ! 1267 !/ ------------------------------------------------------------------- / 1268 REAL :: USTARM, ALPHAM, LEVTAILM 1269 REAL :: CONST1, OMEGA, OMEGAC, LEVTAIL 1270 REAL :: UST, UST0, ZZ0,OMEGACC, CM 1271 REAL :: TAUW, TAUW0 1272 INTEGER, PARAMETER :: JTOT=250 1273 REAL, ALLOCATABLE :: W(:) 1274 REAL :: ZX,ZARG,ZMU,ZLOG,ZBETA 1275 REAL :: Y,YC,DELY 1276 INTEGER :: I, J, K, L 1277 REAL :: X0 1278 ! 1279 USTARM = 5. 1280 ALPHAM = 20.*AALPHA 1281 LEVTAILM = 0.05 1282 DELUST = USTARM/REAL(IUSTAR) 1283 DELALP = ALPHAM/REAL(IALPHA) 1284 DELTAIL = ALPHAM/REAL(ILEVTAIL) 1285 CONST1 = BBETA/KAPPA**2 1286 OMEGAC = TPI*FRMAX 1287 ! 1288 TAUHFT(0:IUSTAR,0:IALPHA)=0. !table initialization 1289 ! 1290 ALLOCATE(W(JTOT)) 1291 W(2:JTOT-1)=1. 1292 W(1)=0.5 1293 W(JTOT)=0.5 1294 X0 = 0.05 1295 ! 1296 DO K=0,IUSTAR 1297 UST0 = MAX(REAL(K)*DELUST,0.000001) 1298 DO L=0,IALPHA 1299 UST=UST0 1300 ZZ0 = UST0**2*(AALPHA+FLOAT(L)*DELALP)/GRAV Page 37 Source Listing TABU_TAUHF2 2014-09-16 16:52 w3src4md.f90 1301 OMEGACC = MAX(OMEGAC,X0*GRAV/UST) 1302 YC = OMEGACC*SQRT(ZZ0/GRAV) 1303 DELY = MAX((1.-YC)/REAL(JTOT),0.) 1304 ! For a given value of UST and ALPHA, 1305 ! the wave-supported stress is integrated all the way 1306 ! to 0.05*g/UST 1307 DO I=0,ILEVTAIL 1308 LEVTAIL=REAL(I)*DELTAIL 1309 TAUHFT(K,L)=0. 1310 TAUHFT2(K,L,I)=0. 1311 TAUW0=UST0**2 1312 TAUW=TAUW0 1313 DO J=1,JTOT 1314 Y = YC+REAL(J-1)*DELY 1315 OMEGA = Y*SQRT(GRAV/ZZ0) 1316 ! This is the deep water phase speed 1317 CM = GRAV/OMEGA 1318 !this is the inverse wave age, shifted by ZZALP (tuning) 1319 ZX = UST0/CM +ZZALP 1320 ZARG = MIN(KAPPA/ZX,20.) 1321 ZMU = MIN(GRAV*ZZ0/CM**2*EXP(ZARG),1.) 1322 ZLOG = MIN(ALOG(ZMU),0.) 1323 ZBETA = CONST1*ZMU*ZLOG**4 1324 ! Power of Y in denominator should be FACHFE-4 1325 TAUHFT(K,L) = TAUHFT(K,L)+W(J)*ZBETA/Y*DELY 1326 ZX = UST/CM +ZZALP 1327 ZARG = MIN(KAPPA/ZX,20.) 1328 ZMU = MIN(GRAV*ZZ0/CM**2*EXP(ZARG),1.) 1329 ZLOG = MIN(ALOG(ZMU),0.) 1330 ZBETA = CONST1*ZMU*ZLOG**4 1331 ! Power of Y in denominator should be FACHFE-4 1332 TAUHFT2(K,L,I) = TAUHFT2(K,L,I)+W(J)*ZBETA*(UST/UST0)**2/Y*DELY 1333 TAUW=TAUW-W(J)*UST**2*ZBETA*LEVTAIL/Y*DELY 1334 UST=SQRT(MAX(TAUW,0.)) 1335 END DO 1336 END DO 1337 END DO 1338 END DO 1339 DEALLOCATE(W) 1340 RETURN 1341 END SUBROUTINE TABU_TAUHF2 Page 38 Source Listing TABU_TAUHF2 2014-09-16 16:52 Entry Points w3src4md.f90 ENTRY POINTS Name w3src4md_mp_tabu_tauhf2_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References AALPHA Local 1245 R(4) 4 scalar PTR 1245,1280,1300 ALOG Func 1322 scalar 1322,1329 ALPHAM Local 1268 R(4) 4 scalar 1280,1283,1284 BBETA Local 1245 R(4) 4 scalar PTR 1245,1285 CM Local 1270 R(4) 4 scalar 1317,1319,1321,1326,1328 CONST1 Local 1269 R(4) 4 scalar 1285,1323,1330 CONSTANTS Module 1244 1244 DELY Local 1275 R(4) 4 scalar 1303,1314,1325,1332,1333 EXP Func 1321 scalar 1321,1328 FACHFE Local 1245 R(4) 4 scalar PTR 1245 FLOAT Func 1300 scalar 1300 FRMAX Dummy 1189 R(4) 4 scalar ARG,IN 1286 GRAV Param 1300 R(4) 4 scalar 1300,1301,1302,1315,1317,1321,1328 I Local 1276 I(4) 4 scalar 1307,1308,1310,1332 J Local 1276 I(4) 4 scalar 1313,1314,1325,1332,1333 JTOT Param 1272 I(4) 4 scalar 1290,1291,1293,1303,1313 K Local 1276 I(4) 4 scalar 1296,1297,1309,1310,1325,1332 KAPPA Param 1285 R(4) 4 scalar 1285,1320,1327 L Local 1276 I(4) 4 scalar 1298,1300,1309,1310,1325,1332 LEVTAIL Local 1269 R(4) 4 scalar 1308,1333 LEVTAILM Local 1268 R(4) 4 scalar 1281 MAX Func 1297 scalar 1297,1301,1303,1334 MIN Func 1320 scalar 1320,1321,1322,1327,1328,1329 OMEGA Local 1269 R(4) 4 scalar 1315,1317 OMEGAC Local 1269 R(4) 4 scalar 1286,1301 OMEGACC Local 1270 R(4) 4 scalar 1301,1302 REAL Func 1282 scalar 1282,1283,1284,1297,1303,1308,1314 SQRT Func 1302 scalar 1302,1315,1334 TABU_TAUHF2 Subr 1189 776 TAUW Local 1271 R(4) 4 scalar 1312,1333,1334 TAUW0 Local 1271 R(4) 4 scalar 1311,1312 TPI Param 1286 R(4) 4 scalar 1286 TTAUWSHELTER Local 1246 R(4) 4 scalar PTR 1246 UST Local 1270 R(4) 4 scalar 1299,1301,1326,1332,1333,1334 UST0 Local 1270 R(4) 4 scalar 1297,1299,1300,1311,1319,1332 USTARM Local 1268 R(4) 4 scalar 1279,1282 W Local 1273 R(4) 4 1 1 ALC 1290,1291,1292,1293,1325,1332,1333 ,1339 W3GDATMD Module 1245 1245 X0 Local 1277 R(4) 4 scalar 1294,1301 XFR Local 1245 R(4) 4 scalar PTR 1245 Y Local 1275 R(4) 4 scalar 1314,1315,1325,1332,1333 YC Local 1275 R(4) 4 scalar 1302,1303,1314 ZARG Local 1274 R(4) 4 scalar 1320,1321,1327,1328 Page 39 Source Listing TABU_TAUHF2 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References ZBETA Local 1274 R(4) 4 scalar 1323,1325,1330,1332,1333 ZLOG Local 1274 R(4) 4 scalar 1322,1323,1329,1330 ZMU Local 1274 R(4) 4 scalar 1321,1322,1323,1328,1329,1330 ZX Local 1274 R(4) 4 scalar 1319,1320,1326,1327 ZZ0 Local 1270 R(4) 4 scalar 1300,1302,1315,1321,1328 ZZ0MAX Local 1246 R(4) 4 scalar PTR 1246 ZZALP Local 1245 R(4) 4 scalar PTR 1245,1319,1326 Page 40 Source Listing TABU_TAUHF2 2014-09-16 16:52 w3src4md.f90 1342 1343 !/ ------------------------------------------------------------------- / 1344 SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN) 1345 !/ 1346 !/ +-----------------------------------+ 1347 !/ | WAVEWATCH III NOAA/NCEP | 1348 !/ | F. Ardhuin | 1349 !/ | FORTRAN 90 | 1350 !/ | Last update 2006/08/14 | 1351 !/ +-----------------------------------+ 1352 !/ 1353 !/ 27-Feb-2004 : Origination in WW3 ( version 2.22-SHOM ) 1354 !/ the resulting table was checked to be identical to the original f77 result 1355 !/ 14-Aug-2006 : Modified following Bidlot ( version 2.22-SHOM ) 1356 !/ 18-Aug-2006 : Ported to version 3.09 1357 !/ 03-Apr-2010 : Adding output of Charnock parameter ( version 3.14-IFREMER ) 1358 ! 1359 ! 1. Purpose : 1360 ! 1361 ! Compute friction velocity based on wind speed U10 1362 ! 1363 ! 2. Method : 1364 ! 1365 ! Computation of u* based on Quasi-linear theory 1366 ! 1367 ! 3. Parameters : 1368 ! 1369 ! Parameter list 1370 ! ---------------------------------------------------------------- 1371 ! U10,TAUW,USTAR,Z0 1372 ! ---------------------------------------------------------------- 1373 ! WINDSPEED Real I 10-m wind speed ... should be NEUTRAL 1374 ! TAUW Real I Wave-supported stress 1375 ! USTAR Real O Friction velocity. 1376 ! Z0 Real O air-side roughness length 1377 ! ---------------------------------------------------------------- 1378 ! 1379 ! 4. Subroutines used : 1380 ! 1381 ! STRACE Service routine. 1382 ! 1383 ! 5. Called by : 1384 ! 1385 ! W3SIN3 Wind input Source term routine. 1386 ! 1387 ! 6. Error messages : 1388 ! 1389 ! 7. Remarks : 1390 ! 1391 ! 8. Structure : 1392 ! 1393 ! See source code. 1394 ! 1395 ! 9. Switches : 1396 ! 1397 ! !/S Enable subroutine tracing. 1398 ! !/T Enable test output. Page 41 Source Listing CALC_USTAR 2014-09-16 16:52 w3src4md.f90 1399 ! 1400 ! 10. Source code : 1401 !-----------------------------------------------------------------------------! 1402 USE CONSTANTS, ONLY: GRAV, KAPPA 1403 USE W3GDATMD, ONLY: ZZWND, AALPHA 1404 IMPLICIT NONE 1405 REAL, intent(in) :: WINDSPEED,TAUW 1406 REAL, intent(out) :: USTAR, Z0, CHARN 1407 ! local variables 1408 real a,b ! constants of parameterisation 1409 real cd ! drag coefficient 1410 REAL SQRTCDM1 1411 REAL X,XI,DELI1,DELI2,XJ,delj1,delj2 1412 REAL UST,DELTOLD,TAUW_LOCAL 1413 INTEGER IND,J 1414 ! 1415 TAUW_LOCAL=MAX(MIN(TAUW,TAUWMAX),0.) 1416 XI = SQRT(TAUW_LOCAL)/DELTAUW 1417 IND = MIN ( ITAUMAX-1, INT(XI)) ! index for stress table 1418 DELI1 = MIN(1.,XI - REAL(IND)) !interpolation coefficient for stress table 1419 DELI2 = 1. - DELI1 1420 XJ = WINDSPEED/DELU 1421 J = MIN ( JUMAX-1, INT(XJ) ) 1422 DELJ1 = MIN(1.,XJ - REAL(J)) 1423 DELJ2 = 1. - DELJ1 1424 USTAR=(TAUT(IND,J)*DELI2+TAUT(IND+1,J )*DELI1)*DELJ2 & 1425 + (TAUT(IND,J+1)*DELI2+TAUT(IND+1,J+1)*DELI1)*DELJ1 1426 ! 1427 ! Determines roughness length 1428 ! 1429 SQRTCDM1 = MIN(WINDSPEED/USTAR,100.0) 1430 Z0 = ZZWND*EXP(-KAPPA*SQRTCDM1) 1431 IF (USTAR.GT.0.001) THEN 1432 CHARN = GRAV*Z0/USTAR**2 1433 ELSE 1434 CHARN = AALPHA 1435 END IF 1436 ! 1437 RETURN 1438 END SUBROUTINE CALC_USTAR Page 42 Source Listing CALC_USTAR 2014-09-16 16:52 Entry Points w3src4md.f90 ENTRY POINTS Name w3src4md_mp_calc_ustar_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References A Local 1408 R(4) 4 scalar AALPHA Local 1403 R(4) 4 scalar PTR 1403,1434 B Local 1408 R(4) 4 scalar CALC_USTAR Subr 1344 276 CD Local 1409 R(4) 4 scalar CHARN Dummy 1344 R(4) 4 scalar ARG,OUT 1432,1434 CONSTANTS Module 1402 1402 DELI1 Local 1411 R(4) 4 scalar 1418,1419,1424,1425 DELI2 Local 1411 R(4) 4 scalar 1419,1424,1425 DELJ1 Local 1411 R(4) 4 scalar 1422,1423,1425 DELJ2 Local 1411 R(4) 4 scalar 1423,1424 DELTOLD Local 1412 R(4) 4 scalar EXP Func 1430 scalar 1430 GRAV Param 1402 R(4) 4 scalar 1402,1432 IND Local 1413 I(4) 4 scalar 1417,1418,1424,1425 INT Func 1417 scalar 1417,1421 J Local 1413 I(4) 4 scalar 1421,1422,1424,1425 KAPPA Param 1402 R(4) 4 scalar 1402,1430 MAX Func 1415 scalar 1415 MIN Func 1415 scalar 1415,1417,1418,1421,1422,1429 REAL Func 1418 scalar 1418,1422 SQRT Func 1416 scalar 1416 SQRTCDM1 Local 1410 R(4) 4 scalar 1429,1430 TAUW Dummy 1344 R(4) 4 scalar ARG,IN 1415 TAUW_LOCAL Local 1412 R(4) 4 scalar 1415,1416 UST Local 1412 R(4) 4 scalar USTAR Dummy 1344 R(4) 4 scalar ARG,OUT 1424,1429,1431,1432 W3GDATMD Module 1403 1403 WINDSPEED Dummy 1344 R(4) 4 scalar ARG,IN 1420,1429 X Local 1411 R(4) 4 scalar XI Local 1411 R(4) 4 scalar 1416,1417,1418 XJ Local 1411 R(4) 4 scalar 1420,1421,1422 Z0 Dummy 1344 R(4) 4 scalar ARG,OUT 1430,1432 ZZWND Local 1403 R(4) 4 scalar PTR 1403,1430 Page 43 Source Listing CALC_USTAR 2014-09-16 16:52 w3src4md.f90 1439 !/ ------------------------------------------------------------------- / 1440 SUBROUTINE W3SDS4 (A, K, CG, USTAR, USDIR, DEPTH, S, & 1441 1442 D, IX, IY, WHITECAP ) 1443 !/ 1444 !/ +-----------------------------------+ 1445 !/ | WAVEWATCH III NOAA/NCEP | 1446 !/ ! F. Ardhuin ! 1447 !/ | FORTRAN 90 | 1448 !/ | Last update : 30-Aug-2010 | 1449 !/ +-----------------------------------+ 1450 !/ 1451 !/ 30-Aug-2010 : Clean up from common ST3-ST4 routine( version 3.14-Ifremer ) 1452 !/ 1453 ! 1. Purpose : 1454 ! 1455 ! Calculate whitecapping source term and diagonal term of derivative. 1456 ! 1457 ! 2. Method : 1458 ! 1459 ! This codes does either one or the other of 1460 ! Ardhuin et al. (JPO 2010) 1461 ! Filipot & Ardhuin (JGR 2012) 1462 ! 1463 ! 3. Parameters : 1464 ! 1465 ! Parameter list 1466 ! ---------------------------------------------------------------- 1467 ! A R.A. I Action density spectrum (1-D). 1468 ! K R.A. I Wavenumber for entire spectrum. *) 1469 ! USTAR Real I Friction velocity. 1470 ! USDIR Real I wind stress direction. 1471 ! DEPTH Real I Water depth. 1472 ! S R.A. O Source term (1-D version). 1473 ! D R.A. O Diagonal term of derivative. *) 1474 ! ---------------------------------------------------------------- 1475 ! *) Stored in 1-D array with dimension NTH*NK 1476 ! 1477 ! 4. Subroutines used : 1478 ! 1479 ! STRACE Subroutine tracing. ( !/S switch ) 1480 ! PRT2DS Print plot of spectrum. ( !/T0 switch ) 1481 ! OUTMAT Print out matrix. ( !/T1 switch ) 1482 ! 1483 ! 5. Called by : 1484 ! 1485 ! W3SRCE Source term integration. 1486 ! W3EXPO Point output program. 1487 ! GXEXPO GrADS point output program. 1488 ! 1489 ! 6. Error messages : 1490 ! 1491 ! 7. Remarks : 1492 ! 1493 ! 8. Structure : 1494 ! 1495 ! See source code. Page 44 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1496 ! 1497 ! 9. Switches : 1498 ! 1499 ! !/S Enable subroutine tracing. 1500 ! !/T Enable general test output. 1501 ! !/T0 2-D print plot of source term. 1502 ! !/T1 Print arrays. 1503 ! 1504 ! 10. Source code : 1505 ! 1506 !/ ------------------------------------------------------------------- / 1507 USE CONSTANTS 1508 USE W3GDATMD, ONLY: NSPEC, NTH, NK, SSDSBR, DDEN, & 1509 SSDSC, & 1510 SIG, SSDSP, ECOS, ESIN, DTH, DSIP, & 1511 SSDSISO, SSDSDTH, SSDSBR2, SSDSBM, & 1512 SSDSBRFDF, SSDSBCK, SSDSBINT, IKTAB, DCKI, & 1513 SATINDICES, SATWEIGHTS, CUMULW, NKHS, NKD, & 1514 NDTAB, QBI 1515 USE W3ODATMD, ONLY: UNDEF, FLOGRD 1516 ! 1517 IMPLICIT NONE 1518 !/ 1519 !/ ------------------------------------------------------------------- / 1520 !/ Parameter list 1521 !/ 1522 REAL, INTENT(IN) :: A(NSPEC), K(NK), CG(NK), & 1523 DEPTH, USTAR, USDIR 1524 INTEGER, INTENT(IN) :: IX, IY 1525 REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC) 1526 REAL, INTENT(OUT) :: WHITECAP(1:4) 1527 !/ 1528 !/ ------------------------------------------------------------------- / 1529 !/ Local parameters 1530 !/ 1531 INTEGER :: IS, IS2, IS0, IS20, IA, J, IKL, ITOT, NTOT, ID, NKL, IO 1532 INTEGER :: IK, ITH, I_INT, IK2, ITH2, IKIND, L, & 1533 IKHS, IKD, SDSNTH, IT 1534 INTEGER :: NSMOOTH(NK) 1535 REAL :: FACTOR, COSWIND, ASUM, SDIAGISO 1536 REAL :: COEF1, COEF2, COEF3, COEF4(NK) 1537 REAL :: ALFAMEAN, KB 1538 REAL :: FACTURB, DTURB, DCUMULATIVE, BREAKFRACTION 1539 REAL :: RENEWALFREQ, EPSR 1540 REAL :: NTIMES(NK), S1(NK), E1(NK) 1541 REAL :: GAM, XT, M1 1542 REAL :: DK(NK), HS(NK), KBAR(NK), DCK(NK) 1543 REAL :: EFDF(NK) ! Energy integrated over a spectral band 1544 INTEGER :: IKSUP(NK) 1545 REAL :: Q1(NK) 1546 REAL :: FACSAT, DKHS 1547 REAL :: BTH0(NK) !saturation spectrum 1548 REAL :: BTH(NSPEC) !saturation spectrum 1549 REAL :: BTH0S(NK) !smoothed saturation spectrum 1550 REAL :: BTHS(NSPEC) !smoothed saturation spectrum 1551 REAL :: W, MICHE, X 1552 REAL :: QB(NK), S2(NK), KD, DC(NK) Page 45 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1553 REAL :: TSTR, TMAX, DT, T, MFT 1554 REAL :: PB(NSPEC), PB2(NSPEC), LAMBDA(NSPEC) 1555 LOGICAL :: MASK(NSPEC) 1556 !/ 1557 !/ ------------------------------------------------------------------- / 1558 !/ 1559 ! 1560 !---------------------------------------------------------------------- 1561 ! 1562 ! 1. Initialization and numerical factors 1563 ! 1564 FACTURB=SSDSC(5)*USTAR**2/GRAV*DAIR/DWAT 1565 BREAKFRACTION=0. 1566 RENEWALFREQ=0. 1567 ! 1568 ! 2. Estimation of spontaneous breaking 1569 ! 1570 IF ( (SSDSBCK-SSDSC(1)).LE.0 ) THEN 1571 ! 1572 ! 2.a Case of a direction-dependent breaking term (TEST441) 1573 ! 1574 S =0. 1575 D =0. 1576 PB =0. 1577 EPSR=SQRT(SSDSBR) 1578 ! 1579 ! 2.a.1 Computes saturation 1580 ! 1581 SDSNTH = MIN(NINT(SSDSDTH/(DTH*RADE)),NTH/2-1) 1582 DO IK=1, NK 1583 FACSAT=SIG(IK)*K(IK)**3*DTH 1584 IS0=(IK-1)*NTH 1585 BTH(IS0+1)=0. 1586 BTH0(IK)=SUM(A(IS0+1:IS0+NTH))*FACSAT 1587 IF (SSDSDTH.GE.180) THEN ! integrates around full circle 1588 BTH(IS0+1:IS0+NTH)=BTH0(IK) 1589 ELSE 1590 DO ITH=1,NTH ! partial integration 1591 IS=ITH+(IK-1)*NTH 1592 BTH(IS)=DOT_PRODUCT(SATWEIGHTS(:,ITH), & 1593 A(IS0+SATINDICES(:,ITH)) )*FACSAT 1594 ! BTH(IS)=SUM( SATWEIGHTS(:,ITH)* & 1595 ! A(IS0+SATINDICES(:,ITH)) )*FACSAT 1596 END DO 1597 IF (SSDSISO.EQ.1) THEN 1598 BTH0(IK)=SUM(A(IS0+1:IS0+NTH))*FACSAT 1599 ELSE 1600 BTH0(IK)=MAXVAL(BTH(IS0+1:IS0+NTH)) 1601 END IF 1602 END IF 1603 END DO 1604 ! 1605 ! Optional smooting of B and B0 over frequencies 1606 ! 1607 IF (SSDSBRFDF.GT.0.AND.SSDSBRFDF.LT.NK/2) THEN 1608 BTH0S(:)=BTH0(:) 1609 BTHS(:)=BTH(:) Page 46 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1610 NSMOOTH(:)=1 1611 DO IK=1, SSDSBRFDF 1612 BTH0S(1+SSDSBRFDF)=BTH0S(1+SSDSBRFDF)+BTH0(IK) 1613 NSMOOTH(1+SSDSBRFDF)=NSMOOTH(1+SSDSBRFDF)+1 1614 DO ITH=1,NTH 1615 IS=ITH+(IK-1)*NTH 1616 BTHS(ITH+SSDSBRFDF*NTH)=BTHS(ITH+SSDSBRFDF*NTH)+BTH(IS) 1617 END DO 1618 END DO 1619 DO IK=2+SSDSBRFDF,1+2*SSDSBRFDF 1620 BTH0S(1+SSDSBRFDF)=BTH0S(1+SSDSBRFDF)+BTH0(IK) 1621 NSMOOTH(1+SSDSBRFDF)=NSMOOTH(1+SSDSBRFDF)+1 1622 DO ITH=1,NTH 1623 IS=ITH+(IK-1)*NTH 1624 BTHS(ITH+SSDSBRFDF*NTH)=BTHS(ITH+SSDSBRFDF*NTH)+BTH(IS) 1625 END DO 1626 END DO 1627 DO IK=SSDSBRFDF,1,-1 1628 BTH0S(IK)=BTH0S(IK+1)-BTH0(IK+SSDSBRFDF+1) 1629 NSMOOTH(IK)=NSMOOTH(IK+1)-1 1630 DO ITH=1,NTH 1631 IS=ITH+(IK-1)*NTH 1632 BTHS(IS)=BTHS(IS+NTH)-BTH(IS+(SSDSBRFDF+1)*NTH) 1633 END DO 1634 END DO 1635 ! 1636 DO IK=2+SSDSBRFDF,NK-SSDSBRFDF 1637 BTH0S(IK)=BTH0S(IK-1)-BTH0(IK-SSDSBRFDF-1)+BTH0(IK+SSDSBRFDF) 1638 NSMOOTH(IK)=NSMOOTH(IK-1) 1639 DO ITH=1,NTH 1640 IS=ITH+(IK-1)*NTH 1641 BTHS(IS)=BTHS(IS-NTH)-BTH(IS-(SSDSBRFDF+1)*NTH)+BTH(IS+(SSDSBRFDF)*NTH) 1642 END DO 1643 END DO 1644 ! 1645 DO IK=NK-SSDSBRFDF+1,NK 1646 BTH0S(IK)=BTH0S(IK-1)-BTH0(IK-SSDSBRFDF) 1647 NSMOOTH(IK)=NSMOOTH(IK-1)-1 1648 DO ITH=1,NTH 1649 IS=ITH+(IK-1)*NTH 1650 BTHS(IS)=BTHS(IS-NTH)-BTH(IS-(SSDSBRFDF+1)*NTH) 1651 END DO 1652 END DO 1653 ! 1654 ! final division by NSMOOTH 1655 ! 1656 BTH0(:)=MAX(0.,BTH0S(:)/NSMOOTH(:)) 1657 DO IK=1,NK 1658 IS0=(IK-1)*NTH 1659 BTH(IS0+1:IS0+NTH)=MAX(0.,BTHS(IS0+1:IS0+NTH)/NSMOOTH(IK)) 1660 END DO 1661 END IF 1662 ! 1663 ! 2.a.2 Computes spontaneous breaking dissipation rate 1664 ! 1665 DO IK=1, NK 1666 ! Page 47 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1667 ! Correction of saturation level for shallow-water kinematics 1668 ! 1669 IF (SSDSBM(0).EQ.1) THEN 1670 MICHE=1. 1671 ELSE 1672 X=TANH(MIN(K(IK)*DEPTH,10.)) 1673 MICHE=(X*(SSDSBM(1)+X*(SSDSBM(2)+X*(SSDSBM(3)+X*SSDSBM(4)))))**2 ! Correction of saturation level for shallow-wa 1673 ter kine 1674 END IF 1675 COEF1=(SSDSBR*MICHE) 1676 ! 1677 ! Computes isotropic part 1678 ! 1679 SDIAGISO = SSDSC(2) * SIG(IK)*SSDSC(6)*(MAX(0.,BTH0(IK)/COEF1-1.))**2 1680 ! 1681 ! Computes anisotropic part and sums isotropic part 1682 ! 1683 COEF2=SSDSC(2) * SIG(IK)*(1-SSDSC(6))/(COEF1*COEF1) 1684 COEF3=-2.*SIG(IK)*K(IK)*FACTURB 1685 D((IK-1)*NTH+1:IK*NTH) = SDIAGISO + & 1686 COEF2*((MAX(0.,BTH((IK-1)*NTH+1:IK*NTH)-COEF1))**SSDSP) 1687 END DO 1688 ! 1689 ! Computes Breaking probability 1690 ! 1691 PB = (MAX(SQRT(BTH)-EPSR,0.))**2 1692 ! 1693 ! Multiplies by 28.16 = 22.0 * 1.6² * 1/2 with 1694 ! 22.0 (Banner & al. 2000, figure 6) 1695 ! 1.6 the coefficient that transforms SQRT(B) to Banner et al. (2000)'s epsilon 1696 ! 1/2 factor to correct overestimation of Banner et al. (2000)'s breaking probability due to zero-crossing analysis 1697 ! 1698 PB = PB * 28.16 1699 !/ 1700 END IF ! End of test for (Ardhuin et al. 2010)'s spontaneous dissipation source term 1701 ! 1702 ! 2.b Computes spontaneous breaking for T500 ////////////// 1703 ! 1704 IF (SSDSBCK.GT.0) THEN ! test for (Filipot et al. 2010)'s disspation source term 1705 E1 =0. 1706 HS =0. 1707 S =0. 1708 D =0. 1709 PB2=0. 1710 ! 1711 ! Computes Wavenumber spectrum E1 integrated over direction and computes dk 1712 ! 1713 DO IK=1, NK 1714 E1(IK)=0. 1715 DO ITH=1,NTH 1716 IS=ITH+(IK-1)*NTH 1717 E1(IK)=E1(IK)+(A(IS)*SIG(IK))*DTH 1718 END DO 1719 DK(IK)=DDEN(IK)/(DTH*SIG(IK)*CG(IK)) 1720 END DO 1721 ! 1722 ! Gets windows indices of IKTAB Page 48 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1723 ! 1724 ID=MIN(NINT(DEPTH),NDTAB) 1725 IF (ID < 1) THEN 1726 ID = 1 1727 ELSE IF(ID > NDTAB) THEN 1728 ID = NDTAB 1729 END IF 1730 ! 1731 ! loop over wave scales 1732 ! 1733 HS=0. 1734 EFDF=0. 1735 KBAR=0. 1736 EFDF=0. 1737 NKL=0. !number of windows 1738 DO IKL=1,NK 1739 IKSUP(IKL)=IKTAB(IKL,ID) 1740 IF (IKSUP(IKL) .LE. NK) THEN 1741 EFDF(IKL) = DOT_PRODUCT(E1(IKL:IKSUP(IKL)-1),DK(IKL:IKSUP(IKL)-1)) 1742 IF (EFDF(IKL) .NE. 0) THEN 1743 KBAR(IKL) = DOT_PRODUCT(K(IKL:IKSUP(IKL)-1)*E1(IKL:IKSUP(IKL)-1), & 1744 DK(IKL:IKSUP(IKL)-1)) / EFDF(IKL) 1745 ELSE 1746 KBAR(IKL)=0. 1747 END IF 1748 ! estimation of Significant wave height of a given scale 1749 HS(IKL) = 4*SQRT(EFDF(IKL)) 1750 NKL = NKL+1 1751 END IF 1752 END DO 1753 ! 1754 ! Computes Dissipation and breaking probability in each scale 1755 ! 1756 DCK=0. 1757 QB =0. 1758 DKHS = KHSMAX/NKHS 1759 DO IKL=1, NKL 1760 IF (HS(IKL) .NE. 0. .AND. KBAR(IKL) .NE. 0.) THEN 1761 ! 1762 ! gets indices for tabulated dissipation DCKI and breaking probability QBI 1763 ! 1764 IKD = FAC_KD2+ANINT(LOG(KBAR(IKL)*DEPTH)/LOG(FAC_KD1)) 1765 IKHS= 1+ANINT(KBAR(IKL)*HS(IKL)/DKHS) 1766 IF (IKD > NKD) THEN ! Deep water 1767 IKD = NKD 1768 ELSE IF (IKD < 1) THEN ! Shallow water 1769 IKD = 1 1770 END IF 1771 IF (IKHS > NKHS) THEN 1772 IKHS = NKHS 1773 ELSE IF (IKHS < 1) THEN 1774 IKHS = 1 1775 END IF 1776 XT = TANH(KBAR(IKL)*DEPTH) 1777 ! 1778 ! Gamma corrected for water depth 1779 ! Page 49 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1780 GAM=1.0314*(XT**3)-1.9958*(XT**2)+1.5522*XT+0.1885 1781 ! 1782 ! Computes the energy dissipated for the scale IKL 1783 ! using DCKI which is tabulated in INSIN4 1784 ! 1785 DCK(IKL)=((KBAR(IKL)**(-2.5))*(KBAR(IKL)/(2*PI)))*DCKI(IKHS,IKD) 1786 ! 1787 ! Get the breaking probability for the scale IKL 1788 ! 1789 QB(IKL) = QBI(IKHS,IKD) ! QBI is tabulated in INSIN4 1790 ELSE 1791 DCK(IKL)=0. 1792 QB(IKL) =0. 1793 END IF 1794 END DO 1795 ! 1796 ! Distributes scale dissipation over the frequency spectrum 1797 ! 1798 S1 = 0. 1799 S2 = 0. 1800 NTIMES = 0. 1801 DO IKL=1, NKL 1802 IF (EFDF(IKL) .GT. 0.) THEN 1803 S1(IKL:IKSUP(IKL)) = S1(IKL:IKSUP(IKL)) + & 1804 DCK(IKL)*E1(IKL:IKSUP(IKL)) / EFDF(IKL) 1805 S2(IKL:IKSUP(IKL)) = S2(IKL:IKSUP(IKL)) + & 1806 QB(IKL) *E1(IKL:IKSUP(IKL)) / EFDF(IKL) 1807 NTIMES(IKL:IKSUP(IKL)) = NTIMES(IKL:IKSUP(IKL)) + 1 1808 END IF 1809 END DO 1810 ! 1811 ! Finish the average 1812 ! 1813 WHERE (NTIMES .NE. 0.) 1814 S1 = S1 / NTIMES 1815 S2 = S2 / NTIMES 1816 ELSEWHERE 1817 S1 = 0. 1818 S2 = 0. 1819 END WHERE 1820 S1 = S1 / SIG ! goes back to action for dissipation source term 1821 ! 1822 ! Makes Isotropic distribution 1823 ! 1824 ASUM = 0. 1825 DO IK = 1, NK 1826 ASUM = (SUM(A(((IK-1)*NTH+1):(IK*NTH)))*DTH) 1827 IF (ASUM.GT.1.E-8) THEN 1828 FORALL (IS=1+(IK-1)*NTH:IK*NTH) D(IS) = S1(IK)/ASUM 1829 FORALL (IS=1+(IK-1)*NTH:IK*NTH) PB2(IS) = S2(IK)/ASUM 1830 ELSE 1831 FORALL (IS=1+(IK-1)*NTH:IK*NTH) D(IS) = 0. 1832 FORALL (IS=1+(IK-1)*NTH:IK*NTH) PB2(IS) = 0. 1833 END IF 1834 IF (PB2(1+(IK-1)*NTH).GT.0.001) THEN 1835 BTH0(IK) = 2.*SSDSBR 1836 ELSE Page 50 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1837 BTH0(IK) = 0. 1838 END IF 1839 END DO 1840 ! 1841 PB = (1-SSDSC(1))*PB2*A + SSDSC(1)*PB 1842 ! 1843 END IF ! END OF TEST ON SSDSBCK 1844 ! 1845 ! 3. Computes Lambda from breaking probability 1846 ! 1847 ! Compute Lambda = PB* l(k,th) 1848 ! with l(k,th)=1/(2*pi²)= the breaking crest density 1849 ! 1850 LAMBDA = PB / (2.*PI**2.) 1851 ! 1852 !/ ------------------------------------------------------------------- / 1853 ! WAVE-TURBULENCE INTERACTION AND CUMULATIVE EFFECT 1854 !/ ------------------------------------------------------------------- / 1855 ! 1856 ! loop over spectrum 1857 ! 1858 DO IK=1, NK 1859 DO ITH=1,NTH 1860 IS=ITH+(IK-1)*NTH 1861 ! 1862 ! Computes cumulative effect from Breaking probability 1863 ! 1864 RENEWALFREQ = 0. 1865 IF (SSDSC(3).NE.0 .AND. IK.GT.DIKCUMUL) THEN 1866 DO IK2=1,IK-DIKCUMUL 1867 IF (BTH0(IK2).GT.SSDSBR) THEN 1868 IS2=(IK2-1)*NTH 1869 RENEWALFREQ=RENEWALFREQ+DOT_PRODUCT(CUMULW(IS2+1:IS2+NTH,IS),LAMBDA(IS2+1:IS2+NTH)) 1870 !DO ITH2=1,NTH 1871 ! IS2=ITH2+(IK2-1)*NTH 1872 ! RENEWALFREQ=RENEWALFREQ+CUMULW(IS2,IS)*LAMBDA(IS2) 1873 ! END DO 1874 END IF 1875 END DO 1876 END IF 1877 ! 1878 ! Computes wave turbulence interaction 1879 ! 1880 COSWIND=(ECOS(IS)*COS(USDIR)+ESIN(IS)*SIN(USDIR)) 1881 DTURB=-2.*SIG(IK)*K(IK)*FACTURB*COSWIND ! Theory -> stress direction 1882 ! 1883 ! Add effects 1884 ! 1885 D(IS) = D(IS) + (SSDSC(3)*RENEWALFREQ+DTURB) 1886 END DO 1887 END DO 1888 ! 1889 !/ ------------------------------------------------------------------- / 1890 ! COMPUTES SOURCES TERM 1891 !/ ------------------------------------------------------------------- / 1892 ! 1893 S = D * A Page 51 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1894 ! 1895 !/ ------------------------------------------------------------------- / 1896 ! COMPUTES WHITECAP PARAMETERS 1897 !/ ------------------------------------------------------------------- / 1898 ! 1899 IF ( .NOT. (FLOGRD(5,7).OR.FLOGRD(5,8) ) ) THEN 1900 RETURN 1901 END IF 1902 !/ 1903 WHITECAP(1:2) = 0. 1904 ! 1905 ! precomputes integration of Lambda over direction 1906 ! times wavelength times a (a=5 in Reul&Chapron JGR 2003) times dk 1907 ! 1908 DO IK=1,NK 1909 COEF4(IK) = SUM(LAMBDA((IK-1)*NTH+1:IK*NTH) * DTH) *(2*PI/K(IK)) * & 1910 SSDSC(7) * DDEN(IK)/(DTH*SIG(IK)*CG(IK)) 1911 ! NB: SSDSC(7) is WHITECAPWIDTH 1912 END DO 1913 !/ 1914 IF ( FLOGRD(5,7) ) THEN 1915 !/ 1916 !/ Computes the Total WhiteCap Coverage (a=5. ; Reul and Chapron, 2003) 1917 !/ 1918 DO IK=1,NK 1919 WHITECAP(1) = WHITECAP(1) + COEF4(IK) * (1-WHITECAP(1)) 1920 END DO 1921 END IF 1922 !/ 1923 IF ( FLOGRD(5,8) ) THEN 1924 !/ 1925 !/ Calculates the Mean Foam Thickness for component K(IK) => Fig.3, Reul and Chapron, 2003 1926 !/ 1927 DO IK=1,NK 1928 ! Duration of active breaking (TAU*) 1929 TSTR = 0.8 * 2*PI/SIG(IK) 1930 ! Time persistence of foam (a=5.) 1931 TMAX = 5. * 2*PI/SIG(IK) 1932 DT = TMAX / 50 1933 MFT = 0. 1934 DO IT = 1, 50 1935 ! integration over time of foam persistance 1936 T = FLOAT(IT) * DT 1937 ! Eq. 5 and 6 of Reul and Chapron, 2003 1938 IF ( T .LT. TSTR ) THEN 1939 MFT = MFT + 0.4 / (K(IK)*TSTR) * T * DT 1940 ELSE 1941 MFT = MFT + 0.4 / K(IK) * EXP(-1*(T-TSTR)/3.8) * DT 1942 END IF 1943 END DO 1944 MFT = MFT / TMAX 1945 ! 1946 ! Computes foam-layer thickness (Reul and Chapron, 2003) 1947 ! 1948 WHITECAP(2) = WHITECAP(2) + COEF4(IK) * MFT 1949 END DO 1950 END IF Page 52 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1951 ! 1952 ! End of output computing 1953 ! 1954 RETURN 1955 ! 1956 ! Formats 1957 ! 1958 !/ 1959 !/ End of W3SDS4 ----------------------------------------------------- / 1960 !/ 1961 END SUBROUTINE W3SDS4 ENTRY POINTS Name w3src4md_mp_w3sds4_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References A Dummy 1440 R(4) 4 1 0 ARG,IN 1586,1593,1598,1717,1826,1841,1893 ALFAMEAN Local 1537 R(4) 4 scalar ANINT Func 1764 scalar 1764,1765 ASUM Local 1535 R(4) 4 scalar 1824,1826,1827,1828,1829 BREAKFRACTION Local 1538 R(4) 4 scalar 1565 BTH Local 1548 R(4) 4 1 0 1585,1588,1592,1600,1609,1616,1624 ,1632,1641,1650,1659,1686,1691 BTH0 Local 1547 R(4) 4 1 0 1586,1588,1598,1600,1608,1612,1620 ,1628,1637,1646,1656,1679,1835,183 7,1867 BTH0S Local 1549 R(4) 4 1 0 1608,1612,1620,1628,1637,1646,1656 BTHS Local 1550 R(4) 4 1 0 1609,1616,1624,1632,1641,1650,1659 CG Dummy 1440 R(4) 4 1 0 ARG,IN 1719,1910 COEF1 Local 1536 R(4) 4 scalar 1675,1679,1683,1686 COEF2 Local 1536 R(4) 4 scalar 1683,1686 COEF3 Local 1536 R(4) 4 scalar 1684 COEF4 Local 1536 R(4) 4 1 0 1909,1919,1948 CONSTANTS Module 1507 1507 COS Func 1880 scalar 1880 COSWIND Local 1535 R(4) 4 scalar 1880,1881 CUMULW Local 1513 R(4) 4 2 1 PTR 1513,1869 D Dummy 1442 R(4) 4 1 0 ARG,OUT 1575,1685,1708,1828,1831,1885,1893 DAIR Param 1564 R(4) 4 scalar 1564 DC Local 1552 R(4) 4 1 0 DCK Local 1542 R(4) 4 1 0 1756,1785,1791,1804 DCKI Local 1512 R(4) 4 2 1 PTR 1512,1785 DCUMULATIVE Local 1538 R(4) 4 scalar DDEN Local 1508 R(4) 4 1 1 PTR 1508,1719,1910 DEPTH Dummy 1440 R(4) 4 scalar ARG,IN 1672,1724,1764,1776 DK Local 1542 R(4) 4 1 0 1719,1741,1744 DKHS Local 1546 R(4) 4 scalar 1758,1765 DOT_PRODUCT Func 1592 scalar 1592,1741,1743,1869 DSIP Local 1510 R(4) 4 1 1 PTR 1510 Page 53 Source Listing W3SDS4 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References DT Local 1553 R(4) 4 scalar 1932,1936,1939,1941 DTH Local 1510 R(4) 4 scalar PTR 1510,1581,1583,1717,1719,1826,1909 ,1910 DTURB Local 1538 R(4) 4 scalar 1881,1885 DWAT Param 1564 R(4) 4 scalar 1564 E1 Local 1540 R(4) 4 1 0 1705,1714,1717,1741,1743,1804,1806 ECOS Local 1510 R(4) 4 1 1 PTR 1510,1880 EFDF Local 1543 R(4) 4 1 0 1734,1736,1741,1742,1744,1749,1802 ,1804,1806 EPSR Local 1539 R(4) 4 scalar 1577,1691 ESIN Local 1510 R(4) 4 1 1 PTR 1510,1880 EXP Func 1941 scalar 1941 FACSAT Local 1546 R(4) 4 scalar 1583,1586,1593,1598 FACTOR Local 1535 R(4) 4 scalar FACTURB Local 1538 R(4) 4 scalar 1564,1684,1881 FLOAT Func 1936 scalar 1936 FLOGRD Local 1515 L(4) 4 2 1 PTR 1515,1899,1914,1923 GAM Local 1541 R(4) 4 scalar 1780 GRAV Param 1564 R(4) 4 scalar 1564 HS Local 1542 R(4) 4 1 0 1706,1733,1749,1760,1765 IA Local 1531 I(4) 4 scalar ID Local 1531 I(4) 4 scalar 1724,1725,1726,1727,1728,1739 IK Local 1532 I(4) 4 scalar 1582,1583,1584,1586,1588,1591,1598 ,1600,1611,1612,1615,1619,1620,162 3,1627,1628,1629,1631,1636,1637,16 38,1640,1645,1646,1647,1649,1657,1 658,1659,1665,1672,1679,1683,1684, 1685,1686,1713,1714,1716,1717,1719 ,1825,1826,1828,1829,1831,1832,183 4,1835,1837,1858,1860,1865,1866,18 81,1908,1909,1910,1918,1919,1927,1 929,1931,1939,1941,1948 IK2 Local 1532 I(4) 4 scalar 1866,1867,1868 IKD Local 1533 I(4) 4 scalar 1764,1766,1767,1768,1769,1785,1789 IKHS Local 1533 I(4) 4 scalar 1765,1771,1772,1773,1774,1785,1789 IKIND Local 1532 I(4) 4 scalar IKL Local 1531 I(4) 4 scalar 1738,1739,1740,1741,1742,1743,1744 ,1746,1749,1759,1760,1764,1765,177 6,1785,1789,1791,1792,1801,1802,18 03,1804,1805,1806,1807 IKSUP Local 1544 I(4) 4 1 0 1739,1740,1741,1743,1744,1803,1804 ,1805,1806,1807 IKTAB Local 1512 I(4) 4 2 1 PTR 1512,1739 IO Local 1531 I(4) 4 scalar IS Local 1531 I(4) 4 scalar 1591,1592,1615,1616,1623,1624,1631 ,1632,1640,1641,1649,1650,1716,171 7,1860,1869,1880,1885 IS Local 1828 I(4) 4 scalar 1828,1829,1831,1832 IS0 Local 1531 I(4) 4 scalar 1584,1585,1586,1588,1593,1598,1600 ,1658,1659 IS2 Local 1531 I(4) 4 scalar 1868,1869 IS20 Local 1531 I(4) 4 scalar IT Local 1533 I(4) 4 scalar 1934,1936 ITH Local 1532 I(4) 4 scalar 1590,1591,1592,1593,1614,1615,1616 ,1622,1623,1624,1630,1631,1639,164 Page 54 Source Listing W3SDS4 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References 0,1648,1649,1715,1716,1859,1860 ITH2 Local 1532 I(4) 4 scalar ITOT Local 1531 I(4) 4 scalar IX Dummy 1442 I(4) 4 scalar ARG,IN IY Dummy 1442 I(4) 4 scalar ARG,IN I_INT Local 1532 I(4) 4 scalar J Local 1531 I(4) 4 scalar K Dummy 1440 R(4) 4 1 0 ARG,IN 1583,1672,1684,1743,1881,1909,1939 ,1941 KB Local 1537 R(4) 4 scalar KBAR Local 1542 R(4) 4 1 0 1735,1743,1746,1760,1764,1765,1776 ,1785 KD Local 1552 R(4) 4 scalar L Local 1532 I(4) 4 scalar LAMBDA Local 1554 R(4) 4 1 0 1850,1869,1909 LOG Func 1764 scalar 1764 M1 Local 1541 R(4) 4 scalar MASK Local 1555 L(4) 4 1 0 MAX Func 1656 scalar 1656,1659,1679,1686,1691 MAXVAL Func 1600 scalar 1600 MFT Local 1553 R(4) 4 scalar 1933,1939,1941,1944,1948 MICHE Local 1551 R(4) 4 scalar 1670,1673,1675 MIN Func 1581 scalar 1581,1672,1724 NDTAB Param 1514 I(4) 4 scalar 1514,1724,1727,1728 NINT Func 1581 scalar 1581,1724 NK Local 1508 I(4) 4 scalar PTR 1508,1522,1534,1536,1540,1542,1543 ,1544,1545,1547,1549,1552,1582,160 7,1636,1645,1657,1665,1713,1738,17 40,1825,1858,1908,1918,1927 NKD Param 1513 I(4) 4 scalar 1513,1766,1767 NKHS Param 1513 I(4) 4 scalar 1513,1758,1771,1772 NKL Local 1531 I(4) 4 scalar 1737,1750,1759,1801 NSMOOTH Local 1534 I(4) 4 1 0 1610,1613,1621,1629,1638,1647,1656 ,1659 NSPEC Local 1508 I(4) 4 scalar PTR 1508,1522,1525,1548,1550,1554,1555 NTH Local 1508 I(4) 4 scalar PTR 1508,1581,1584,1586,1588,1590,1591 ,1598,1600,1614,1615,1616,1622,162 3,1624,1630,1631,1632,1639,1640,16 41,1648,1649,1650,1658,1659,1685,1 686,1715,1716,1826,1828,1829,1831, 1832,1834,1859,1860,1868,1869,1909 NTIMES Local 1540 R(4) 4 1 0 1800,1807,1813,1814,1815 NTOT Local 1531 I(4) 4 scalar PB Local 1554 R(4) 4 1 0 1576,1691,1698,1841,1850 PB2 Local 1554 R(4) 4 1 0 1709,1829,1832,1834,1841 PI Param 1785 R(4) 4 scalar 1785,1850,1909,1929,1931 Q1 Local 1545 R(4) 4 1 0 QB Local 1552 R(4) 4 1 0 1757,1789,1792,1806 QBI Local 1514 R(4) 4 2 1 PTR 1514,1789 RADE Param 1581 R(4) 4 scalar 1581 RENEWALFREQ Local 1539 R(4) 4 scalar 1566,1864,1869,1885 S Dummy 1440 R(4) 4 1 0 ARG,OUT 1574,1707,1893 S1 Local 1540 R(4) 4 1 0 1798,1803,1814,1817,1820,1828 S2 Local 1552 R(4) 4 1 0 1799,1805,1815,1818,1829 SATINDICES Local 1513 I(4) 4 2 1 PTR 1513,1593 Page 55 Source Listing W3SDS4 2014-09-16 16:52 Symbol Table w3src4md.f90 Name Object Declared Type Bytes Dimen Elements Attributes References SATWEIGHTS Local 1513 R(4) 4 2 1 PTR 1513,1592 SDIAGISO Local 1535 R(4) 4 scalar 1679,1685 SDSNTH Local 1533 I(4) 4 scalar 1581 SIG Local 1510 R(4) 4 1 1 PTR 1510,1583,1679,1683,1684,1717,1719 ,1820,1881,1910,1929,1931 SIN Func 1880 scalar 1880 SQRT Func 1577 scalar 1577,1691,1749 SSDSBCK Local 1512 R(4) 4 scalar PTR 1512,1570,1704 SSDSBINT Local 1512 R(4) 4 scalar PTR 1512 SSDSBM Local 1511 R(4) 4 1 1 PTR 1511,1669,1673 SSDSBR Local 1508 R(4) 4 scalar PTR 1508,1577,1675,1835,1867 SSDSBR2 Local 1511 R(4) 4 scalar PTR 1511 SSDSBRFDF Local 1512 I(4) 4 scalar PTR 1512,1607,1611,1612,1613,1616,1619 ,1620,1621,1624,1627,1628,1632,163 6,1637,1641,1645,1646,1650 SSDSC Local 1509 R(4) 4 1 1 PTR 1509,1564,1570,1679,1683,1841,1865 ,1885,1910 SSDSDTH Local 1511 R(4) 4 scalar PTR 1511,1581,1587 SSDSISO Local 1511 I(4) 4 scalar PTR 1511,1597 SSDSP Local 1510 R(4) 4 scalar PTR 1510,1686 SUM Func 1586 scalar 1586,1598,1826,1909 T Local 1553 R(4) 4 scalar 1936,1938,1939,1941 TANH Func 1672 scalar 1672,1776 TMAX Local 1553 R(4) 4 scalar 1931,1932,1944 TSTR Local 1553 R(4) 4 scalar 1929,1938,1939,1941 UNDEF Local 1515 R(4) 4 scalar 1515 USDIR Dummy 1440 R(4) 4 scalar ARG,IN 1880 USTAR Dummy 1440 R(4) 4 scalar ARG,IN 1564 W Local 1551 R(4) 4 scalar W3GDATMD Module 1508 1508 W3ODATMD Module 1515 1515 W3SDS4 Subr 1440 WHITECAP Dummy 1442 R(4) 4 1 4 ARG,OUT 1903,1919,1948 X Local 1551 R(4) 4 scalar 1672,1673 XT Local 1541 R(4) 4 scalar 1776,1780 Page 56 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 1962 1963 1964 END MODULE W3SRC4MD SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References BETATB Local 92 R(4) 4 2 16842 PRIV 92 DDRAG Local 91 R(4) 4 scalar PRIV 91 DRAGMX Param 89 R(4) 4 scalar PRIV 89 DSIGA Local 91 R(4) 4 scalar PRIV 91 FIRST Local 94 L(4) 4 scalar PRIV 94 G Param 81 R(4) 4 scalar HA Local 64 R(4) 4 2 75 HA2 Local 64 R(4) 4 2 75 HD Local 64 R(4) 4 2 75 NH Local 63 I(4) 4 1 3 NHMAX Param 62 I(4) 4 scalar 63,64 NRDRAG Param 87 I(4) 4 scalar PRIV 87,92 NRSIGA Param 86 I(4) 4 scalar PRIV 86,92 PI Param 81 R(4) 4 scalar SIGAMX Param 88 R(4) 4 scalar PRIV 88 THO Local 63 I(4) 4 3 150 W3SRC4MD Module 2 XSTRESS Local 77 R(4) 4 1 1 ALC YSTRESS Local 77 R(4) 4 1 1 ALC Page 57 Source Listing W3SDS4 2014-09-16 16:52 Subprograms/Common Blocks w3src4md.f90 SUBPROGRAMS/COMMON BLOCKS Name Object Declared Type Bytes Dimen Elements Attributes References CALC_USTAR Subr 1344 276 INSIN4 Subr 677 TABU_STRESS Subr 931 772 TABU_TAUHF Subr 1053 773 TABU_TAUHF2 Subr 1189 776 W3SDS4 Subr 1440 W3SIN4 Subr 292 W3SPR4 Subr 98 W3SRC4MD Module 2 COMPILER OPTIONS BEING USED -align nocommons -align nodcommons -align noqcommons -align records -align nosequence -align norec1byte -align norec2byte -align norec4byte -align norec8byte -align norec16byte -altparam -assume accuracy_sensitive -assume nobscc -assume nobuffered_io -assume byterecl -assume nocc_omp -assume nocstring -assume nodummy_aliases -assume nofpe_summary -assume noieee_fpe_flags -assume nominus0 -assume noold_boz -assume old_unit_star -assume old_ldout_format -assume noold_logical_ldio -assume old_maxminloc -assume old_xor -assume protect_constants -assume noprotect_parens -assume split_common -assume source_include -assume nostd_intent_in -assume nostd_mod_proc_name -assume norealloc_lhs -assume underscore -assume no2underscores -auto no -auto_scalar no -bintext -ccdefault default -check noargs -check noarg_temp_created -check nobounds -check noformat -check nooutput_conversion -check nooverflow -check nopointers -check power -check noshape -check nounderflow -check nouninitialized -coarray-num-procs 0 no -coarray-config-file -convert big_endian -cross_reference -D __INTEL_COMPILER=1210 -D __unix__ -D __unix -D __linux__ -D __linux -D __gnu_linux__ -D unix -D linux -D __ELF__ -D __x86_64 -D __x86_64__ -D _MT -D __INTEL_COMPILER_BUILD_DATE=20120612 -D __pentium4 -D __pentium4__ -D __tune_pentium4__ -D __SSE2__ -D __SSE3__ -D __SSSE3__ -D __SSE4_1__ -D __SSE4_2__ -D __SSE__ -D __MMX__ Page 58 Source Listing W3SDS4 2014-09-16 16:52 w3src4md.f90 -D __AVX__ -double_size 64 no -d_lines no -Qdyncom -error_limit 30 no -f66 no -f77rtl no -fast -fpscomp nofilesfromcmd -fpscomp nogeneral -fpscomp noioformat -fpscomp noldio_spacing -fpscomp nologicals no -fpconstant -fpe3 -fprm nearest no -ftz -fp_model noprecise -fp_model fast -fp_model nostrict -fp_model nosource -fp_model nodouble -fp_model noextended -fp_model novery_fast -fp_model noexcept -fp_model nono_except -heap_arrays 0 no -threadprivate_compat -free -g0 -iface nomixed_str_len_arg -iface nono_mixed_str_len_arg no -intconstant -integer_size 32 no -mixed_str_len_arg no -module -names lowercase no -noinclude -O2 no -pad_source -real_size 32 no -recursive -reentrancy none no -sharable_localsaves -vec=simd -show noinclude -show map -show options no -syntax_only no -threadcom no -U no -vms -w noall -w nonone -w alignments -w noargument_checking -w nodeclarations -w general -w noignore_bounds -w noignore_loc -w nointerfaces -w notruncated_source -w uncalled -w uninitialized -w nounused -w usage -includepath : /gpfs/gp1/usrx/local/intel/composer_xe_2011_sp1.11.339/compiler/include/,.f,./.f,/usrx/local/intel/composerxe/mkl/include/.f, /usrx/local/intel/composerxe/tbb/include/.f,/gpfs/gp1/usrx/local/intel/composer_xe_2011_sp1.11.339/compiler/include/intel64/.f, /gpfs/gp1/usrx/local/intel/composer_xe_2011_sp1.11.339/compiler/include/.f,/usr/local/include/.f,/usr/lib/gcc/x86_64-redhat-linux/4.4.7/include/.f, /usr/include/.f,/usr/include/.f -list filename : w3src4md.lst -o filename : none COMPILER: Intel(R) Fortran 12.1-2100