WAVEWATCH III  beta 0.0.1
Functions/Subroutines | Variables
w3src4md Module Reference

The 'SHOM/Ifremer' source terms based on P.A.E.M. More...

Functions/Subroutines

subroutine w3spr4 (A, CG, WN, EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, U, UDIR, ifdef W3_FLX5
 Calculate mean wave parameters for the use in the source term routines. More...
 
subroutine w3sin4 (A, CG, K, U, USTAR, DRAT, AS, USDIR, Z0, CD, TAUWX, TAUWY, TAUWNX, TAUWNY, S, D, LLWS, IX, IY, BRLAMBDA)
 Calculate diagonal and input source term for WAM4+ approach. More...
 
subroutine insin4 (FLTABS)
 Initialization for source term routine. More...
 
subroutine tabu_stress
 To generate friction velocity table TAUT(TAUW,U10)=SQRT(TAU). More...
 
subroutine tabu_tauhf (SIGMAX)
 Tabulation of the high-frequency wave-supported stress. More...
 
subroutine tabu_tauhf2 (SIGMAX)
 Tabulation of the high-frequency wave-supported stress as a function of ustar, alpha (modified Charnock), and tail energy level. More...
 
subroutine calc_ustar (WINDSPEED, TAUW, USTAR, Z0, CHARN)
 Compute friction velocity based on wind speed U10. More...
 
subroutine w3sds4 (A, K, CG, USTAR, USDIR, DEPTH, DAIR, SRHS, DDIAG, IX, IY, BRLAMBDA, WHITECAP, DLWMEAN)
 Calculate whitecapping source term and diagonal term of derivative. More...
 

Variables

integer, parameter itaumax =200
 
integer, parameter jumax =200
 
integer, parameter iustar =100
 
integer, parameter ialpha =200
 
integer, parameter ilevtail =50
 
real, dimension(:,:), allocatable taut
 
real, dimension(:,:), allocatable tauhft
 
real, dimension(:,:,:), allocatable tauhft2
 
real delust
 
real delalp
 
real deltauw
 
real delu
 
real deltail
 
real, parameter umax = 50.
 
real, parameter tauwmax = 2.2361
 
integer dikcumul
 
integer, parameter nkhi =100
 
integer, parameter fac_kd2 =1000
 
real, parameter fac_kd1 =1.01
 
real, parameter khsmax =2.
 
real, parameter khmax =2.
 
real, parameter kdmax =200000.
 

Detailed Description

The 'SHOM/Ifremer' source terms based on P.A.E.M.

Janssen's wind input and dissipation functions by Ardhuin et al. (2009,2010) and Filipot & Ardhuin (2010) The wind input is converted from the original WAM codes, courtesy of P.A.E.M. Janssen and J. Bidlot

Author
F. Ardhuin
Date
13-Nov-2013
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.

Function/Subroutine Documentation

◆ calc_ustar()

subroutine w3src4md::calc_ustar ( real, intent(in)  WINDSPEED,
real, intent(in)  TAUW,
real, intent(out)  USTAR,
real, intent(out)  Z0,
real, intent(out)  CHARN 
)

Compute friction velocity based on wind speed U10.

Computation of u* based on Quasi-linear theory.

Parameters
[in]WINDSPEED10-m wind speed – should be NEUTRAL.
[in]TAUWWave-supported stress.
[out]USTARFriction velocity.
[out]Z0Air-side roughness length.
[out]CHARNCharnock.
Author
F. Ardhuin
Date
14-Aug-2006

Definition at line 1839 of file w3src4md.F90.

1839  !/
1840  !/ +-----------------------------------+
1841  !/ | WAVEWATCH III NOAA/NCEP |
1842  !/ | F. Ardhuin |
1843  !/ | FORTRAN 90 |
1844  !/ | Last update 2006/08/14 |
1845  !/ +-----------------------------------+
1846  !/
1847  !/ 27-Feb-2004 : Origination in WW3 ( version 2.22-SHOM )
1848  !/ the resulting table was checked to be identical to the original f77 result
1849  !/ 14-Aug-2006 : Modified following Bidlot ( version 2.22-SHOM )
1850  !/ 18-Aug-2006 : Ported to version 3.09
1851  !/ 03-Apr-2010 : Adding output of Charnock parameter ( version 3.14-IFREMER )
1852  !/ 03-May-2024 : Optional functional form of ( version 7.15 )
1853  !/ Charnock coefficient and surface drag (UK Met Office).
1854  !
1855  ! 1. Purpose :
1856  !
1857  ! Compute friction velocity based on wind speed U10
1858  !
1859  ! 2. Method :
1860  !
1861  ! Computation of u* based on Quasi-linear theory
1862  ! uses Charnock relation with modified roughness Z1=Z0/SQRT(1-TAUW/TAU)
1863  !
1864  ! 3. Parameters :
1865  !
1866  ! Parameter list
1867  ! ----------------------------------------------------------------
1868  ! U10,TAUW,USTAR,Z0
1869  ! ----------------------------------------------------------------
1870  ! WINDSPEED Real I 10-m wind speed ... should be NEUTRAL
1871  ! TAUW Real I Wave-supported stress
1872  ! USTAR Real O Friction velocity.
1873  ! Z0 Real O air-side roughness length
1874  ! ----------------------------------------------------------------
1875  !
1876  ! 4. Subroutines used :
1877  !
1878  ! STRACE Service routine.
1879  !
1880  ! 5. Called by :
1881  !
1882  ! W3SIN3 Wind input Source term routine.
1883  !
1884  ! 6. Error messages :
1885  !
1886  ! 7. Remarks :
1887  !
1888  ! 8. Structure :
1889  !
1890  ! See source code.
1891  !
1892  ! 9. Switches :
1893  !
1894  ! !/S Enable subroutine tracing.
1895  ! !/T Enable test output.
1896  !
1897  ! 10. Source code :
1898  !-----------------------------------------------------------------------------!
1899  USE constants, ONLY: grav, kappa, nu_air
1900  USE w3gdatmd, ONLY: zzwnd, aalpha, zz0max, sintailpar, capchnk
1901 #ifdef W3_T
1902  USE w3odatmd, ONLY: ndst
1903 #endif
1904  IMPLICIT NONE
1905  REAL, intent(in) :: WINDSPEED,TAUW
1906  REAL, intent(out) :: USTAR, Z0, CHARN
1907  ! local variables
1908  REAL :: SQRTCDM1
1909  REAL :: XI,DELI1,DELI2,XJ,delj1,delj2 ! used for table version
1910  INTEGER :: IND,J
1911  REAL :: TAUW_LOCAL
1912  REAL :: TAUOLD,CDRAG,WCD,USTOLD,X,UST,ZZ0,ZNU,ZZ00,F,DELF
1913  REAL :: CHATH, XMIN ! used for reduction of high winds
1914  INTEGER, PARAMETER :: NITER=10
1915  REAL , PARAMETER :: XM=0.50, eps1=0.00001
1916  INTEGER :: ITER
1917  ! VARIABLE. TYPE. PURPOSE.
1918  ! *XM* REAL POWER OF TAUW/TAU IN ROUGHNESS LENGTH.
1919  ! *XNU* REAL KINEMATIC VISCOSITY OF AIR.
1920  ! *NITER* INTEGER NUMBER OF ITERATIONS TO OBTAIN TOTAL STRESS
1921  ! *EPS1* REAL SMALL NUMBER TO MAKE SURE THAT A SOLUTION
1922  ! IS OBTAINED IN ITERATION WITH TAU>TAUW.
1923 
1924  chath = aalpha
1925  !
1926  IF (sintailpar(1).GT.0.5) THEN
1927  tauw_local=max(min(tauw,tauwmax),0.)
1928  xi = sqrt(tauw_local)/deltauw
1929  ind = min( itaumax-1, int(xi)) ! index for stress table
1930  deli1 = min(1.,xi - real(ind)) !interpolation coefficient for stress table
1931  deli2 = 1. - deli1
1932  xj = windspeed/delu
1933  j = min( jumax-1, int(xj) )
1934  delj1 = min(1.,xj - real(j))
1935  delj2 = 1. - delj1
1936  ustar=(taut(ind,j)*deli2+taut(ind+1,j )*deli1)*delj2 &
1937  + (taut(ind,j+1)*deli2+taut(ind+1,j+1)*deli1)*delj1
1938  ELSE
1939  IF (capchnk(1).EQ.1.) THEN
1940  ! Computation of sea surface roughness and charnock coefficient based
1941  ! on Donelan (2018). Determines minimum charnock; reduction for winds
1942  ! above a particular threshold
1943  chath = capchnk(2) + 0.5 * (capchnk(3) - capchnk(2)) * (1 &
1944  - tanh((windspeed-capchnk(4))/capchnk(5)))
1945  xmin = 0.15 * (capchnk(3)-chath)
1946  ELSE
1947  xmin = 0.
1948  END IF
1949 
1950  ! This max is for comparison ... to be removed later
1951  ! TAUW_LOCAL=MAX(MIN(TAUW,TAUWMAX),0.)
1952  tauw_local=tauw
1953  cdrag = 0.0012875
1954  wcd = sqrt(cdrag)
1955  ustold = windspeed*wcd
1956  tauold = max(ustold**2, tauw_local+eps1)
1957  ! Newton method to solve for ustar in U=ustar*log(Z/Z0)
1958  DO iter=1,niter
1959  x = max(tauw_local/tauold, xmin)
1960  ust = sqrt(tauold)
1961  zz00 = chath*tauold/grav
1962  IF (zz0max.NE.0) zz00=min(zz00,zz0max)
1963  ! Corrects roughness ZZ00 for quasi-linear effect
1964  zz0 = zz00/(1.-x)**xm
1965  znu = 0.11*nu_air/max(ust,1e-6)
1966  zz0 = sintailpar(5)*znu+zz0
1967  f = ust-kappa*windspeed/(alog(zzwnd/zz0))
1968  delf= 1.-kappa*windspeed/(alog(zzwnd/zz0))**2*2./ust &
1969  *(1.-(xm+1)*x)/(1.-x)
1970  ust = ust-f/delf
1971  tauold= max(ust**2., tauw_local+eps1)
1972  END DO
1973  ustar=ust
1974  END IF
1975  !
1976  ! Determines roughness length
1977  !
1978  IF (ustar.GT.0.001) THEN
1979  sqrtcdm1 = min(windspeed/ustar,100.0)
1980  z0 = zzwnd*exp(-kappa*sqrtcdm1)
1981  charn = grav*z0/ustar**2
1982  ELSE
1983  IF (ustar.GT.0) THEN
1984  sqrtcdm1 = min(windspeed/ustar,100.0)
1985  z0 = zzwnd*exp(-kappa*sqrtcdm1)
1986  ELSE
1987  z0 = chath*0.001*0.001/grav
1988  END IF
1989  charn = chath
1990  END IF
1991  IF(capchnk(1) .EQ. 1) THEN
1992  ! Problem with large values of CHARN for low winds
1993  charn = min( 0.09 , charn )
1994  IF(charn.LT.chath) charn = chath
1995  ENDIF
1996 
1997  ! WRITE(6,*) 'CALC_USTAR:',WINDSPEED,TAUW,AALPHA,CHARN,Z0,USTAR
1998  !
1999  RETURN

References w3gdatmd::aalpha, w3gdatmd::capchnk, deltauw, delu, constants::grav, itaumax, jumax, constants::kappa, w3odatmd::ndst, constants::nu_air, w3gdatmd::sintailpar, taut, tauwmax, w3gdatmd::zz0max, and w3gdatmd::zzwnd.

Referenced by w3spr4().

◆ insin4()

subroutine w3src4md::insin4 ( logical, intent(in)  FLTABS)

Initialization for source term routine.

Parameters
[in]FLTABS
Author
F. Ardhuin
Date
30-Aug-2010

Definition at line 1012 of file w3src4md.F90.

1012  !/
1013  !/ +-----------------------------------+
1014  !/ | WAVEWATCH III NOAA/NCEP |
1015  !/ | SHOM |
1016  !/ | F. Ardhuin |
1017  !/ | FORTRAN 90 |
1018  !/ | Last update : 30-Aug-2010 |
1019  !/ +-----------------------------------+
1020  !/
1021  !/ 30-Aug-2010 : Origination. ( version 3.14-Ifremer )
1022  !
1023  ! 1. Purpose :
1024  !
1025  ! Initialization for source term routine.
1026  !
1027  ! 2. Method :
1028  !
1029  ! 3. Parameters :
1030  !
1031  ! ----------------------------------------------------------------
1032  ! FLTABS Logical
1033  ! ----------------------------------------------------------------
1034  !
1035  ! 4. Subroutines used :
1036  !
1037  ! Name Type Module Description
1038  ! ----------------------------------------------------------------
1039  ! STRACE Subr. W3SERVMD Subroutine tracing.
1040  ! ----------------------------------------------------------------
1041  !
1042  ! 5. Called by :
1043  !
1044  ! Name Type Module Description
1045  ! ----------------------------------------------------------------
1046  ! W3SIN4 Subr. W3SRC3MD Corresponding source term.
1047  ! ----------------------------------------------------------------
1048  !
1049  ! 6. Error messages :
1050  !
1051  ! None.
1052  !
1053  ! 7. Remarks :
1054  !
1055  ! 8. Structure :
1056  !
1057  ! See source code.
1058  !
1059  ! 9. Switches :
1060  !
1061  ! !/S Enable subroutine tracing.
1062  !
1063  ! 10. Source code :
1064  !
1065  !/ ------------------------------------------------------------------- /
1066  USE constants, ONLY: tpiinv, rade, grav
1067  USE w3odatmd, ONLY: ndse
1068  USE w3servmd, ONLY: extcde
1069  USE w3dispmd, ONLY: wavnu2
1070  USE w3gdatmd, ONLY: sig, dsip, nk, nth, ttauwshelter, &
1071  ssdsdth, ssdscos, th, dth, xfr, ecos, esin, &
1072  ssdsc, ssdsbrf1, ssdsbck, ssdsbint, ssdspbk, &
1073  ssdsabk, ssdshck, iktab, dcki, satindices, &
1074  satweights, cumulw, nkhs, nkd, ndtab, qbi, &
1075  sintailpar
1076 #ifdef W3_S
1077  USE w3servmd, ONLY: strace
1078 #endif
1079  !/
1080  IMPLICIT NONE
1081  !/
1082  !/ ------------------------------------------------------------------- /
1083  !/ Parameter list
1084  !/
1085  LOGICAL, INTENT(IN) :: FLTABS
1086  !/
1087  !/ ------------------------------------------------------------------- /
1088  !/
1089  INTEGER SDSNTH, ITH, I_INT, J_INT, IK, IK2, ITH2 , IS, IS2
1090  INTEGER IKL, ID, ICON, IKD, IKHS, IKH, TOTO
1091  REAL C, C2
1092  REAL DIFF1, DIFF2, BINF, BSUP, CGG, PROF
1093  REAL KIK, DHS, KD, KHS, KH, XT, GAM, DKH, PR, W, EPS
1094  REAL DKD
1095  REAL, DIMENSION(:,:) , ALLOCATABLE :: SIGTAB
1096  REAL, DIMENSION(:,:) , ALLOCATABLE :: K1, K2
1097  !/
1098  !/ ------------------------------------------------------------------- /
1099  !/ Local parameters
1100  !/
1101 #ifdef W3_S
1102  INTEGER, SAVE :: IENT = 0
1103 #endif
1104  !/
1105  !/ ------------------------------------------------------------------- /
1106  !/
1107 #ifdef W3_S
1108  CALL strace (ient, 'INSIN4')
1109 #endif
1110  !
1111  ! 1. Initializations ------------------------------------------------ *
1112  !
1113  !
1114  ! These precomputed tables are written in mod_def.ww3
1115  !
1116  IF (sintailpar(1).GT.0.5) THEN
1117  IF (.NOT. ALLOCATED(taut)) ALLOCATE(taut(0:itaumax,0:jumax))
1118  IF (.NOT. ALLOCATED(tauhft)) ALLOCATE(tauhft(0:iustar,0:ialpha))
1119  IF (fltabs) THEN
1120  CALL tabu_stress
1121  CALL tabu_tauhf(sig(nk) ) !tabulate high-frequency stress: 2D table
1122  END IF
1123  IF (ttauwshelter.GT.0) THEN
1124  IF (.NOT. ALLOCATED(tauhft2)) ALLOCATE(tauhft2(0:iustar,0:ialpha,0:ilevtail))
1125  IF (fltabs) CALL tabu_tauhf2(sig(nk) ) !tabulate high-frequency stress: 3D table
1126  END IF
1127  END IF
1128  !
1129  ! 2. SPONTANEOUS BREAKING
1130  ! 2.a Precomputes the indices for integrating the spectrum to get saturation (TEST 4xx )
1131  !
1132  IF (ssdsdth.LT.180) THEN
1133  sdsnth = min(nint(ssdsdth/(dth*rade)),nth/2-1)
1134  satindices(:,:)=1
1135  satweights(:,:)=0.
1136  DO ith=1,nth
1137  DO i_int=ith-sdsnth, ith+sdsnth
1138  j_int=i_int
1139  IF (i_int.LT.1) j_int=i_int+nth
1140  IF (i_int.GT.nth) j_int=i_int-nth
1141  satindices(i_int-(ith-sdsnth)+1,ith)=j_int
1142  satweights(i_int-(ith-sdsnth)+1,ith)= &
1143  cos(th(ith)-th(j_int))**ssdscos
1144  END DO
1145  END DO
1146  ELSE
1147  satindices(:,:)=1
1148  satweights(:,:)=1.
1149  END IF
1150  !/ ------------------------------------------------------------------- /
1151  !
1152  ! Precomputes QBI and DCKI (TEST 500)
1153  !
1154  IF (ssdsbck.GT.0) THEN
1155  !
1156  ! Precomputes the indices for integrating the spectrum over frequency bandwidth
1157  !
1158  binf=(1-ssdsbint) ! Banner et al 2002: Hp=4*sqrt(int_0.7^1.3fp E df), SSDSBINT=0.3
1159  bsup=(1+ssdsbint)
1160  kik=0.
1161  !
1162  ! High frequency tail for convolution calculation
1163  !
1164  ALLOCATE(k1(nk,ndtab))
1165  ALLOCATE(k2(nk,ndtab))
1166  ALLOCATE(sigtab(nk,ndtab))
1167 
1168  sigtab=0. !contains frequency for upper windows boundaries
1169  iktab=0 ! contains indices for upper windows boundaries
1170 
1171  DO id=1,ndtab
1172  toto=0
1173  prof=real(id)
1174  DO ikl=1,nk ! last window starts at IK=NK
1175  CALL wavnu2(sig(ikl), prof, kik, cgg, 1e-7, 15, icon)
1176  k1(ikl,id)=kik ! wavenumber lower boundary (is directly related to the frequency indices, IK)
1177  k2(ikl,id)=((bsup/binf)**2.)*k1(ikl,id)! wavenumber upper boundary
1178  sigtab(ikl,id)=sqrt(grav*k2(ikl,id)*tanh(k2(ikl,id)*id)) ! corresponding frequency upper boundary
1179  IF(sigtab(ikl,id) .LE. sig(1)) THEN
1180  iktab(ikl,id)=1
1181  END IF
1182  IF(sigtab(ikl,id) .GT. sig(nk)) THEN
1183  iktab(ikl,id)=nk+toto ! in w3sds4 only windows with IKSUP<=NK will be kept
1184  toto=1
1185  END IF
1186  DO ik=1,nk-1
1187  diff1=0.
1188  diff2=0.
1189  IF(sig(ik)<sigtab(ikl,id) .AND. sig(ik+1)>=sigtab(ikl,id)) THEN
1190  diff1=sigtab(ikl,id)-sig(ik) ! seeks the indices of the upper boundary
1191  diff2=sig(ik+1)-sigtab(ikl,id)! the indices of lower boudary = IK
1192  IF (diff1<diff2) THEN
1193  iktab(ikl,id)=ik
1194  ELSE
1195  iktab(ikl,id)=ik+1
1196  END IF
1197  END IF
1198  END DO
1199  END DO
1200  END DO
1201  !
1202  ! Tabulates DCKI and QBI
1203  !
1204  dhs=khsmax/nkhs ! max value of KHS=KHSMAX
1205  dkh=khmax/nkhi ! max value of KH=KHMAX
1206  dkd=kdmax/nkd
1207  ALLOCATE(dcki(nkhs,nkd))
1208  ALLOCATE(qbi(nkhs,nkd))
1209  dcki=0.
1210  qbi =0.
1211  DO ikd=1,nkd
1212  khs=0.
1213  kd=(fac_kd1**(ikd-fac_kd2))
1214  xt=tanh(kd)
1215  gam=1.0314*(xt**3)-1.9958*(xt**2)+1.5522*xt+0.1885
1216  gam=gam/2.15
1217  DO ikhs=1,nkhs ! max value of KHS=1.
1218  kh=0.
1219  khs=khs+dhs
1220  DO ikh=1,nkhi
1221  kh=kh+dkh
1222  pr=(4.*kh/(khs**2.))*exp(-(2*((kh/khs)**2.)))
1223  ! W=1.5*(((KHS)/(SQRT(2.)*GAM*XT))**2.)*(1-exp(-(((KH)/(GAM*XT))**4.))) !CK2002 parameterization
1224  w=ssdsabk*(((khs)/(sqrt(2.)*gam*xt))**2.)*(1-exp(-(((kh)/(gam*xt))**ssdspbk)))
1225  eps=-((((ssdsbck/(xt**ssdshck))*kh)**3.)/4)*sqrt(grav/xt)
1226  dcki(ikhs, ikd)= dcki(ikhs, ikd)+pr*w*eps*dkh
1227  qbi(ikhs, ikd) = qbi(ikhs, ikd) +pr*w* dkh
1228  END DO
1229  END DO
1230  END DO
1231 
1232  WHERE ( qbi .GT. 1. )
1233  qbi = 1.
1234  END WHERE
1235 
1236  DEALLOCATE(k1,k2)
1237  DEALLOCATE(sigtab)
1238  ELSE
1239  iktab(:,:)=1
1240  dcki(:,:) =0.
1241  qbi(:,:) =0.
1242  END IF
1243  !
1244  !/ ------------------------------------------------------------------- /
1245  ! CUMULATIVE EFFECT
1246  !/ ------------------------------------------------------------------- /
1247  !
1248  ! Precomputes the weights for the cumulative effect (TEST 441 and 500)
1249  !
1250  dikcumul = 0
1251  IF (ssdsc(3).LT.0.) THEN
1252  ! DIKCUMUL is the integer difference in frequency bands
1253  ! between the "large breakers" and short "wiped-out waves"
1254  dikcumul = nint(ssdsbrf1/(xfr-1.))
1255  ! WRITE(6,*) 'INSIN4b:',DIKCUMUL
1256  cumulw(:,:)=0.
1257  DO ik=1,nk
1258  c = grav/sig(ik) ! Valid in deep water only
1259  !C = SIG(IK)/K(IK) ! Valid in all water depth ???
1260  DO ith=1,nth
1261  is=ith+(ik-1)*nth
1262  DO ik2=1,ik-dikcumul
1263  c2 = grav/sig(ik2) ! Valid in deep water only
1264  !C2 = SIG(IK2)/K(IK2) ! Valid in all water depth ???
1265  DO ith2=1,nth
1266  is2=ith2+(ik2-1)*nth
1267  cumulw(is2,is)=sqrt(c**2+c2**2-2*c*c2*ecos(1+abs(ith2-ith))) & ! = deltaC
1268  *dsip(ik2)/(0.5*c2) * dth ! = dk*dtheta (Valid in deep water only)
1269  END DO
1270  END DO
1271  END DO
1272  END DO
1273  ELSE
1274  cumulw(:,:)=0.
1275  END IF
1276  !/
1277  !/ End of INSIN4 ----------------------------------------------------- /
1278  !/

References w3gdatmd::cumulw, w3gdatmd::dcki, dikcumul, w3gdatmd::dsip, w3gdatmd::dth, w3gdatmd::ecos, w3gdatmd::esin, w3servmd::extcde(), fac_kd1, fac_kd2, constants::grav, ialpha, w3gdatmd::iktab, ilevtail, itaumax, iustar, jumax, kdmax, khmax, khsmax, w3odatmd::ndse, w3gdatmd::ndtab, w3gdatmd::nk, w3gdatmd::nkd, nkhi, w3gdatmd::nkhs, w3gdatmd::nth, w3gdatmd::qbi, constants::rade, w3gdatmd::satindices, w3gdatmd::satweights, w3gdatmd::sig, w3gdatmd::sintailpar, w3gdatmd::ssdsabk, w3gdatmd::ssdsbck, w3gdatmd::ssdsbint, w3gdatmd::ssdsbrf1, w3gdatmd::ssdsc, w3gdatmd::ssdscos, w3gdatmd::ssdsdth, w3gdatmd::ssdshck, w3gdatmd::ssdspbk, w3servmd::strace(), tabu_stress(), tabu_tauhf(), tabu_tauhf2(), tauhft, tauhft2, taut, w3gdatmd::th, constants::tpiinv, w3gdatmd::ttauwshelter, w3dispmd::wavnu2(), and w3gdatmd::xfr.

Referenced by w3iogrmd::w3iogr().

◆ tabu_stress()

subroutine w3src4md::tabu_stress

To generate friction velocity table TAUT(TAUW,U10)=SQRT(TAU).

Author
F. Ardhuin
Date
17-Oct-2007

Definition at line 1289 of file w3src4md.F90.

1289  !/
1290  !/ +-----------------------------------+
1291  !/ | WAVEWATCH III NOAA/NCEP |
1292  !/ | F. Ardhuin |
1293  !/ | FORTRAN 90 |
1294  !/ | Last update : 17-Oct-2007 |
1295  !/ +-----------------------------------+
1296  !/
1297  !/ 23-Jun-2006 : Origination. ( version 3.13 )
1298  !/ adapted from WAM, original:P.A.E.M. JANSSEN KNMI AUGUST 1990
1299  !/ adapted version (subr. STRESS): J. BIDLOT ECMWF OCTOBER 2004
1300  !/ Table values were checkes against the original f90 result and found to
1301  !/ be identical (at least at 0.001 m/s accuracy)
1302  !/
1303  ! 1. Purpose :
1304  ! TO GENERATE friction velocity table TAUT(TAUW,U10)=SQRT(TAU).
1305  ! METHOD.
1306  ! A STEADY STATE WIND PROFILE IS ASSUMED.
1307  ! THE WIND STRESS IS COMPUTED USING THE ROUGHNESSLENGTH
1308  ! Z1=Z0/SQRT(1-TAUW/TAU)
1309  ! WHERE Z0 IS THE CHARNOCK RELATION , TAUW IS THE WAVE-
1310  ! INDUCED STRESS AND TAU IS THE TOTAL STRESS.
1311  ! WE SEARCH FOR STEADY-STATE SOLUTIONS FOR WHICH TAUW/TAU < 1.
1312  ! FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990.
1313  !
1314  ! Initialization for source term routine.
1315  !
1316  ! 2. Method :
1317  !
1318  ! 3. Parameters :
1319  !
1320  ! Parameter list
1321  ! ----------------------------------------------------------------
1322  ! ----------------------------------------------------------------
1323  !
1324  ! 4. Subroutines used :
1325  !
1326  ! Name Type Module Description
1327  ! ----------------------------------------------------------------
1328  ! STRACE Subr. W3SERVMD Subroutine tracing.
1329  ! ----------------------------------------------------------------
1330  !
1331  ! 5. Called by :
1332  !
1333  ! Name Type Module Description
1334  ! ----------------------------------------------------------------
1335  ! W3SIN3 Subr. W3SRC3MD Corresponding source term.
1336  ! ----------------------------------------------------------------
1337  !
1338  ! 6. Error messages :
1339  !
1340  ! None.
1341  !
1342  ! 7. Remarks :
1343  !
1344  ! 8. Structure :
1345  !
1346  ! See source code.
1347  !
1348  ! 9. Switches :
1349  !
1350  ! !/S Enable subroutine tracing.
1351  !
1352  ! 10. Source code :
1353  !
1354  !/ ------------------------------------------------------------------- /
1355  USE constants, ONLY: kappa, grav
1356  USE w3gdatmd, ONLY: zzwnd, aalpha, zz0max
1357  IMPLICIT NONE
1358  INTEGER, PARAMETER :: NITER=10
1359  REAL , PARAMETER :: XM=0.50, eps1=0.00001
1360  ! VARIABLE. TYPE. PURPOSE.
1361  ! *XM* REAL POWER OF TAUW/TAU IN ROUGHNESS LENGTH.
1362  ! *XNU* REAL KINEMATIC VISCOSITY OF AIR.
1363  ! *NITER* INTEGER NUMBER OF ITERATIONS TO OBTAIN TOTAL STRESS
1364  ! *EPS1* REAL SMALL NUMBER TO MAKE SURE THAT A SOLUTION
1365  ! IS OBTAINED IN ITERATION WITH TAU>TAUW.
1366  ! ----------------------------------------------------------------------
1367  INTEGER I,J,ITER
1368  REAL ZTAUW,UTOP,CDRAG,WCD,USTOLD,TAUOLD
1369  REAL X,UST,ZZ0,F,DELF,ZZ00
1370  !
1371  !
1372  delu = umax/float(jumax)
1373  deltauw = tauwmax/float(itaumax)
1374  DO i=0,itaumax
1375  ztauw = (real(i)*deltauw)**2
1376  DO j=0,jumax
1377  utop = float(j)*delu
1378  cdrag = 0.0012875
1379  wcd = sqrt(cdrag)
1380  ustold = utop*wcd
1381  tauold = max(ustold**2, ztauw+eps1)
1382  DO iter=1,niter
1383  x = ztauw/tauold
1384  ust = sqrt(tauold)
1385  zz00=aalpha*tauold/grav
1386  IF (zz0max.NE.0) zz00=min(zz00,zz0max)
1387  ! Corrects roughness ZZ00 for quasi-linear effect
1388  zz0 = zz00/(1.-x)**xm
1389  !ZNU = 0.1*nu_air/UST ! This was removed by Bidlot in 1996
1390  !ZZ0 = MAX(ZNU,ZZ0)
1391  f = ust-kappa*utop/(alog(zzwnd/zz0))
1392  delf= 1.-kappa*utop/(alog(zzwnd/zz0))**2*2./ust &
1393  *(1.-(xm+1)*x)/(1.-x)
1394  ust = ust-f/delf
1395  tauold= max(ust**2., ztauw+eps1)
1396  END DO
1397  taut(i,j) = sqrt(tauold)
1398  END DO
1399  END DO
1400  i=itaumax
1401  j=jumax
1402  !
1403  ! Force zero wind to have zero stress (Bidlot 1996)
1404  !
1405  DO i=0,itaumax
1406  taut(i,0)=0.0
1407  END DO
1408  RETURN

References w3gdatmd::aalpha, deltauw, delu, constants::grav, itaumax, jumax, constants::kappa, taut, tauwmax, umax, w3gdatmd::zz0max, and w3gdatmd::zzwnd.

Referenced by insin4().

◆ tabu_tauhf()

subroutine w3src4md::tabu_tauhf ( real, intent(in)  SIGMAX)

Tabulation of the high-frequency wave-supported stress.

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.

Parameters
[in]SIGMAXMaximum frequency * TPI.
Author
F. Ardhuin
Date
14-Aug-2006

Definition at line 1426 of file w3src4md.F90.

1426  !/
1427  !/ +-----------------------------------+
1428  !/ | WAVEWATCH III NOAA/NCEP |
1429  !/ | F. Ardhuin |
1430  !/ | FORTRAN 90 |
1431  !/ | Last update 2006/08/14 |
1432  !/ +-----------------------------------+
1433  !/
1434  !/ 27-Feb-2004 : Origination in WW3 ( version 2.22.SHOM )
1435  !/ the resulting table was checked to be identical to the original f77 result
1436  !/ 14-Aug-2006 : Modified following Bidlot ( version 2.22.SHOM )
1437  !/ 18-Aug-2006 : Ported to version 3.09
1438  !
1439  ! 1. Purpose :
1440  !
1441  ! Tabulation of the high-frequency wave-supported stress
1442  !
1443  ! 2. Method :
1444  !
1445  ! SEE REFERENCE FOR WAVE STRESS CALCULATION.
1446  ! FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990.
1447  ! See tech. Memo ECMWF 03 december 2003 by Bidlot & Janssen
1448  !
1449  ! 3. Parameters :
1450  !
1451  ! Parameter list
1452  ! ----------------------------------------------------------------
1453  ! SIGMAX Real I maximum frequency * TPI
1454  ! ----------------------------------------------------------------
1455  !
1456  ! 4. Subroutines used :
1457  !
1458  ! STRACE Service routine.
1459  !
1460  ! 5. Called by :
1461  !
1462  ! W3SIN3 Wind input Source term routine.
1463  !
1464  ! 6. Error messages :
1465  !
1466  ! 7. Remarks :
1467  !
1468  ! 8. Structure :
1469  !
1470  ! See source code.
1471  !
1472  ! 9. Switches :
1473  !
1474  ! !/S Enable subroutine tracing.
1475  ! !/T Enable test output.
1476  !
1477  ! 10. Source code :
1478  !
1479  !/ ------------------------------------------------------------------- /
1480  USE constants, ONLY: kappa, grav
1481 #ifdef W3_S
1482  USE w3servmd, ONLY: strace
1483 #endif
1484  USE w3gdatmd, ONLY: aalpha, bbeta, zzalp, fachfe, zz0max
1485 #ifdef W3_T
1486  USE w3odatmd, ONLY: ndst
1487 #endif
1488  !
1489  IMPLICIT NONE
1490  !/
1491  !/ ------------------------------------------------------------------- /
1492  !/ Parameter list
1493  !/
1494  REAL, intent(in) :: SIGMAX ! maximum frequency
1495  !/
1496  !/ ------------------------------------------------------------------- /
1497  !/ Local parameters
1498  !/
1499  ! USTARM R.A. Maximum friction velocity
1500  ! ALPHAM R.A. Maximum Charnock Coefficient
1501  ! WLV R.A. Water levels.
1502  ! UA R.A. Absolute wind speeds.
1503  ! UD R.A. Absolute wind direction.
1504  ! U10 R.A. Wind speed used.
1505  ! U10D R.A. Wind direction used.
1506  ! 10. Source code :
1507  !
1508  !/ ------------------------------------------------------------------- /
1509 #ifdef W3_S
1510  INTEGER, SAVE :: IENT = 0
1511 #endif
1512  REAL :: USTARM, ALPHAM
1513  REAL :: CONST1, OMEGA, OMEGAC
1514  REAL :: UST, ZZ0,OMEGACC, CM
1515  INTEGER, PARAMETER :: JTOT=250
1516  REAL, ALLOCATABLE :: W(:)
1517  REAL :: ZX,ZARG,ZMU,ZLOG,ZZ00,ZBETA
1518  REAL :: Y,YC,DELY
1519  INTEGER :: J,K,L
1520  REAL :: X0
1521  !
1522 #ifdef W3_S
1523  CALL strace (ient, 'TABU_HF')
1524 #endif
1525  !
1526  ustarm = 5.
1527  alpham = 20.*aalpha
1528  delust = ustarm/real(iustar)
1529  delalp = alpham/real(ialpha)
1530  const1 = bbeta/kappa**2
1531  omegac = sigmax
1532  !
1533  tauhft(0:iustar,0:ialpha)=0. !table initialization
1534  !
1535  ALLOCATE(w(jtot))
1536  w(2:jtot-1)=1.
1537  w(1)=0.5
1538  w(jtot)=0.5
1539  x0 = 0.05
1540  !
1541  DO l=0,ialpha
1542  DO k=0,iustar
1543  ust = max(real(k)*delust,0.000001)
1544  zz00 = ust**2*aalpha/grav
1545  IF (zz0max.NE.0) zz00=min(zz00,zz0max)
1546  zz0 = zz00*(1+float(l)*delalp/aalpha)
1547  omegacc = max(omegac,x0*grav/ust)
1548  yc = omegacc*sqrt(zz0/grav)
1549  dely = max((1.-yc)/real(jtot),0.)
1550  ! For a given value of UST and ALPHA,
1551  ! the wave-supported stress is integrated all the way
1552  ! to 0.05*g/UST
1553  DO j=1,jtot
1554  y = yc+real(j-1)*dely
1555  omega = y*sqrt(grav/zz0)
1556  ! This is the deep water phase speed
1557  cm = grav/omega
1558  !this is the inverse wave age, shifted by ZZALP (tuning)
1559  zx = ust/cm +zzalp
1560  zarg = min(kappa/zx,20.)
1561  zmu = min(grav*zz0/cm**2*exp(zarg),1.)
1562  zlog = min(alog(zmu),0.)
1563  zbeta = const1*zmu*zlog**4
1564  ! Power of Y in denominator should be FACHFE-4 tail applied here
1565  tauhft(k,l) = tauhft(k,l)+w(j)*zbeta/y*dely
1566  END DO
1567 #ifdef W3_T
1568  WRITE (ndst,9000) l,k,aalpha+float(l)*delalp,ust,tauhft(k,l)
1569 #endif
1570  END DO
1571  END DO
1572  DEALLOCATE(w)
1573  RETURN
1574 #ifdef W3_T
1575 9000 FORMAT ('TABU_HF, L, K, ALPHA, UST, TAUHFT(K,L) :',(2i4,3f8.3))
1576 #endif

References w3gdatmd::aalpha, w3gdatmd::bbeta, delalp, delust, w3gdatmd::fachfe, constants::grav, ialpha, iustar, constants::kappa, w3odatmd::ndst, w3servmd::strace(), tauhft, w3gdatmd::zz0max, and w3gdatmd::zzalp.

Referenced by insin4().

◆ tabu_tauhf2()

subroutine w3src4md::tabu_tauhf2 ( real, intent(in)  SIGMAX)

Tabulation of the high-frequency wave-supported stress as a function of ustar, alpha (modified Charnock), and tail energy level.

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

Parameters
[in]SIGMAXMaximum frequency*TPI.
Author
F. Ardhuin
Date
24-Jan-2013

Definition at line 1596 of file w3src4md.F90.

1596  !/
1597  !/ +-----------------------------------+
1598  !/ | WAVEWATCH III NOAA/NCEP |
1599  !/ | F. Ardhuin |
1600  !/ | FORTRAN 90 |
1601  !/ | Last update 2006/08/14 |
1602  !/ | Last update 2013/01/24 |
1603  !/ +-----------------------------------+
1604  !/
1605  !/ 15-May-2007 : Origination in WW3 ( version 3.10.SHOM )
1606  !/ 24-Jan-2013 : Allows to read in table ( version 4.08 )
1607  !
1608  ! 1. Purpose :
1609  !
1610  ! Tabulation of the high-frequency wave-supported stress as a function of
1611  ! ustar, alpha (modified Charnock), and tail energy level
1612  !
1613  ! 2. Method :
1614  !
1615  ! SEE REFERENCE FOR WAVE STRESS CALCULATION.
1616  ! FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990.
1617  ! See tech. Memo ECMWF 03 december 2003 by Bidlot & Janssen
1618  !
1619  ! 3. Parameters :
1620  !
1621  ! Parameter list
1622  ! ----------------------------------------------------------------
1623  ! SIGMAX Real I maximum frequency*TPI
1624  ! ----------------------------------------------------------------
1625  !
1626  ! 4. Subroutines used :
1627  !
1628  ! STRACE Service routine.
1629  !
1630  ! 5. Called by :
1631  !
1632  ! W3SIN3 Wind input Source term routine.
1633  !
1634  ! 6. Error messages :
1635  !
1636  ! 7. Remarks :
1637  !
1638  ! 8. Structure :
1639  !
1640  ! See source code.
1641  !
1642  ! 9. Switches :
1643  !
1644  ! !/S Enable subroutine tracing.
1645  ! !/T Enable test output.
1646  !
1647  ! 10. Source code :
1648  !
1649  !/ ------------------------------------------------------------------- /
1650  USE constants, ONLY: kappa, grav, file_endian
1651 #ifdef W3_S
1652  USE w3servmd, ONLY: strace
1653 #endif
1654  USE w3gdatmd, ONLY: aalpha, bbeta, zzalp, fachfe, &
1655  ttauwshelter, zz0max
1656  USE w3odatmd, ONLY: ndse
1657 #ifdef W3_T
1658  USE w3odatmd, ONLY: ndst
1659 #endif
1660  !
1661  IMPLICIT NONE
1662  !/
1663  !/ ------------------------------------------------------------------- /
1664  !/ Parameter list
1665  !/
1666  REAL, intent(in) :: SIGMAX ! maximum frequency * TPI
1667  !/
1668  !/ ------------------------------------------------------------------- /
1669  !/ Local parameters
1670  !/
1671  ! USTARM R.A. Maximum friction velocity
1672  ! ALPHAM R.A. Maximum Charnock Coefficient
1673  ! WLV R.A. Water levels.
1674  ! UA R.A. Absolute wind speeds.
1675  ! UD R.A. Absolute wind direction.
1676  ! U10 R.A. Wind speed used.
1677  ! U10D R.A. Wind direction used.
1678  ! 10. Source code :
1679  !
1680  !/ ------------------------------------------------------------------- /
1681 #ifdef W3_S
1682  INTEGER, SAVE :: IENT = 0
1683 #endif
1684  REAL :: USTARM, ALPHAM, LEVTAILM
1685  REAL :: CONST1, OMEGA, OMEGAC, LEVTAIL
1686  REAL :: UST, UST0, ZZ0,OMEGACC, CM
1687  REAL :: TAUW, TAUW0
1688  INTEGER, PARAMETER :: JTOT=250
1689  REAL, ALLOCATABLE :: W(:)
1690  REAL :: ZX,ZARG,ZMU,ZLOG,ZBETA
1691  REAL :: Y,YC,DELY
1692  INTEGER :: I, J, K, L
1693  REAL :: X0, INSIGMAX, INAALPHA, INBBETA, INZZALP, INKAPPA, INGRAV
1694  INTEGER :: INIUSTAR, INIALPHA, INILEVTAIL, IERR
1695  CHARACTER(160) :: FNAMETAB
1696  LOGICAL :: NOFILE
1697  CHARACTER(LEN=10), PARAMETER :: VERGRD = '2018-06-08'
1698  CHARACTER(LEN=35), PARAMETER :: IDSTR = 'WAVEWATCH III ST4 TABLE FOR STRESS '
1699  CHARACTER(LEN=10) :: VERTST=' '
1700  CHARACTER(LEN=35) :: IDTST=' '
1701  !
1702 #ifdef W3_S
1703  CALL strace (ient, 'TABU_HF')
1704 #endif
1705  !
1706  fnametab='ST4TABUHF2.bin'
1707  nofile=.true.
1708  OPEN (993,file=fnametab,form='UNFORMATTED', convert=file_endian,iostat=ierr,status='OLD')
1709  IF (ierr.EQ.0) THEN
1710  READ(993,iostat=ierr) idtst, vertst, insigmax, inaalpha, inbbeta, iniustar, &
1711  inialpha, inilevtail, inzzalp, inkappa, ingrav
1712  IF (vertst.EQ.vergrd.AND.idtst.EQ.idstr.AND.ierr.EQ.0 &
1713  .AND.insigmax.EQ.sigmax.AND.inaalpha.EQ.aalpha.AND.inbbeta.EQ.bbeta) THEN
1714  IF (iniustar.EQ.iustar.AND.inialpha.EQ.ialpha.AND.inilevtail.EQ.ilevtail.AND. &
1715  inzzalp.EQ.zzalp.AND.ingrav.EQ.grav.AND.inkappa.EQ.kappa) THEN
1716  nofile=.false.
1717  ELSE
1718  CLOSE(993)
1719  END IF
1720  END IF
1721  END IF
1722  !
1723  ustarm = 5.
1724  alpham = 20.*aalpha
1725  levtailm = 0.05
1726  delust = ustarm/real(iustar)
1727  delalp = alpham/real(ialpha)
1728  deltail = alpham/real(ilevtail)
1729  const1 = bbeta/kappa**2
1730  omegac = sigmax
1731 800 CONTINUE
1732  IF ( nofile ) THEN
1733  WRITE(ndse,*) 'Filling 3D look-up table for SIN4. please wait'
1734  WRITE(ndse,*) idstr, vergrd, sigmax, aalpha, bbeta, iustar, ialpha, &
1735  ilevtail, zzalp, kappa, grav
1736  !
1737  tauhft(0:iustar,0:ialpha)=0. !table initialization
1738  !
1739  ALLOCATE(w(jtot))
1740  w(2:jtot-1)=1.
1741  w(1)=0.5
1742  w(jtot)=0.5
1743  x0 = 0.05
1744  !
1745  DO k=0,iustar
1746  ust0 = max(real(k)*delust,0.000001)
1747  DO l=0,ialpha
1748  ust=ust0
1749  zz0 = ust0**2*(aalpha+float(l)*delalp)/grav
1750  omegacc = max(omegac,x0*grav/ust)
1751  yc = omegacc*sqrt(zz0/grav)
1752  dely = max((1.-yc)/real(jtot),0.)
1753  ! For a given value of UST and ALPHA,
1754  ! the wave-supported stress is integrated all the way
1755  ! to 0.05*g/UST
1756  DO i=0,ilevtail
1757  levtail=real(i)*deltail
1758  tauhft(k,l)=0.
1759  tauhft2(k,l,i)=0.
1760  tauw0=ust0**2
1761  tauw=tauw0
1762  DO j=1,jtot
1763  y = yc+real(j-1)*dely
1764  omega = y*sqrt(grav/zz0)
1765  ! This is the deep water phase speed
1766  cm = grav/omega
1767  !this is the inverse wave age, shifted by ZZALP (tuning)
1768  zx = ust0/cm +zzalp
1769  zarg = min(kappa/zx,20.)
1770  zmu = min(grav*zz0/cm**2*exp(zarg),1.)
1771  zlog = min(alog(zmu),0.)
1772  zbeta = const1*zmu*zlog**4
1773  ! Power of Y in denominator should be FACHFE-4
1774  tauhft(k,l) = tauhft(k,l)+w(j)*zbeta/y*dely
1775  zx = ust/cm +zzalp
1776  zarg = min(kappa/zx,20.)
1777  zmu = min(grav*zz0/cm**2*exp(zarg),1.)
1778  zlog = min(alog(zmu),0.)
1779  zbeta = const1*zmu*zlog**4
1780  ! Power of Y in denominator should be FACHFE-4
1781  tauhft2(k,l,i) = tauhft2(k,l,i)+w(j)*zbeta*(ust/ust0)**2/y*dely
1782  tauw=tauw-w(j)*ust**2*zbeta*levtail/y*dely
1783  ust=sqrt(max(tauw,0.))
1784  END DO
1785 #ifdef W3_T
1786  WRITE (ndst,9000) k,l,i,ust0,aalpha+float(l)*delalp,levtail,tauhft2(k,l,i)
1787 #endif
1788  END DO
1789  END DO
1790  END DO
1791  DEALLOCATE(w)
1792  OPEN (993,file=fnametab,form='UNFORMATTED', convert=file_endian,iostat=ierr,status='UNKNOWN')
1793  WRITE(993) idstr, vergrd, sigmax, aalpha, bbeta, iustar, ialpha, ilevtail, zzalp, kappa, grav
1794  WRITE(993) tauhft(0:iustar,0:ialpha)
1795  WRITE(993) tauhft2
1796  CLOSE(993)
1797  !DO K=0,IUSTAR
1798  ! DO L=0,IALPHA
1799  ! DO I=0,ILEVTAIL
1800  ! WRITE(995,*) K,L,I,MAX(REAL(K)*DELUST,0.000001),AALPHA+FLOAT(L)*DELALP,REAL(I)*DELTAIL,TAUHFT(K,L),TAUHFT2(K,L,I)
1801  ! END DO
1802  ! END DO
1803  ! END DO
1804  !
1805  ELSE
1806  WRITE(ndse,*) 'Reading 3D look-up table for SIN4 from file.'
1807  READ(993,err=2000,iostat=ierr ) tauhft(0:iustar,0:ialpha)
1808  READ(993,err=2000,iostat=ierr ) tauhft2
1809  CLOSE(993)
1810  END IF
1811  !
1812  GOTO 2001
1813 2000 nofile=.true.
1814  GOTO 800
1815 2001 CONTINUE
1816  RETURN
1817 #ifdef W3_T
1818 9000 FORMAT (' TEST TABU_HFT2, K, L, I, UST, ALPHA, LEVTAIL, TAUHFT2(K,L,I) :',(3i4,4f10.5))
1819 #endif

References w3gdatmd::aalpha, w3gdatmd::bbeta, delalp, deltail, delust, w3gdatmd::fachfe, file(), constants::file_endian, constants::grav, ialpha, ilevtail, iustar, constants::kappa, w3odatmd::ndse, w3odatmd::ndst, w3servmd::strace(), tauhft, tauhft2, w3gdatmd::ttauwshelter, w3gdatmd::zz0max, and w3gdatmd::zzalp.

Referenced by insin4().

◆ w3sds4()

subroutine w3src4md::w3sds4 ( real, dimension(nspec), intent(in)  A,
real, dimension(nk), intent(in)  K,
real, dimension(nk), intent(in)  CG,
real, intent(in)  USTAR,
real, intent(in)  USDIR,
real, intent(in)  DEPTH,
real, intent(in)  DAIR,
real, dimension(nspec), intent(out)  SRHS,
real, dimension(nspec), intent(out)  DDIAG,
integer, intent(in), optional  IX,
integer, intent(in), optional  IY,
real, dimension(nspec), intent(out)  BRLAMBDA,
real, dimension(1:4), intent(out)  WHITECAP,
real, intent(in)  DLWMEAN 
)

Calculate whitecapping source term and diagonal term of derivative.

This codes does either one or the other of Ardhuin et al. (JPO 2010) Filipot & Ardhuin (JGR 2012) Romero (GRL 2009) the choice depends on SDSBCHOICE

Parameters
[in]AAction density spectrum (1-D).
[in]KWavenumber for entire spectrum.
[in]CGGroup velocity.
[in]USTARFriction velocity.
[in]USDIRWind stress direction.
[in]DEPTHWater depth.
[in]DAIRAir density.
[out]SRHS
[out]DDIAG
[in]IXGrid Index.
[in]IYGrid Index.
[out]BRLAMBDAPhillips' Lambdas.
[out]WHITECAP
[in]DLWMEAN
Author
F. Ardhuin
F. Leckler
L. Romero
Date
13-Aug-2021

Definition at line 2034 of file w3src4md.F90.

2034  !/
2035  !/ +-----------------------------------+
2036  !/ | WAVEWATCH III NOAA/NCEP |
2037  !/ ! F. Ardhuin, F. Leckler, L. Romero !
2038  !/ | FORTRAN 90 |
2039  !/ | Last update : 13-Aug-2021 |
2040  !/ +-----------------------------------+
2041  !/
2042  !/ 30-Aug-2010 : Clean up from common ST3-ST4 routine( version 3.14-Ifremer )
2043  !/ 23-Jan-2012 : Add output of lambdas to be used in SIN
2044  !/ 13-Nov-2013 : Reduced frequency range with IG1 switch
2045  !/ 06-Jun-2018 : Add optional DEBUGSRC ( version 6.04 )
2046  !/ 22-Feb-2020 : Option to use Romero (GRL 2019) ( version 7.06 )
2047  !/ 13-Aug-2021 : Consider DAIR a variable ( version 7.14 )
2048  !/ 01-Mar-2023 : Clean up of SDS4 ( version 7.xx )
2049  !/
2050  ! 1. Purpose :
2051  !
2052  ! Calculate wave dissipation source term and diagonal term of derivative.
2053  !
2054  ! 2. Method :
2055  !
2056  ! This codes does either one or the other of
2057  ! Ardhuin et al. (JPO 2010)
2058  ! Filipot & Ardhuin (JGR 2012)
2059  ! Romero (GRL 2009)
2060  ! the choice depends on SDSBCHOICE
2061  !
2062  ! 3. Parameters :
2063  !
2064  ! Parameter list
2065  ! ----------------------------------------------------------------
2066  ! IX, IY Int I Grid Index
2067  ! A R.A. I Action density spectrum (1-D).
2068  ! K R.A. I Wavenumber for entire spectrum. *)
2069  ! USTAR Real I Friction velocity.
2070  ! USDIR Real I wind stress direction.
2071  ! DEPTH Real I Water depth.
2072  ! DAIR Real I Air density
2073  ! S R.A. O Source term (1-D version).
2074  ! D R.A. O Diagonal term of derivative. *)
2075  ! BRLAMBDA R.A. O Phillips' Lambdas
2076  ! ----------------------------------------------------------------
2077  ! *) Stored in 1-D array with dimension NTH*NK
2078  !
2079  ! 4. Subroutines used :
2080  !
2081  ! STRACE Subroutine tracing. ( !/S switch )
2082  ! PRT2DS Print plot of spectrum. ( !/T0 switch )
2083  ! OUTMAT Print out matrix. ( !/T1 switch )
2084  !
2085  ! 5. Called by :
2086  !
2087  ! W3SRCE Source term integration.
2088  ! W3EXPO Point output program.
2089  ! GXEXPO GrADS point output program.
2090  !
2091  ! 6. Error messages :
2092  !
2093  ! 7. Remarks :
2094  !
2095  ! 8. Structure :
2096  !
2097  ! See source code.
2098  !
2099  ! 9. Switches :
2100  !
2101  ! !/S Enable subroutine tracing.
2102  ! !/T Enable general test output.
2103  ! !/T0 2-D print plot of source term.
2104  ! !/T1 Print arrays.
2105  !
2106  ! 10. Source code :
2107  !
2108  !/ ------------------------------------------------------------------- /
2109  USE constants,ONLY: grav, dwat, pi, tpi, rade, debug_node
2110  USE w3gdatmd, ONLY: nspec, nth, nk, ssdsbr, ssdsbt, dden, &
2111  ssdsc, ec2, es2, esc, &
2112  sig, ssdsp, ecos, esin, dth, aairgb, &
2113  ssdsiso, ssdsdth, ssdsbm, aaircmin, &
2114  ssdsbrfdf, ssdsbck, iktab, dcki, &
2115  satindices, satweights, cumulw, nkhs, nkd, &
2116  ndtab, qbi, dsip, ssdsbrf1,xfr
2117 #ifdef W3_IG1
2118  USE w3gdatmd, ONLY: igpars
2119 #endif
2120  USE w3odatmd, ONLY: flogrd
2121 #ifdef W3_S
2122  USE w3servmd, ONLY: strace
2123 #endif
2124 #ifdef W3_T
2125  USE w3odatmd, ONLY: ndst
2126 #endif
2127 #ifdef W3_T0
2128  USE w3odatmd, ONLY: ndst
2129  USE w3arrymd, ONLY: prt2ds
2130 #endif
2131 #ifdef W3_T1
2132  USE w3arrymd, ONLY: outmat
2133 #endif
2134  !
2135  IMPLICIT NONE
2136  !/
2137  !/ ------------------------------------------------------------------- /
2138  !/ Parameter list
2139  !/
2140  INTEGER, OPTIONAL, INTENT(IN) :: IX, IY
2141  REAL, INTENT(IN) :: A(NSPEC), K(NK), CG(NK), &
2142  DEPTH, DAIR, USTAR, USDIR, DLWMEAN
2143  REAL, INTENT(OUT) :: SRHS(NSPEC), DDIAG(NSPEC), BRLAMBDA(NSPEC)
2144  REAL, INTENT(OUT) :: WHITECAP(1:4)
2145  !/
2146  !/ ------------------------------------------------------------------- /
2147  !/ Local parameters
2148  !/
2149  INTEGER :: IS, IS2, IS0, IKL, IKC, ID, NKL
2150 #ifdef W3_S
2151  INTEGER, SAVE :: IENT = 0
2152 #endif
2153  INTEGER :: IK, IK1, ITH, IK2, JTH, ITH2, &
2154  IKHS, IKD, SDSNTH, IT, IKM, NKM
2155  INTEGER :: NSMOOTH(NK)
2156  REAL :: C, C2, CUMULWISO, COSWIND, ASUM, SDIAGISO
2157  REAL :: COEF1, COEF2, COEF4(NK), &
2158  COEF5(NK)
2159 
2160  REAL :: FACTURB, FACTURB2, DTURB, DVISC, DIAG2, BREAKFRACTION
2161  REAL :: RENEWALFREQ, EPSR
2162  REAL :: S1(NK), E1(NK)
2163  INTEGER :: NTIMES(NK)
2164  REAL :: GAM, XT
2165  REAL :: DK(NK), HS(NK), KBAR(NK), DCK(NK)
2166  REAL :: EFDF(NK) ! Energy integrated over a spectral band
2167  INTEGER :: IKSUP(NK)
2168  REAL :: FACSAT, DKHS, FACSTRAINB, FACSTRAINL
2169  REAL :: BTH0(NK) !saturation spectrum
2170  REAL :: BTH(NSPEC) !saturation spectrum
2171  REAL :: MSSSUM(NK,5), FACHF
2172  REAL :: MSSLONG
2173  REAL :: MSSPCS, MSSPC2, MSSPS2, MSSP, MSSD, MSSTH
2174  REAL :: MICHE, X, KLOC
2175 #ifdef W3_T0
2176  REAL :: DOUT(NK,NTH)
2177 #endif
2178  REAL :: QB(NK), S2(NK)
2179  REAL :: TSTR, TMAX, DT, T, MFT, DIRFORCUM
2180  REAL :: PB(NSPEC), PB2(NSPEC), BRM12(NK), BTOVER
2181  REAL :: KO, LMODULATION(NTH)
2182  !/
2183  !/ ------------------------------------------------------------------- /
2184  !/
2185 #ifdef W3_S
2186  CALL strace (ient, 'W3SDS4')
2187 #endif
2188  !
2189  !
2190  !----------------------------------------------------------------------
2191  !
2192  ! 0. Pre-Initialization to zero out arrays. All arrays should be reset
2193  ! within the computation, but these are helping with some bugs
2194  ! found in certain compilers
2195  nsmooth=0
2196  s1=0.; e1=0.
2197  ntimes=0;iksup=0
2198  dk=0.; hs=0.; kbar=0.; dck=0.; efdf=0.
2199  bth0=0.; bth=0.; ddiag=0.; srhs=0.; pb=0.
2200  msssum(:,:)=0.
2201 #ifdef W3_T0
2202  dout=0.
2203 #endif
2204  qb=0.; s2=0.;pb=0.; pb2=0.
2205  brm12(:)=0.
2206  !
2207  ! 1. Initialization and numerical factors
2208  !
2209  facturb=ssdsc(5)*ustar**2/grav*dair/dwat
2210  dikcumul = nint(ssdsbrf1/(xfr-1.))
2211  breakfraction=0.
2212  renewalfreq=0.
2213  ik1=1
2214 #ifdef W3_IG1
2215  ik1=nint(igpars(5))+1
2216 #endif
2217  !
2218  ! 1.b MSS parameters used for Modulation factors for lambda (Romero )
2219  !
2220  IF (ssdsc(8).GT.0.OR.ssdsc(11).GT.0.OR.ssdsc(18).GT.0) THEN
2221  DO ik=1,nk
2222  mssp = 0.
2223  msspc2 = 0.
2224  mssps2 = 0.
2225  msspcs = 0.
2226  !
2227  ! Sums the contributions to the directional MSS for all angles
2228  !
2229  DO ith=1,nth
2230  is=ith+(ik-1)*nth
2231  msslong = k(ik)**ssdsc(20) * a(is) * dden(ik) / cg(ik) ! contribution to MSS
2232  msspc2 = msspc2 +msslong*ec2(ith)
2233  mssps2 = mssps2 +msslong*es2(ith)
2234  msspcs = msspcs +msslong*esc(ith)
2235  mssp = mssp +msslong
2236  END DO
2237  msssum(ik:nk,1) = msssum(ik:nk,1) +mssp
2238  msssum(ik:nk,3) = msssum(ik:nk,3) +msspc2
2239  msssum(ik:nk,4) = msssum(ik:nk,4) +mssps2
2240  msssum(ik:nk,5) = msssum(ik:nk,5) +msspcs
2241  !
2242  ! Direction of long wave mss summed up to IK
2243  !
2244  mssd=0.5*(atan2(2*msssum(ik,5),msssum(ik,3)-msssum(ik,4)))
2245  IF (mssd.LT.0) mssd = mssd + pi
2246  msssum(ik,2) = mssd
2247  END DO
2248  END IF ! SSDSC(8).GT.0) THEN
2249  !
2250  ! 2. Estimation of spontaneous breaking from local saturation
2251  !
2252  !############################################################################################"
2253  SELECT CASE (nint(ssdsc(1)))
2254  CASE (1)
2255  !
2256  ! 2.a Case of a direction-dependent breaking term following Ardhuin et al. 2010
2257  !
2258  epsr = sqrt(ssdsbr)
2259  !
2260  ! 2.a.1 Computes saturation
2261  !
2262  bth(:) = 0.
2263 
2264  DO ik=ik1, nk
2265 
2266  facsat=sig(ik)*k(ik)**3*dth
2267  is0=(ik-1)*nth
2268  bth(is0+1)=0.
2269  asum = sum(a(is0+1:is0+nth))
2270  bth0(ik)=asum*facsat
2271  !
2272  IF (ssdsdth.GE.180) THEN ! integrates around full circle
2273  bth(is0+1:is0+nth)=bth0(ik)
2274  ELSE
2275  DO ith=1,nth ! partial integration
2276  is=ith+(ik-1)*nth
2277  bth(is)=dot_product(satweights(:,ith), a(is0+satindices(:,ith)) ) &
2278  *facsat
2279  END DO
2280 
2281  bth0(ik)=maxval(bth(is0+1:is0+nth))
2282  END IF
2283  !
2284  END DO !IK=NK
2285  !
2286  ! 2.a.2 Computes spontaneous breaking dissipation rate
2287  !
2288  DO ik=ik1, nk
2289  !
2290  ! Correction of saturation level for shallow-water kinematics
2291  !
2292  IF (ssdsbm(0).EQ.1) THEN
2293  miche=1.
2294  ELSE
2295  x=tanh(min(k(ik)*depth,10.))
2296  ! Correction of saturation threshold for shallow-water kinematics
2297  miche=(x*(ssdsbm(1)+x*(ssdsbm(2)+x*(ssdsbm(3)+x*ssdsbm(4)))))**2
2298  END IF
2299  coef1=(ssdsbr*miche)
2300  !
2301  ! Computes isotropic part
2302  !
2303  sdiagiso = ssdsc(2) * sig(ik)*ssdsc(6)*(max(0.,bth0(ik)/coef1-1.))**2
2304  !
2305  ! Computes anisotropic part and sums isotropic part
2306  !
2307  coef2=ssdsc(2) * sig(ik)*(1-ssdsc(6))/(coef1*coef1)
2308  ddiag((ik-1)*nth+1:ik*nth) = sdiagiso + &
2309  coef2*((max(0.,bth((ik-1)*nth+1:ik*nth)-coef1))**ssdsp)
2310  END DO
2311  !
2312  ! Computes Breaking probability
2313  !
2314  pb = (max(sqrt(bth)-epsr,0.))**2
2315  !
2316  ! Multiplies by 28.16 = 22.0 * 1.6² * 1/2 with
2317  ! 22.0 (Banner & al. 2000, figure 6)
2318  ! 1.6 the coefficient that transforms SQRT(B) to Banner et al. (2000)'s epsilon
2319  ! 1/2 factor to correct overestimation of Banner et al. (2000)'s breaking probability due to zero-crossing analysis
2320  !
2321  pb = pb * 28.16
2322  ! Compute Lambda = PB* l(k,th)
2323  ! with l(k,th)=1/(2*pi²)= the breaking crest density
2324  brlambda = pb / (2.*pi**2.)
2325  srhs = ddiag * a
2326 
2327  !############################################################################################"
2328  CASE(2)
2329  !
2330  ! 2.b Computes spontaneous breaking for T500 (Filipot et al. JGR 2010)
2331  !
2332  e1 = 0.
2333  hs = 0.
2334  srhs = 0.
2335  ddiag = 0.
2336  pb2 = 0.
2337  !
2338  ! Computes Wavenumber spectrum E1 integrated over direction and computes dk
2339  !
2340  DO ik=ik1, nk
2341  e1(ik)=0.
2342  DO ith=1,nth
2343  is=ith+(ik-1)*nth
2344  e1(ik)=e1(ik)+(a(is)*sig(ik))*dth
2345  END DO
2346  dk(ik)=dden(ik)/(dth*sig(ik)*cg(ik))
2347  END DO
2348  !
2349  ! Gets windows indices of IKTAB
2350  !
2351  id=min(nint(depth),ndtab)
2352  IF (id < 1) THEN
2353  id = 1
2354  ELSE IF(id > ndtab) THEN
2355  id = ndtab
2356  END IF
2357  !
2358  ! loop over wave scales
2359  !
2360  hs=0.
2361  efdf=0.
2362  kbar=0.
2363  nkl=0. !number of windows
2364  DO ikl=1,nk
2365  iksup(ikl)=iktab(ikl,id)
2366  IF (iksup(ikl) .LE. nk) THEN
2367  efdf(ikl) = dot_product(e1(ikl:iksup(ikl)-1),dk(ikl:iksup(ikl)-1))
2368  IF (efdf(ikl) .NE. 0) THEN
2369  kbar(ikl) = dot_product(k(ikl:iksup(ikl)-1)*e1(ikl:iksup(ikl)-1), &
2370  dk(ikl:iksup(ikl)-1)) / efdf(ikl)
2371  ELSE
2372  kbar(ikl)=0.
2373  END IF
2374  ! estimation of Significant wave height of a given scale
2375  hs(ikl) = 4*sqrt(efdf(ikl))
2376  nkl = nkl+1
2377  END IF
2378  END DO
2379  !
2380  ! Computes Dissipation and breaking probability in each scale
2381  !
2382  dck=0.
2383  qb =0.
2384  dkhs = khsmax/nkhs
2385  DO ikl=1, nkl
2386  IF (hs(ikl) .NE. 0. .AND. kbar(ikl) .NE. 0.) THEN
2387  ! gets indices for tabulated dissipation DCKI and breaking probability QBI
2388  !
2389  ikd = fac_kd2+anint(log(kbar(ikl)*depth)/log(fac_kd1))
2390  ikhs= 1+anint(kbar(ikl)*hs(ikl)/dkhs)
2391  IF (ikd > nkd) THEN ! Deep water
2392  ikd = nkd
2393  ELSE IF (ikd < 1) THEN ! Shallow water
2394  ikd = 1
2395  END IF
2396  IF (ikhs > nkhs) THEN
2397  ikhs = nkhs
2398  ELSE IF (ikhs < 1) THEN
2399  ikhs = 1
2400  END IF
2401  xt = tanh(kbar(ikl)*depth)
2402  !
2403  ! Gamma corrected for water depth
2404  !
2405  gam=1.0314*(xt**3)-1.9958*(xt**2)+1.5522*xt+0.1885
2406  !
2407  ! Computes the energy dissipated for the scale IKL
2408  ! using DCKI which is tabulated in INSIN4
2409  !
2410  dck(ikl)=((kbar(ikl)**(-2.5))*(kbar(ikl)/(2*pi)))*dcki(ikhs,ikd)
2411  !
2412  ! Get the breaking probability for the scale IKL
2413  !
2414  qb(ikl) = qbi(ikhs,ikd) ! QBI is tabulated in INSIN4
2415  ELSE
2416  dck(ikl)=0.
2417  qb(ikl) =0.
2418  END IF
2419  END DO
2420  !
2421  ! Distributes scale dissipation over the frequency spectrum
2422  !
2423  s1 = 0.
2424  s2 = 0.
2425  ntimes = 0
2426  DO ikl=1, nkl
2427  IF (efdf(ikl) .GT. 0.) THEN
2428  s1(ikl:iksup(ikl)) = s1(ikl:iksup(ikl)) + &
2429  dck(ikl)*e1(ikl:iksup(ikl)) / efdf(ikl)
2430  s2(ikl:iksup(ikl)) = s2(ikl:iksup(ikl)) + &
2431  qb(ikl) *e1(ikl:iksup(ikl)) / efdf(ikl)
2432  ntimes(ikl:iksup(ikl)) = ntimes(ikl:iksup(ikl)) + 1
2433  END IF
2434  END DO
2435  !
2436  ! Finish the average
2437  !
2438  WHERE (ntimes .GT. 0)
2439  s1 = s1 / ntimes
2440  s2 = s2 / ntimes
2441  ELSEWHERE
2442  s1 = 0.
2443  s2 = 0.
2444  END WHERE
2445  ! goes back to action for dissipation source term
2446  s1(1:nk) = s1(1:nk) / sig(1:nk)
2447  !
2448  ! Makes Isotropic distribution
2449  !
2450  asum = 0.
2451  DO ik = 1, nk
2452  asum = (sum(a(((ik-1)*nth+1):(ik*nth)))*dth)
2453  IF (asum.GT.1.e-8) THEN
2454  FORALL (is=1+(ik-1)*nth:ik*nth) ddiag(is) = s1(ik)/asum
2455  FORALL (is=1+(ik-1)*nth:ik*nth) pb2(is) = s2(ik)/asum
2456  ELSE
2457  FORALL (is=1+(ik-1)*nth:ik*nth) ddiag(is) = 0.
2458  FORALL (is=1+(ik-1)*nth:ik*nth) pb2(is) = 0.
2459  END IF
2460  IF (pb2(1+(ik-1)*nth).GT.0.001) THEN
2461  bth0(ik) = 2.*ssdsbr
2462  ELSE
2463  bth0(ik) = 0.
2464  END IF
2465  END DO
2466  !
2467  pb = (1-ssdsc(1))*pb2*a + ssdsc(1)*pb
2468  ! Compute Lambda = PB* l(k,th)
2469  ! with l(k,th)=1/(2*pi²)= the breaking crest density
2470  brlambda = pb / (2.*pi**2.)
2471  srhs = ddiag * a
2472  !############################################################################################"
2473  CASE(3)
2474  !
2475  ! 2c Romero (GRL 2019)
2476  !
2477  ! directional saturation I
2478  ! integrate in azimuth
2479  ko=(grav/(1e-6+ustar**2))/(28./ssdsc(16))**2
2480  DO ik=1,nk
2481  is0=(ik-1)*nth
2482  kloc=k(ik)**(2-ssdsc(20)) ! local wavenumber factor, if mss not used.
2483  bth(1:nth)=max(a(is0+1:is0+nth)*sig(ik)*k(ik)**3,.00000000000001)
2484  !
2485  dirforcum=dlwmean
2486  IF (ssdsc(11).GT.0) dirforcum=msssum(ik,2)
2487 
2488  c=sig(ik)/k(ik)
2489  bth0(ik)=sum(bth(1:nth)*dth)
2490  IF (ssdsc(18).GT.0) THEN ! Applies modulation factor on Lambda
2491  DO ith=1,nth
2492  facstrainl=1.+ssdsc(18)*((msssum(ik,1)*kloc)**ssdsc(14) * & ! Romero
2493  (ecos(ith)*cos(dirforcum)+esin(ith)*sin(dirforcum))**2)
2494  lmodulation(ith)= facstrainl**ssdsc(19)
2495  END DO
2496  ELSE
2497  lmodulation(:)= 1.
2498  END IF
2499 
2500  ! Lambda
2501  brlambda(is0+1:is0+nth)=ssdsc(9)*exp(-ssdsbr/bth(1:nth)) &
2502  *( 1.0+ssdsc(13)*max(1.,(k(ik)/ko))**ssdsc(15) ) &
2503  /(ssdsc(13)+1)*lmodulation(1:nth)
2504  ! Breaking strength : generalisation of Duncan's b parameter
2505  btover = sqrt(bth0(ik))-sqrt(ssdsbt)
2506  brm12(ik)=ssdsc(2)*(max(0.,btover))**(2.5)/sig(ik) ! not function of direction
2507  ! For consistency set BRLAMBDA set to zero if b is zero
2508  brlambda(is0+1:is0+nth)= max(0.,sign(brlambda(is0+1:is0+nth),btover))
2509  ! Source term / sig2 (action dissipation)
2510  srhs(is0+1:is0+nth)= brm12(ik)/grav**2*brlambda(is0+1:is0+nth)*c**5
2511  ! diagonal
2512  ddiag(is0+1:is0+nth) = srhs(is0+1:is0+nth)*ssdsbr/max(1.e-20,bth(1:nth))/max(1e-20,a(is0+1:is0+nth)) !
2513  END DO
2514  ! Breaking probability (Is actually the breaking rate)
2515  pb = brlambda *c
2516  !
2517  END SELECT
2518  !############################################################################################"
2519  !
2520  !
2521  !/ ------------------------------------------------------------------- /
2522  ! WAVE-TURBULENCE INTERACTION AND CUMULATIVE EFFECT
2523  !/ ------------------------------------------------------------------- /
2524  !
2525  !
2526  ! loop over spectrum
2527  !
2528  IF ( (ssdsc(3).NE.0.) .OR. (ssdsc(5).NE.0.) .OR. (ssdsc(21).NE.0.) ) THEN
2529  DO ik=ik1, nk
2530  renewalfreq = 0.
2531  facturb2=-2.*sig(ik)*k(ik)*facturb
2532  dvisc=-4.*ssdsc(21)*k(ik)*k(ik)
2533  c = sig(ik)/k(ik) ! phase speed
2534  !
2535  IF (ssdsc(3).GT.0 .AND. ik.GT.dikcumul) THEN
2536  ! this is the cheap isotropic version
2537  DO ik2=ik1,ik-dikcumul
2538  c2 = sig(ik2)/k(ik2)
2539  is2=(ik2-1)*nth
2540  cumulwiso=abs(c2-c)*dsip(ik2)/(0.5*c2) * dth
2541  renewalfreq=renewalfreq-cumulwiso*sum(brlambda(is2+1:is2+nth))
2542  END DO
2543  END IF
2544 
2545  DO ith=1,nth
2546  is=ith+(ik-1)*nth
2547  !
2548  ! Computes cumulative effect from Breaking probability
2549  !
2550  IF (ssdsc(3).LT.0 .AND. ik.GT.dikcumul) THEN
2551  renewalfreq = 0.
2552  ! this is the expensive and largely useless version
2553  DO ik2=ik1,ik-dikcumul
2554  IF (bth0(ik2).GT.ssdsbr) THEN
2555  is2=(ik2-1)*nth
2556  renewalfreq=renewalfreq+dot_product(cumulw(is2+1:is2+nth,is),brlambda(is2+1:is2+nth))
2557  END IF
2558  END DO
2559  END IF
2560  !
2561  ! Computes wave turbulence interaction
2562  !
2563  coswind=(ecos(ith)*cos(usdir)+esin(ith)*sin(usdir))
2564  dturb=facturb2*max(0.,coswind) ! Theory -> stress direction
2565  !
2566  ! Add effects
2567  !
2568  diag2 = (ssdsc(3)*renewalfreq+dturb+dvisc)
2569  ddiag(is) = ddiag(is) + diag2
2570  srhs(is) = srhs(is) + a(is)* diag2
2571  END DO
2572  END DO
2573  END IF
2574  !
2575  ! COMPUTES WHITECAP PARAMETERS
2576  !
2577  IF ( .NOT. (flogrd(5,7).OR.flogrd(5,8) ) ) THEN
2578  RETURN
2579  END IF
2580  !
2581  whitecap(1:4) = 0.
2582  !
2583  ! precomputes integration of Lambda over direction
2584  ! times wavelength times a (a=5 in Reul&Chapron JGR 2003) times dk
2585  !
2586  DO ik=1,min(floor(aaircmin),nk)
2587  c=sig(ik)/k(ik)
2588  is0=(ik-1)*nth
2589  coef4(ik) = c*c*sum(brlambda(is0+1:is0+nth)) &
2590  *2.*pi/grav*ssdsc(7) * dden(ik)/(sig(ik)*cg(ik))
2591  coef5(ik) = c**3*sum(brlambda(is0+1:is0+nth) &
2592  *brm12(ik)) &
2593  *aairgb/grav * dden(ik)/(sig(ik)*cg(ik))
2594  ! COEF4(IK) = SUM(BRLAMBDA((IK-1)*NTH+1:IK*NTH) * DTH) *(2*PI/K(IK)) * &
2595  ! SSDSC(7) * DDEN(IK)/(DTH*SIG(IK)*CG(IK))
2596  ! NB: SSDSC(7) is WHITECAPWIDTH
2597  END DO
2598  ! Need to extrapolate above NK if necessary ... to be added later.
2599  DO ik=min(floor(aaircmin),nk),nk
2600  coef4(ik)=0.
2601  coef5(ik)=0.
2602  END DO
2603 
2604  !/
2605  IF ( flogrd(5,7) ) THEN
2606  !
2607  ! Computes the Total WhiteCap Coverage (a=5. ; Reul and Chapron, 2003)
2608  !
2609  DO ik=ik1,min(floor(aaircmin),nk)
2610  whitecap(1) = whitecap(1) + coef4(ik) * (1-whitecap(1))
2611  whitecap(4) = whitecap(4) + coef5(ik)
2612  END DO
2613  END IF
2614  !/
2615  IF ( flogrd(5,8) ) THEN
2616  !
2617  ! Calculates the Mean Foam Thickness for component K(IK) => Fig.3, Reul and Chapron, 2003
2618  ! ( Copied from ST4 - not yet tested/validated with Romero 2019 (Lambda model)
2619  !
2620  DO ik=ik1,nk
2621  ! Duration of active breaking (TAU*)
2622  tstr = 0.8 * 2*pi/sig(ik)
2623  ! Time persistence of foam (a=5.)
2624  tmax = 5. * 2*pi/sig(ik)
2625  dt = tmax / 50
2626  mft = 0.
2627  DO it = 1, 50
2628  ! integration over time of foam persistance
2629  t = float(it) * dt
2630  ! Eq. 5 and 6 of Reul and Chapron, 2003
2631  IF ( t .LT. tstr ) THEN
2632  mft = mft + 0.4 / (k(ik)*tstr) * t * dt
2633  ELSE
2634  mft = mft + 0.4 / k(ik) * exp(-1*(t-tstr)/3.8) * dt
2635  END IF
2636  END DO
2637  mft = mft / tmax
2638  !
2639  ! Computes foam-layer thickness (Reul and Chapron, 2003)
2640  !
2641  whitecap(2) = whitecap(2) + coef4(ik) * mft
2642  END DO
2643  END IF
2644  !
2645  ! End of output computing
2646  !
2647  RETURN
2648  !
2649  ! Formats
2650  !
2651  !/
2652  !/ End of W3SDS4 ----------------------------------------------------- /
2653  !/

References w3gdatmd::aaircmin, w3gdatmd::aairgb, w3gdatmd::cumulw, w3gdatmd::dcki, w3gdatmd::dden, constants::debug_node, dikcumul, w3gdatmd::dsip, w3gdatmd::dth, constants::dwat, w3gdatmd::ec2, w3gdatmd::ecos, w3gdatmd::es2, w3gdatmd::esc, w3gdatmd::esin, fac_kd1, fac_kd2, w3odatmd::flogrd, constants::grav, w3gdatmd::igpars, w3gdatmd::iktab, khsmax, w3odatmd::ndst, w3gdatmd::ndtab, w3gdatmd::nk, w3gdatmd::nkd, w3gdatmd::nkhs, w3gdatmd::nspec, w3gdatmd::nth, w3arrymd::outmat(), constants::pi, w3arrymd::prt2ds(), w3gdatmd::qbi, constants::rade, w3gdatmd::satindices, w3gdatmd::satweights, w3gdatmd::sig, w3gdatmd::ssdsbck, w3gdatmd::ssdsbm, w3gdatmd::ssdsbr, w3gdatmd::ssdsbrf1, w3gdatmd::ssdsbrfdf, w3gdatmd::ssdsbt, w3gdatmd::ssdsc, w3gdatmd::ssdsdth, w3gdatmd::ssdsiso, w3gdatmd::ssdsp, w3servmd::strace(), constants::tpi, and w3gdatmd::xfr.

Referenced by gxexpo(), w3exnc(), w3expo(), and w3srcemd::w3srce().

◆ w3sin4()

subroutine w3src4md::w3sin4 ( real, dimension(nspec), intent(in)  A,
real, dimension(nk), intent(in)  CG,
real, dimension(nspec), intent(in)  K,
real, intent(in)  U,
real, intent(in)  USTAR,
real, intent(in)  DRAT,
real, intent(in)  AS,
real, intent(in)  USDIR,
real, intent(in)  Z0,
real, intent(in)  CD,
real, intent(out)  TAUWX,
real, intent(out)  TAUWY,
real, intent(out)  TAUWNX,
real, intent(out)  TAUWNY,
real, dimension(nspec), intent(out)  S,
real, dimension(nspec), intent(out)  D,
logical, dimension(nspec), intent(out)  LLWS,
integer, intent(in)  IX,
integer, intent(in)  IY,
real, dimension(nspec), intent(in)  BRLAMBDA 
)

Calculate diagonal and input source term for WAM4+ approach. 

       WAM-4     : Janssen et al.
       WAM-"4.5" : gustiness effect (Cavaleri et al. )
       SAT       : high-frequency input reduction for balance with
                   saturation dissipation (Ardhuin et al., 2008)
       SWELL     : negative wind input (Ardhuin et al. 2008)
Parameters
[in]AAction density spectrum (1-D).
[in]CGGroup speed.
[in]KWavenumber for entire spectrum.
[in]UWind speed.
[in]USTARFriction velocity.
[in]DRATAir/water density ratio.
[in]ASAir-sea temperature difference.
[in]USDIRWind stress direction.
[in]Z0Air-sea roughness length.
[in]CDWind drag coefficient.
[out]TAUWXComponent of the wave-supported stress.
[out]TAUWYComponent of the wave-supported stress.
[out]TAUWNXComponent of the negative wave-supported stress.
[out]TAUWNYComponent of the negative wave-supported stress.
[out]SSource term (1-D version).
[out]DDiagonal term of derivative.
[out]LLWS
[in]IX
[in]IY
[in]BRLAMBDA
Author
F. Ardhuin
H. L. Tolman
Date
05-Dec-2013

Definition at line 426 of file w3src4md.F90.

426  !/
427  !/ +-----------------------------------+
428  !/ | WAVEWATCH III SHOM |
429  !/ ! F. Ardhuin !
430  !/ | H. L. Tolman |
431  !/ | FORTRAN 90 |
432  !/ | Last update : 05-Dec-2013 |
433  !/ +-----------------------------------+
434  !/
435  !/ 09-Oct-2007 : Origination. ( version 3.13 )
436  !/ 24-Jan-2013 : Adding breaking-related input ( version 4.16 )
437  !/ 05-Dec-2013 : Cleaning up the ICE input ( version 4.16 )
438  !/
439  ! 1. Purpose :
440  !
441  ! Calculate diagonal and input source term for WAM4+ approach.
442  !
443  ! 2. Method :
444  !
445  ! WAM-4 : Janssen et al.
446  ! WAM-"4.5" : gustiness effect (Cavaleri et al. )
447  ! SAT : high-frequency input reduction for balance with
448  ! saturation dissipation (Ardhuin et al., 2008)
449  ! SWELL : negative wind input (Ardhuin et al. 2008)
450  !
451  ! 3. Parameters :
452  !
453  ! Parameter list
454  ! ----------------------------------------------------------------
455  ! A R.A. I Action density spectrum (1-D).
456  ! CG R.A. I Group speed *)
457  ! K R.A. I Wavenumber for entire spectrum. *)
458  ! U Real I WIND SPEED
459  ! USTAR Real I Friction velocity.
460  ! DRAT Real I Air/water density ratio.
461  ! AS Real I Air-sea temperature difference
462  ! USDIR Real I wind stress direction
463  ! Z0 Real I Air-side roughness lengh.
464  ! CD Real I Wind drag coefficient.
465  ! USDIR Real I Direction of friction velocity
466  ! TAUWX-Y Real I Components of the wave-supported stress.
467  ! TAUWNX Real I Component of the negative wave-supported stress.
468  ! TAUWNY Real I Component of the negative wave-supported stress.
469  ! S R.A. O Source term (1-D version).
470  ! D R.A. O Diagonal term of derivative. *)
471  ! ----------------------------------------------------------------
472  ! *) Stored as 1-D array with dimension NTH*NK
473  !
474  ! 4. Subroutines used :
475  !
476  ! STRACE Subroutine tracing. ( !/S switch )
477  ! PRT2DS Print plot of spectrum. ( !/T0 switch )
478  ! OUTMAT Print out matrix. ( !/T1 switch )
479  !
480  ! 5. Called by :
481  !
482  ! W3SRCE Source term integration.
483  ! W3EXPO Point output program.
484  ! GXEXPO GrADS point output program.
485  !
486  ! 6. Error messages :
487  !
488  ! 7. Remarks :
489  !
490  ! 8. Structure :
491  !
492  ! See source code.
493  !
494  ! 9. Switches :
495  !
496  ! !/S Enable subroutine tracing.
497  ! !/T Enable general test output.
498  ! !/T0 2-D print plot of source term.
499  ! !/T1 Print arrays.
500  !
501  ! 10. Source code :
502  !
503  !/ ------------------------------------------------------------------- /
504  USE constants, ONLY: grav,nu_air,kappa,tpi,fwtable,sizefwtable, &
505 #ifdef w3_t
506  rade, &
507 #endif
508  delab,abmin
509  USE w3gdatmd, ONLY: nk, nth, nspec, dden, sig, sig2, th, &
510  esin, ecos, ec2, zzwnd, aalpha, bbeta, zzalp,&
511  ttauwshelter, sswellf, dden2, dth, ssinthp, &
512  zz0rat, ssinbr, sintailpar
513 #ifdef W3_S
514  USE w3servmd, ONLY: strace
515 #endif
516 #ifdef W3_T
517  USE w3odatmd, ONLY: ndst
518 #endif
519 #ifdef W3_T0
520  USE w3odatmd, ONLY: ndst
521 #endif
522  USE w3odatmd, ONLY: iaproc
523 #ifdef W3_T0
524  USE w3arrymd, ONLY: prt2ds
525 #endif
526 #ifdef W3_T1
527  USE w3arrymd, ONLY: outmat
528 #endif
529  !
530  IMPLICIT NONE
531  !/
532  !/ ------------------------------------------------------------------- /
533  !/ Parameter list
534  !/
535  REAL, INTENT(IN) :: A(NSPEC), BRLAMBDA(NSPEC)
536  REAL, INTENT(IN) :: CG(NK), K(NSPEC),Z0,U, CD
537  REAL, INTENT(IN) :: USTAR, USDIR, AS, DRAT
538  REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC), TAUWX, TAUWY, TAUWNX, TAUWNY
539  LOGICAL, INTENT(OUT) :: LLWS(NSPEC)
540  INTEGER, INTENT(IN) :: IX, IY
541  !/
542  !/ ------------------------------------------------------------------- /
543  !/ Local parameters
544  !/
545  INTEGER :: IS,IK,ITH
546 #ifdef W3_S
547  INTEGER, SAVE :: IENT = 0
548 #endif
549  REAL :: FACLN1, FACLN2, LAMBDA
550  REAL :: COSU, SINU, TAUX, TAUY, USDIRP, USTP
551  REAL :: TAUPX, TAUPY, UST2, TAUW, TAUWB
552  REAL , PARAMETER :: EPS1 = 0.00001, eps2 = 0.000001
553  REAL :: Usigma !standard deviation of U due to gustiness
554  REAL :: USTARsigma !standard deviation of USTAR due to gustiness
555  REAL :: CM,UCN,ZCN, &
556  Z0VISC, Z0NOZ, EB, &
557  EBX, EBY, AORB, AORB1, FW, UORB, TH2, &
558  RE, FU, FUD, SWELLCOEFV, SWELLCOEFT
559  REAL :: PTURB, PVISC, SMOOTH
560  REAL XI,DELI1,DELI2
561  REAL XJ,DELJ1,DELJ2
562  REAL XK,DELK1,DELK2
563  REAL :: CONST, CONST0, CONST2, TAU1, TAU1NT, ZINF, TENSK
564  REAL X,ZARG,ZLOG,UST
565  REAL :: COSWIND, XSTRESS, YSTRESS, TAUHF
566  REAL TEMP, TEMP2
567  INTEGER IND,J,I,ISTAB
568  REAL DSTAB(3,NSPEC), DVISC, DTURB
569  REAL STRESSSTAB(3,2),STRESSSTABN(3,2)
570  !
571  INTEGER, PARAMETER :: JTOT=50
572  REAL , PARAMETER :: KM=363.,cmm=0.2325 ! K and C at phase speed minimum in rad/m
573  REAL :: OMEGACC, OMEGA, ZZ0, ZX, ZBETA, USTR, TAUR, &
574  CONST1, LEVTAIL0, X0, Y, DELY, YC, ZMU, &
575  LEVTAIL, CGTAIL, ALPHAM, FM, ALPHAT, FMEAN
576 
577  REAL, ALLOCATABLE :: W(:)
578 #ifdef W3_T0
579  REAL :: DOUT(NK,NTH)
580 #endif
581  !/
582  !/ ------------------------------------------------------------------- /
583  !/
584 #ifdef W3_S
585  CALL strace (ient, 'W3SIN4')
586 #endif
587  !
588 #ifdef W3_T
589  WRITE (ndst,9000) bbeta, ustar, usdir*rade
590 #endif
591  !
592  ! 1. Preparations
593  !
594  !JDM: Initializing values to zero, they shouldn't be used unless
595  !set in another place, but seems to solve some bugs with certain
596  !compilers.
597  dstab =0.
598  stressstab =0.
599  stressstabn =0.
600  !
601  ! Coupling coefficient times density ratio DRAT
602  !
603  const1=bbeta/kappa**2 ! needed for the tail
604  const0=const1*drat ! needed for the resolved spectrum
605  !
606  ! 1.a estimation of surface roughness parameters
607  !
608  z0visc = 0.1*nu_air/max(ustar,0.0001)
609  z0noz = max(z0visc,zz0rat*z0)
610  facln1 = u / log(zzwnd/z0noz)
611  facln2 = log(z0noz)
612  !
613  ! 1.b estimation of surface orbital velocity and displacement
614  !
615  uorb=0.
616  aorb=0.
617 
618  DO ik=1, nk
619  eb = 0.
620  ebx = 0.
621  eby = 0.
622  DO ith=1, nth
623  is=ith+(ik-1)*nth
624  eb = eb + a(is)
625  END DO
626  !
627  ! At this point UORB and AORB are the variances of the orbital velocity and surface elevation
628  !
629  uorb = uorb + eb *sig(ik)**2 * dden(ik) / cg(ik)
630  aorb = aorb + eb * dden(ik) / cg(ik) !correct for deep water only
631  END DO
632  ! FMEAN = SQRT((UORB+1E-6)/(AORB+1E-6))
633  uorb = 2*sqrt(uorb) ! significant orbital amplitude
634  aorb1 = 2*aorb**(1-0.5*sswellf(6)) ! half the significant wave height ... if SWELLF(6)=1
635  re = 4*uorb*aorb1 / nu_air ! Reynolds number
636  !
637  ! Defines the swell dissipation based on the "Reynolds number"
638  !
639  IF (sswellf(4).GT.0) THEN
640  IF (sswellf(7).GT.0.) THEN
641  smooth = 0.5*tanh((re-sswellf(4))/sswellf(7))
642  pturb=(0.5+smooth)
643  pvisc=(0.5-smooth)
644  ELSE
645  IF (re.LE.sswellf(4)) THEN
646  pturb = 0.
647  pvisc = 1.
648  ELSE
649  pturb = 1.
650  pvisc = 0.
651  END IF
652  END IF
653  ELSE
654  pturb=1.
655  pvisc=1.
656  END IF
657 
658  !
659  IF (sswellf(2).EQ.0) THEN
660  fw=max(abs(sswellf(3)),0.)
661  fu=0.
662  fud=0.
663  ELSE
664  fu=abs(sswellf(3))
665  fud=sswellf(2)
666  aorb=2*sqrt(aorb)
667  xi=(alog10(max(aorb/z0noz,3.))-abmin)/delab
668  ind = min(sizefwtable-1, int(xi))
669  deli1= min(1. ,xi-float(ind))
670  deli2= 1. - deli1
671  fw =fwtable(ind)*deli2+fwtable(ind+1)*deli1
672  END IF
673  !
674  ! 2. Diagonal
675  !
676  ! Here AS is the air-sea temperature difference in degrees. Expression given by
677  ! Abdalla & Cavaleri, JGR 2002 for Usigma. For USTARsigma ... I do not see where
678  ! I got it from, maybe just made up from drag law ...
679  !
680 #ifdef W3_STAB3
681  IF ( isnan(as) ) THEN
682  ! AS is typically NaN on land and can propagate into the domain by interpolation
683  usigma = 0.
684  ELSE
685  usigma = max(0.,-0.025*as)
686  END IF
687  ustarsigma=(1.0+u/(10.+u))*usigma
688 #endif
689 #ifdef W3_T
690  WRITE (ndst,9003) as, usigma, ustarsigma, u
691 #endif
692  ust=ustar
693  istab=3
694 #ifdef W3_STAB3
695  DO istab=1,2
696  IF (istab.EQ.1) ust=ustar*(1.-ustarsigma)
697  IF (istab.EQ.2) ust=ustar*(1.+ustarsigma)
698 #endif
699  taux = ust**2* cos(usdir)
700  tauy = ust**2* sin(usdir)
701 #ifdef W3_T
702  WRITE (ndst,9001) istab, taux, tauy, ust
703 #endif
704  !
705  ! Loop over the resolved part of the spectrum
706  !
707  stressstab(istab,:)=0.
708  stressstabn(istab,:)=0.
709  !
710  DO ik=1, nk
711  taupx=taux-abs(ttauwshelter)*stressstab(istab,1)
712  taupy=tauy-abs(ttauwshelter)*stressstab(istab,2)
713  ! With MIN and MAX the bug should disappear.... but where did it come from?
714  ustp=min((taupx**2+taupy**2)**0.25,max(ust,0.3))
715  usdirp=atan2(taupy,taupx)
716  cosu = cos(usdirp)
717  sinu = sin(usdirp)
718  is=1+(ik-1)*nth
719  cm=k(is)/sig2(is) !inverse of phase speed
720  ucn=ustp*cm+zzalp !this is the inverse wave age
721  ! the stress is the real stress (N/m^2) divided by
722  ! rho_a, and thus comparable to USTAR**2
723  ! it is the integral of rho_w g Sin/C /rho_a
724  ! (air-> waves momentum flux)
725  const2=dden2(is)/cg(ik) & !Jacobian to get energy in band
726  *grav/(sig(ik)/k(is)*drat) ! coefficient to get momentum
727  const=sig2(is)*const0
728  ! CM parameter is 1 / C_phi
729  ! Z0 corresponds to Z0+Z1 of the Janssen eq. 14
730  zcn=alog(k(is)*z0)
731  !
732  ! precomputes swell factors
733  !
734  swellcoefv=-sswellf(5)*drat*2*k(is)*sqrt(2*nu_air*sig2(is))
735  swellcoeft=-drat*sswellf(1)*16*sig2(is)**2/grav
736  !
737  DO ith=1,nth
738  is=ith+(ik-1)*nth
739  coswind=(ecos(is)*cosu+esin(is)*sinu)
740  IF (coswind.GT.0.01) THEN
741  x=coswind*ucn
742  ! this ZARG term is the argument of the exponential
743  ! in Janssen 1991 eq. 16.
744  zarg=kappa/x
745  ! ZLOG is ALOG(MU) where MU is defined by Janssen 1991 eq. 15
746  ! MU=
747  zlog=zcn+zarg
748 
749  IF (zlog.LT.0.) THEN
750  ! The source term Sp is beta * omega * X**2
751  ! as given by Janssen 1991 eq. 19
752  ! Note that this is slightly diffent from ECWAM code CY45R2 where ZLOG is replaced by ??
753  dstab(istab,is) = const*exp(zlog)*zlog**4*ucn*ucn*coswind**ssinthp
754 
755  ! Below is an example with breaking probability feeding back to the input...
756  !DSTAB(ISTAB,IS) = CONST*EXP(ZLOG)*ZLOG**4 &
757  ! *UCN*UCN*COSWIND**SSINTHP *(1+BRLAMBDA(IS)*20*SSINBR)
758  llws(is)=.true.
759  ELSE
760  dstab(istab,is) = 0.
761  llws(is)=.false.
762  END IF
763  !
764  ! Added for consistency with ECWAM implsch.F
765  !
766  IF (28.*cm*ustar*coswind.GE.1) THEN
767  llws(is)=.true.
768  END IF
769  ELSE ! (COSWIND.LE.0.01)
770  dstab(istab,is) = 0.
771  llws(is)=.false.
772  END IF
773  !
774  IF ((sswellf(1).NE.0.AND.dstab(istab,is).LT.1e-7*sig2(is)) &
775  .OR.sswellf(3).GT.0) THEN
776  !
777  dvisc=swellcoefv
778  dturb=swellcoeft*(fw*uorb+(fu+fud*coswind)*ustp)
779  !
780  dstab(istab,is) = dstab(istab,is) + pturb*dturb + pvisc*dvisc
781  END IF
782  !
783  ! Sums up the wave-supported stress
784  !
785  ! Wave direction is "direction to"
786  ! therefore there is a PLUS sign for the stress
787  temp2=const2*dstab(istab,is)*a(is)
788  IF (dstab(istab,is).LT.0) THEN
789  stressstabn(istab,1)=stressstabn(istab,1)+temp2*ecos(is)
790  stressstabn(istab,2)=stressstabn(istab,2)+temp2*esin(is)
791  ELSE
792  stressstab(istab,1)=stressstab(istab,1)+temp2*ecos(is)
793  stressstab(istab,2)=stressstab(istab,2)+temp2*esin(is)
794  END IF
795  END DO
796  END DO
797  !
798  d(:)=dstab(3,:)
799  xstress=stressstab(3,1)
800  ystress=stressstab(3,2)
801  tauwnx =stressstabn(3,1)
802  tauwny =stressstabn(3,2)
803 #ifdef W3_STAB3
804  END DO
805  d(:)=0.5*(dstab(1,:)+dstab(2,:))
806  xstress=0.5*(stressstab(1,1)+stressstab(2,1))
807  ystress=0.5*(stressstab(1,2)+stressstab(2,2))
808  tauwnx=0.5*(stressstabn(1,1)+stressstabn(2,1))
809  tauwny=0.5*(stressstabn(1,2)+stressstabn(2,2))
810 #endif
811 #ifdef W3_T
812  WRITE (ndst,9002) sum(d), sum(a), xstress, ystress, tauwnx, tauwny
813 #endif
814  s = d * a
815  !
816  ! ... Test output of arrays
817  !
818 #ifdef W3_T0
819  DO ik=1, nk
820  DO ith=1, nth
821  dout(ik,ith) = d(ith+(ik-1)*nth)
822  END DO
823  END DO
824  CALL prt2ds (ndst, nk, nk, nth, dout, sig(1), ' ', 1., &
825  0.0, 0.001, 'Diag Sin', ' ', 'NONAME')
826 #endif
827  !
828 #ifdef W3_T1
829  CALL outmat (ndst, d, nth, nth, nk, 'diag Sin')
830 #endif
831  !
832  taupx=taux-abs(ttauwshelter)*xstress
833  taupy=tauy-abs(ttauwshelter)*ystress
834  ustp=(taupx**2+taupy**2)**0.25
835  usdirp=atan2(taupy,taupx)
836 
837  ust=ustp
838  !
839  ! Computes HF tail
840  !
841  ! Computes the high-frequency contribution
842  ! the difference in spectal density (kx,ky) to (f,theta)
843  ! is integrated in this modified CONST0
844  const0=dth*sig(nk)**5/((grav**2)*tpi) &
845  *tpi*sig(nk) / cg(nk) !conversion WAM (E(f,theta) to WW3 A(k,theta)
846  temp=0.
847  DO ith=1,nth
848  is=ith+(nk-1)*nth
849  coswind=(ecos(is)*cosu+esin(is)*sinu)
850  temp=temp+a(is)*(max(coswind,0.))**3
851  END DO
852  !
853  levtail0= const0*temp ! LEVTAIL is sum over theta of A(k,theta)*cos^3(theta-wind)*DTH*SIG^5/(g^2*2pi)*2*pi*SIG/CG
854  ! which is the same as sum of E(f,theta)*cos^3(theta-wind)*DTH*SIG^5/(g^2*2pi)
855  ! reminder: sum of E(f,theta)*DTH*SIG^5/(g^2*2pi) is 2*k^3*E(k)
856 !
857 ! Computation of stress supported by tail: uses table if SINTAILPAR(1)=1 , correspoding to SINTABLE = 1
858 !
859  IF (sintailpar(1).LT.0.5) THEN
860  ALLOCATE(w(jtot))
861  w(2:jtot-1)=1.
862  w(1)=0.5
863  w(jtot)=0.5
864  x0 = 0.05
865  !
866  ustr= ust
867  zz0=z0
868  omegacc = max(sig(nk),x0*grav/ust)
869  yc = omegacc*sqrt(zz0/grav)
870 
871  ! DELY = MAX((1.-YC)/REAL(JTOT),0.)
872  ! Changed integration variable from Y to LOG(Y) and to log(K)
873  !ZINF = LOG(YC)
874  !DELY = MAX((1.-ZINF)/REAL(JTOT),0.)
875  zinf = log(sig(nk)**2/grav)
876  dely = (log(tpi/0.005)-zinf)/real(jtot)
877 
878  taur=ust**2
879  tau1=0.
880 
881  ! Integration loop over the tail wavenumbers or frequencies ...
882  DO j=1,jtot
883  !Y = YC+REAL(J-1)*DELY
884  !OMEGA = Y*SQRT(GRAV/ZZ0)
885  !OMEGA = SQRT(GRAV*Y)
886  ! This is the deep water phase speed... No surface tension !!
887  !CM = GRAV/OMEGA
888  ! With this form, Y is the wavenumber in the tail;
889  y= exp(zinf+real(j-1)*dely)
890  tensk =1+(y/km)**2
891  omega = sqrt(grav*y*tensk)
892  cm = sqrt(grav*tensk/y)
893  cgtail = 0.5*(3*(y/km)**2+1)*sqrt(grav/(y*tensk))
894  !this is the inverse wave age, shifted by ZZALP (tuning)
895  zx = ustr/cm +zzalp
896  zarg = min(kappa/zx,20.)
897  ! ZMU corresponds to EXP(ZCN)
898  zmu = min(grav*zz0/cm**2*exp(zarg),1.)
899  zlog = min(alog(zmu),0.)
900  zbeta = const1*zmu*zlog**4
901  !
902  ! Optional addition of capillary wave peak if SINTAIL2=1
903  !
904  IF (sintailpar(3).GT.0) THEN
905  IF (ustr.LT.cm) THEN
906  alpham=max(0.,0.01*(1.+alog(ustr/cm)))
907  ELSE
908  alpham=0.01*(1+3.*alog(ustr/cm))
909  END IF
910  fm=exp(-0.25*(y/km-1)**2)
911 
912  alphat=alpham*(cmm/cm)*fm ! equivalent to 2*Bh in Elfouhaily et al.
913  levtail=levtail0*0.5*(1-tanh((y-20)/5))+sintailpar(3)*0.5*(1+tanh((y-20)/5))*alphat
914  ELSE
915  levtail=levtail0
916  END IF
917  ! WRITE(991,*) 'TAIL??',SINTAILPAR(3),LEVTAIL0,LEVTAIL,ALPHAT,Y,Y/KM,OMEGA/(TPI)
918 
919  !TAU1=TAU1+W(J)*ZBETA*(USTR/UST)**2/Y*DELY ! integration over LOG(Y)
920  tau1=tau1+w(j)*zbeta*ustr**2*levtail*dely*cgtail/cm ! integration over LOG(K)
921 
922  ! NB: the factor ABS(TTAUWSHELTER) was forgotten in the TAUHFT2 table
923  !TAUR=TAUR-W(J)*ABS(TTAUWSHELTER)*USTR**2*ZBETA*LEVTAIL/Y*DELY
924  !TAUR=TAUR-W(J)*USTR**2*ZBETA*LEVTAIL*DELY ! integration over LOG(Y)
925  taur=taur-w(j)*sintailpar(2)*ustr**2*zbeta*levtail*dely*cgtail/cm ! DK/K*CG/C = D OMEGA / OMEGA
926  ustr=sqrt(max(taur,0.))
927  END DO
928  DEALLOCATE(w)
929  tau1nt=tau1
930  tauhf = tau1
931  !
932  ! In this case, uses tables for high frequency contribution to TAUW.
933  !
934  ELSE
935  ! finds the values in the tabulated stress TAUHFT
936  xi=ust/delust
937  ind = max(1,min(iustar-1, int(xi)))
938  deli1= max(min(1. ,xi-float(ind)),0.)
939  deli2= 1. - deli1
940  xj=max(0.,(grav*z0/max(ust,0.00001)**2-aalpha) / delalp)
941  j = max(1 ,min(ialpha-1, int(xj)))
942  delj1= max(0.,min(1. , xj-float(j)))
943  delj2=1. - delj1
944  IF (ttauwshelter.GT.0) THEN
945  xk = levtail0/ deltail
946  i = min(ilevtail-1, int(xk))
947  delk1= min(1. ,xk-float(i))
948  delk2=1. - delk1
949  tau1 =((tauhft2(ind,j,i)*deli2+tauhft2(ind+1,j,i)*deli1 )*delj2 &
950  +(tauhft2(ind,j+1,i)*deli2+tauhft2(ind+1,j+1,i)*deli1)*delj1)*delk2 &
951  +((tauhft2(ind,j,i+1)*deli2+tauhft2(ind+1,j,i+1)*deli1 )*delj2 &
952  +(tauhft2(ind,j+1,i+1)*deli2+tauhft2(ind+1,j+1,i+1)*deli1)*delj1)*delk1
953  ELSE
954  tau1 =(tauhft(ind,j)*deli2+tauhft(ind+1,j)*deli1 )*delj2 &
955  +(tauhft(ind,j+1)*deli2+tauhft(ind+1,j+1)*deli1)*delj1
956  END IF
957  !
958  tauhf = levtail0*ust**2*tau1
959  END IF ! End of test on use of table
960 
961  tauwx = xstress+tauhf*cos(usdirp)
962  tauwy = ystress+tauhf*sin(usdirp)
963  !
964  ! Reduces tail effect to make sure that wave-supported stress
965  ! is less than total stress, this is borrowed from ECWAM Stresso.F
966  !
967  tauw = sqrt(tauwx**2+tauwy**2)
968  ust2 = max(ustar,eps2)**2
969  tauwb = min(tauw,max(ust2-eps1,eps2**2))
970  IF (tauwb.LT.tauw) THEN
971  tauwx=tauwx*tauwb/tauw
972  tauwy=tauwy*tauwb/tauw
973  END IF
974  !
975  RETURN
976  !
977  ! Formats
978  !
979 #ifdef W3_T
980 9000 FORMAT (' TEST W3SIN4 : COMMON FACT.: ',3e10.3)
981 9001 FORMAT (' TEST W3SIN4 : ISTAB :',i2/ &
982  ' TAUX :',e12.3/ &
983  ' TAUY :',e12.3/ &
984  ' UST :',e12.3)
985 9002 FORMAT (' TEST W3SIN4 : SUM(D) :',e12.3/ &
986  ' SUM(A) :',e12.3/ &
987  ' STRESSX :',e12.3/ &
988  ' STRESSY :',e12.3/ &
989  ' TAUWNX :',e12.3/ &
990  ' TAUWNY :',e12.3)
991 9003 FORMAT (' TEST W3SIN4 : AS :',f8.4/ &
992  ' Usigma :',e12.3/ &
993  ' USTARsigma :',e12.3/ &
994  ' U :',e12.3)
995 
996 #endif
997  !/
998  !/ End of W3SIN4 ----------------------------------------------------- /
999  !/

References w3gdatmd::aalpha, w3gdatmd::bbeta, w3gdatmd::dden, w3gdatmd::dden2, delalp, deltail, delust, w3gdatmd::dth, w3gdatmd::ec2, w3gdatmd::ecos, w3gdatmd::esin, constants::fwtable, constants::grav, ialpha, w3odatmd::iaproc, ilevtail, iustar, constants::kappa, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, constants::nu_air, w3arrymd::outmat(), w3arrymd::prt2ds(), w3gdatmd::sig, w3gdatmd::sig2, w3gdatmd::sintailpar, constants::sizefwtable, w3gdatmd::ssinbr, w3gdatmd::ssinthp, w3gdatmd::sswellf, w3servmd::strace(), tauhft, tauhft2, w3gdatmd::th, constants::tpi, w3gdatmd::ttauwshelter, w3gdatmd::zz0rat, w3gdatmd::zzalp, and w3gdatmd::zzwnd.

Referenced by gxexpo(), w3exnc(), w3expo(), and w3srcemd::w3srce().

◆ w3spr4()

subroutine w3src4md::w3spr4 ( real, dimension(nth,nk), intent(in)  A,
real, dimension(nk), intent(in)  CG,
real, dimension(nk), intent(in)  WN,
real, intent(out)  EMEAN,
real, intent(out)  FMEAN,
real, intent(out)  FMEAN1,
real, intent(out)  WNMEAN,
real, intent(out)  AMAX,
real, intent(in)  U,
real, intent(in)  UDIR,
  ifdef,
  W3_FLX5 
)

Calculate mean wave parameters for the use in the source term routines.

Parameters
[in]AAction density spectrum.
[in]CGGroup velocities.
[in]WNWavenumbers.
[out]EMEANEnergy.
[out]FMEANMean frequency for determination of tail.
[out]FMEAN1Mean frequency (fm0,-1) used for reflection.
[out]WNMEANMean wavenumber.
[out]AMAXMaximum of action spectrum.
[in]UWind speed.
[in]UDIRWind direction.
[in]TAUAAtm total stress.
[in]TAUADIRAtm total stress direction.
[in]DAIRAir density.
[in,out]USTARFriction velocity.
[in,out]USDIRWind stress direction.
[in]TAUWXComponent of wave-supported stress.
[in]TAUWYComponent of wave-supported stress.
[out]CDDrag coefficient at wind level ZWND.
[out]Z0Corresponding z0.
[out]CHARNCorresponding Charnock coefficient.
[in]LLWSWind sea true/false array for each component.
[out]FMEANWSMean frequency of wind sea, used for tail.
[out]DLWMEANMean Long wave direction (L. Romero 2019).
Author
F. Ardhuin
H. L. Tolman
Date
22-Feb-2020

Definition at line 145 of file w3src4md.F90.

145  taua, tauadir, dair, &
146 #endif
147  ustar, usdir, &
148  tauwx, tauwy, cd, z0, charn, llws, fmeanws, dlwmean)
149  !/
150  !/ +-----------------------------------+
151  !/ | WAVEWATCH III SHOM |
152  !/ ! F. Ardhuin !
153  !/ | H. L. Tolman |
154  !/ | FORTRAN 90 |
155  !/ | Last update : 22-Feb-2020 |
156  !/ +-----------------------------------+
157  !/
158  !/ 03-Oct-2007 : Origination. ( version 3.13 )
159  !/ 13-Jun-2011 : Adds f_m0,-1 as FMEAN in the outout ( version 4.04 )
160  !/ 08-Jun-2018 : use STRACE and FLUSH ( version 6.04 )
161  !/ 22-Feb-2020 : Merge Romero (2019) and cleanup ( version 7.06 )
162  !/ 22-Jun-2021 : Add FLX5 to use stresses with the ST( version 7.14 )
163  !/
164  ! 1. Purpose :
165  !
166  ! Calculate mean wave parameters for the use in the source term
167  ! routines.
168  !
169  ! 2. Method :
170  !
171  ! See source term routines.
172  !
173  ! 3. Parameters :
174  !
175  ! Parameter list
176  ! ----------------------------------------------------------------
177  ! A R.A. I Action density spectrum.
178  ! CG R.A. I Group velocities.
179  ! WN R.A. I Wavenumbers.
180  ! EMEAN Real O Energy
181  ! FMEAN1 Real O Mean frequency (fm0,-1) used for reflection
182  ! FMEAN Real O Mean frequency for determination of tail
183  ! WNMEAN Real O Mean wavenumber.
184  ! AMAX Real O Maximum of action spectrum.
185  ! U Real I Wind speed.
186  ! UDIR Real I Wind direction.
187  ! TAUA Real I Atm. total stress. ( /!FLX5 )
188  ! TAUADIR Real I Atm. total stress direction. ( /!FLX5 )
189  ! DAIR Real I Air density. ( /!FLX5 )
190  ! USTAR Real I/O Friction velocity.
191  ! USDIR Real I/O wind stress direction.
192  ! TAUWX-Y Real I Components of wave-supported stress.
193  ! CD Real O Drag coefficient at wind level ZWND.
194  ! Z0 Real O Corresponding z0.
195  ! CHARN Real O Corresponding Charnock coefficient
196  ! LLWS L.A. I Wind sea true/false array for each component
197  ! FMEANWS Real O Mean frequency of wind sea, used for tail
198  ! DLWMEAN Real O Mean Long wave direction (L. Romero 2019)
199  ! ----------------------------------------------------------------
200  !
201  ! 4. Subroutines used :
202  !
203  ! STRACE Service routine.
204  !
205  ! 5. Called by :
206  !
207  ! W3SRCE Source term integration routine.
208  ! W3OUTP Point output program.
209  ! GXEXPO GrADS point output program.
210  !
211  ! 6. Error messages :
212  !
213  ! 7. Remarks :
214  !
215  ! 8. Structure :
216  !
217  ! See source code.
218  !
219  ! 9. Switches :
220  !
221  ! !/S Enable subroutine tracing.
222  ! !/T Enable test output.
223  ! !/FLX5 Direct use of stress from atmoshperic model.
224  !
225  ! 10. Source code :
226  !
227  !/ ------------------------------------------------------------------- /
228  USE w3odatmd, ONLY: iaproc
229  USE constants, ONLY: tpiinv, grav, nu_air
230  USE w3gdatmd, ONLY: nk, nth, nspec, sig, dth, dden, wwnmeanp, &
231  wwnmeanptail, fte, ftf, sstxftf, sstxftwn,&
232  sstxftftail, sswellf, esin, ecos, aaircmin, &
233  aairgb, aalpha, zzwnd, ssdsc
234 #ifdef W3_S
235  USE w3servmd, ONLY: strace
236 #endif
237 #ifdef W3_T
238  USE w3odatmd, ONLY: ndst
239  USE w3odatmd, ONLY: ndst
240 #endif
241  !
242 #ifdef W3_FLX5
243  USE w3flx5md
244 #endif
245  IMPLICIT NONE
246  !/
247  !/ ------------------------------------------------------------------- /
248  !/ Parameter list
249  !/
250  REAL, INTENT(IN) :: A(NTH,NK), CG(NK), WN(NK), U, UDIR
251 #ifdef W3_FLX5
252  REAL, INTENT(IN) :: TAUA, TAUADIR, DAIR
253 #endif
254  REAL, INTENT(IN) :: TAUWX, TAUWY
255  LOGICAL, INTENT(IN) :: LLWS(NSPEC)
256  REAL, INTENT(INOUT) :: USTAR ,USDIR
257  REAL, INTENT(OUT) :: EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, &
258  CD, Z0, CHARN, FMEANWS, DLWMEAN
259  !/
260  !/ ------------------------------------------------------------------- /
261  !/ Local parameters
262  !/
263  INTEGER :: IS, IK, ITH
264 #ifdef W3_S
265  INTEGER, SAVE :: IENT = 0
266 #endif
267 
268  REAL :: TAUW, EBAND, EMEANWS,UNZ, &
269  EB(NK),EB2(NK),ELCS, ELSN, SIGFAC
270  !/
271  !/ ------------------------------------------------------------------- /
272  !/
273 #ifdef W3_S
274  CALL strace (ient, 'W3SPR4')
275 #endif
276  !
277  unz = max( 0.01 , u )
278  ustar = max( 0.0001 , ustar )
279  !
280  emean = 0.
281  emeanws= 0.
282  fmeanws= 0.
283  fmean = 0.
284  fmean1 = 0.
285  wnmean = 0.
286  amax = 0.
287  dlwmean =0.
288  elcs =0.
289  elsn =0.
290  !
291  ! 1. Integral over directions and maximum --------------------------- *
292  !
293  DO ik=1, nk
294  eb(ik) = 0.
295  eb2(ik) = 0.
296  sigfac=sig(ik)**ssdsc(12) * dden(ik) / cg(ik)
297  DO ith=1, nth
298  is=ith+(ik-1)*nth
299  eb(ik) = eb(ik) + a(ith,ik)
300  elcs = elcs + a(ith,ik)*ecos(is)*sigfac
301  elsn = elsn + a(ith,ik)*esin(is)*sigfac
302  IF (llws(is)) eb2(ik) = eb2(ik) + a(ith,ik)
303  amax = max( amax , a(ith,ik) )
304  END DO
305  END DO
306  !
307  dlwmean=atan2(elsn,elcs)
308  !
309  ! 2. Integrate over directions -------------------------------------- *
310  !
311  DO ik=1, nk
312  eb(ik) = eb(ik) * dden(ik) / cg(ik)
313  eb2(ik) = eb2(ik) * dden(ik) / cg(ik)
314  emean = emean + eb(ik)
315  fmean = fmean + eb(ik) /sig(ik)
316  fmean1 = fmean1 + eb(ik) *(sig(ik)**(2.*wwnmeanptail))
317  wnmean = wnmean + eb(ik) *(wn(ik)**wwnmeanp)
318  emeanws = emeanws+ eb2(ik)
319  fmeanws = fmeanws+ eb2(ik)*(sig(ik)**(2.*wwnmeanptail))
320  END DO
321  !
322  ! 3. Add tail beyond discrete spectrum and get mean pars ------------ *
323  ! ( DTH * SIG absorbed in FTxx )
324  !
325  eband = eb(nk) / dden(nk)
326  emean = emean + eband * fte
327  fmean = fmean + eband * ftf
328  fmean1 = fmean1 + eband * sstxftftail
329  wnmean = wnmean + eband * sstxftwn
330  eband = eb2(nk) / dden(nk)
331  emeanws = emeanws + eband * fte
332  fmeanws = fmeanws + eband * sstxftftail
333  !
334  ! 4. Final processing
335  !
336  fmean = tpiinv * emean / max( 1.e-7 , fmean )
337  IF (fmean1.LT.1.e-7) THEN
338  fmean1=tpiinv * sig(nk)
339  ELSE
340  fmean1 = tpiinv *( max( 1.e-7 , fmean1 ) &
341  / max( 1.e-7 , emean ))**(1/(2.*wwnmeanptail))
342  ENDIF
343  wnmean = ( max( 1.e-7 , wnmean ) &
344  / max( 1.e-7 , emean ) )**(1/wwnmeanp)
345  IF (fmeanws.LT.1.e-7.OR.emeanws.LT.1.e-7) THEN
346  fmeanws=tpiinv * sig(nk)
347  ELSE
348  fmeanws = tpiinv *( max( 1.e-7 , fmeanws ) &
349  / max( 1.e-7 , emeanws ))**(1/(2.*wwnmeanptail))
350  END IF
351 
352  !
353  ! 5. Cd and z0 ----------------------------------------------- *
354  !
355  tauw = sqrt(tauwx**2+tauwy**2)
356  !
357 #ifdef W3_FLX5
358  CALL w3flx5 ( zzwnd, u, udir, taua, tauadir, dair, &
359  ustar, usdir, z0, cd, charn )
360 #else
361  CALL calc_ustar(u,tauw,ustar,z0,charn)
362  unz = max( 0.01 , u )
363  cd = (ustar/unz)**2
364  usdir = udir
365 #endif
366  !
367  ! 6. Final test output ---------------------------------------------- *
368  !
369 #ifdef W3_T
370  WRITE (ndst,9060) emean, wnmean, tpiinv, ustar, cd, z0
371 #endif
372  !
373  RETURN
374  !
375  ! Formats
376  !
377 #ifdef W3_T
378 9060 FORMAT (' TEST W3SPR4 : E,WN MN :',f8.3,f8.4/ &
379  ' TPIINV, USTAR, CD, Z0 :',f8.3,f7.2,1x,2f9.5)
380 #endif
381  !/
382  !/ End of W3SPR4 ----------------------------------------------------- /
383  !/

References w3gdatmd::aaircmin, w3gdatmd::aairgb, w3gdatmd::aalpha, calc_ustar(), w3gdatmd::dden, w3gdatmd::dth, w3gdatmd::ecos, w3gdatmd::esin, w3gdatmd::fte, w3gdatmd::ftf, constants::grav, w3odatmd::iaproc, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, constants::nu_air, w3gdatmd::sig, w3gdatmd::ssdsc, w3gdatmd::sstxftf, w3gdatmd::sstxftftail, w3gdatmd::sstxftwn, w3gdatmd::sswellf, w3servmd::strace(), constants::tpiinv, w3flx5md::w3flx5(), w3gdatmd::wwnmeanp, w3gdatmd::wwnmeanptail, and w3gdatmd::zzwnd.

Referenced by wmesmfmd::fieldindex(), gxexpo(), pdlib_w3profsmd::pdlib_jacobi_gauss_seidel_block(), w3exnc(), w3expo(), and w3srcemd::w3srce().

Variable Documentation

◆ delalp

real w3src4md::delalp

Definition at line 96 of file w3src4md.F90.

Referenced by tabu_tauhf(), tabu_tauhf2(), w3iogrmd::w3iogr(), and w3sin4().

◆ deltail

real w3src4md::deltail

Definition at line 98 of file w3src4md.F90.

98  REAL :: DELTAIL

Referenced by tabu_tauhf2(), w3iogrmd::w3iogr(), and w3sin4().

◆ deltauw

real w3src4md::deltauw

Definition at line 96 of file w3src4md.F90.

Referenced by calc_ustar(), tabu_stress(), and w3iogrmd::w3iogr().

◆ delu

real w3src4md::delu

Definition at line 96 of file w3src4md.F90.

Referenced by calc_ustar(), tabu_stress(), and w3iogrmd::w3iogr().

◆ delust

real w3src4md::delust

Definition at line 96 of file w3src4md.F90.

96  REAL :: DELUST, DELALP,DELTAUW, DELU

Referenced by tabu_tauhf(), tabu_tauhf2(), w3iogrmd::w3iogr(), and w3sin4().

◆ dikcumul

integer w3src4md::dikcumul

Definition at line 101 of file w3src4md.F90.

101  INTEGER :: DIKCUMUL

Referenced by insin4(), w3iogrmd::w3iogr(), and w3sds4().

◆ fac_kd1

real, parameter w3src4md::fac_kd1 =1.01

Definition at line 104 of file w3src4md.F90.

104  REAL, PARAMETER :: FAC_KD1=1.01, khsmax=2., khmax=2.

Referenced by insin4(), and w3sds4().

◆ fac_kd2

integer, parameter w3src4md::fac_kd2 =1000

Definition at line 103 of file w3src4md.F90.

Referenced by insin4(), and w3sds4().

◆ ialpha

integer, parameter w3src4md::ialpha =200

Definition at line 93 of file w3src4md.F90.

Referenced by insin4(), tabu_tauhf(), tabu_tauhf2(), and w3sin4().

◆ ilevtail

integer, parameter w3src4md::ilevtail =50

Definition at line 93 of file w3src4md.F90.

Referenced by insin4(), tabu_tauhf2(), and w3sin4().

◆ itaumax

integer, parameter w3src4md::itaumax =200

Definition at line 92 of file w3src4md.F90.

92  INTEGER, PARAMETER :: ITAUMAX=200,jumax=200

Referenced by calc_ustar(), insin4(), and tabu_stress().

◆ iustar

integer, parameter w3src4md::iustar =100

Definition at line 93 of file w3src4md.F90.

93  INTEGER, PARAMETER :: IUSTAR=100,ialpha=200, ilevtail=50

Referenced by insin4(), tabu_tauhf(), tabu_tauhf2(), and w3sin4().

◆ jumax

integer, parameter w3src4md::jumax =200

Definition at line 92 of file w3src4md.F90.

Referenced by calc_ustar(), insin4(), and tabu_stress().

◆ kdmax

real, parameter w3src4md::kdmax =200000.

Definition at line 105 of file w3src4md.F90.

105  REAL, PARAMETER ::KDMAX=200000.

Referenced by insin4().

◆ khmax

real, parameter w3src4md::khmax =2.

Definition at line 104 of file w3src4md.F90.

Referenced by insin4().

◆ khsmax

real, parameter w3src4md::khsmax =2.

Definition at line 104 of file w3src4md.F90.

Referenced by insin4(), and w3sds4().

◆ nkhi

integer, parameter w3src4md::nkhi =100

Definition at line 103 of file w3src4md.F90.

103  INTEGER, PARAMETER :: NKHI=100, fac_kd2=1000

Referenced by insin4().

◆ tauhft

real, dimension(:,:), allocatable w3src4md::tauhft

Definition at line 95 of file w3src4md.F90.

Referenced by insin4(), tabu_tauhf(), tabu_tauhf2(), w3iogrmd::w3iogr(), and w3sin4().

◆ tauhft2

real, dimension(:,:,:), allocatable w3src4md::tauhft2

Definition at line 95 of file w3src4md.F90.

Referenced by insin4(), tabu_tauhf2(), w3iogrmd::w3iogr(), and w3sin4().

◆ taut

real, dimension(:,:), allocatable w3src4md::taut

Definition at line 95 of file w3src4md.F90.

95  REAL, ALLOCATABLE :: TAUT(:,:),TAUHFT(:,:),TAUHFT2(:,:,:)

Referenced by calc_ustar(), insin4(), tabu_stress(), and w3iogrmd::w3iogr().

◆ tauwmax

real, parameter w3src4md::tauwmax = 2.2361

Definition at line 100 of file w3src4md.F90.

100  REAL, PARAMETER :: TAUWMAX = 2.2361 !SQRT(5.)

Referenced by calc_ustar(), and tabu_stress().

◆ umax

real, parameter w3src4md::umax = 50.

Definition at line 99 of file w3src4md.F90.

99  REAL, PARAMETER :: UMAX = 50.

Referenced by tabu_stress().

w3gdatmd::nk
integer, pointer nk
Definition: w3gdatmd.F90:1230
w3gdatmd::esc
real, dimension(:), pointer esc
Definition: w3gdatmd.F90:1234
constants::pi
real, parameter pi
PI Value of Pi.
Definition: constants.F90:71
w3gdatmd::ndtab
integer, parameter ndtab
Definition: w3gdatmd.F90:637
w3gdatmd::ssinbr
real, pointer ssinbr
Definition: w3gdatmd.F90:1324
w3gdatmd::dth
real, pointer dth
Definition: w3gdatmd.F90:1232
w3gdatmd::sstxftftail
real, pointer sstxftftail
Definition: w3gdatmd.F90:1340
w3gdatmd::nspec
integer, pointer nspec
Definition: w3gdatmd.F90:1230
w3gdatmd::ssdsbrfdf
integer, pointer ssdsbrfdf
Definition: w3gdatmd.F90:1321
w3gdatmd::wwnmeanp
real, pointer wwnmeanp
Definition: w3gdatmd.F90:1311
w3flx5md::w3flx5
subroutine w3flx5(ZWND, U10, U10D, TAUA, TAUADIR, RHOAIR, UST, USTD, Z0, CD, CHARN)
Unified process to obtain friction velocity and drag when stresses are an input (from atmospheric mod...
Definition: w3flx5md.F90:119
w3gdatmd::ssdsbck
real, pointer ssdsbck
Definition: w3gdatmd.F90:1324
w3gdatmd::sstxftwn
real, pointer sstxftwn
Definition: w3gdatmd.F90:1311
constants::nu_air
real, parameter nu_air
NU_AIR Kinematic viscosity of air (m2/s).
Definition: constants.F90:64
w3gdatmd::iktab
integer, dimension(:,:), pointer iktab
Definition: w3gdatmd.F90:1321
w3gdatmd::ssdsc
real, dimension(:), pointer ssdsc
Definition: w3gdatmd.F90:1324
w3src4md::delust
real delust
Definition: w3src4md.F90:96
w3gdatmd::cumulw
real, dimension(:,:), pointer cumulw
Definition: w3gdatmd.F90:1323
w3gdatmd::sig
real, dimension(:), pointer sig
Definition: w3gdatmd.F90:1234
w3odatmd::iaproc
integer, pointer iaproc
Definition: w3odatmd.F90:457
w3gdatmd::wwnmeanptail
real, pointer wwnmeanptail
Definition: w3gdatmd.F90:1340
w3gdatmd::satweights
real, dimension(:,:), pointer satweights
Definition: w3gdatmd.F90:1323
w3gdatmd::fachfe
real, pointer fachfe
Definition: w3gdatmd.F90:1232
w3gdatmd::dcki
real, dimension(:,:), pointer dcki
Definition: w3gdatmd.F90:1323
w3gdatmd::ecos
real, dimension(:), pointer ecos
Definition: w3gdatmd.F90:1234
constants::rade
real, parameter rade
RADE Conversion factor from radians to degrees.
Definition: constants.F90:76
w3gdatmd::ssdsbrf1
real, pointer ssdsbrf1
Definition: w3gdatmd.F90:1324
w3gdatmd::dden2
real, dimension(:), pointer dden2
Definition: w3gdatmd.F90:1234
w3gdatmd::ssdsabk
real, pointer ssdsabk
Definition: w3gdatmd.F90:1324
w3gdatmd::zzwnd
real, pointer zzwnd
Definition: w3gdatmd.F90:1311
w3gdatmd::dsip
real, dimension(:), pointer dsip
Definition: w3gdatmd.F90:1234
w3gdatmd::ssdsbt
real, pointer ssdsbt
Definition: w3gdatmd.F90:1311
w3arrymd::outmat
subroutine outmat(NDS, A, MX, NX, NY, MNAME)
Definition: w3arrymd.F90:988
w3gdatmd::th
real, dimension(:), pointer th
Definition: w3gdatmd.F90:1234
w3src4md::taut
real, dimension(:,:), allocatable taut
Definition: w3src4md.F90:95
w3src4md::tauhft2
real, dimension(:,:,:), allocatable tauhft2
Definition: w3src4md.F90:95
w3gdatmd::zzalp
real, pointer zzalp
Definition: w3gdatmd.F90:1311
w3gdatmd::zz0rat
real, pointer zz0rat
Definition: w3gdatmd.F90:1311
scrip_timers::status
character(len=8), dimension(max_timers), save status
Definition: scrip_timers.f:63
w3odatmd::ndse
integer, pointer ndse
Definition: w3odatmd.F90:456
w3flx5md
Unified process to obtain friction velocity and drag when stresses are an input (from atmospheric mod...
Definition: w3flx5md.F90:28
w3gdatmd::nkd
integer, parameter nkd
Definition: w3gdatmd.F90:636
w3gdatmd::es2
real, dimension(:), pointer es2
Definition: w3gdatmd.F90:1234
w3gdatmd::sstxftf
real, pointer sstxftf
Definition: w3gdatmd.F90:1311
constants::sizefwtable
integer, parameter sizefwtable
SIZEFWTABLE.
Definition: constants.F90:91
w3gdatmd::esin
real, dimension(:), pointer esin
Definition: w3gdatmd.F90:1234
w3gdatmd::ssdspbk
real, pointer ssdspbk
Definition: w3gdatmd.F90:1324
w3gdatmd::ssdsp
real, pointer ssdsp
Definition: w3gdatmd.F90:1311
w3dispmd::wavnu2
subroutine wavnu2(W, H, K, CG, EPS, NMAX, ICON)
Definition: w3dispmd.F90:204
w3gdatmd::aaircmin
real, pointer aaircmin
Definition: w3gdatmd.F90:1183
w3gdatmd::bbeta
real, pointer bbeta
Definition: w3gdatmd.F90:1311
constants::kappa
real, parameter kappa
KAPPA von Karman's constant (N.D.).
Definition: constants.F90:69
w3src4md::deltail
real deltail
Definition: w3src4md.F90:98
w3src4md::delalp
real delalp
Definition: w3src4md.F90:96
w3servmd
Definition: w3servmd.F90:3
w3odatmd::flogrd
logical, dimension(:,:), pointer flogrd
Definition: w3odatmd.F90:478
w3gdatmd::ssdsbm
real, dimension(:), pointer ssdsbm
Definition: w3gdatmd.F90:1311
constants::tpiinv
real, parameter tpiinv
TPIINV Inverse of 2*Pi.
Definition: constants.F90:74
w3gdatmd::nth
integer, pointer nth
Definition: w3gdatmd.F90:1230
w3odatmd
Definition: w3odatmd.F90:3
w3gdatmd::ssdsbint
real, pointer ssdsbint
Definition: w3gdatmd.F90:1324
w3gdatmd::satindices
integer, dimension(:,:), pointer satindices
Definition: w3gdatmd.F90:1321
constants::dwat
real, parameter dwat
DWAT Density of water (kg/m3).
Definition: constants.F90:62
w3src4md::dikcumul
integer dikcumul
Definition: w3src4md.F90:101
file
file(STRINGS ${CMAKE_BINARY_DIR}/switch switch_strings) separate_arguments(switches UNIX_COMMAND $
Definition: CMakeLists.txt:3
w3gdatmd::sig2
real, dimension(:), pointer sig2
Definition: w3gdatmd.F90:1234
w3gdatmd::ttauwshelter
real, pointer ttauwshelter
Definition: w3gdatmd.F90:1311
constants::tpi
real, parameter tpi
TPI 2*Pi.
Definition: constants.F90:72
w3servmd::strace
subroutine strace(IENT, SNAME)
Definition: w3servmd.F90:148
w3snl4md::diag2
real, dimension(:,:), allocatable diag2
Definition: w3snl4md.F90:196
constants::fwtable
real, dimension(0:sizefwtable) fwtable
FWTABLE.
Definition: constants.F90:92
w3gdatmd::sintailpar
real, dimension(:), pointer sintailpar
Definition: w3gdatmd.F90:1324
w3gdatmd::nkhs
integer, parameter nkhs
Definition: w3gdatmd.F90:636
w3arrymd
Definition: w3arrymd.F90:3
w3gdatmd::xfr
real, pointer xfr
Definition: w3gdatmd.F90:1232
w3gdatmd::ssinthp
real, pointer ssinthp
Definition: w3gdatmd.F90:1311
w3gdatmd::fte
real, pointer fte
Definition: w3gdatmd.F90:1232
w3gdatmd::ssdsiso
integer, pointer ssdsiso
Definition: w3gdatmd.F90:1321
w3src4md::tauhft
real, dimension(:,:), allocatable tauhft
Definition: w3src4md.F90:95
w3gdatmd::aairgb
real, pointer aairgb
Definition: w3gdatmd.F90:1183
w3src4md::deltauw
real deltauw
Definition: w3src4md.F90:96
w3gdatmd::qbi
real, dimension(:,:), pointer qbi
Definition: w3gdatmd.F90:1323
w3odatmd::ndst
integer, pointer ndst
Definition: w3odatmd.F90:456
constants::debug_node
integer, parameter debug_node
DEBUG_NODE Node number used for debugging.
Definition: constants.F90:99
constants
Define some much-used constants for global use (all defined as PARAMETER).
Definition: constants.F90:20
w3gdatmd::dden
real, dimension(:), pointer dden
Definition: w3gdatmd.F90:1234
w3gdatmd
Definition: w3gdatmd.F90:16
w3gdatmd::sswellf
real, dimension(:), pointer sswellf
Definition: w3gdatmd.F90:1311
w3gdatmd::ssdscos
real, pointer ssdscos
Definition: w3gdatmd.F90:1311
w3gdatmd::igpars
real, dimension(:), pointer igpars
Definition: w3gdatmd.F90:1142
constants::file_endian
character(*), parameter file_endian
FILE_ENDIAN Filled by preprocessor with 'big_endian', 'little_endian', or 'native'.
Definition: constants.F90:86
w3servmd::extcde
subroutine extcde(IEXIT, UNIT, MSG, FILE, LINE, COMM)
Definition: w3servmd.F90:736
w3arrymd::prt2ds
subroutine prt2ds(NDS, NFR0, NFR, NTH, E, FR, UFR, FACSP, FSC, RRCUT, PRVAR, PRUNIT, PNTNME)
Definition: w3arrymd.F90:1943
w3gdatmd::aalpha
real, pointer aalpha
Definition: w3gdatmd.F90:1311
w3gdatmd::capchnk
real, dimension(:), pointer capchnk
Definition: w3gdatmd.F90:1324
w3gdatmd::ftf
real, pointer ftf
Definition: w3gdatmd.F90:1232
w3gdatmd::ssdsdth
real, pointer ssdsdth
Definition: w3gdatmd.F90:1311
w3src4md::delu
real delu
Definition: w3src4md.F90:96
w3gdatmd::ssdshck
real, pointer ssdshck
Definition: w3gdatmd.F90:1324
w3gdatmd::ec2
real, dimension(:), pointer ec2
Definition: w3gdatmd.F90:1234
w3gdatmd::zz0max
real, pointer zz0max
Definition: w3gdatmd.F90:1311
w3dispmd
Definition: w3dispmd.F90:3
constants::grav
real, parameter grav
GRAV Acc.
Definition: constants.F90:61
w3gdatmd::ssdsbr
real, pointer ssdsbr
Definition: w3gdatmd.F90:1311