WW3/w3fld1md_8F90_source.html

!/ ------------------------------------------------------------------- /

Module w3fld1md

  !/

  !/                  +-----------------------------------+

  !/                  | WAVEWATCH III      NOAA/NCEP/NOPP |

  !/                  |           B. G. Reichl            |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         22-Mar-2021 |

  !/                  +-----------------------------------+

  !/

  !/    01-Jul-2013 : Origination.                        ( version 4.xx )

  !/    18-Mar-2015 : Clean-up/prepare for distribution   ( version 5.12 )

  !/    15-Jan-2016 : Updates                             ( version 5.12 )

  !/                                                      ( B. G. Reichl )

  !/    27-Jul-2016 : Added Charnock output  (J.Meixner)  ( version 5.12 )

  !/    22-Jun-2018 : updated SIG2WN subroutine (X.Chen)  ( version 6.06 )

  !/                  modified the range of wind profile computation;

  !/                  corrected direction of the shortest waves

  !/    22-Mar-2021 : Consider DAIR a variable            ( version 7.13 )

  !/

  !/

  !/    Copyright 2009 National Weather Service (NWS),

  !/       National Oceanic and Atmospheric Administration.  All rights

  !/       reserved.  WAVEWATCH III is a trademark of the NWS.

  !/       No unauthorized use without permission.

  !/

  !  1. Purpose :

  !

  !     This Module contains routines to compute the wind stress vector

  !     from the wave spectrum, the wind vector, and the lower atmospheric

  !     stability (the form included here is for neutral conditions, but

  !     the structure needed to include stability is included in comments).

  !     The stress calculated via this subroutine is

  !     intended for coupling to serve as the boundary condition

  !     between the ocean and atmosphere, and (for now)

  !     and has no impact on the wave spectrum calculated.

  !     The calculation in w3fld1 is based on the method

  !     presented in Reichl, Hara, and Ginis (2014), "Sea State Dependence

  !     of the Wind Stress under Hurricane Winds."

  !

  !  2. Variables and types :

  !

  !     Not applicable.

  !

  !  3. Subroutines and functions :

  !

  !      Name       Type  Scope    Description

  !     ----------------------------------------------------------------

  !      W3FLD1     Subr. Public   Reichl et al. 2014 stress calculation

  !      INFLD1     Subr. Public   Corresponding initialization routine.

  !      APPENDTAIL Subr. Public   Modification of tail for calculation

  !      SIG2WN     Subr. Public   Depth-dependent dispersion relation

  !      WND2Z0M    Subr. Public   Wind to roughness length

  !     ----------------------------------------------------------------

  !

  !  4. Subroutines and functions used :

  !

  !      Name      Type  Module   Description

  !     ----------------------------------------------------------------

  !      STRACE    Subr. W3SERVMD Subroutine tracing.

  !     ----------------------------------------------------------------

  !

  !  5. Remarks :

  !

  !  6. Switches :

  !

  !     !/S  Enable subroutine tracing.

  !     !/

  !

  !  7. Source code :

  !/

  !/ ------------------------------------------------------------------- /

  !/

  !

  PUBLIC

  ! Tail_Choice: Chose the method to determine the level of the tail

  INTEGER, SAVE :: tail_choice

#ifdef W3_OMPG

  !$omp threadprivate(Tail_Choice)

#endif

  REAL, SAVE    :: tail_level !if Tail_Choice=0, tail is constant

  REAL, SAVE    :: tail_transition_ratio1! freq/fpi where tail

  !  adjustment begins

  REAL, SAVE    :: tail_transition_ratio2! freq/fpi where tail

  !  adjustment ends

#ifdef W3_OMPG

  !$omp threadprivate(Tail_Level)

  !$omp threadprivate(Tail_transition_ratio1,Tail_transition_ratio2)

#endif

  !/

CONTAINS

  !/ ------------------------------------------------------------------- /

  SUBROUTINE w3fld1( ASPC, FPI, WNDX,WNDY, ZWND,               &

       DEPTH, RIB, DAIR, UST, USTD, Z0,               &

       TAUNUX, TAUNUY, CHARN)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III      NOAA/NCEP/NOPP |

    !/                  |           B. G. Reichl            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         22-Mar-2021 |

    !/                  +-----------------------------------+

    !/

    !/    01-Jul-2013 : Origination.                        ( version 4.xx )

    !/    18-Mar-2015 : Prepare for submission              ( version 5.12 )

    !/    22-Mar-2021 : Consider DAIR a variable            ( version 7.13 )

    !/

    !  1. Purpose :

    !

    !     Diagnostic stress vector calculation from wave spectrum, lower

    !     atmosphere stability, and wind vector (at some given height).

    !     The height of wind vector is assumed to be within the constant

    !     stress layer.  These parameterizations are meant to be performed

    !     at wind speeds > 10 m/s, and may not converge for extremely young

    !     seas (i.e. starting from flat sea conditions).

    !

    !  2. Method :

    !     See Reichl et al. (2014).

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       ASPC    Real   I   1-D Wave action spectrum.

    !       FPI     Real   I   Peak input frequency.

    !       WNDX    Real   I   X-dir wind (assumed referenced to current)

    !       WNDY    Real   I   Y-dir wind (assumed referenced to current)

    !       ZWND    Real   I   Wind height.

    !       DEPTH   Real   I   Water depth.

    !       RIB     REAL   I   Bulk Richardson in lower atmosphere

    !                          (for determining stability in ABL to get

    !                          10 m neutral wind)

    !       DAIR    REAL   I   Air density

    !       TAUNUX  Real   0   X-dir viscous stress (guessed from prev.)

    !       TAUNUY  Real   0   Y-dir viscous stress (guessed from prev.)

    !       UST     Real   O   Friction velocity.

    !       USTD    Real   O   Direction of friction velocity.

    !       Z0      Real   O   Surface roughness length

    !       CHARN   Real   O,optional    Charnock parameter

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3ASIM    Subr. W3ASIMMD Air-sea interface module.

    !      W3EXPO    Subr.   N/A    Point output post-processor.

    !      GXEXPO    Subr.   N/A    GrADS point output post-processor.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE constants, ONLY: grav, dwat, tpi, pi, kappa

    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dth, xfr, th

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    REAL, INTENT(IN)        :: ASPC(NSPEC), WNDX, WNDY,  &

         ZWND, DEPTH, RIB, DAIR, FPI

    REAL, INTENT(OUT)       :: UST, USTD, Z0

    REAL, INTENT(OUT), OPTIONAL :: CHARN

    REAL, INTENT(INOUT)     :: TAUNUX, TAUNUY

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    REAL, PARAMETER         ::  NU=0.105/10000.0

    REAL, PARAMETER         ::  DELTA=0.03

    ! Commonly used parameters

    REAL                    ::  wnd_in_mag, wnd_in_dir

    !For Calculating Tail

    REAL                    ::  KMAX, KTAILA, KTAILB, KTAILC

    REAL                    ::  SAT, z01, z02, u10

    LOGICAL                 ::  ITERFLAG

    INTEGER                 ::  COUNT

    !For Iterations

    REAL                    ::  DTX, DTY, iter_thresh, &

         USTSM, Z0SM, Z1

    !For stress calculation

    REAL                    ::  WAGE, CBETA, BP, CD,       &

         USTRB, ANGDIF, USTAR, ZNU, &

         TAUT, TAUX, TAUY, BETAG, TAUDIR, &

         TAUDIRB

    !For wind profile calculation

    REAL                    ::  UPROFV, VPROFV

    !For wind profile iteration

    REAL                    ::  WND_1X, WND_1Y, &

         WND_2X, WND_2Y, &

         WND_3X, WND_3Y, &

         DIFU10XX, DIFU10YX, DIFU10XY, DIFU10YY, &

         FD_A, FD_B, FD_C, FD_D, &

         DWNDX, DWNDY, &

         APAR, CH,UITV, VITV,USTL,&

         CK

    !For adding stability to wind profile

    REAL                    ::  WND_TOP, ANG_TOP, WND_PA, WND_PE,   &

         WND_PEx, WND_PEy, WND_PAx, WND_PAy, &

         CDM

    INTEGER                 ::  NKT, K, T, Z2, ITER, ZI, ZII, &

         I, CTR, ITERATION, KA1, KA2, &

         KA3, KB

    ! For defining extended spectrum with appended tail.

    REAL, ALLOCATABLE, DIMENSION(:)   :: WN, DWN, CP,SIG2

    REAL, ALLOCATABLE, DIMENSION(:,:) :: SPC2

    REAL, ALLOCATABLE, DIMENSION(:)   :: TLTN, TLTE, TAUD, &

         TLTND, &

         TLTED, ZOFK, UPROF, VPROF, &

         FTILDE, UP1, VP1, UP, VP, &

         TLTNA, TLTEA

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    LOGICAL                 :: FSFL1,FSFL2, CRIT1, CRIT2

    LOGICAL                 :: IT_FLAG1, IT_FLAG2

    LOGICAL, SAVE           :: FIRST = .true.

#ifdef W3_OMPG

    !$omp threadprivate( FIRST)

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3FLD1')

#endif

    !

    ! 0.  Initializations ------------------------------------------------ *

    !

    !     **********************************************************

    !     ***    The initialization routine should include all   ***

    !     *** initialization, including reading data from files. ***

    !     **********************************************************

    !

    IF ( first ) THEN

      CALL infld

      first  = .false.

    END IF

    wnd_in_mag = sqrt( wndx**2 + wndy**2 )

    wnd_in_dir = atan2(wndy, wndx)

    !----------------------------------------------------------+

    ! Assume wind input is neutral 10 m wind.  If wind input   +

    ! is not 10 m, tail level will need to be calculated based +

    ! on esimation of 10 m wind.                               +

    !----------------------------------------------------------+

    u10 = wnd_in_mag

    ! - Get tail level

    if (tail_choice.eq.0) then

      sat=tail_level

    elseif (tail_choice.eq.1) then

      CALL wnd2sat(u10,sat)

    endif

    !

    ! 1.  Attach Tail ---------------------------------------------------- *

    !

    ! If the depth remains constant, the allocation could be limited to the

    !   first time step.  Since this code is designed for coupled

    !   implementation where the water level can change, I keep it the

    !   allocation on every time step.  When computational efficiency is

    !   important, this process may be rethought.

    !

    ! i. Find maximum wavenumber of input spectrum

    call sig2wn(sig(nk),depth,kmax)

    nkt = nk

    ! ii. Find additional wavenumber bins to extended to cm scale waves

    DO WHILE ( kmax .LT. 366.0 )

      nkt = nkt + 1

      kmax = ( kmax * xfr**2 )

    enddo!K<366

    ! iii. Allocate new "extended" spectrum

    ALLOCATE( wn(nkt), dwn(nkt), cp(nkt), sig2(nkt),spc2(nkt,nth), &

         tltn(nkt), tlte(nkt), taud(nkt), &

         tltnd(nkt), tlted(nkt), zofk(nkt), uprof(nkt+1),&

         vprof(nkt+1), ftilde(nkt), up1(nkt+1),vp1(nkt+1), &

         up(nkt+1), vp(nkt+1), tltna(nkt),tltea(nkt))

    !

    ! 1a. Build Discrete Wavenumbers for defining extended spectrum on---- *

    !

    !i. Copy existing sig to extended sig2, calculate phase speed.

    DO k = 1, nk !existing spectrum

      call sig2wn(sig(k),depth,wn(k))

      cp(k) = ( sig(k) / wn(k) )

      sig2(k) = sig(k)

    enddo!K

    !ii. Calculate extended sig2 and phase speed.

    DO k = ( nk + 1 ), ( nkt) !extension

      sig2(k) = sig2(k-1) *xfr

      call sig2wn(sig2(k),depth,wn(k))

      cp(k) = sig2(k) / wn(k)

    enddo!K

    !iii. Calculate dk's for integrations.

    DO k = 1, nkt-1

      dwn(k) = wn(k+1) - wn(k)

    ENDDO

    dwn(nkt) = wn(nkt)*xfr**2 - wn(nkt)

    !

    ! 1b. Attach initial tail--------------------------------------------- *

    !

    !i. Convert action spectrum to variance spectrum

    !   SPC(k,theta) = A(k,theta) * sig(k)

    ! This could be redone for computational efficiency

    i=0

    DO k=1, nk

      DO t=1, nth

        i = i + 1

        spc2(k,t) = aspc(i)  * sig(k)

      enddo!T

    enddo!K

    !ii. Extend k^-3 tail to extended spectrum

    DO k=nk+1, nkt

      DO t=1, nth

        spc2(k,t)=spc2(nk,t)*wn(nk)**3.0/wn(k)**(3.0)

      enddo!T

    enddo!K

    !

    ! 1c. Calculate transitions for new (constant saturation ) tail ------ *

    !

    !

    !i. Find wavenumber for beginning spc level transition to tail

    call sig2wn (fpi*tpi*tail_transition_ratio1,depth,ktaila )

    !ii. Find wavenumber for ending spc level transition to tail

    call sig2wn (fpi*tpi*tail_transition_ratio2,depth,ktailb )

    !iii. Find wavenumber for ending spc direction transition to tail

    ktailc= ktailb * 2.0

    !iv. Find corresponding indices of wavenumber transitions

    ka1 = 2     ! Do not modify 1st wavenumber bin

    DO WHILE ( ( ktaila .GE. wn(ka1) ) .AND. (ka1 .LT. nkt-6) )

      ka1 = ka1 + 1

    ENDDO

    ka2 = ka1+2

    DO WHILE ( ( ktailb .GE. wn(ka2) ) .AND. (ka2 .LT. nkt-4) )

      ka2 = ka2 + 1

    ENDDO

    ka3 = ka2+2

    DO WHILE ( ( ktailc .GE. wn(ka3)) .AND. (ka3 .LT. nkt-2) )

      ka3 = ka3 + 1

    ENDDO

    !v. Call subroutine to perform actually tail truncation

    ! only if there is some energy in spectrum

    CALL appendtail(spc2, wn, nkt, ka1, ka2, ka3,&

         wnd_in_dir, sat)

    ! Spectrum is now set for stress integration

    !

    ! 2.  Prepare for iterative calculation of wave-form stress----------- *

    !

    dtx = 0.00005

    dty = 0.00005

    iter_thresh = 0.001

    !

    ! 2a. Calculate initial guess for viscous stress from smooth-law------ *

    ! (Would be preferable to use prev. step)

    !

    z0sm = 0.001  !Guess

    it_flag1 = .true.

    iteration = 0

    DO WHILE( it_flag1 )

      iteration = iteration + 1

      z1 = z0sm

      ustsm = kappa * wnd_in_mag / ( log( zwnd / z1 ) )

      z0sm = 0.132 * nu / ustsm

      IF ( (abs( z0sm - z1 ) .LT. 10.0**(-6)) .OR.&

           ( iteration .GT. 5 )) THEN

        it_flag1 = .false.

      ENDIF

    ENDDO


    iteration = 1

    ! Guessed values of viscous stress

    taunux = ustsm**2 * dair * wndx / wnd_in_mag

    taunuy = ustsm**2 * dair * wndy / wnd_in_mag

    !

    ! 3.  Enter iterative calculation of wave form/skin stress----------  *

    !

    it_flag1 = .true.

    DO WHILE (it_flag1)

      DO iter=1, 3 !3 loops for TAUNU iteration

        z2 = nkt

        ! First : TAUNUX + DX

        IF (iter .EQ. 1) THEN

          taunux = taunux + dtx

          ! Second : TAUNUY + DY

        ELSEIF (iter .EQ. 2) THEN

          taunux = taunux - dtx

          taunuy = taunuy + dty

          ! Third : unmodified

        ELSEIF (iter .EQ. 3) THEN

          taunuy = taunuy - dty

        ENDIF

        ! Near surface turbulent stress = taunu

        tltn(1) = taunuy

        tlte(1) = taunux


        CALL appendtail(spc2, wn, nkt, ka1, ka2, ka3,&

             atan2(taunuy,taunux), sat)

        !|---------------------------------------------------------------------|

        !|-----Calculate first guess at growth rate and local turbulent stress-|

        !|-----for integration as a function of wavedirection------------------|

        !|---------------------------------------------------------------------|

        DO zi = 2, nkt

          ustl=0.0

          tltnd(zi)=0.0

          tlted(zi)=0.0

          z2 = z2 - 1

          ! Use value of prev. wavenumber/height

          taud(zi) = atan2( tltn(zi-1), tlte(zi-1))

          ustl = sqrt( sqrt( tltn(zi-1)**2 + tlte(zi-1)**2 )/ dair )

          DO t = 1, nth

            angdif=taud(zi)-th(t) !stress/wave angle

            IF ( cos( angdif ) .GE. 0.0 ) THEN !Waves aligned

              wage = cp(z2) / (ustl)

              ! First, waves much slower than wind.

              IF ( wage .LT. 10. ) THEN

                cbeta = 25.0

                ! Transition from waves slower than wind to faster

              ELSEIF ( ( wage .GE. 10.0 ) .AND. &

                   ( wage .LE. 25.0 ) ) THEN

                cbeta = 10.0 + 15.0 * cos( pi * ( wage - 10.0 ) &

                     / 15.0 )

                ! Waves faster than wind

              ELSEIF ( wage .GT. 25.0 ) THEN

                cbeta = -5.0

              ENDIF

              ! Waves opposing wind

            ELSE

              cbeta = -25.0

            ENDIF

            !Integrate turbulent stress

            tltnd(zi) =tltnd(zi)+( sin( th(t) ) * cos( angdif )**2)&

                 * cbeta * spc2(z2,t) * &

                 sqrt( tlte(zi-1)**2 + tltn(zi-1)**2.0 ) &

                 * ( wn(z2)**2.0 )*dth

            tlted(zi) = tlted(zi)+(cos( th(t) ) * cos( angdif )**2)&

                 * cbeta * spc2(z2,t) * &

                 sqrt( tlte(zi-1)**2 + tltn(zi-1)**2.0 ) &

                 * ( wn(z2)**2.0 )*dth

          ENDDO

          !|---------------------------------------------------------------------|

          !|-----Complete the integrations---------------------------------------|

          !|---------------------------------------------------------------------|

          IF (zi .EQ. 2) THEN

            !First turbulent stress bin above taunu

            tltna(zi) = tltnd(zi) * dwn(z2) * 0.5

            tltea(zi) = tlted(zi) * dwn(z2) * 0.5

          ELSE

            tltna(zi)=(tltnd(zi)+tltnd(zi-1))*0.5*dwn(z2)

            tltea(zi)=(tlted(zi)+tlted(zi-1))*0.5*dwn(z2)

          ENDIF

          tltn(zi)=tltn(zi-1)+tltna(zi)

          tlte(zi)=tlte(zi-1)+tltea(zi)

        ENDDO

        tauy=tltn(nkt)

        taux=tlte(nkt)

        ! This is the first guess at the stress.

        !|---------------------------------------------------------------------|

        !|----Iterate til convergence------------------------------------------|

        !|---------------------------------------------------------------------|

        ustrb=sqrt(sqrt(tauy**2.0+taux**2.0)/dair)

        taudirb=atan2(tauy,taux)

        it_flag2 = .true.

        ctr=1

        DO WHILE ( (it_flag2) .AND. ( ctr .LT. 10 ) )

          z2=nkt+1

          DO zi=1, nkt

            z2=z2-1

            ustl=0.0

            tlted(zi)=0.0

            tltnd(zi)=0.0

            ftilde(z2)=0.0

            taud(zi) = atan2(tltn(zi),tlte(zi))

            ustl = sqrt(sqrt(tltn(zi)**2+tlte(zi)**2)/dair)

            DO t=1, nth

              betag=0.0

              angdif = taud(zi)-th(t)

              IF ( cos( angdif ) .GE. 0.0 ) THEN

                wage = cp(z2)  / (ustl)

                IF ( wage .LT. 10 ) THEN

                  cbeta = 25.0

                ELSEIF ( ( wage .GE. 10.0 ) .AND. &

                     ( wage .LE. 25.0 ) ) THEN

                  cbeta = 10.0 + 15.0 * cos( pi * ( wage - 10.0 ) &

                       / 15.0 )

                ELSEIF ( wage .GT. 25.0 ) THEN

                  cbeta = -5.0

                ENDIF

              ELSE

                cbeta = -25.0

              ENDIF

              bp = sqrt( (cos( th(t) ) * cos( angdif )**2.0)**2.0 &

                   + (sin( th(t) ) * cos( angdif )**2.0)**2.0 )

              betag=bp*cbeta*sqrt(tlte(zi)**2.0+tltn(zi)**2.0) &

                   /(dwat)*sig2(z2)/cp(z2)**2

              ftilde(z2) = ftilde(z2) + betag * dwat * grav &

                   * spc2(z2,t) * dth

              tltnd(zi) =tltnd(zi)+ (sin( th(t) ) * cos( angdif )**2.0)&

                   * cbeta * spc2(z2,t) * sqrt( &

                   tlte(zi)**2.0 + tltn(zi)**2.0 ) * &

                   ( wn(z2)**2.0 )*dth

              tlted(zi) = tlted(zi)+(cos( th(t) ) * cos( angdif )**2.0)&

                   * cbeta * spc2(z2,t) * sqrt( &

                   tlte(zi)**2.0 + tltn(zi)**2.0 ) * &

                   ( wn(z2)**2.0 )*dth

            ENDDO

            IF (zi .EQ. 1) THEN

              tltna(zi)=tltnd(zi)*dwn(z2)*0.5

              tltea(zi)=tlted(zi)*dwn(z2)*0.5

            ELSE

              tltna(zi)=(tltnd(zi)+tltnd(zi-1))*0.5*dwn(z2)

              tltea(zi)=(tlted(zi)+tlted(zi-1))*0.5*dwn(z2)

            ENDIF

            IF (zi.GT.1) then

              tltn(zi)=tltn(zi-1)+tltna(zi)

              tlte(zi)=tlte(zi-1)+tltea(zi)

            else

              tltn(zi)=taunuy+tltna(zi)

              tlte(zi)=taunux+tltea(zi)

            endif

          ENDDO

          tauy=tltn(nkt) !by NKT full stress is entirely

          taux=tlte(nkt) !from turbulent stress

          taut=sqrt(tauy**2.0+taux**2.0)

          ustar=sqrt(sqrt(tauy**2.0+taux**2.0)/dair)

          taudir=atan2(tauy, taux)

          ! Note: add another criterion (stress direction) for iteration.

          crit1=(abs(ustar-ustrb)*100.0)/((ustar+ustrb)*0.5) .GT. 0.1

          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

    !

  END SUBROUTINE w3fld1

  !/ ------------------------------------------------------------------- /

  SUBROUTINE infld

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           B. G. Reichl            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         15-Jan-2016 |

    !/                  +-----------------------------------+

    !/

    !/    15-Jan-2016 : Origination.                        ( version 5.12 )

    !/

    !  1. Purpose :

    !

    !     Initialization for w3fld1 (also used by w3fld2)

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3FLDX    Subr. W3FLDXMD Corresponding source term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE w3odatmd, ONLY: ndse

    USE w3gdatmd, ONLY: tail_id, tail_lev, tail_tran1, tail_tran2

    USE w3servmd, ONLY: extcde

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'INFLD')

#endif

    !

    ! 1.  .... ----------------------------------------------------------- *

    !

    tail_choice=tail_id

    tail_level=tail_lev

    tail_transition_ratio1 = tail_tran1

    tail_transition_ratio2 = tail_tran2


    !

    RETURN

    !

    ! Formats

    !


    !/

    !/ End of INFLD1 ----------------------------------------------------- /

    !/

  END SUBROUTINE infld

  !/

  !/ ------------------------------------------------------------------- /

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

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           B. G. Reichl            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         15-Jan-2016 |

    !/                  +-----------------------------------+

    !/

    !/    15-Jan-2016 : Origination.                        ( version 5.12 )

    !/

    !  1. Purpose :

    !

    !     Set tail for stress calculation.

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3FLD1    Subr. W3FLD1MD Corresponding source term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE constants, ONLY: tpi, pi

    USE w3gdatmd, ONLY: nth, th, dth

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    INTEGER, INTENT(IN) :: NKT, KA1, KA2, KA3

    REAL, INTENT(IN)    :: WN2(NKT), WNDDIR,SAT

    REAL, INTENT(INOUT)   :: INSPC(NKT,NTH)

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    REAL                :: BT(NKT), IC, ANGLE2, ANG(NKT),&

         NORMSPC(NTH), AVG, ANGDIF, M, MAXANG, &

         MAXAN, MINAN

    INTEGER             :: MAI, I, K, T

    REAL, ALLOCATABLE, DIMENSION(:)  :: ANGLE1

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'APPENDTAIL')

#endif

    !

    ! 1.  .... ----------------------------------------------------------- *

    !

    !|###############################################################|

    !|##1. Get the level of the saturation spectrum in transition

    !|##   region A

    !|###############################################################|

    !-------------------------------------------

    ! 1a, get saturation level at KA1 (1.25xFPI)

    !-------------------------------------------

    bt(ka1) = 0

    ang = 0.0

    DO t=1, nth

      bt(ka1)=bt(ka1)+inspc(ka1,t)*wn2(ka1)**3.0*dth

    ENDDO

    !-----------------------------------------------

    ! 1b, Set saturation level at KA2 (3xFPI) to SAT

    !-----------------------------------------------

    bt(ka2) = sat

    !-------------------------------------------------------------

    ! 1c, Find slope of saturation spectrum in transition region A

    !-------------------------------------------------------------

    m = ( bt(ka2) - bt(ka1) ) / ( wn2(ka2) - wn2(ka1) )

    !----------------------------------------------------------------

    ! 1d, Find intercept of saturation spectrum in transition region

    !     A

    !----------------------------------------------------------------

    ic = bt(ka1) - m * wn2(ka1)

    !------------------------------------------------------

    ! 1e, Calculate saturation level for all wavenumbers in

    !     transition region A

    !------------------------------------------------------

    DO k=ka1,ka2

      bt(k)=m*wn2(k)+ic

    ENDDO

    !|###############################################################|

    !|##2. Determine the directionality at each wavenumber in

    !|##   transition region B

    !|###############################################################|

    !-----------------------------------------------

    ! 2a, Find angle of spectral peak at KA2 (3xFPI)

    !-----------------------------------------------

    maxang = 0.0

    DO t=1, nth

      IF (inspc(ka2,t) .GT. maxang) THEN

        maxang=inspc(ka2,t)

      ENDIF

    ENDDO

    !-------------------------------

    ! 2b, Check if peak spans 2 bins

    !-------------------------------

    !MAI = total number of angles of peak (if it spans more than 1)

    mai = 0

    DO t=1, nth

      IF (maxang .EQ. inspc(ka2,t)) THEN

        mai = mai+1

      ENDIF

    ENDDO

    !ANGLE1 = angles that correspond to peak (array)

    mai = max(1,mai)

    ALLOCATE(angle1(mai))

    !-----------------------------------------------------

    ! 2c, If peak spans 2 or more bins it must be averaged

    !-----------------------------------------------------

    k=1

    DO t=1, nth

      IF (maxang .EQ. inspc(ka2,t)) THEN

        angle1(k) = th(t)

        k=k+1

      ENDIF

    ENDDO

    DO k=1, mai

      DO WHILE (angle1(k) .LT. 0.0)

        angle1(k) = angle1(k) + tpi

      ENDDO

      DO WHILE (angle1(k) .GE. tpi)

        angle1(k) = angle1(k) - tpi

      ENDDO

    ENDDO

    IF (mai .GT. 1) THEN

      maxan = angle1(1)

      minan = angle1(1)

      DO i=2, mai

        IF (maxan .LT. angle1(i) )THEN

          maxan = angle1(i)

        ENDIF

        IF (minan .GT. angle1(i) )THEN

          minan = angle1(i)

        ENDIF

      ENDDO

      !------------------------------------------------------

      !  Need to distinguish if mean cross the origin (0/2pi)

      !------------------------------------------------------

      IF (maxan-minan .GT. pi) THEN

        DO i=1, mai

          IF (maxan - angle1(i) .GT. pi) THEN

            angle1(i) = angle1(i) + tpi

          ENDIF

        ENDDO

        angle2=sum(angle1)/max(real(mai),1.)

      ELSE

        angle2=sum(angle1)/max(real(mai),1.)

      ENDIF

    ELSE

      angle2=angle1(1)

    ENDIF

    DO WHILE (angle2 .LT. 0.0)

      angle2 = angle2 + tpi

    ENDDO

    DO WHILE (angle2 .GE. tpi)

      angle2 = angle2 - tpi

    ENDDO

    !

    !---------------------------------------------------

    ! This deals with angles that are less than 90

    !---------------------------------------------------

    if (cos(angle2-wnddir) .ge. 0.) then  !Less than 90

      m=asin(sin(wnddir-angle2))/(wn2(ka3)-wn2(ka2))

      ic=angle2

      do k=ka2, ka3

        ang(k)=ic +m*(wn2(k)-wn2(ka2))

      enddo

    else

      !----------------------------------------------------

      !  This deals with angels that turn clockwise

      !----------------------------------------------------

      if (sin(wnddir-angle2).GE.0) then

        m=acos(cos(wnddir-angle2))/(wn2(ka3)-wn2(ka2))

        ic=angle2

        do k=ka2, ka3

          ang(k)=ic+m*(wn2(k)-wn2(ka2))

        enddo

      else

        !-----------------------------------------------------

        !  This deals with angels that cross counter-clockwise

        !-----------------------------------------------------

        m=acos(cos(wnddir-angle2))/(wn2(ka3)-wn2(ka2))

        ic=angle2

        do k=ka2, ka3

          ang(k)=ic-m*(wn2(k)-wn2(ka2))

        enddo

      endif

    endif

    !----------------------------------------------

    ! Region A, Saturation level decreased linearly

    ! while direction is maintained

    !----------------------------------------------

    DO k=ka1, ka2-1

      avg=sum(inspc(k,:))/max(real(nth),1.)

      DO t=1,nth

        if (avg /= 0.0) then

          inspc(k,t)=bt(k)*inspc(k,t)/tpi/(wn2(k)**3.0)/avg

        else

          inspc(k,t) = 0.0

        end if

      ENDDO

    ENDDO

    !-----------------------------------------------------------

    ! Region B, Saturation level left flat while spectrum turned

    ! to direction of wind

    !-----------------------------------------------------------

    DO k = ka2, ka3

      DO t=1, nth

        angdif=th(t)-ang(k)

        IF (cos(angdif) .GT. 0.0) THEN

          normspc(t) = cos(angdif)**2.0

        ELSE

          normspc(t)=0.0

        ENDIF

      ENDDO

      avg=sum(normspc)/max(real(nth),1.)

      DO t=1, nth

        if (avg /= 0.0) then

          inspc(k,t) = sat * normspc(t)/tpi/(wn2(k)**3.0)/avg

        else

          inspc(k,t) = 0.0

        end if

      ENDDO

    ENDDO

    DO t=1, nth

      angdif=th(t)-wnddir

      IF (cos(angdif) .GT. 0.0) THEN

        normspc(t) = cos(angdif)**2.0

      ELSE

        normspc(t) = 0.0

      ENDIF

    ENDDO

    avg=sum(normspc)/max(real(nth),1.)!1./4.

    DO k=ka3+1, nkt

      DO t=1, nth

        if (avg /= 0.0) then

          inspc(k,t)=normspc(t)*(sat)/tpi/(wn2(k)**3.0)/avg

        else

          inspc(k,t) = 0.0

        end if

      ENDDO

    ENDDO

    DEALLOCATE(angle1)

    !

    ! Formats

    !

    !/

    !/ End of APPENDTAIL ----------------------------------------------------- /

    !/

    RETURN

    !

  END SUBROUTINE appendtail

  !/ ------------------------------------------------------------------- /

  !/

  !/ ------------------------------------------------------------------- /

  SUBROUTINE sig2wn(SIG,DEPTH,WN)

    !/ ------------------------------------------------------------------- /

    !Author: Brandon Reichl (GSO/URI)

    !Origination  : 2013

    !Update       : March - 18 - 2015

    !             : June -22 -2018  (XYC)

    !Puropse      : Convert from angular frequency to wavenumber

    !               using full gravity wave dispersion relation

    !               if tanh(kh)<0.99, otherwise uses deep-water

    !               approximation.

    !NOTE: May be a better version internal to WW3 that can replace this.

    !      Improved by using newton's method for iteration.(2018)

    !/ ------------------------------------------------------------------- /

    !/

    use constants, only: grav

    !/

    implicit none

    !/

    REAL,INTENT(IN)    :: SIG,DEPTH

    REAL,INTENT(OUT)   :: WN

    !/

    real    :: wn1,wn2 !,sig1,sig2,dsigdk

    real    :: fk, fk_slp

    integer :: i

    logical :: SWITCH

    !/ ------------------------------------------------------------------- /

    wn1=sig**2/grav

    switch=.true.

    !/ Updated code with Newton's method by XYC:

    if (tanh(wn1*depth) .LT. 0.99) then

      do while (switch)

        fk=grav*wn1*tanh(wn1*depth) - sig**2

        fk_slp = grav*tanh(wn1*depth) + grav*wn1*depth/(cosh(wn1*depth))**2


        wn2=wn1 - fk/fk_slp


        if (abs(wn2-wn1)/wn1 .LT. 0.0001 ) then

          switch = .false.

        else

          wn1=wn2

        endif

      enddo

    else

      wn2=wn1

    endif

    wn=wn2

    !/ END of update

    !/

    !/ Previous code by BR:

    !/ ------------------------------------------------------------------- /

    !        wn1=sig**2/GRAV

    !        wn2=wn1+0.00001

    !        SWITCH=.true.

    !/ ------------------------------------------------------------------- /

    !        if (tanh(wn1*depth).LT.0.99) then

    !           do i=1,5

    !              if (SWITCH) then

    !                 sig1=sqrt(GRAV*wn1*tanh(wn1*depth))

    !                 sig2=sqrt(GRAV*wn2*tanh(wn2*depth))

    !                 if (sig1.lt.sig*.99999.or.sig1.gt.sig*1.00001) then

    !                    dsigdk=(sig2-sig1)/(wn2-wn1)

    !                    WN1=WN1+(SIG2-SIG1)/dsigdk

    !                    wn2=wn1+wn1*0.00001

    !                 else

    !                    SWITCH = .FALSE.

    !                 endif

    !              endif

    !           enddo

    !        endif

    !/

    !        WN=WN1

    !/

    RETURN

  END SUBROUTINE sig2wn

  !/ ------------------------------------------------------------------- /

  !/

  !/ ------------------------------------------------------------------- /

  SUBROUTINE wnd2z0m( W10M , ZNOTM )

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           B. G. Reichl            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         04-Aug-2016 |

    !/                  +-----------------------------------+

    !/

    !/    09-Apr-2014 : Last Update.                        ( version 5.12 )

    !/    15-Aug-2016 : Updated for 2016 HWRF z0            ( J. Meixner   )

    !/

    !  1. Purpose :

    !

    !     Get bulk momentum z0 from 10-m wind.

    !          Bulk stress corresponds to 2015 GFDL Hurricane model

    !          Not published yet, but contact Brandon Reichl or

    !          Isaac Ginis (Univ. of Rhode Island) for further info

    !

    !  2. Method :

    !          This has now been updated for 2016 HWRF z0 using routines

    !          from HWRF  znot_m_v1, Biju Thomas, 02/07/2014

    !           and       znot_wind10m Weiguo Wang, 02/24/2016

    !

    !  3. Parameters :

    !       Name  Unit  Type      Description

    !     ----------------------------------------------------------------

    !       W10M   m/s  input    10 m neutral wind [m/s]

    !       ZNOTM  m    output   Roughness scale for momentum

    !     ----------------------------------------------------------------

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3FLD1    Subr. W3FLD1MD Corresponding source term.

    !      W3FLD2    Subr. W3FLD2MD Corresponding source term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    REAL, INTENT(IN) :: W10M

    REAL, INTENT(OUT):: ZNOTM


    !Parameters from znot_m_v1

    REAL, PARAMETER :: bs0 = -8.367276172397277e-12

    REAL, PARAMETER :: bs1 = 1.7398510865876079e-09

    REAL, PARAMETER :: bs2 = -1.331896578363359e-07

    REAL, PARAMETER :: bs3 = 4.507055294438727e-06

    REAL, PARAMETER :: bs4 = -6.508676881906914e-05

    REAL, PARAMETER :: bs5 = 0.00044745137674732834

    REAL, PARAMETER :: bs6 = -0.0010745704660847233


    REAL, PARAMETER :: cf0 = 2.1151080765239772e-13

    REAL, PARAMETER :: cf1 = -3.2260663894433345e-11

    REAL, PARAMETER :: cf2 = -3.329705958751961e-10

    REAL, PARAMETER :: cf3 = 1.7648562021709124e-07

    REAL, PARAMETER :: cf4 = 7.107636825694182e-06

    REAL, PARAMETER :: cf5 = -0.0013914681964973246

    REAL, PARAMETER :: cf6 = 0.0406766967657759


    !Variables from znot_wind10m

    REAL            :: Z10, U10,AAA,TMP


#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

    CALL strace (ient, 'WND2Z0M')

#endif


    !Values as set in znot_wind10m

    z10=10.0

    u10=w10m

    if (u10 > 85.0) u10=85.0


    !Calculation of z0 as in znot_m_v1

    IF ( u10 .LE. 5.0 ) THEN

      znotm = (0.0185 / 9.8*(7.59e-4*u10**2+2.46e-2*u10)**2)

    ELSEIF (u10 .GT. 5.0 .AND. u10 .LT. 10.0) THEN

      znotm =.00000235*(u10**2 - 25 ) + 3.805129199617346e-05

    ELSEIF ( u10 .GE. 10.0  .AND. u10 .LT. 60.0) THEN

      znotm = bs6 + bs5*u10 + bs4*u10**2 + bs3*u10**3 + bs2*u10**4 +&

           bs1*u10**5 + bs0*u10**6

    ELSE

      znotm = cf6 + cf5*u10 + cf4*u10**2 + cf3*u10**3 + cf2*u10**4 +&

           cf1*u10**5 + cf0*u10**6

    END IF


    !Modifications as in znot_wind10m for icoef_sf=4


    !for wind<20, cd similar to icoef=2 at 10m, then reduced

    tmp=0.4*0.4/(alog(10.0/znotm))**2   ! cd at zlev

    aaa=0.75

    IF (u10 < 20) THEN

      aaa=0.99

    ELSEIF(u10 < 45.0) THEN

      aaa=0.99+(u10-20)*(0.75-0.99)/(45.0-20.0)

    END IF

    znotm=z10/exp( sqrt(0.4*0.4/(tmp*aaa)) )


  END SUBROUTINE wnd2z0m

  !/ ------------------------------------------------------------------- /

  SUBROUTINE wnd2sat(WND10,SAT)

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           B. G. Reichl            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         04-Aug-2016 |

    !/                  +-----------------------------------+

    !/

    !/    15-Jan-2016 : Origination.                        ( version 5.12 )

    !/    04-Aug-2016 : Updated for 2016 HWRF CD/U10 curve  ( J. Meixner   )

    !/

    !  1. Purpose :

    !

    !     Gives level of saturation spectrum to produce

    !         equivalent Cd as in wnd2z0m (for neutral 10m wind)

    !         tuned for method of Reichl et al. 2014

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !Input: WND10 - 10 m neutral wind [m/s]

    !Output: SAT  - Level 1-d saturation spectrum in tail [non-dim]

    !     ----------------------------------------------------------------

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3FLD1    Subr. W3FLD1MD Corresponding source term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    REAL, INTENT(IN) :: WND10

    REAL, INTENT(OUT) :: SAT

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

    CALL strace (ient, 'WND2SAT')

#endif

    !/ Old HWRF 2015 and ST2

    !        SAT =0.000000000001237 * WND10**6 +&

    !             -0.000000000364155 * WND10**5 +&

    !             0.000000037435015 * WND10**4 +&

    !             -0.000001424719473 * WND10**3 +&

    !             0.000000471570975 * WND10**2 +&

    !             0.000778467452178 * WND10**1 +&

    !             0.002962335390055

    !

    !     SAT values based on

    !     HWRF 2016 CD curve, created with  fetch limited cases ST4 physics

    IF (wnd10<20.0) THEN

      sat = -0.000018541921682*wnd10**2  &

           +0.000751515452434*wnd10     &

           +0.002466529381004

    ELSEIF (wnd10<45) THEN

      sat = -0.000000009060349*wnd10**4  &

           +0.000001276678367*wnd10**3  &

           -0.000068274393789*wnd10**2  &

           +0.001418180888868*wnd10     &

           +0.000262277682984

    ELSE

      sat = -0.000155976275073*wnd10     &

           +0.012027763023184

    ENDIF


    sat = min(max(sat,0.002),0.014)

  END SUBROUTINE wnd2sat

  !

  !/ End of module W3FLD1MD -------------------------------------------- /

  !/

END MODULE w3fld1md