!/ ------------------------------------------------------------------- /
      Module W3FLD1MD
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III      NOAA/NCEP/NOPP |
!/                  |           B. G. Reichl            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         18-Mar-2015 |
!/                  +-----------------------------------+
!/
!/    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
!/                  
!/
!/    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
      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
!/
      CONTAINS
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3FLD1( ASPC, FPI, WNDX,WNDY, ZWND,               & 
                    DEPTH, RIB, UST, USTD, Z0, TAUNUX,TAUNUY, CHARN)
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III      NOAA/NCEP/NOPP |
!/                  |           B. G. Reichl            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         18-Mar-2015 |
!/                  +-----------------------------------+
!/
!/    01-Jul-2013 : Origination.                        ( version 4.xx )
!/    18-Mar-2015 : Prepare for submission              ( version 5.12 )
!/
!  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)
!       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, DAIR, TPI, PI, KAPPA
      USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, DTH, XFR, TH
      USE W3ODATMD, ONLY: NDSE
      USE W3SERVMD, ONLY: EXTCDE
!/S      USE W3SERVMD, ONLY: STRACE
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      REAL, INTENT(IN)        :: ASPC(NSPEC), WNDX, WNDY,  & 
                                 ZWND, DEPTH, RIB, 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
!/S      INTEGER, SAVE           :: IENT = 0
      LOGICAL                 :: FSFL1,FSFL2, CRIT1, CRIT2
      LOGICAL                 :: IT_FLAG1, IT_FLAG2
      LOGICAL, SAVE           :: FIRST = .TRUE.
!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'W3FLD1')
!
! 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
!/S      USE W3SERVMD, ONLY: STRACE
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S      INTEGER, SAVE           :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'INFLD')
!
! 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
!/S      USE W3SERVMD, ONLY: STRACE
!/
      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
!/
!/S      INTEGER, SAVE           :: IENT = 0
      REAL                :: BT(NKT), IC, ANGLE2, ANG(NKT),&
                             NORMSPC(NTH), AVG, ANGDIF, M, MAXANG, &
                             MAXAN, MINAN
      INTEGER             :: MAI, I, K, T
      REAL, ALLOCATABLE, DIMENSION(:)  :: ANGLE1
!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'APPENDTAIL')
!
! 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
            INSPC(K,T)=BT(K)*INSPC(K,T)/TPI/(WN2(K)**3.0)/AVG
         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
          INSPC(K,T) = SAT * NORMSPC(T)/TPI/(WN2(K)**3.0)/AVG
        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
         INSPC(K,T)=NORMSPC(T)*(SAT)/TPI/(WN2(K)**3.0)/AVG
        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 :
!
!/ ------------------------------------------------------------------- /
!/S      USE W3SERVMD, ONLY: STRACE
!/
      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

!/S      INTEGER, SAVE           :: IENT = 0
!/S      CALL STRACE (IENT, 'WND2Z0M')

      !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 :
!
!/ ------------------------------------------------------------------- /
!/S      USE W3SERVMD, ONLY: STRACE
!/
        IMPLICIT NONE
!/
        REAL, INTENT(IN) :: WND10
        REAL, INTENT(OUT) :: SAT
!/S      INTEGER, SAVE           :: IENT = 0
!
!/S      CALL STRACE (IENT, 'WND2SAT')
!/ 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