Functions/Subroutines
subroutine	w3fld1 (ASPC, FPI, WNDX, WNDY, ZWND, DEPTH, RIB, DAIR, UST, USTD, Z0, TAUNUX, TAUNUY, CHARN)

subroutine	infld

subroutine	appendtail (INSPC, WN2, NKT, KA1, KA2, KA3, WNDDIR, SAT)

subroutine	sig2wn (SIG, DEPTH, WN)

subroutine	wnd2z0m (W10M, ZNOTM)

subroutine	wnd2sat (WND10, SAT)

Variables
integer, save	tail_choice

real, save	tail_level

real, save	tail_transition_ratio1

real, save	tail_transition_ratio2

Function/Subroutine Documentation

◆ appendtail()

subroutine w3fld1md::appendtail	(	real, dimension(nkt,nth), intent(inout)	INSPC,
		real, dimension(nkt), intent(in)	WN2,
		integer, intent(in)	NKT,
		integer, intent(in)	KA1,
		integer, intent(in)	KA2,
		integer, intent(in)	KA3,
		real, intent(in)	WNDDIR,
		real, intent(in)	SAT
	)

Definition at line 886 of file w3fld1md.F90.

     !/
     !/                  +-----------------------------------+
     !/                  | 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
     !

References w3gdatmd::dth, w3servmd::extcde(), w3odatmd::ndse, w3gdatmd::nth, constants::pi, w3servmd::strace(), w3gdatmd::th, and constants::tpi.

Referenced by w3fld1(), and w3fld2md::w3fld2().

◆ infld()

subroutine w3fld1md::infld

Definition at line 787 of file w3fld1md.F90.

     !/
     !/                  +-----------------------------------+
     !/                  | 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 ----------------------------------------------------- /
     !/

References w3servmd::extcde(), w3odatmd::ndse, w3servmd::strace(), tail_choice, w3gdatmd::tail_id, w3gdatmd::tail_lev, tail_level, w3gdatmd::tail_tran1, w3gdatmd::tail_tran2, tail_transition_ratio1, and tail_transition_ratio2.

Referenced by w3fld1(), and w3fld2md::w3fld2().

◆ sig2wn()

subroutine w3fld1md::sig2wn	(	real, intent(in)	SIG,
		real, intent(in)	DEPTH,
		real, intent(out)	WN
	)

Definition at line 1184 of file w3fld1md.F90.

     !/ ------------------------------------------------------------------- /
     !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

References constants::grav.

Referenced by w3fld1(), and w3fld2md::w3fld2().

◆ w3fld1()

subroutine w3fld1md::w3fld1	(	real, dimension(nspec), intent(in)	ASPC,
		real, intent(in)	FPI,
		real, intent(in)	WNDX,
		real, intent(in)	WNDY,
		real, intent(in)	ZWND,
		real, intent(in)	DEPTH,
		real, intent(in)	RIB,
		real, intent(in)	DAIR,
		real, intent(out)	UST,
		real, intent(out)	USTD,
		real, intent(out)	Z0,
		real, intent(inout)	TAUNUX,
		real, intent(inout)	TAUNUY,
		real, intent(out), optional	CHARN
	)

Definition at line 96 of file w3fld1md.F90.

     !/
     !/                  +-----------------------------------+
     !/                  | 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
           IF ((taudir+taudirb).NE.0.) THEN
               crit2=(abs(taudir-taudirb)*100.0/(taudir+taudirb)*0.5) .GT. 0.1
           ELSE
               crit2=.true.
           ENDIF
           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
     !

References appendtail(), w3gdatmd::dth, constants::dwat, w3servmd::extcde(), constants::grav, infld(), constants::kappa, w3odatmd::ndse, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, constants::pi, w3gdatmd::sig, sig2wn(), w3servmd::strace(), tail_choice, tail_level, tail_transition_ratio1, tail_transition_ratio2, w3gdatmd::th, constants::tpi, wnd2sat(), wnd2z0m(), and w3gdatmd::xfr.

Referenced by w3srcemd::w3srce().

◆ wnd2sat()

subroutine w3fld1md::wnd2sat	(	real, intent(in)	WND10,
		real, intent(out)	SAT
	)

Definition at line 1387 of file w3fld1md.F90.

     !/                  +-----------------------------------+
     !/                  | 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)

References w3servmd::strace().

Referenced by w3fld1().

◆ wnd2z0m()

subroutine w3fld1md::wnd2z0m	(	real, intent(in)	W10M,
		real, intent(out)	ZNOTM
	)

Definition at line 1261 of file w3fld1md.F90.

     !/                  +-----------------------------------+
     !/                  | 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)) )
  