!/ ------------------------------------------------------------------- /
Module W3FLD1MD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP/NOPP |
!/ | B. G. Reichl |
!/ | FORTRAN 90 |
!/ | Last update : 22-Mar-2021 |
!/ +-----------------------------------+
!/
!/ 01-Jul-2013 : Origination. ( version 4.xx )
!/ 18-Mar-2015 : Clean-up/prepare for distribution ( version 5.12 )
!/ 15-Jan-2016 : Updates ( version 5.12 )
!/ ( B. G. Reichl )
!/ 27-Jul-2016 : Added Charnock output (J.Meixner) ( version 5.12 )
!/ 22-Jun-2018 : updated SIG2WN subroutine (X.Chen) ( version 6.06 )
!/ modified the range of wind profile computation;
!/ corrected direction of the shortest waves
!/ 22-Mar-2021 : Consider DAIR a variable ( version 7.13 )
!/
!/
!/ Copyright 2009 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.
!/
! 1. Purpose :
!
! This Module contains routines to compute the wind stress vector
! from the wave spectrum, the wind vector, and the lower atmospheric
! stability (the form included here is for neutral conditions, but
! the structure needed to include stability is included in comments).
! The stress calculated via this subroutine is
! intended for coupling to serve as the boundary condition
! between the ocean and atmosphere, and (for now)
! and has no impact on the wave spectrum calculated.
! The calculation in w3fld1 is based on the method
! presented in Reichl, Hara, and Ginis (2014), "Sea State Dependence
! of the Wind Stress under Hurricane Winds."
!
! 2. Variables and types :
!
! Not applicable.
!
! 3. Subroutines and functions :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! W3FLD1 Subr. Public Reichl et al. 2014 stress calculation
! INFLD1 Subr. Public Corresponding initialization routine.
! APPENDTAIL Subr. Public Modification of tail for calculation
! SIG2WN Subr. Public Depth-dependent dispersion relation
! WND2Z0M Subr. Public Wind to roughness length
! ----------------------------------------------------------------
!
! 4. Subroutines and functions used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Remarks :
!
! 6. Switches :
!
! !/S Enable subroutine tracing.
! !/
!
! 7. Source code :
!/
!/ ------------------------------------------------------------------- /
!/
!
PUBLIC
! Tail_Choice: Chose the method to determine the level of the tail
INTEGER, SAVE :: Tail_Choice
#ifdef W3_OMPG
!$omp threadprivate(Tail_Choice)
#endif
REAL, SAVE :: Tail_Level !if Tail_Choice=0, tail is constant
REAL, SAVE :: Tail_transition_ratio1! freq/fpi where tail
! adjustment begins
REAL, SAVE :: Tail_transition_ratio2! freq/fpi where tail
! adjustment ends
#ifdef W3_OMPG
!$omp threadprivate(Tail_Level)
!$omp threadprivate(Tail_transition_ratio1,Tail_transition_ratio2)
#endif
!/
CONTAINS
!/ ------------------------------------------------------------------- /
SUBROUTINE W3FLD1( ASPC, FPI, WNDX,WNDY, ZWND, &
DEPTH, RIB, DAIR, UST, USTD, Z0, &
TAUNUX, TAUNUY, CHARN)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP/NOPP |
!/ | B. G. Reichl |
!/ | FORTRAN 90 |
!/ | Last update : 22-Mar-2021 |
!/ +-----------------------------------+
!/
!/ 01-Jul-2013 : Origination. ( version 4.xx )
!/ 18-Mar-2015 : Prepare for submission ( version 5.12 )
!/ 22-Mar-2021 : Consider DAIR a variable ( version 7.13 )
!/
! 1. Purpose :
!
! Diagnostic stress vector calculation from wave spectrum, lower
! atmosphere stability, and wind vector (at some given height).
! The height of wind vector is assumed to be within the constant
! stress layer. These parameterizations are meant to be performed
! at wind speeds > 10 m/s, and may not converge for extremely young
! seas (i.e. starting from flat sea conditions).
!
! 2. Method :
! See Reichl et al. (2014).
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ASPC Real I 1-D Wave action spectrum.
! FPI Real I Peak input frequency.
! WNDX Real I X-dir wind (assumed referenced to current)
! WNDY Real I Y-dir wind (assumed referenced to current)
! ZWND Real I Wind height.
! DEPTH Real I Water depth.
! RIB REAL I Bulk Richardson in lower atmosphere
! (for determining stability in ABL to get
! 10 m neutral wind)
! DAIR REAL I Air density
! TAUNUX Real 0 X-dir viscous stress (guessed from prev.)
! TAUNUY Real 0 Y-dir viscous stress (guessed from prev.)
! UST Real O Friction velocity.
! USTD Real O Direction of friction velocity.
! Z0 Real O Surface roughness length
! CHARN Real O,optional Charnock parameter
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3ASIM Subr. W3ASIMMD Air-sea interface module.
! W3EXPO Subr. N/A Point output post-processor.
! GXEXPO Subr. N/A GrADS point output post-processor.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS, ONLY: GRAV, DWAT, TPI, PI, KAPPA
USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, DTH, XFR, TH
USE W3ODATMD, ONLY: NDSE
USE W3SERVMD, ONLY: EXTCDE
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: ASPC(NSPEC), WNDX, WNDY, &
ZWND, DEPTH, RIB, DAIR, FPI
REAL, INTENT(OUT) :: UST, USTD, Z0
REAL, INTENT(OUT), OPTIONAL :: CHARN
REAL, INTENT(INOUT) :: TAUNUX, TAUNUY
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
REAL, PARAMETER :: NU=0.105/10000.0
REAL, PARAMETER :: DELTA=0.03
! Commonly used parameters
REAL :: wnd_in_mag, wnd_in_dir
!For Calculating Tail
REAL :: KMAX, KTAILA, KTAILB, KTAILC
REAL :: SAT, z01, z02, u10
LOGICAL :: ITERFLAG
INTEGER :: COUNT
!For Iterations
REAL :: DTX, DTY, iter_thresh, &
USTSM, Z0SM, Z1
!For stress calculation
REAL :: WAGE, CBETA, BP, CD, &
USTRB, ANGDIF, USTAR, ZNU, &
TAUT, TAUX, TAUY, BETAG, TAUDIR, &
TAUDIRB
!For wind profile calculation
REAL :: UPROFV, VPROFV
!For wind profile iteration
REAL :: WND_1X, WND_1Y, &
WND_2X, WND_2Y, &
WND_3X, WND_3Y, &
DIFU10XX, DIFU10YX, DIFU10XY, DIFU10YY, &
FD_A, FD_B, FD_C, FD_D, &
DWNDX, DWNDY, &
APAR, CH,UITV, VITV,USTL,&
CK
!For adding stability to wind profile
REAL :: WND_TOP, ANG_TOP, WND_PA, WND_PE, &
WND_PEx, WND_PEy, WND_PAx, WND_PAy, &
CDM
INTEGER :: NKT, K, T, Z2, ITER, ZI, ZII, &
I, CTR, ITERATION, KA1, KA2, &
KA3, KB
! For defining extended spectrum with appended tail.
REAL, ALLOCATABLE, DIMENSION(:) :: WN, DWN, CP,SIG2
REAL, ALLOCATABLE, DIMENSION(:,:) :: SPC2
REAL, ALLOCATABLE, DIMENSION(:) :: TLTN, TLTE, TAUD, &
TLTND, &
TLTED, ZOFK, UPROF, VPROF, &
FTILDE, UP1, VP1, UP, VP, &
TLTNA, TLTEA
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
LOGICAL :: FSFL1,FSFL2, CRIT1, CRIT2
LOGICAL :: IT_FLAG1, IT_FLAG2
LOGICAL, SAVE :: FIRST = .TRUE.
#ifdef W3_OMPG
!$omp threadprivate( FIRST)
#endif
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'W3FLD1')
#endif
!
! 0. Initializations ------------------------------------------------ *
!
! **********************************************************
! *** The initialization routine should include all ***
! *** initialization, including reading data from files. ***
! **********************************************************
!
IF ( FIRST ) THEN
CALL INFLD
FIRST = .FALSE.
END IF
wnd_in_mag = sqrt( wndx**2 + wndy**2 )
wnd_in_dir = atan2(wndy, wndx)
!----------------------------------------------------------+
! Assume wind input is neutral 10 m wind. If wind input +
! is not 10 m, tail level will need to be calculated based +
! on esimation of 10 m wind. +
!----------------------------------------------------------+
u10 = wnd_in_mag
! - Get tail level
if (Tail_Choice.eq.0) then
SAT=Tail_Level
elseif (Tail_Choice.eq.1) then
CALL WND2SAT(U10,SAT)
endif
!
! 1. Attach Tail ---------------------------------------------------- *
!
! If the depth remains constant, the allocation could be limited to the
! first time step. Since this code is designed for coupled
! implementation where the water level can change, I keep it the
! allocation on every time step. When computational efficiency is
! important, this process may be rethought.
!
! i. Find maximum wavenumber of input spectrum
call sig2wn(sig(nk),depth,kmax)
NKT = NK
! ii. Find additional wavenumber bins to extended to cm scale waves
DO WHILE ( KMAX .LT. 366.0 )
NKT = NKT + 1
KMAX = ( KMAX * XFR**2 )
ENDDO!K<366
! iii. Allocate new "extended" spectrum
ALLOCATE( WN(NKT), DWN(NKT), CP(NKT), SIG2(NKT),SPC2(NKT,NTH), &
TLTN(NKT), TLTE(NKT), TAUD(NKT), &
TLTND(NKT), TLTED(NKT), ZOFK(NKT), UPROF(NKT+1),&
VPROF(NKT+1), FTILDE(NKT), UP1(NKT+1),VP1(NKT+1), &
UP(NKT+1), VP(NKT+1), TLTNA(NKT),TLTEA(NKT))
!
! 1a. Build Discrete Wavenumbers for defining extended spectrum on---- *
!
!i. Copy existing sig to extended sig2, calculate phase speed.
DO K = 1, NK !existing spectrum
call sig2wn(sig(k),depth,wn(k))
CP(K) = ( SIG(K) / WN(K) )
sig2(k) = sig(k)
ENDDO!K
!ii. Calculate extended sig2 and phase speed.
DO K = ( NK + 1 ), ( NKT) !extension
sig2(k) = sig2(k-1) *XFR
call sig2wn(sig2(k),depth,wn(k))
CP(K) = SIG2(K) / WN(K)
ENDDO!K
!iii. Calculate dk's for integrations.
DO K = 1, NKT-1
DWN(K) = WN(K+1) - WN(K)
ENDDO
DWN(NKT) = WN(NKT)*XFR**2 - WN(NKT)
!
! 1b. Attach initial tail--------------------------------------------- *
!
!i. Convert action spectrum to variance spectrum
! SPC(k,theta) = A(k,theta) * sig(k)
! This could be redone for computational efficiency
I=0
DO K=1, NK
DO T=1, NTH
I = I + 1
SPC2(K,T) = ASPC(I) * SIG(K)
ENDDO!T
ENDDO!K
!ii. Extend k^-3 tail to extended spectrum
DO K=NK+1, NKT
DO T=1, NTH
SPC2(K,T)=SPC2(NK,T)*WN(NK)**3.0/WN(K)**(3.0)
ENDDO!T
ENDDO!K
!
! 1c. Calculate transitions for new (constant saturation ) tail ------ *
!
!
!i. Find wavenumber for beginning spc level transition to tail
call sig2wn (FPI*TPI*tail_transition_ratio1,depth,ktaila )
!ii. Find wavenumber for ending spc level transition to tail
call sig2wn (FPI*TPI*tail_transition_ratio2,depth,ktailb )
!iii. Find wavenumber for ending spc direction transition to tail
KTAILC= KTAILB * 2.0
!iv. Find corresponding indices of wavenumber transitions
KA1 = 2 ! Do not modify 1st wavenumber bin
DO WHILE ( ( KTAILA .GE. WN(KA1) ) .AND. (KA1 .LT. NKT-6) )
KA1 = KA1 + 1
ENDDO
KA2 = KA1+2
DO WHILE ( ( KTAILB .GE. WN(KA2) ) .AND. (KA2 .LT. NKT-4) )
KA2 = KA2 + 1
ENDDO
KA3 = KA2+2
DO WHILE ( ( KTAILC .GE. WN(KA3)) .AND. (KA3 .LT. NKT-2) )
KA3 = KA3 + 1
ENDDO
!v. Call subroutine to perform actually tail truncation
! only if there is some energy in spectrum
CALL APPENDTAIL(SPC2, WN, NKT, KA1, KA2, KA3,&
wnd_in_dir, SAT)
! Spectrum is now set for stress integration
!
! 2. Prepare for iterative calculation of wave-form stress----------- *
!
DTX = 0.00005
DTY = 0.00005
iter_thresh = 0.001
!
! 2a. Calculate initial guess for viscous stress from smooth-law------ *
! (Would be preferable to use prev. step)
!
Z0SM = 0.001 !Guess
IT_FLAG1 = .true.
ITERATION = 0
DO WHILE( IT_FLAG1 )
ITERATION = ITERATION + 1
Z1 = Z0SM
USTSM = KAPPA * wnd_in_mag / ( LOG( ZWND / Z1 ) )
Z0SM = 0.132 * NU / USTSM
IF ( (ABS( Z0SM - Z1 ) .LT. 10.0**(-6)) .OR.&
( ITERATION .GT. 5 )) THEN
IT_FLAG1 = .false.
ENDIF
ENDDO
ITERATION = 1
! Guessed values of viscous stress
TAUNUX = USTSM**2 * DAIR * wndx / wnd_in_mag
TAUNUY = USTSM**2 * DAIR * wndy / wnd_in_mag
!
! 3. Enter iterative calculation of wave form/skin stress---------- *
!
IT_FLAG1 = .true.
DO WHILE (IT_FLAG1)
DO ITER=1, 3 !3 loops for TAUNU iteration
Z2 = NKT
! First : TAUNUX + DX
IF (ITER .EQ. 1) THEN
TAUNUX = TAUNUX + DTX
! Second : TAUNUY + DY
ELSEIF (ITER .EQ. 2) THEN
TAUNUX = TAUNUX - DTX
TAUNUY = TAUNUY + DTY
! Third : unmodified
ELSEIF (ITER .EQ. 3) THEN
TAUNUY = TAUNUY - DTY
ENDIF
! Near surface turbulent stress = taunu
TLTN(1) = TAUNUY
TLTE(1) = TAUNUX
CALL APPENDTAIL(SPC2, WN, NKT, KA1, KA2, KA3,&
atan2(TAUNUY,TAUNUX), SAT)
!|---------------------------------------------------------------------|
!|-----Calculate first guess at growth rate and local turbulent stress-|
!|-----for integration as a function of wavedirection------------------|
!|---------------------------------------------------------------------|
DO ZI = 2, NKT
USTL=0.0
TLTND(zi)=0.0
TLTED(zi)=0.0
Z2 = Z2 - 1
! Use value of prev. wavenumber/height
TAUD(ZI) = ATAN2( TLTN(ZI-1), TLTE(ZI-1))
USTL = SQRT( SQRT( TLTN(ZI-1)**2 + TLTE(ZI-1)**2 )/ DAIR )
DO T = 1, NTH
ANGDIF=TAUD(ZI)-TH(T) !stress/wave angle
IF ( COS( ANGDIF ) .GE. 0.0 ) THEN !Waves aligned
WAGE = CP(Z2) / (USTL)
! First, waves much slower than wind.
IF ( WAGE .LT. 10. ) THEN
CBETA = 25.0
! Transition from waves slower than wind to faster
ELSEIF ( ( WAGE .GE. 10.0 ) .AND. &
( WAGE .LE. 25.0 ) ) THEN
CBETA = 10.0 + 15.0 * COS( PI * ( WAGE - 10.0 ) &
/ 15.0 )
! Waves faster than wind
ELSEIF ( WAGE .GT. 25.0 ) THEN
CBETA = -5.0
ENDIF
! Waves opposing wind
ELSE
CBETA = -25.0
ENDIF
!Integrate turbulent stress
TLTND(ZI) =TLTND(ZI)+( SIN( TH(T) ) * COS( ANGDIF )**2)&
* CBETA * SPC2(Z2,T) * &
SQRT( TLTE(ZI-1)**2 + TLTN(ZI-1)**2.0 ) &
* ( WN(Z2)**2.0 )*DTH
TLTED(ZI) = TLTED(ZI)+(COS( TH(T) ) * COS( ANGDIF )**2)&
* CBETA * SPC2(Z2,T) * &
SQRT( TLTE(ZI-1)**2 + TLTN(ZI-1)**2.0 ) &
* ( WN(Z2)**2.0 )*DTH
ENDDO
!|---------------------------------------------------------------------|
!|-----Complete the integrations---------------------------------------|
!|---------------------------------------------------------------------|
IF (ZI .EQ. 2) THEN
!First turbulent stress bin above taunu
TLTNA(ZI) = TLTND(ZI) * DWN(Z2) * 0.5
TLTEA(ZI) = TLTED(ZI) * DWN(Z2) * 0.5
ELSE
TLTNA(ZI)=(TLTND(ZI)+TLTND(ZI-1))*0.5*DWN(Z2)
TLTEA(ZI)=(TLTED(ZI)+TLTED(ZI-1))*0.5*DWN(Z2)
ENDIF
TLTN(ZI)=TLTN(ZI-1)+TLTNA(ZI)
TLTE(ZI)=TLTE(ZI-1)+TLTEA(ZI)
ENDDO
TAUY=TLTN(NKT)
TAUX=TLTE(NKT)
! This is the first guess at the stress.
!|---------------------------------------------------------------------|
!|----Iterate til convergence------------------------------------------|
!|---------------------------------------------------------------------|
USTRB=SQRT(SQRT(TAUY**2.0+TAUX**2.0)/DAIR)
TAUDIRB=atan2(TAUY,TAUX)
IT_FLAG2 = .TRUE.
CTR=1
DO WHILE ( (IT_FLAG2) .AND. ( CTR .LT. 10 ) )
Z2=NKT+1
DO ZI=1, NKT
Z2=Z2-1
USTL=0.0
TLTED(zi)=0.0
TLTND(zi)=0.0
FTILDE(z2)=0.0
TAUD(ZI) = ATAN2(TLTN(ZI),TLTE(ZI))
USTL = SQRT(SQRT(TLTN(ZI)**2+TLTE(ZI)**2)/DAIR)
DO T=1, NTH
BETAG=0.0
ANGDIF = TAUD(ZI)-TH(T)
IF ( COS( ANGDIF ) .GE. 0.0 ) THEN
WAGE = CP(Z2) / (USTL)
IF ( WAGE .LT. 10 ) THEN
CBETA = 25.0
ELSEIF ( ( WAGE .GE. 10.0 ) .AND. &
( WAGE .LE. 25.0 ) ) THEN
CBETA = 10.0 + 15.0 * COS( PI * ( WAGE - 10.0 ) &
/ 15.0 )
ELSEIF ( WAGE .GT. 25.0 ) THEN
CBETA = -5.0
ENDIF
ELSE
CBETA = -25.0
ENDIF
BP = SQRT( (COS( TH(T) ) * COS( ANGDIF )**2.0)**2.0 &
+ (SIN( TH(T) ) * COS( ANGDIF )**2.0)**2.0 )
BETAG=BP*CBETA*SQRT(TLTE(ZI)**2.0+TLTN(ZI)**2.0) &
/(DWAT)*SIG2(Z2)/CP(Z2)**2
FTILDE(Z2) = FTILDE(Z2) + BETAG * DWAT * GRAV &
* SPC2(Z2,T) * DTH
TLTND(zi) =tltnd(zi)+ (SIN( TH(T) ) * COS( ANGDIF )**2.0)&
* CBETA * SPC2(Z2,T) * SQRT( &
TLTE(ZI)**2.0 + TLTN(ZI)**2.0 ) * &
( WN(Z2)**2.0 )*dth
TLTED(zi) = tlted(zi)+(COS( TH(T) ) * COS( ANGDIF )**2.0)&
* CBETA * SPC2(Z2,T) * SQRT( &
TLTE(ZI)**2.0 + TLTN(ZI)**2.0 ) * &
( WN(Z2)**2.0 )*dth
ENDDO
IF (ZI .EQ. 1) THEN
TLTNA(ZI)=TLTND(ZI)*DWN(Z2)*0.5
TLTEA(ZI)=TLTED(ZI)*DWN(Z2)*0.5
ELSE
TLTNA(ZI)=(TLTND(ZI)+TLTND(ZI-1))*0.5*DWN(Z2)
TLTEA(ZI)=(TLTED(ZI)+TLTED(ZI-1))*0.5*DWN(Z2)
ENDIF
IF (ZI.GT.1) then
TLTN(ZI)=TLTN(ZI-1)+TLTNA(ZI)
TLTE(ZI)=TLTE(ZI-1)+TLTEA(ZI)
else
TLTN(ZI)=TAUNUY+TLTNA(ZI)
TLTE(ZI)=TAUNUX+TLTEA(ZI)
endif
ENDDO
TAUY=TLTN(NKT) !by NKT full stress is entirely
TAUX=TLTE(NKT) !from turbulent stress
TAUT=SQRT(TAUY**2.0+TAUX**2.0)
USTAR=SQRT(SQRT(TAUY**2.0+TAUX**2.0)/DAIR)
TAUDIR=atan2(TAUY, TAUX)
! Note: add another criterion (stress direction) for iteration.
CRIT1=(ABS(USTAR-USTRB)*100.0)/((USTAR+USTRB)*0.5) .GT. 0.1
CRIT2=(ABS(TAUDIR-TAUDIRB)*100.0/(TAUDIR+TAUDIRB)*0.5) .GT. 0.1
IF (CRIT1 .OR. CRIT2) THEN
! IF ((ABS(USTAR-USTRB)*100.0)/((USTAR+USTRB)*0.5) .GT. 0.1) THEN
USTRB=USTAR
TAUDIRB=TAUDIR
CTR=CTR+1
ELSE
IT_FLAG2 = .FALSE.
ENDIF
ENDDO
! Note: search for the top of WBL from top to bottom (avoid problems
! caused by for very long swell)
KB=NKT
DO WHILE(((TLTN(KB)**2+TLTE(KB)**2)/(TAUX**2+TAUY**2)).GT. &
.99)
KB=KB-1
ENDDO
KB=KB+1
!|---------------------------------------------------------------------|
!|----Now begin work on wind profile-----------------------------------|
!|---------------------------------------------------------------------|
DO I=1,NKT
ZOFK(I)=DELTA/WN(I)
ENDDO
ZNU=0.1 * 1.45E-5 / SQRT(SQRT(TAUNUX**2.0+TAUNUY**2.0)/DAIR)
UPROF(1:NKT+1)=0.0
VPROF(1:NKT+1)=0.0
UPROFV=0.0
VPROFV=0.0
ZI=1
Z2=NKT
UP1(ZI) = ( ( ( WN(Z2)**2 / DELTA ) * FTILDE(z2) ) + &
( DAIR / ( ZOFK(Z2) * KAPPA ) ) * ( SQRT( &
TLTN(ZI)**2 + TLTE(ZI)**2 ) / DAIR )**(3/2) ) &
* ( TLTE(ZI) ) / ( TLTE(ZI) * TAUX &
+ TLTN(ZI) * TAUY )
VP1(ZI) = ( ( ( WN(Z2)**2 / DELTA ) * FTILDE(z2) ) + &
( DAIR / ( ZOFK(Z2) * KAPPA ) ) * ( SQRT ( &
TLTN(ZI)**2 + TLTE(ZI)**2 ) / DAIR )**(3/2) ) &
* ( TLTN(ZI) ) / ( TLTE(ZI) * TAUX &
+ TLTN(ZI) * TAUY )
UP(ZI) = UP1(ZI)
VP(ZI) = VP1(ZI)
UPROF(ZI) = DAIR / KAPPA * ( SQRT( TAUNUX**2.0 + TAUNUY**2.0 ) &
/ DAIR )**(1.5) * ( TAUNUX / ( TAUX * &
TAUNUX + TAUY * TAUNUY ) ) * LOG( &
ZOFK(Z2) / ZNU )
VPROF(ZI) = DAIR / KAPPA * ( SQRT( TAUNUX**2.0 + TAUNUY**2.0 ) &
/ DAIR )**(1.5) * ( TAUNUY / ( TAUX * &
TAUNUX + TAUY * TAUNUY ) ) * LOG( &
ZOFK(Z2) / ZNU )
!Noted: wind profile computed till the inner layer height of the longest
!wave, not just to the top of wave boundary layer (previous)
DO ZI=2, NKT
Z2 = Z2 - 1
UP1(ZI) = ( ( ( WN(Z2)**2.0 / DELTA ) * FTILDE(Z2) ) + &
( DAIR / ( ZOFK(Z2) * KAPPA ) ) * ( SQRT( &
TLTN(ZI)**2.0 + TLTE(ZI)**2.0 ) / DAIR )**(1.5) ) &
* ( TLTE(ZI) ) / ( TLTE(ZI) * TAUX + &
TLTN(ZI) * TAUY )
VP1(ZI) = ( ( ( WN(Z2)**2.0 / DELTA ) * FTILDE(Z2) ) + &
( DAIR / ( ZOFK(Z2) * KAPPA ) ) * ( SQRT( &
TLTN(ZI)**2.0 + TLTE(ZI)**2.0 ) / DAIR )**(1.5) ) &
* ( TLTN(ZI) ) / ( TLTE(ZI) * TAUX + &
TLTN(ZI) * TAUY )
UP(ZI) = UP1(ZI) * 0.5 + UP1(ZI-1) * 0.5
VP(ZI) = VP1(ZI) * 0.5 + VP1(ZI-1) * 0.5
UPROF(ZI) = UPROF(ZI-1) + UP(ZI) * ( ZOFK(Z2) - ZOFK(Z2+1) )
VPROF(ZI) = VPROF(ZI-1) + VP(ZI) * ( ZOFK(Z2) - ZOFK(Z2+1) )
ENDDO
!|---------------------------------------------------------------------|
!|----Iteration completion/checks--------------------------------------|
!|---------------------------------------------------------------------|
!ZI = ( KB + 1 )
! Now solving for 'ZWND' height wind
UPROF(NKT+1) = UPROF(NKT) + ( SQRT( SQRT( TAUY**2.0 + &
TAUX**2.0 ) / DAIR ) ) / KAPPA * TAUX &
/ SQRT( TAUY**2.0 +TAUX**2.0 ) * LOG( ZWND &
/ ZOFK(Z2) )
VPROF(NKT+1) = VPROF(NKT) + ( SQRT( SQRT( TAUY**2.0 + &
TAUX**2.0 ) / DAIR ) ) / KAPPA * TAUY &
/ SQRT( TAUY**2.0 +TAUX**2.0 ) * LOG( ZWND &
/ ZOFK(Z2) )
IF (ITER .EQ. 3) THEN
WND_1X = UPROF(NKT+1)
WND_1Y = VPROF(NKT+1)
ELSEIF (ITER .EQ. 2) THEN
WND_2X = UPROF(NKT+1)
WND_2Y = VPROF(NKT+1)
ELSEIF (ITER .EQ. 1) THEN
WND_3X = UPROF(NKT+1)
WND_3Y = VPROF(NKT+1)
ENDIF
!-------------------------------------+
! Guide for adding stability effects +
!-------------------------------------+
!Get Wind at top of wave boundary layer
! WND_TOP=SQRT(UPROF(KB)**2+VPROF(KB)**2)
! Get Wind Angle at top of wave boundary layer
! ANG_TOP=ATAN2(VPROF(KB),UPROF(KB))
! Stress and direction
! USTD = ATAN2(TAUY,TAUX)
! UST = SQRT( SQRT( TAUX**2 + TAUY**2 ) / DAIR)
! Calclate along (PA) and across (PE) wind components
! WND_PA=WND_TOP*COS(ANG_TOP-USTD)
! WND_PE=WND_TOP*SIN(ANG_TOP-USTD)
! Calculate cartesian aligned wind
! WND_PAx=WND_PA*cos(ustd)
! WND_PAy=WND_PA*sin(USTd)
!Calculate cartesion across wind
! WND_PEx=WND_PE*cos(ustd+pi/2.)
! WND_PEy=WND_PE*sin(ustd+pi/2.)
!----------------------------------------------------+
! If a non-neutral profile is used the effective z0 +
! should be computed. This z0 can then be used +
! with stability information to derive a Cd, which +
! can be used to project the along-stress wind to +
! the given height. +
! i.e.: Assume neutral inside WBL calculate Z0 +
! Z0=ZOFK(Z2)*EXP(-WND_PA*kappa/UST) +
! WND_PA=UST/SQRT(CDM) +
!----------------------------------------------------+
! WND_PAx=WND_PA*cos(ustd)
! WND_PAy=WND_PA*sin(USTd)
! IF (ITER .EQ. 3) THEN
! WND_1X = WND_PAx+WND_PEx
! WND_1Y = WND_PAy+WND_PEy
! ELSEIF (ITER .EQ. 2) THEN
! WND_2X = WND_PAx+WND_PEx
! WND_2Y = WND_PAy+WND_PEy
! ELSEIF (ITER .EQ. 1) THEN
! WND_3X = WND_PAx+WND_PEx
! WND_3Y = WND_PAy+WND_PEy
! ENDIF
ENDDO
ITERATION = ITERATION + 1
DIFU10XX = WND_3X - WND_1X
DIFU10YX = WND_3Y - WND_1Y
DIFU10XY = WND_2X - WND_1X
DIFU10YY = WND_2Y - WND_1Y
FD_A = DIFU10XX / DTX
FD_B = DIFU10XY / DTY
FD_C = DIFU10YX / DTX
FD_D = DIFU10YY / DTY
DWNDX = - WNDX + WND_1X
DWNDY = - WNDY + WND_1Y
UITV = ABS( DWNDX )
VITV = ABS( DWNDY )
CH = SQRT( UITV**2.0 + VITV**2.0 )
IF (CH .GT. 15.) THEN
APAR = 0.5 / ( FD_A * FD_D - FD_B * FD_C )
ELSE
APAR = 1.0 / ( FD_A * FD_D - FD_B * FD_C )
ENDIF
CK=4.
IF (((VITV/MAX(ABS(WNDY),CK) .GT. iter_thresh) .OR. &
(UITV/MAX(ABS(WNDX),CK) .GT. iter_thresh)) .AND. &
(ITERATION .LT. 2)) THEN
TAUNUX = TAUNUX - APAR * ( FD_D * DWNDX - FD_B * DWNDY )
TAUNUY = TAUNUY - APAR * ( -FD_C * DWNDX +FD_A * DWNDY )
ELSEIF (((VITV/MAX(ABS(WNDY),CK) .GT. iter_thresh) .OR. &
(UITV/MAX(ABS(WNDX),CK) .GT. iter_thresh)) .AND. &
(ITERATION .LT. 24)) THEN
iter_thresh = 0.001
TAUNUX = TAUNUX - APAR * ( FD_D * DWNDX - FD_B * DWNDY )
TAUNUY = TAUNUY - APAR * ( -FD_C * DWNDX +FD_A * DWNDY )
ELSEIF (((VITV/MAX(ABS(WNDY),CK) .GT. iter_thresh) .OR. &
(UITV/MAX(ABS(WNDX),CK) .GT. iter_thresh)) .AND. &
(ITERATION .LT. 26)) THEN
iter_thresh = 0.01
TAUNUX = TAUNUX - APAR * ( FD_D * DWNDX - FD_B * DWNDY )
TAUNUY = TAUNUY - APAR * ( -FD_C * DWNDX +FD_A * DWNDY )
ELSEIF (((VITV/MAX(ABS(WNDY),CK) .GT. iter_thresh) .OR. &
(UITV/MAX(ABS(WNDX),CK) .GT. iter_thresh)) .AND. &
(ITERATION .LT. 30)) THEN
iter_thresh = 0.05
TAUNUX = TAUNUX - APAR * ( FD_D * DWNDX - FD_B * DWNDY )
TAUNUY = TAUNUY - APAR * ( -FD_C * DWNDX +FD_A * DWNDY )
ELSEIF (ITERATION .GE. 30) THEN
write(*,*)'Attn: W3FLD1 not converged.'
write(*,*)' Wind (X/Y): ',WNDX,WNDY
IT_FLAG1 = .FALSE.
UST=-999
TAUNUX=0.
TAUNUY=0.
ELSEIF (((VITV/MAX(ABS(WNDY),CK) .LT. iter_thresh) .AND.&
(UITV/MAX(ABS(WNDX),CK) .LT. iter_thresh)) .AND. &
(ITERATION .GE. 2)) THEN
IT_FLAG1 = .FALSE.
ENDIF
! if taunu iteration is unstable try to reset with new guess...
if (.not.(cos(wnd_in_dir-atan2(taunuy,taunux)).ge.0.0)) then
TAUNUX = USTSM**2 * DAIR * wndx / wnd_in_mag*.95
TAUNUY = USTSM**2 * DAIR * wndy / wnd_in_mag*.95
endif
ENDDO
!|---------------------------------------------------------------------|
!|----Finish-----------------------------------------------------------|
!|---------------------------------------------------------------------|
USTD = ATAN2(TAUY,TAUX)
UST = SQRT( SQRT( TAUX**2 + TAUY**2 ) / DAIR)
! Get Z0 from aligned wind
WND_PA=wnd_in_mag*COS(wnd_in_dir-USTD)
Z0 = ZWND/exp(wnd_pa*kappa/ust)
CD = UST**2 / wnd_in_mag**2
IF (PRESENT(CHARN)) THEN
CHARN = 0.01/SQRT(SQRT( TAUNUX**2 + TAUNUY**2 )/(UST**2))
ENDIF
FSFL1=.not.((CD .LT. 0.01).AND.(CD .GT. 0.0001))
FSFL2=.not.(cos(wnd_in_dir-ustd).GT.0.9)
IF (FSFL1 .or. FSFL2) THEN
!Fail safe to bulk
write(*,*)'Attn: W3FLD1 failed, will output bulk...'
CALL wnd2z0m(wnd_in_mag,z0)
UST = wnd_in_mag*kappa/log(zwnd/z0)
USTD = wnd_in_dir
CD = UST**2 / wnd_in_mag**2
ENDIF
DEALLOCATE(WN, DWN, CP,SIG2, SPC2, TLTN, TLTE, TAUD, &
TLTND, TLTED, ZOFK, UPROF, &
VPROF, FTILDE, UP1, VP1, UP, VP, TLTNA, TLTEA)
!/ End of W3FLD1 ----------------------------------------------------- /
!/
RETURN
!
END SUBROUTINE W3FLD1
!/ ------------------------------------------------------------------- /
SUBROUTINE INFLD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | B. G. Reichl |
!/ | FORTRAN 90 |
!/ | Last update : 15-Jan-2016 |
!/ +-----------------------------------+
!/
!/ 15-Jan-2016 : Origination. ( version 5.12 )
!/
! 1. Purpose :
!
! Initialization for w3fld1 (also used by w3fld2)
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3FLDX Subr. W3FLDXMD Corresponding source term.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE W3ODATMD, ONLY: NDSE
USE W3GDATMD, ONLY: TAIL_ID, TAIL_LEV, TAIL_TRAN1, TAIL_TRAN2
USE W3SERVMD, ONLY: EXTCDE
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'INFLD')
#endif
!
! 1. .... ----------------------------------------------------------- *
!
Tail_Choice=Tail_ID
Tail_Level=TAIL_Lev
Tail_transition_ratio1 = TAIL_TRAN1
Tail_transition_ratio2 = TAIL_TRAN2
!
RETURN
!
! Formats
!
!/
!/ End of INFLD1 ----------------------------------------------------- /
!/
END SUBROUTINE INFLD
!/
!/ ------------------------------------------------------------------- /
SUBROUTINE APPENDTAIL(INSPC, WN2, NKT, KA1, KA2, KA3, WNDDIR,SAT)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | B. G. Reichl |
!/ | FORTRAN 90 |
!/ | Last update : 15-Jan-2016 |
!/ +-----------------------------------+
!/
!/ 15-Jan-2016 : Origination. ( version 5.12 )
!/
! 1. Purpose :
!
! Set tail for stress calculation.
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3FLD1 Subr. W3FLD1MD Corresponding source term.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS, ONLY: TPI, PI
USE W3GDATMD, ONLY: NTH, TH, DTH
USE W3ODATMD, ONLY: NDSE
USE W3SERVMD, ONLY: EXTCDE
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
INTEGER, INTENT(IN) :: NKT, KA1, KA2, KA3
REAL, INTENT(IN) :: WN2(NKT), WNDDIR,SAT
REAL, INTENT(INOUT) :: INSPC(NKT,NTH)
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
REAL :: BT(NKT), IC, ANGLE2, ANG(NKT),&
NORMSPC(NTH), AVG, ANGDIF, M, MAXANG, &
MAXAN, MINAN
INTEGER :: MAI, I, K, T
REAL, ALLOCATABLE, DIMENSION(:) :: ANGLE1
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'APPENDTAIL')
#endif
!
! 1. .... ----------------------------------------------------------- *
!
!|###############################################################|
!|##1. Get the level of the saturation spectrum in transition
!|## region A
!|###############################################################|
!-------------------------------------------
! 1a, get saturation level at KA1 (1.25xFPI)
!-------------------------------------------
BT(KA1) = 0
ANG = 0.0
DO T=1, NTH
BT(KA1)=BT(KA1)+INSPC(KA1,T)*WN2(KA1)**3.0*DTH
ENDDO
!-----------------------------------------------
! 1b, Set saturation level at KA2 (3xFPI) to SAT
!-----------------------------------------------
BT(KA2) = SAT
!-------------------------------------------------------------
! 1c, Find slope of saturation spectrum in transition region A
!-------------------------------------------------------------
M = ( BT(KA2) - BT(KA1) ) / ( WN2(KA2) - WN2(KA1) )
!----------------------------------------------------------------
! 1d, Find intercept of saturation spectrum in transition region
! A
!----------------------------------------------------------------
IC = BT(KA1) - M * WN2(KA1)
!------------------------------------------------------
! 1e, Calculate saturation level for all wavenumbers in
! transition region A
!------------------------------------------------------
DO K=KA1,KA2
BT(K)=M*WN2(K)+IC
ENDDO
!|###############################################################|
!|##2. Determine the directionality at each wavenumber in
!|## transition region B
!|###############################################################|
!-----------------------------------------------
! 2a, Find angle of spectral peak at KA2 (3xFPI)
!-----------------------------------------------
MAXANG = 0.0
DO T=1, NTH
IF (INSPC(KA2,T) .GT. MAXANG) THEN
MAXANG=INSPC(KA2,T)
ENDIF
ENDDO
!-------------------------------
! 2b, Check if peak spans 2 bins
!-------------------------------
!MAI = total number of angles of peak (if it spans more than 1)
MAI = 0
DO T=1, NTH
IF (MAXANG .EQ. INSPC(KA2,T)) THEN
MAI = MAI+1
ENDIF
ENDDO
!ANGLE1 = angles that correspond to peak (array)
MAI = MAX(1,MAI)
ALLOCATE(ANGLE1(MAI))
!-----------------------------------------------------
! 2c, If peak spans 2 or more bins it must be averaged
!-----------------------------------------------------
K=1
DO T=1, NTH
IF (MAXANG .EQ. INSPC(KA2,T)) THEN
ANGLE1(K) = TH(T)
K=K+1
ENDIF
ENDDO
DO K=1, MAI
DO WHILE (ANGLE1(K) .LT. 0.0)
ANGLE1(K) = ANGLE1(K) + TPI
ENDDO
DO WHILE (ANGLE1(K) .GE. TPI)
ANGLE1(K) = ANGLE1(K) - TPI
ENDDO
ENDDO
IF (MAI .GT. 1) THEN
MAXAN = ANGLE1(1)
MINAN = ANGLE1(1)
DO I=2, MAI
IF (MAXAN .LT. ANGLE1(I) )THEN
MAXAN = ANGLE1(I)
ENDIF
IF (MINAN .GT. ANGLE1(I) )THEN
MINAN = ANGLE1(I)
ENDIF
ENDDO
!------------------------------------------------------
! Need to distinguish if mean cross the origin (0/2pi)
!------------------------------------------------------
IF (MAXAN-MINAN .GT. PI) THEN
DO I=1, MAI
IF (MAXAN - ANGLE1(I) .GT. PI) THEN
ANGLE1(I) = ANGLE1(I) + TPI
ENDIF
ENDDO
ANGLE2=SUM(ANGLE1)/MAX(REAL(MAI),1.)
ELSE
ANGLE2=SUM(ANGLE1)/MAX(REAL(MAI),1.)
ENDIF
ELSE
ANGLE2=ANGLE1(1)
ENDIF
DO WHILE (ANGLE2 .LT. 0.0)
ANGLE2 = ANGLE2 + TPI
ENDDO
DO WHILE (ANGLE2 .GE. TPI)
ANGLE2 = ANGLE2 - TPI
ENDDO
!
!---------------------------------------------------
! This deals with angles that are less than 90
!---------------------------------------------------
if (cos(angle2-wnddir) .ge. 0.) then !Less than 90
m=asin(sin(wnddir-angle2))/(wn2(ka3)-wn2(ka2))
ic=angle2
do k=ka2, ka3
ang(k)=ic +m*(wn2(k)-wn2(ka2))
enddo
else
!----------------------------------------------------
! This deals with angels that turn clockwise
!----------------------------------------------------
if (sin(wnddir-angle2).GE.0) then
m=acos(cos(wnddir-angle2))/(wn2(ka3)-wn2(ka2))
ic=angle2
do k=ka2, ka3
ang(k)=ic+m*(wn2(k)-wn2(ka2))
enddo
else
!-----------------------------------------------------
! This deals with angels that cross counter-clockwise
!-----------------------------------------------------
m=acos(cos(wnddir-angle2))/(wn2(ka3)-wn2(ka2))
ic=angle2
do k=ka2, ka3
ang(k)=ic-m*(wn2(k)-wn2(ka2))
enddo
endif
endif
!----------------------------------------------
! Region A, Saturation level decreased linearly
! while direction is maintained
!----------------------------------------------
DO K=KA1, KA2-1
AVG=SUM(INSPC(K,:))/MAX(REAL(NTH),1.)
DO T=1,NTH
if (avg /= 0.0) then
INSPC(K,T)=BT(K)*INSPC(K,T)/TPI/(WN2(K)**3.0)/AVG
else
inspc(k,t) = 0.0
end if
ENDDO
ENDDO
!-----------------------------------------------------------
! Region B, Saturation level left flat while spectrum turned
! to direction of wind
!-----------------------------------------------------------
DO K = KA2, KA3
DO T=1, NTH
angdif=th(t)-ang(k)
IF (COS(ANGDIF) .GT. 0.0) THEN
NORMSPC(T) = COS(ANGDIF)**2.0
ELSE
NORMSPC(T)=0.0
ENDIF
ENDDO
AVG=SUM(NORMSPC)/MAX(REAL(NTH),1.)
DO T=1, NTH
if (avg /= 0.0) then
INSPC(K,T) = SAT * NORMSPC(T)/TPI/(WN2(K)**3.0)/AVG
else
inspc(k,t) = 0.0
end if
ENDDO
ENDDO
DO T=1, NTH
angdif=th(t)-wnddir
IF (COS(ANGDIF) .GT. 0.0) THEN
NORMSPC(T) = COS(ANGDIF)**2.0
ELSE
NORMSPC(T) = 0.0
ENDIF
ENDDO
AVG=SUM(NORMSPC)/MAX(REAL(NTH),1.)!1./4.
DO K=KA3+1, NKT
DO T=1, NTH
if (avg /= 0.0) then
INSPC(K,T)=NORMSPC(T)*(SAT)/TPI/(WN2(K)**3.0)/AVG
else
inspc(k,t) = 0.0
end if
ENDDO
ENDDO
DEALLOCATE(ANGLE1)
!
! Formats
!
!/
!/ End of APPENDTAIL ----------------------------------------------------- /
!/
RETURN
!
END SUBROUTINE APPENDTAIL
!/ ------------------------------------------------------------------- /
!/
!/ ------------------------------------------------------------------- /
SUBROUTINE SIG2WN(SIG,DEPTH,WN)
!/ ------------------------------------------------------------------- /
!Author: Brandon Reichl (GSO/URI)
!Origination : 2013
!Update : March - 18 - 2015
! : June -22 -2018 (XYC)
!Puropse : Convert from angular frequency to wavenumber
! using full gravity wave dispersion relation
! if tanh(kh)<0.99, otherwise uses deep-water
! approximation.
!NOTE: May be a better version internal to WW3 that can replace this.
! Improved by using newton's method for iteration.(2018)
!/ ------------------------------------------------------------------- /
!/
use constants, only: GRAV
!/
implicit none
!/
REAL,INTENT(IN) :: SIG,DEPTH
REAL,INTENT(OUT) :: WN
!/
real :: wn1,wn2 !,sig1,sig2,dsigdk
real :: fk, fk_slp
integer :: i
logical :: SWITCH
!/ ------------------------------------------------------------------- /
wn1=sig**2/GRAV
SWITCH=.true.
!/ Updated code with Newton's method by XYC:
if (tanh(wn1*depth) .LT. 0.99) then
do while (SWITCH)
fk=grav*wn1*tanh(wn1*depth) - sig**2
fk_slp = grav*tanh(wn1*depth) + grav*wn1*depth/(cosh(wn1*depth))**2
wn2=wn1 - fk/fk_slp
if (abs(wn2-wn1)/wn1 .LT. 0.0001 ) then
SWITCH = .FALSE.
else
wn1=wn2
endif
enddo
else
wn2=wn1
endif
WN=WN2
!/ END of update
!/
!/ Previous code by BR:
!/ ------------------------------------------------------------------- /
! wn1=sig**2/GRAV
! wn2=wn1+0.00001
! SWITCH=.true.
!/ ------------------------------------------------------------------- /
! if (tanh(wn1*depth).LT.0.99) then
! do i=1,5
! if (SWITCH) then
! sig1=sqrt(GRAV*wn1*tanh(wn1*depth))
! sig2=sqrt(GRAV*wn2*tanh(wn2*depth))
! if (sig1.lt.sig*.99999.or.sig1.gt.sig*1.00001) then
! dsigdk=(sig2-sig1)/(wn2-wn1)
! WN1=WN1+(SIG2-SIG1)/dsigdk
! wn2=wn1+wn1*0.00001
! else
! SWITCH = .FALSE.
! endif
! endif
! enddo
! endif
!/
! WN=WN1
!/
RETURN
END SUBROUTINE SIG2WN
!/ ------------------------------------------------------------------- /
!/
!/ ------------------------------------------------------------------- /
SUBROUTINE WND2Z0M( W10M , ZNOTM )
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | B. G. Reichl |
!/ | FORTRAN 90 |
!/ | Last update : 04-Aug-2016 |
!/ +-----------------------------------+
!/
!/ 09-Apr-2014 : Last Update. ( version 5.12 )
!/ 15-Aug-2016 : Updated for 2016 HWRF z0 ( J. Meixner )
!/
! 1. Purpose :
!
! Get bulk momentum z0 from 10-m wind.
! Bulk stress corresponds to 2015 GFDL Hurricane model
! Not published yet, but contact Brandon Reichl or
! Isaac Ginis (Univ. of Rhode Island) for further info
!
! 2. Method :
! This has now been updated for 2016 HWRF z0 using routines
! from HWRF znot_m_v1, Biju Thomas, 02/07/2014
! and znot_wind10m Weiguo Wang, 02/24/2016
!
! 3. Parameters :
! Name Unit Type Description
! ----------------------------------------------------------------
! W10M m/s input 10 m neutral wind [m/s]
! ZNOTM m output Roughness scale for momentum
! ----------------------------------------------------------------
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3FLD1 Subr. W3FLD1MD Corresponding source term.
! W3FLD2 Subr. W3FLD2MD Corresponding source term.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
REAL, INTENT(IN) :: W10M
REAL, INTENT(OUT):: ZNOTM
!Parameters from znot_m_v1
REAL, PARAMETER :: bs0 = -8.367276172397277e-12
REAL, PARAMETER :: bs1 = 1.7398510865876079e-09
REAL, PARAMETER :: bs2 = -1.331896578363359e-07
REAL, PARAMETER :: bs3 = 4.507055294438727e-06
REAL, PARAMETER :: bs4 = -6.508676881906914e-05
REAL, PARAMETER :: bs5 = 0.00044745137674732834
REAL, PARAMETER :: bs6 = -0.0010745704660847233
REAL, PARAMETER :: cf0 = 2.1151080765239772e-13
REAL, PARAMETER :: cf1 = -3.2260663894433345e-11
REAL, PARAMETER :: cf2 = -3.329705958751961e-10
REAL, PARAMETER :: cf3 = 1.7648562021709124e-07
REAL, PARAMETER :: cf4 = 7.107636825694182e-06
REAL, PARAMETER :: cf5 = -0.0013914681964973246
REAL, PARAMETER :: cf6 = 0.0406766967657759
!Variables from znot_wind10m
REAL :: Z10, U10,AAA,TMP
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
CALL STRACE (IENT, 'WND2Z0M')
#endif
!Values as set in znot_wind10m
Z10=10.0
U10=W10M
if (U10 > 85.0) U10=85.0
!Calculation of z0 as in znot_m_v1
IF ( U10 .LE. 5.0 ) THEN
ZNOTM = (0.0185 / 9.8*(7.59e-4*U10**2+2.46e-2*U10)**2)
ELSEIF (U10 .GT. 5.0 .AND. U10 .LT. 10.0) THEN
ZNOTM =.00000235*(U10**2 - 25 ) + 3.805129199617346e-05
ELSEIF ( U10 .GE. 10.0 .AND. U10 .LT. 60.0) THEN
ZNOTM = bs6 + bs5*U10 + bs4*U10**2 + bs3*U10**3 + bs2*U10**4 +&
bs1*U10**5 + bs0*U10**6
ELSE
ZNOTM = cf6 + cf5*U10 + cf4*U10**2 + cf3*U10**3 + cf2*U10**4 +&
cf1*U10**5 + cf0*U10**6
END IF
!Modifications as in znot_wind10m for icoef_sf=4
!for wind<20, cd similar to icoef=2 at 10m, then reduced
TMP=0.4*0.4/(ALOG(10.0/ZNOTM))**2 ! cd at zlev
AAA=0.75
IF (U10 < 20) THEN
AAA=0.99
ELSEIF(U10 < 45.0) THEN
AAA=0.99+(U10-20)*(0.75-0.99)/(45.0-20.0)
END IF
ZNOTM=Z10/EXP( SQRT(0.4*0.4/(TMP*AAA)) )
END SUBROUTINE WND2Z0M
!/ ------------------------------------------------------------------- /
SUBROUTINE WND2SAT(WND10,SAT)
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | B. G. Reichl |
!/ | FORTRAN 90 |
!/ | Last update : 04-Aug-2016 |
!/ +-----------------------------------+
!/
!/ 15-Jan-2016 : Origination. ( version 5.12 )
!/ 04-Aug-2016 : Updated for 2016 HWRF CD/U10 curve ( J. Meixner )
!/
! 1. Purpose :
!
! Gives level of saturation spectrum to produce
! equivalent Cd as in wnd2z0m (for neutral 10m wind)
! tuned for method of Reichl et al. 2014
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
!Input: WND10 - 10 m neutral wind [m/s]
!Output: SAT - Level 1-d saturation spectrum in tail [non-dim]
! ----------------------------------------------------------------
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3FLD1 Subr. W3FLD1MD Corresponding source term.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
REAL, INTENT(IN) :: WND10
REAL, INTENT(OUT) :: SAT
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
CALL STRACE (IENT, 'WND2SAT')
#endif
!/ Old HWRF 2015 and ST2
! SAT =0.000000000001237 * WND10**6 +&
! -0.000000000364155 * WND10**5 +&
! 0.000000037435015 * WND10**4 +&
! -0.000001424719473 * WND10**3 +&
! 0.000000471570975 * WND10**2 +&
! 0.000778467452178 * WND10**1 +&
! 0.002962335390055
!
! SAT values based on
! HWRF 2016 CD curve, created with fetch limited cases ST4 physics
IF (WND10<20.0) THEN
SAT = -0.000018541921682*WND10**2 &
+0.000751515452434*WND10 &
+0.002466529381004
ELSEIF (WND10<45) THEN
SAT = -0.000000009060349*WND10**4 &
+0.000001276678367*WND10**3 &
-0.000068274393789*WND10**2 &
+0.001418180888868*WND10 &
+0.000262277682984
ELSE
SAT = -0.000155976275073*WND10 &
+0.012027763023184
ENDIF
SAT = min(max(SAT,0.002),0.014)
END SUBROUTINE WND2SAT
!
!/ End of module W3FLD1MD -------------------------------------------- /
!/
END MODULE W3FLD1MD