!> @file !> @brief The 'WAM4+' source terms based on P.A.E.M. Janssen's work, with !> extensions by him and by J.-R. Bidlot. !> !> @author F. Ardhuin !> @author H. L. Tolman !> @date 02-Sep-2012 !> #include "w3macros.h" !/ ------------------------------------------------------------------- / !> !> @brief The 'WAM4+' source terms based on P.A.E.M. Janssen's work, with !> extensions by him and by J.-R. Bidlot. !> !> @details Converted from the original WAM codes by F. Ardhuin, !> with further extensions to adapt to a saturation-based breaking !> and observation-based swell dissipation. !> !> @author F. Ardhuin !> @author H. L. Tolman !> @date 02-Sep-2012 !> !> @copyright Copyright 2009-2022 National Weather Service (NWS), !> National Oceanic and Atmospheric Administration. All rights !> reserved. WAVEWATCH III is a trademark of the NWS. !> No unauthorized use without permission. !> MODULE W3SRC3MD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III SHOM | !/ ! F. Ardhuin ! !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 02-Sep-2012 | !/ +-----------------------------------+ !/ !/ 09-Oct-2007 : Origination. ( version 3.13 ) !/ 10-Oct-2010 : Adding Janssen-style swell damping ( version 3.14) !/ 02-Sep-2012 : Clean up test output T, T0, T1 ( version 4.07 ) !/ ! 1. Purpose : ! ! The 'WAM4+' source terms based on P.A.E.M. Janssen's work, with ! extensions by him and by J.-R. Bidlot. Converted from the original ! WAM codes by F. Ardhuin, with further extensions to adapt to a ! saturation-based breaking and observation-based swell dissipation ! ! ! 2. Variables and types : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 3. Subroutines and functions : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3SPR3 Subr. Public Mean parameters from spectrum. ! W3SIN3 Subr. Public WAM4+ input source term. ! INSIN3 Subr. Public Corresponding initialization routine. ! TABU_STRESS, TABU_TAUHF, TABU_TAUHF2 ! Subr. Public Populate various tables. ! CALC_USTAR ! Subr. Public Compute stresses. ! W3SDS3 Subr. Public User supplied dissipation. ! ---------------------------------------------------------------- ! ! 4. Subroutines and functions used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 5. Remarks : ! ! 6. Switches : ! ! 7. Source code : !/ !/ ------------------------------------------------------------------- / !/ USE CONSTANTS, ONLY: KAPPA, nu_air IMPLICIT NONE PUBLIC !/ !/ Public variables !/ !air kinematic viscosity (used in WAM) INTEGER, PARAMETER :: ITAUMAX=200,JUMAX=200 INTEGER, PARAMETER :: IUSTAR=100,IALPHA=200, ILEVTAIL=50 REAL :: TAUT(0:ITAUMAX,0:JUMAX), DELTAUW, DELU ! Table for H.F. stress as a function of 2 variables REAL :: TAUHFT(0:IUSTAR,0:IALPHA), DELUST, DELALP REAL, PARAMETER :: UMAX = 50. REAL, PARAMETER :: TAUWMAX = 2.2361 !SQRT(5.) !/ CONTAINS !/ ------------------------------------------------------------------- / !> !> @brief Calculate mean wave parameters for the use in the source term !> routines. !> !> @param[in] A Action density spectrum. !> @param[in] CG Group velocities. !> @param[in] WN Wavenumbers. !> @param[out] EMEAN Energy. !> @param[out] FMEAN Mean frequency for determination of tail. !> @param[out] FMEANS Mean frequency for dissipation source term. !> @param[out] WNMEAN Mean wavenumber. !> @param[out] AMAX Maximum of action spectrum. !> @param[in] U Wind speed. !> @param[in] UDIR Wind direction. !> @param[inout] USTAR Friction velocity !> @param[inout] USDIR Wind stress direction. !> @param[in] TAUWX Component of wave-supported stress. !> @param[in] TAUWY Component of wave-supported stress. !> @param[out] CD Drag coefficient at wind level ZWND. !> @param[out] Z0 Corresponding z0. !> @param[out] CHARN !> @param[in] LLWS Wind sea true/false array for each component. !> @param[out] FMEANWS Mean frequency of wind sea, used for tail. !> !> @author F. Ardhuin !> @author H. L. Tolman !> @date 17-Oct-2007 !> SUBROUTINE W3SPR3 (A, CG, WN, EMEAN, FMEAN, FMEANS, WNMEAN, & AMAX, U, UDIR, USTAR, USDIR, TAUWX, TAUWY, CD, Z0,& CHARN, LLWS, FMEANWS) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III SHOM | !/ ! F. Ardhuin ! !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 17-Oct-2007 | !/ +-----------------------------------+ !/ !/ 03-Oct-2007 : Origination. ( version 3.13 ) !/ ! 1. Purpose : ! ! Calculate mean wave parameters for the use in the source term ! routines. ! ! 2. Method : ! ! See source term routines. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! A R.A. I Action density spectrum. ! CG R.A. I Group velocities. ! WN R.A. I Wavenumbers. ! EMEAN Real O Energy ! FMEAN Real O Mean frequency for determination of tail ! FMEANS Real O Mean frequency for dissipation source term ! WNMEAN Real O Mean wavenumber. ! AMAX Real O Maximum of action spectrum. ! U Real I Wind speed. ! UDIR Real I Wind direction. ! USTAR Real I/O Friction velocity. ! USDIR Real I/O wind stress direction. ! TAUWX-Y Real I Components of wave-supported stress. ! CD Real O Drag coefficient at wind level ZWND. ! Z0 Real O Corresponding z0. ! LLWS L.A. I Wind sea true/false array for each component ! FMEANWS Real O Mean frequency of wind sea, used for tail ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! STRACE Service routine. ! ! 5. Called by : ! ! W3SRCE Source term integration routine. ! W3OUTP Point output program. ! GXEXPO GrADS point output program. ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable test output. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS, ONLY: TPIINV USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, DTH, DDEN, WWNMEANP, & WWNMEANPTAIL, FTE, FTF, SSTXFTF, SSTXFTWN,& SSTXFTFTAIL, SSWELLF #ifdef W3_T USE W3ODATMD, ONLY: NDST #endif #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL, INTENT(IN) :: A(NTH,NK), CG(NK), WN(NK), U, UDIR REAL, INTENT(IN) :: TAUWX, TAUWY LOGICAL, INTENT(IN) :: LLWS(NSPEC) REAL, INTENT(INOUT) :: USTAR ,USDIR REAL, INTENT(OUT) :: EMEAN, FMEAN, FMEANS, WNMEAN, AMAX, & CD, Z0, CHARN, FMEANWS !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IS, IK, ITH #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif REAL :: TAUW, EBAND, EMEANWS, UNZ, & EB(NK),EB2(NK),ALFA(NK) !/ !/ ------------------------------------------------------------------- / !/ #ifdef W3_S CALL STRACE (IENT, 'W3SPR3') #endif ! UNZ = MAX ( 0.01 , U ) USTAR = MAX ( 0.0001 , USTAR ) ! EMEAN = 0. EMEANWS= 0. FMEANWS= 0. FMEAN = 0. FMEANS = 0. WNMEAN = 0. AMAX = 0. ! ! 1. Integral over directions and maximum --------------------------- * ! DO IK=1, NK EB(IK) = 0. EB2(IK) = 0. DO ITH=1, NTH IS=ITH+(IK-1)*NTH EB(IK) = EB(IK) + A(ITH,IK) IF (LLWS(IS)) EB2(IK) = EB2(IK) + A(ITH,IK) AMAX = MAX ( AMAX , A(ITH,IK) ) END DO END DO ! ! 2. Integrate over directions -------------------------------------- * ! DO IK=1, NK ALFA(IK) = 2. * DTH * SIG(IK) * EB(IK) * WN(IK)**3 EB(IK) = EB(IK) * DDEN(IK) / CG(IK) EB2(IK) = EB2(IK) * DDEN(IK) / CG(IK) EMEAN = EMEAN + EB(IK) FMEAN = FMEAN + EB(IK) *(SIG(IK)**(2.*WWNMEANPTAIL)) FMEANS = FMEANS + EB(IK) *(SIG(IK)**(2.*WWNMEANP)) WNMEAN = WNMEAN + EB(IK) *(WN(IK)**WWNMEANP) EMEANWS = EMEANWS+ EB2(IK) FMEANWS = FMEANWS+ EB2(IK)*(SIG(IK)**(2.*WWNMEANPTAIL)) END DO ! ! 3. Add tail beyond discrete spectrum and get mean pars ------------ * ! ( DTH * SIG absorbed in FTxx ) ! EBAND = EB(NK) / DDEN(NK) EMEAN = EMEAN + EBAND * FTE FMEAN = FMEAN + EBAND * SSTXFTFTAIL FMEANS = FMEANS + EBAND * SSTXFTF WNMEAN = WNMEAN + EBAND * SSTXFTWN EBAND = EB2(NK) / DDEN(NK) EMEANWS = EMEANWS + EBAND * FTE FMEANWS = FMEANWS + EBAND * SSTXFTFTAIL ! ! 4. Final processing ! IF (FMEAN.LT.1.E-7) THEN FMEAN=TPIINV * SIG(NK) ELSE FMEAN = TPIINV *( MAX ( 1.E-7 , FMEAN ) & / MAX ( 1.E-7 , EMEAN ))**(1/(2.*WWNMEANPTAIL)) END IF IF (FMEANS.LT.1.E-7) THEN FMEANS=TPIINV * SIG(NK) ELSE FMEANS = TPIINV *( MAX ( 1.E-7 , FMEANS ) & / MAX ( 1.E-7 , EMEAN ))**(1/(2.*WWNMEANP)) END IF WNMEAN = ( MAX ( 1.E-7 , WNMEAN ) & / MAX ( 1.E-7 , EMEAN ) )**(1/WWNMEANP) IF (FMEANWS.LT.1.E-7.OR.EMEANWS.LT.1.E-7) THEN FMEANWS=TPIINV * SIG(NK) ELSE FMEANWS = TPIINV *( MAX ( 1.E-7 , FMEANWS ) & / MAX ( 1.E-7 , EMEANWS ))**(1/(2.*WWNMEANPTAIL)) END IF ! ! 5. Cd and z0 ----------------------------------------------- * ! TAUW = SQRT(TAUWX**2+TAUWY**2) Z0=0. CALL CALC_USTAR(U,TAUW,USTAR,Z0,CHARN) UNZ = MAX ( 0.01 , U ) CD = (USTAR/UNZ)**2 USDIR = UDIR ! ! 6. Final test output ---------------------------------------------- * ! #ifdef W3_T WRITE (NDST,9060) EMEAN, WNMEAN, TPIINV, USTAR, CD, Z0 #endif ! RETURN ! ! Formats ! #ifdef W3_T 9060 FORMAT (' TEST W3SPR3 : E,WN MN :',F8.3,F8.4/ & ' TPIINV, USTAR, CD, Z0 :',F8.3,F7.2,1X,2F9.5) #endif !/ !/ End of W3SPR3 ----------------------------------------------------- / !/ END SUBROUTINE W3SPR3 !/ ------------------------------------------------------------------- / !> !> @brief Calculate diagonal and input source term for WAM4+ approach. !> !> @verbatim !> WAM-4 : Janssen et al. !> WAM-"4.5" : gustiness effect (Cavaleri et al. ) !> SWELLF: damping coefficient (=1) for Janssen (2004) theory !> @endverbatim !> !> @param[in] A Action density spectrum (1-D). !> @param[in] CG Group speed. !> @param[in] K Wavenumber for entire spectrum. !> @param[in] U Wind speed. !> @param[in] USTAR Friction velocity. !> @param[in] DRAT Air/water density ratio. !> @param[in] AS Air-sea temperature difference. !> @param[in] USDIR Wind stress direction. !> @param[in] Z0 Air-sea roughness length. !> @param[in] CD Wind drag coefficient. !> @param[out] TAUWX Component of the wave-supported stress. !> @param[out] TAUWY Component of the wave-supported stress. !> @param[out] TAUWNX Component of the negative wave-supported stress. !> @param[out] TAUWNY Component of the negative wave-supported stress. !> @param[in] ICE Sea ice fraction. !> @param[out] S Source term (1-D version). !> @param[out] D Diagonal term of derivative. !> @param[out] LLWS Wind sea true/false array for each component. !> @param[in] IX !> @param[in] IY !> !> @author F. Ardhuin !> @author H. L. Tolman !> @date 02-Sep-2012 !> SUBROUTINE W3SIN3 (A, CG, K, U, USTAR, DRAT, AS, USDIR, Z0, CD, & TAUWX, TAUWY, TAUWNX, TAUWNY, ICE, S, D, LLWS, IX, IY) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III SHOM | !/ ! F. Ardhuin ! !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 02-Sep-2012 | !/ +-----------------------------------+ !/ !/ 09-Oct-2007 : Origination. ( version 3.13 ) !/ 16-May-2010 : Adding sea ice ( version 3.14_Ifremer ) !/ 02-Sep-2012 : Clean up test output T, T0, T1 ( version 4.07 ) !/ ! 1. Purpose : ! ! Calculate diagonal and input source term for WAM4+ approach. ! ! 2. Method : ! ! WAM-4 : Janssen et al. ! WAM-"4.5" : gustiness effect (Cavaleri et al. ) ! SWELLF: damping coefficient (=1) for Janssen (2004) theory ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! A R.A. I Action density spectrum (1-D). ! CG R.A. I Group speed *) ! K R.A. I Wavenumber for entire spectrum. *) ! U Real I WIND SPEED ! USTAR Real I Friction velocity. ! DRAT Real I Air/water density ratio. ! AS Real I Air-sea temperature difference ! USDIR Real I wind stress direction ! Z0 Real I Air-side roughness lengh. ! CD Real I Wind drag coefficient. ! USDIR Real I Direction of friction velocity ! TAUWX-Y Real I Components of the wave-supported stress. ! TAUWNX Real I Component of the negative wave-supported stress. ! TAUWNY Real I Component of the negative wave-supported stress. ! ICE Real I Sea ice fraction. !/Stefan: ICE is DUMMY argument; remove later. ! S R.A. O Source term (1-D version). ! D R.A. O Diagonal term of derivative. *) ! LLWS L.A. O Wind sea true/false array for each component ! ---------------------------------------------------------------- ! *) Stored as 1-D array with dimension NTH*NK ! ! 4. Subroutines used : ! ! STRACE Subroutine tracing. ( !/S switch ) ! PRT2DS Print plot of spectrum. ( !/T0 switch ) ! OUTMAT Print out matrix. ( !/T1 switch ) ! ! 5. Called by : ! ! W3SRCE Source term integration. ! W3EXPO Point output program. ! GXEXPO GrADS point output program. ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable general test output. ! !/T0 2-D print plot of source term. ! !/T1 Print arrays. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS, ONLY: GRAV, TPI #ifdef W3_T USE CONSTANTS, ONLY: RADE #endif USE W3GDATMD, ONLY: NK, NTH, NSPEC, XFR, DDEN, SIG, SIG2, TH, & ESIN, ECOS, EC2, ZZWND, AALPHA, BBETA, ZZALP,& SSWELLF, & DDEN2, DTH, SSINTHP,ZZ0RAT #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif #ifdef W3_T USE W3ODATMD, ONLY: NDST #endif #ifdef W3_T0 USE W3ODATMD, ONLY: NDST #endif #ifdef W3_T1 USE W3ODATMD, ONLY: NDST #endif #ifdef W3_T0 USE W3ARRYMD, ONLY: PRT2DS #endif #ifdef W3_T1 USE W3ARRYMD, ONLY: OUTMAT #endif ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL, INTENT(IN) :: A(NSPEC) REAL, INTENT(IN) :: CG(NK), K(NSPEC),Z0,U, CD REAL, INTENT(IN) :: USTAR, USDIR, AS, DRAT, ICE !/Stefan: ICE is DUMMY argument; remove later. REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC), TAUWX, TAUWY, TAUWNX, TAUWNY LOGICAL, INTENT(OUT) :: LLWS(NSPEC) INTEGER, INTENT(IN) :: IX, IY !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IS,IK,ITH #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif REAL :: COSU, SINU, TAUX, TAUY REAL :: UST2, TAUW, TAUWB REAL , PARAMETER :: EPS1 = 0.00001, EPS2 = 0.000001 #ifdef W3_STAB3 REAL :: Usigma !standard deviation of U due to gustiness REAL :: USTARsigma !standard deviation of USTAR due to gustiness #endif REAL :: CM,UCO,UCN,ZCN, & Z0VISC REAL XI,DELI1,DELI2 REAL XJ,DELJ1,DELJ2 REAL :: CONST, CONST0, CONST2, CONST3, TAU1 REAL :: X,ZARG,ZLOG,ZBETA,UST REAL COSWIND,XSTRESS,YSTRESS,TAUHF REAL TEMP, TEMP2 INTEGER IND,J,ISTAB REAL DSTAB(3,NSPEC) REAL STRESSSTAB(3,2),STRESSSTABN(3,2) #ifdef W3_T0 REAL :: DOUT(NK,NTH) #endif !/ !/ ------------------------------------------------------------------- / !/ #ifdef W3_S CALL STRACE (IENT, 'W3SIN3') #endif ! #ifdef W3_T WRITE (NDST,9000) BBETA, USTAR, USDIR*RADE #endif ! ! 1. Preparations ! ! JDM: initializing arrays (shouldn't change answers) DSTAB=0.; STRESSSTAB=0. ;STRESSSTABN=0. ! ! 1.a estimation of surface roughness parameters ! Z0VISC = 0.1*nu_air/MAX(USTAR,0.0001) ! ! 2. Diagonal ! ! Here AS is the air-sea temperature difference in degrees. Expression given by ! Abdalla & Cavaleri, JGR 2002 for Usigma. For USTARsigma ... I do not see where ! I got it from, maybe just made up from drag law ... ! #ifdef W3_STAB3 Usigma=MAX(0.,-0.025*AS) USTARsigma=(1.0+U/(10.+U))*Usigma #endif UST=USTAR ISTAB=3 #ifdef W3_STAB3 DO ISTAB=1,2 IF (ISTAB.EQ.1) UST=USTAR*(1.-USTARsigma) IF (ISTAB.EQ.2) UST=USTAR*(1.+USTARsigma) #endif TAUX = UST**2* COS(USDIR) TAUY = UST**2* SIN(USDIR) ! ! Loop over the resolved part of the spectrum ! STRESSSTAB(ISTAB,:)=0. STRESSSTABN(ISTAB,:)=0. CONST0=BBETA*DRAT/(kappa**2) ! SSWELLF(1) is IDAMPING in ECMWAM CONST3 = SSWELLF(1)*2.*KAPPA*DRAT ! COSU = COS(USDIR) SINU = SIN(USDIR) DO IK=1, NK IS=1+(IK-1)*NTH CM=K(IS)/SIG2(IS) !inverse of phase speed UCN=UST*CM+ZZALP !this is the inverse wave age ! ! the stress is the real stress (N/m^2) divided by ! rho_a, and thus comparable to USTAR**2 ! it is the integral of rho_w g Sin/C /rho_a ! (air-> waves momentum flux) ! CONST2=DDEN2(IS)/CG(IK) & !Jacobian to get energy in band *GRAV/(SIG(IK)/K(IS)*DRAT) ! coefficient to get momentum CONST=SIG2(IS)*CONST0 ! this CM parameter is 1 / C_phi ! this is the "correct" shallow-water expression ! here Z0 corresponds to Z0+Z1 of the Janssen eq. 14 ZCN=ALOG(K(IS)*Z0) ! commented below is the original WAM version (OK for deep water) ! ZCN=ALOG(G*Z0b(I)*CM(I)**2) DO ITH=1,NTH IS=ITH+(IK-1)*NTH COSWIND=(ECOS(IS)*COSU+ESIN(IS)*SINU) IF (COSWIND.GT.0.01) THEN X=COSWIND*UCN ! this ZARG term is the argument of the exponential ! in Janssen 1991 eq. 16. ZARG=KAPPA/X ! ZLOG is ALOG(MU) where MU is defined by Janssen 1991 eq. 15 ZLOG=ZCN+ZARG ZBETA = CONST3*(COSWIND+KAPPA/(ZCN*UST*CM))*UCN**2 IF (ZLOG.LT.0.) THEN ! The source term Sp is beta * omega * X**2 ! as given by Janssen 1991 eq. 19 DSTAB(ISTAB,IS) = CONST*EXP(ZLOG)*ZLOG**4*UCN**2*COSWIND**SSINTHP LLWS(IS)=.TRUE. ELSE DSTAB(ISTAB,IS) = ZBETA LLWS(IS)=.TRUE. END IF ELSE ZBETA = CONST3*(COSWIND+KAPPA/(ZCN*UST*CM))*UCN**2 DSTAB(ISTAB,IS) = ZBETA LLWS(IS)=.FALSE. END IF ! TEMP2=CONST2*DSTAB(ISTAB,IS)*A(IS) IF (DSTAB(ISTAB,IS).LT.0) THEN STRESSSTABN(ISTAB,1)=STRESSSTABN(ISTAB,1)+TEMP2*ECOS(IS) STRESSSTABN(ISTAB,2)=STRESSSTABN(ISTAB,2)+TEMP2*ESIN(IS) END IF STRESSSTAB(ISTAB,1)=STRESSSTAB(ISTAB,1)+TEMP2*ECOS(IS) STRESSSTAB(ISTAB,2)=STRESSSTAB(ISTAB,2)+TEMP2*ESIN(IS) END DO END DO ! D(:)=DSTAB(3,:) XSTRESS=STRESSSTAB (3,1) YSTRESS=STRESSSTAB (3,2) TAUWNX =STRESSSTABN(3,1) TAUWNY =STRESSSTABN(3,2) ! WRITE(995,'(A,11G14.5)') 'NEGSTRESS: ',TAUWNX,TAUWNY,FW*UORB**3 #ifdef W3_STAB3 END DO D(:)=0.5*(DSTAB(1,:)+DSTAB(2,:)) XSTRESS=0.5*(STRESSSTAB(1,1)+STRESSSTAB(2,1)) YSTRESS=0.5*(STRESSSTAB(1,2)+STRESSSTAB(2,2)) TAUWNX=0.5*(STRESSSTABN(1,1)+STRESSSTABN(2,1)) TAUWNY=0.5*(STRESSSTABN(1,2)+STRESSSTABN(2,2)) #endif S = D * A ! ! ... Test output of arrays ! #ifdef W3_T0 DO IK=1, NK DO ITH=1, NTH DOUT(IK,ITH) = D(ITH+(IK-1)*NTH) END DO END DO CALL PRT2DS (NDST, NK, NK, NTH, DOUT, SIG(1:NK), ' ', 1., & 0.0, 0.001, 'Diag Sin', ' ', 'NONAME') #endif ! #ifdef W3_T1 CALL OUTMAT (NDST, D, NTH, NTH, NK, 'diag Sin') #endif ! ! Computes the high-frequency contribution ! the difference in spectal density (kx,ky) to (f,theta) ! is integrated in this modified CONST0 CONST0=DTH*SIG(NK)**5/((GRAV**2)*tpi) & *TPI*SIG(NK) / CG(NK) !conversion WAM (E(f,theta) to WW3 A(k,theta) TEMP=0. DO ITH=1,NTH IS=ITH+(NK-1)*NTH COSWIND=(ECOS(IS)*COSU+ESIN(IS)*SINU) TEMP=TEMP+A(IS)*(MAX(COSWIND,0.))**3 END DO ! ! finds the values in the tabulated stress TAUHFT ! XI=UST/DELUST IND = MAX(1,MIN (IUSTAR-1, INT(XI))) DELI1= MAX(MIN (1. ,XI-FLOAT(IND)),0.) DELI2= 1. - DELI1 XJ=MAX(0.,(GRAV*Z0/MAX(UST,0.00001)**2-AALPHA) / DELALP) J = MAX(1 ,MIN (IALPHA-1, INT(XJ))) DELJ1= MAX(0.,MIN (1. , XJ-FLOAT(J))) DELJ2=1. - DELJ1 TAU1 =(TAUHFT(IND,J)*DELI2+TAUHFT(IND+1,J)*DELI1 )*DELJ2 & +(TAUHFT(IND,J+1)*DELI2+TAUHFT(IND+1,J+1)*DELI1)*DELJ1 TAUHF = CONST0*TEMP*UST**2*TAU1 TAUWX = XSTRESS+TAUHF*COS(USDIR) TAUWY = YSTRESS+TAUHF*SIN(USDIR) ! ! Reduces tail effect to make sure that wave-supported stress ! is less than total stress, this is borrowed from ECWAM Stresso.F ! TAUW = SQRT(TAUWX**2+TAUWY**2) UST2 = MAX(USTAR,EPS2)**2 TAUWB = MIN(TAUW,MAX(UST2-EPS1,EPS2**2)) IF (TAUWB.LT.TAUW) THEN TAUWX=TAUWX*TAUWB/TAUW TAUWY=TAUWY*TAUWB/TAUW END IF RETURN ! ! Formats ! #ifdef W3_T 9000 FORMAT (' TEST W3SIN3 : COMMON FACT.: ',3E10.3) #endif !/ !/ End of W3SIN3 ----------------------------------------------------- / !/ END SUBROUTINE W3SIN3 !/ ------------------------------------------------------------------- / !> !> @brief Initialization for source term routine. !> !> @author F. Ardhuin !> @date 23-Jul-2009 !> SUBROUTINE INSIN3 !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | SHOM | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 23-Jul-2009 | !/ +-----------------------------------+ !/ !/ 23-Jun-2006 : Origination. ( version 3.09 ) !/ 23-Jul-2007 : Cleaning up convolutions ( version 3.14-SHOM) ! ! 1. Purpose : ! ! Initialization for source term routine. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3SIN3 Subr. W3SRC3MD Corresponding source term. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS, ONLY: TPIINV USE W3GDATMD, ONLY: SIG, NK #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif !/ IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif !/ !/ ------------------------------------------------------------------- / !/ #ifdef W3_S CALL STRACE (IENT, 'INSIN3') #endif ! ! 1. .... ----------------------------------------------------------- * ! CALL TABU_STRESS CALL TABU_TAUHF(SIG(NK) * TPIINV) !tabulate high-frequency stress !/ !/ End of INSIN3 ----------------------------------------------------- / !/ END SUBROUTINE INSIN3 ! ---------------------------------------------------------------------- !> !> @brief To generate the friction velocity table, TAUT(TAUW,U10)=SQRT(TAU). !> !> @author F. Ardhuin !> @date 17-Oct-2007 !> SUBROUTINE TABU_STRESS !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 17-Oct-2007 | !/ +-----------------------------------+ !/ !/ 23-Jun-2006 : Origination. ( version 3.13 ) !/ adapted from WAM, original:P.A.E.M. JANSSEN KNMI AUGUST 1990 !/ adapted version (subr. STRESS): J. BIDLOT ECMWF OCTOBER 2004 !/ Table values were checkes against the original f90 result and found to !/ be identical (at least at 0.001 m/s accuracy) !/ ! 1. Purpose : ! TO GENERATE friction velocity table TAUT(TAUW,U10)=SQRT(TAU). ! METHOD. ! A STEADY STATE WIND PROFILE IS ASSUMED. ! THE WIND STRESS IS COMPUTED USING THE ROUGHNESSLENGTH ! Z1=Z0/SQRT(1-TAUW/TAU) ! WHERE Z0 IS THE CHARNOCK RELATION , TAUW IS THE WAVE- ! INDUCED STRESS AND TAU IS THE TOTAL STRESS. ! WE SEARCH FOR STEADY-STATE SOLUTIONS FOR WHICH TAUW/TAU < 1. ! FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990. ! ! Initialization for source term routine. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3SIN3 Subr. W3SRC3MD Corresponding source term. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS, ONLY: GRAV USE W3GDATMD, ONLY: ZZWND, AALPHA, ZZ0MAX IMPLICIT NONE INTEGER, PARAMETER :: NITER=10 REAL , PARAMETER :: XM=0.50, EPS1=0.00001 ! VARIABLE. TYPE. PURPOSE. ! *XM* REAL POWER OF TAUW/TAU IN ROUGHNESS LENGTH. ! *XNU* REAL KINEMATIC VISCOSITY OF AIR. ! *NITER* INTEGER NUMBER OF ITERATIONS TO OBTAIN TOTAL STRESS ! *EPS1* REAL SMALL NUMBER TO MAKE SURE THAT A SOLUTION ! IS OBTAINED IN ITERATION WITH TAU>TAUW. ! ---------------------------------------------------------------------- INTEGER I,J,ITER REAL ZTAUW,UTOP,CDRAG,WCD,USTOLD,TAUOLD REAL X,UST,ZZ0,ZNU,F,DELF,ZZ00 ! ! DELU = UMAX/FLOAT(JUMAX) DELTAUW = TAUWMAX/FLOAT(ITAUMAX) DO I=0,ITAUMAX ZTAUW = (REAL(I)*DELTAUW)**2 DO J=0,JUMAX UTOP = FLOAT(J)*DELU CDRAG = 0.0012875 WCD = SQRT(CDRAG) USTOLD = UTOP*WCD TAUOLD = MAX(USTOLD**2, ZTAUW+EPS1) DO ITER=1,NITER X = ZTAUW/TAUOLD UST = SQRT(TAUOLD) ZZ00=AALPHA*TAUOLD/GRAV IF (ZZ0MAX.NE.0) ZZ00=MIN(ZZ00,ZZ0MAX) ! Corrects roughness ZZ00 for quasi-linear effect ZZ0 = ZZ00/(1.-X)**XM !ZNU = 0.1*nu_air/UST ! This was removed by Bidlot in 1996 !ZZ0 = MAX(ZNU,ZZ0) F = UST-KAPPA*UTOP/(ALOG(ZZWND/ZZ0)) DELF= 1.-KAPPA*UTOP/(ALOG(ZZWND/ZZ0))**2*2./UST & *(1.-(XM+1)*X)/(1.-X) UST = UST-F/DELF TAUOLD= MAX(UST**2., ZTAUW+EPS1) END DO TAUT(I,J) = SQRT(TAUOLD) END DO END DO I=ITAUMAX J=JUMAX ! ! Force zero wind to have zero stress (Bidlot 1996) ! DO I=0,ITAUMAX TAUT(I,0)=0.0 END DO RETURN END SUBROUTINE TABU_STRESS !/ ------------------------------------------------------------------- / !> !> @brief Tabulation of the high-frequency wave-supported stress. !> !> @details SEE REFERENCE FOR WAVE STRESS CALCULATION. !> FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990. !> See tech. Memo ECMWF 03 december 2003 by Bidlot & Janssen. !> !> @param[in] FRMAX Maximum frequency. !> !> @author F. Ardhuin !> @date 14-Aug-2006 !> SUBROUTINE TABU_TAUHF(FRMAX) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update 2006/08/14 | !/ +-----------------------------------+ !/ !/ 27-Feb-2004 : Origination in WW3 ( version 2.22.SHOM ) !/ the resulting table was checked to be identical to the original f77 result !/ 14-Aug-2006 : Modified following Bidlot ( version 2.22.SHOM ) !/ 18-Aug-2006 : Ported to version 3.09 ! ! 1. Purpose : ! ! Tabulation of the high-frequency wave-supported stress ! ! 2. Method : ! ! SEE REFERENCE FOR WAVE STRESS CALCULATION. ! FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990. ! See tech. Memo ECMWF 03 december 2003 by Bidlot & Janssen ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! FRMAX Real I maximum frequency. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! STRACE Service routine. ! ! 5. Called by : ! ! W3SIN3 Wind input Source term routine. ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable test output. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS, ONLY: GRAV, TPI USE W3GDATMD, ONLY: AALPHA, BBETA, ZZALP, XFR, FACHFE, ZZ0MAX #ifdef W3_T USE W3ODATMD, ONLY: NDST #endif #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL, intent(in) :: FRMAX ! maximum frequency !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ ! USTARM R.A. Maximum friction velocity ! ALPHAM R.A. Maximum Charnock Coefficient ! WLV R.A. Water levels. ! UA R.A. Absolute wind speeds. ! UD R.A. Absolute wind direction. ! U10 R.A. Wind speed used. ! U10D R.A. Wind direction used. ! 10. Source code : ! !/ ------------------------------------------------------------------- / REAL :: USTARM, ALPHAM REAL :: CONST1, OMEGA, OMEGAC REAL :: UST, ZZ0,OMEGACC, CM INTEGER, PARAMETER :: JTOT=250 REAL, ALLOCATABLE :: W(:) REAL :: ZX,ZARG,ZMU,ZLOG,ZZ00,ZBETA REAL :: Y,YC,DELY INTEGER :: I,J,K,L REAL :: X0 #ifdef W3_S INTEGER, SAVE :: IENT = 0 CALL STRACE (IENT, 'TABU_HF') #endif ! USTARM = 5. ALPHAM = 20.*AALPHA DELUST = USTARM/REAL(IUSTAR) DELALP = ALPHAM/REAL(IALPHA) CONST1 = BBETA/KAPPA**2 OMEGAC = TPI*FRMAX ! TAUHFT(0:IUSTAR,0:IALPHA)=0. !table initialization ! ALLOCATE(W(JTOT)) W(2:JTOT-1)=1. W(1)=0.5 W(JTOT)=0.5 X0 = 0.05 ! DO L=0,IALPHA DO K=0,IUSTAR UST = MAX(REAL(K)*DELUST,0.000001) ZZ00 = UST**2*AALPHA/GRAV IF (ZZ0MAX.NE.0) ZZ00=MIN(ZZ00,ZZ0MAX) ZZ0 = ZZ00*(1+FLOAT(L)*DELALP/AALPHA) OMEGACC = MAX(OMEGAC,X0*GRAV/UST) YC = OMEGACC*SQRT(ZZ0/GRAV) DELY = MAX((1.-YC)/REAL(JTOT),0.) ! For a given value of UST and ALPHA, ! the wave-supported stress is integrated all the way ! to 0.05*g/UST DO J=1,JTOT Y = YC+REAL(J-1)*DELY OMEGA = Y*SQRT(GRAV/ZZ0) ! This is the deep water phase speed CM = GRAV/OMEGA !this is the inverse wave age, shifted by ZZALP (tuning) ZX = UST/CM +ZZALP ZARG = MIN(KAPPA/ZX,20.) ZMU = MIN(GRAV*ZZ0/CM**2*EXP(ZARG),1.) ZLOG = MIN(ALOG(ZMU),0.) ZBETA = CONST1*ZMU*ZLOG**4 ! Power of Y in denominator should be FACHFE-4 TAUHFT(K,L) = TAUHFT(K,L)+W(J)*ZBETA/Y*DELY END DO !IF (MOD(K,5).EQ.0.AND.MOD(L,5).EQ.0) & !WRITE(102,'(2I4,3G16.8)') L,K,UST,AALPHA+FLOAT(L)*DELALP,TAUHFT(K,L) #ifdef W3_T WRITE (NDST,9000) L,K,AALPHA+FLOAT(L)*DELALP,UST,TAUHFT(K,L) #endif END DO END DO DEALLOCATE(W) ! DO L=0,IALPHA ! DO K=0,IUSTAR !WRITE(101,'(A,2I4,G16.8)') 'L,K,TAUHFT(K,L):',L,K,TAUHFT(K,L) ! END DO ! END DO !WRITE(101,*) 'TAUHFT:',FRMAX,BBETA,AALPHA,CONST1,OMEGAC,TPI !WRITE(101,'(20G16.8)') TAUHFT RETURN #ifdef W3_T 9000 FORMAT (' TABU_HF, L, K :',(2I4,3F8.3)/) #endif END SUBROUTINE TABU_TAUHF !/ ------------------------------------------------------------------- / !/ ------------------------------------------------------------------- / !> !> @brief Compute friction velocity based on wind speed U10. !> !> @details Computation of u* based on Quasi-linear theory. !> !> @param[in] WINDSPEED 10-m wind speed ... should be NEUTRAL !> @param[in] TAUW Wave-supported stress. !> @param[out] USTAR Friction velocity. !> @param[out] Z0 Air-side roughness length. !> @param[out] CHARN !> !> @author F. Ardhuin !> @date 14-Aug-2006 !> SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update 2006/08/14 | !/ +-----------------------------------+ !/ !/ 27-Feb-2004 : Origination in WW3 ( version 2.22-SHOM ) !/ the resulting table was checked to be identical to the original f77 result !/ 14-Aug-2006 : Modified following Bidlot ( version 2.22-SHOM ) !/ 18-Aug-2006 : Ported to version 3.09 !/ 03-Apr-2010 : Adding output of Charnock parameter ( version 3.14-IFREMER ) ! ! 1. Purpose : ! ! Compute friction velocity based on wind speed U10 ! ! 2. Method : ! ! Computation of u* based on Quasi-linear theory ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! U10,TAUW,USTAR,Z0 ! ---------------------------------------------------------------- ! WINDSPEED Real I 10-m wind speed ... should be NEUTRAL ! TAUW Real I Wave-supported stress ! USTAR Real O Friction velocity. ! Z0 Real O air-side roughness length ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! STRACE Service routine. ! ! 5. Called by : ! ! W3SIN3 Wind input Source term routine. ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable test output. ! ! 10. Source code : !-----------------------------------------------------------------------------! USE CONSTANTS, ONLY: GRAV USE W3GDATMD, ONLY: ZZWND, AALPHA IMPLICIT NONE REAL, intent(in) :: WINDSPEED,TAUW REAL, intent(out) :: USTAR, Z0, CHARN ! local variables REAL SQRTCDM1 REAL X,XI,DELI1,DELI2,XJ,delj1,delj2 REAL UST,DELTOLD,TAUW_LOCAL INTEGER IND,J ! TAUW_LOCAL=MAX(MIN(TAUW,TAUWMAX),0.) XI = SQRT(TAUW_LOCAL)/DELTAUW IND = MIN ( ITAUMAX-1, INT(XI)) ! index for stress table DELI1 = MIN(1.,XI - REAL(IND)) !interpolation coefficient for stress table DELI2 = 1. - DELI1 XJ = WINDSPEED/DELU J = MIN ( JUMAX-1, INT(XJ) ) DELJ1 = MIN(1.,XJ - REAL(J)) DELJ2 = 1. - DELJ1 USTAR=(TAUT(IND,J)*DELI2+TAUT(IND+1,J )*DELI1)*DELJ2 & + (TAUT(IND,J+1)*DELI2+TAUT(IND+1,J+1)*DELI1)*DELJ1 ! ! Determines roughness length ! SQRTCDM1 = MIN(WINDSPEED/USTAR,100.0) Z0 = ZZWND*EXP(-KAPPA*SQRTCDM1) IF (USTAR.GT.0.001) THEN CHARN = GRAV*Z0/USTAR**2 ELSE CHARN = AALPHA END IF ! RETURN END SUBROUTINE CALC_USTAR !/ ------------------------------------------------------------------- / !> !> @brief Calculate whitecapping source term and diagonal term of derivative. !> !> @details WAM-Cycle 4 and following. !> The last update (09-May-2005) follows the redefinition of !> the mean wavenumber as in Bidlot et al. (2005). !> !> @param[in] A Action density spectrum (1-D). !> @param[in] K Wavenumber for entire spectrum. !> @param[in] CG !> @param[in] EMEAN Mean wave energy. !> @param[in] FMEAN Mean wave frequency. !> @param[in] WNMEAN Mean wavenumber. !> @param[in] USTAR Friction velocity. !> @param[in] USDIR Wind stress direction. !> @param[in] DEPTH Water depth. !> @param[out] S Source term (1-D version). !> @param[out] D Diagonal term of derivative. !> @param[in] IX !> @param[in] IY !> !> @author F. Ardhuin !> @date 23-Jul-2009 !> SUBROUTINE W3SDS3 (A, K, CG, EMEAN, FMEAN, WNMEAN, USTAR, USDIR, & DEPTH, S, D, IX, IY) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ ! F. Ardhuin ! !/ | FORTRAN 90 | !/ | Last update : 23-Jul-2009 | !/ +-----------------------------------+ !/ !/ 05-Dec-1996 : Final FORTRAN 77 ( version 1.18 ) !/ 08-Dec-1999 : Upgrade to FORTRAN 90 ( version 2.00 ) !/ 14-Aug-2006 : Generic WAM4+ dissipation term ( version 2.22SHOM ) !/ 23-Jul-2009 : Addition of Filipot &al convolution ( version 3.14-SHOM ) !/ ! 1. Purpose : ! ! Calculate whitecapping source term and diagonal term of derivative. ! ! 2. Method : ! ! WAM-Cycle 4 and following. ! The last update (09-May-2005) follows the redefinition of ! the mean wavenumber as in Bidlot et al. (2005). ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! A R.A. I Action density spectrum (1-D). ! K R.A. I Wavenumber for entire spectrum. *) ! EMEAN Real I Mean wave energy. ! FMEAN Real I Mean wave frequency. ! WNMEAN Real I Mean wavenumber. ! USTAR Real I Friction velocity. ! USDIR Real I wind stress direction. ! DEPTH Real I Water depth. ! S R.A. O Source term (1-D version). ! D R.A. O Diagonal term of derivative. *) ! ---------------------------------------------------------------- ! *) Stored in 1-D array with dimension NTH*NK ! ! 4. Subroutines used : ! ! STRACE Subroutine tracing. ( !/S switch ) ! PRT2DS Print plot of spectrum. ( !/T0 switch ) ! OUTMAT Print out matrix. ( !/T1 switch ) ! ! 5. Called by : ! ! W3SRCE Source term integration. ! W3EXPO Point output program. ! GXEXPO GrADS point output program. ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable general test output. ! !/T0 2-D print plot of source term. ! !/T1 Print arrays. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS, ONLY: GRAV, TPI USE W3GDATMD, ONLY: NSPEC, NTH, NK, DDELTA1, DDELTA2, & #ifdef W3_T0 SIG, & #endif SSDSC1 #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif #ifdef W3_T USE W3ODATMD, ONLY: NDST #endif #ifdef W3_T1 USE W3ODATMD, ONLY: NDST #endif #ifdef W3_T0 USE W3ARRYMD, ONLY: PRT2DS #endif #ifdef W3_T1 USE W3ARRYMD, ONLY: OUTMAT #endif ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL, INTENT(IN) :: A(NSPEC), K(NK), CG(NK), & DEPTH, USTAR, USDIR, EMEAN, FMEAN, WNMEAN INTEGER, INTENT(IN) :: IX, IY REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IS, IK, ITH #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif REAL :: FACTOR, FACTOR2 REAL :: ALFAMEAN, WNMEAN2 #ifdef W3_T0 REAL :: DOUT(NK,NTH) #endif !/ !/ ------------------------------------------------------------------- / !/ #ifdef W3_S CALL STRACE (IENT, 'W3SDS3') #endif ! ! 0. Pre-initialization of arrays, should be set before being used ! but this is helping with bit reproducibility D=0. ! 1. Common factor ! WNMEAN2 = MAX( 1.E-10 , WNMEAN ) ALFAMEAN=WNMEAN**2*EMEAN FACTOR = SSDSC1 * TPI*FMEAN * ALFAMEAN**2 ! #ifdef W3_T WRITE (NDST,9000) SSDSC1, FMEAN, WNMEAN, EMEAN, FACTOR #endif ! !---------------------------------------------------------------------- ! ! 2. Source term ! DO IK=1, NK ! ! Original WAM4/WAM4+ dissipation term ! FACTOR2=FACTOR*(DDELTA1*K(IK)/WNMEAN2 + DDELTA2*(K(IK)/WNMEAN2)**2) DO ITH=1,NTH IS=ITH+(IK-1)*NTH D(IS)= FACTOR2 END DO END DO ! S = D * A ! ! ... Test output of arrays ! #ifdef W3_T0 DO IK=1, NK DO ITH=1, NTH DOUT(IK,ITH) = D(ITH+(IK-1)*NTH) END DO END DO CALL PRT2DS (NDST, NK, NK, NTH, DOUT, SIG(1:NK), ' ', 1., & 0.0, 0.001, 'Diag Sds', ' ', 'NONAME') #endif ! #ifdef W3_T1 CALL OUTMAT (NDST, D, NTH, NTH, NK, 'diag Sds') #endif ! RETURN ! ! Formats ! #ifdef W3_T 9000 FORMAT (' TEST W3SDS3 : COMMON FACT.: ',5E10.3) #endif !/ !/ End of W3SDS3 ----------------------------------------------------- / !/ END SUBROUTINE W3SDS3 END MODULE W3SRC3MD