w3updtmd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
 
MODULE w3updtmd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         22-Mar-2021 |
  !/                  +-----------------------------------+
  !/
  !/    21-Jan-2000 : Origination.                        ( version 2.00 )
  !/    24-Jan-2001 : Flat grid version.                  ( version 2.06 )
  !/    02-Apr-2001 : Adding sub-grid obstacles.          ( version 2.10 )
  !/    18-May-2001 : Clean up and bug fixes.             ( version 2.11 )
  !/    11-Jan-2002 : Sub-grid ice.                       ( version 2.15 )
  !/    30-Apr-2002 : Water level fixes.                  ( version 2.20 )
  !/    13-Nov-2002 : Add stress vector.                  ( version 3.00 )
  !/    26-Dec-2002 : Moving grid wind correction.        ( version 3.02 )
  !/    15-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    15-Jul-2005 : Adding MAPST2.                      ( version 3.07 )
  !/    07-Sep-2005 : Upgrading W3UBPT.                   ( version 3.08 )
  !/    04-Jul-2006 : Consolidate stress arrays.          ( version 3.09 )
  !/    11-Jan-2007 : Clean-up W3UTRN boundary points.    ( version 3.10 )
  !/    11-May-2007 : Adding NTPROC/NAPROC separation.    ( version 3.11 )
  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
  !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )
  !/                  (W. E. Rogers & T. J. Campbell, NRL)
  !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
  !/                  (W. E. Rogers & T. J. Campbell, NRL)
  !/    17-Aug-2010 : ABPI0-N(:,0) init. bug fix.         ( version 3.14 )
  !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
  !/                  specify index closure for a grid.   ( version 3.14 )
  !/                  (T. J. Campbell, NRL)
  !/    05-Apr-2011 : Place holder for XGR in UNGTYPE     ( version 4.04 )
  !/                  (A. Roland/F. Ardhuin)
  !/    13-Mar-2012 : Add initialization of UST on re-    ( version 4.07 )
  !/                  activation of grid point.
  !/    06-Jun-2012 : Porting bugfixes from 3.14 to 4.07  ( version 4.07 )
  !/    12-Jun-2012 : Add /RTD option or rotated grid option.
  !/                  (Jian-Guo Li)                       ( version 4.07 )
  !/    26-Sep-2012 : Adding update from tidal analysis   ( version 4.08 )
  !/                  (F. Ardhuin)
  !/    16-Sep-2013 : Add Arctic part for SMC grid.       ( version 4.11 )
  !/    11-Nov-2013 : SMC and rotated grid incorporated in the main
  !/                  trunk                               ( version 4.13 )
  !/    13-Nov-2013 : Moved reflection from ww3_grid.ftn  ( version 4.13 )
  !/    27-May-2014 : Ading OMPG parallelizations dir,    ( version 5.02 )
  !/    08-May-2014 : Implement tripolar grid for first order propagation
  !/                  scheme                              ( version 5.03 )
  !/                  (W. E. Rogers, NRL)
  !/    27-Aug-2015 : New function to update variables    ( version 5.08 )
  !/                  ICEF and ICEDMAX at the first time step
  !/                  and add ICEH initialization in W3UICE.
  !/    13-Jan-2016 : Changed initial value of ICEDMAX    ( version 5.08 )
  !/    26-Mar-2018 : Sea-point only Wnd/Cur input.  JGLi ( version 6.04 )
  !/    07-Oct-2019 : RTD option with standard lat-lon
  !/                  grid when nesting to rotated grid   ( version 7.11 )
  !/    22-Mar-2021 : Add W3UTAU, W3URHO routines         ( version 7.13 )
  !/    06-May-2021 : Use ARCTC option for SMC grid. JGLi ( version 7.13 )
  !/
  !/    Copyright 2009-2014 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 :
  !
  !     Bundles all input updating routines for WAVEWATCH III.
  !
  !  2. Variables and types :
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3UCUR    Subr. Public   Update current fields.
  !      W3UWND    Subr. Public   Update wind fields.
  !      W3UTAU    Subr. Public   Update atmospheric momentum fields.
  !      W3UINI    Subr. Public   Update initial conditions.
  !      W3UBPT    Subr. Public   Update boundary conditions.
  !      W3UICE    Subr. Public   Update ice concentrations.
  !      W3ULEV    Subr. Public   Update water levels.
  !      W3URHO    Subr. Public   Update air density.
  !      W3UTRN    Subr. Public   Update cell boundary transparancies.
  !      W3DZXY    Subr. Public   Calculate derivatives of a field.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      DSEC21    Func. W3TIMEMD Difference in time.
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      EXTCDE    Subr. W3SERVMD Exit program with error code.
  !      PRTBLK    Subr. W3ARRYMD Print plot output.
  !      PRT2DS    Subr. W3ARRYMD Print plot output.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !       !/SHRD   Switch for shared / distributed memory architecture.
  !       !/DIST   Id.
  !
  !       !/OMPG   OpenMP compiler directives.
  !
  !       !/CRT0   No current interpolation.
  !       !/CRT1   Linear current interpolation.
  !       !/CRT2   Quasi-quadratic current interpolation.
  !
  !       !/WNT0   No wind/momentum interpolation.
  !       !/WNT1   Linear wind/momentum interpolation.
  !       !/WNT2   Energy conservation in wind/momentum interpolation.
  !
  !       !/RWND   Use wind speeds relative to currents.
  !
  !       !/STAB2  Calculate effective wind speed factor for stability
  !                to be used with !/ST2.
  !
  !       !/S      Enable subroutine tracing.
  !       !/Tn     Test output
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  USE constants
  USE w3odatmd, ONLY: ndse, ndst, naproc, iaproc, naperr
#ifdef W3_S
  USE w3servmd, ONLY : strace
#endif
  USE w3timemd, ONLY: dsec21
  !/
  !/ ------------------------------------------------------------------- /
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE w3ucur ( FLFRST )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         15-Dec-2004 |
    !/                  +-----------------------------------+
    !/
    !/    09-Dec-1996 : Final FORTRAN 77                    ( version 1.18 )
    !/    20-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    15-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    27-Aug-2015 : Rename DT0,DTT by DT0T,DT0N         ( version 5.10 )
    !/    23-Mar-2016 : SMC grid Arctic part adjustment.    ( version 5.18 )
    !/    26-Mar-2018 : Sea-point only current on SMC grid. ( version 6.02 )
    !/
    !  1. Purpose :
    !
    !     Interpolate the current field to the present time.
    !
    !  2. Method :
    !
    !     Linear interpolation of speed and direction, with optionally
    !     a correction to get approximate quadratic interpolation of speed
    !     only.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       FLFRST  Log.  I   Flag for first pass through routine.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - Only currents at sea points are considered.
    !     - Time ranges checked in W3WAVE.
    !     - Currents are stored by components to save on the use of
    !       SIN and COS functions. The actual interpolations, however
    !       are by absolute value and direction.
    !
    !  8. Structure :
    !
    !     --------------------------------------
    !      1.  Prepare auxiliary arrays.
    !      2.  Calculate interpolation factors.
    !      3.  Get actual winds.
    !     --------------------------------------
    !
    !  9. Switches :
    !
    !     !/CRT0  No current interpolation.
    !     !/CRT1  Linear current interpolation.
    !     !/CRT2  Quasi-quadratic current interpolation.
    !
    !     !/S     Enable subroutine tracing.
    !     !/T     Test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsf
#ifdef W3_SMC
    USE w3gdatmd, ONLY: narc, nglo, angarc
    USE w3gdatmd, ONLY: fswnd, arctc
#endif
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cx, cy, ca0, cai, cd0, cdi
    USE w3idatmd, ONLY: tc0, cx0, cy0, tcn, cxn, cyn
#ifdef W3_TIDE
    USE w3gdatmd, ONLY: ygrd
    USE w3timemd
    USE w3idatmd, ONLY: flcurtide, cxtide, cytide, ntide
    USE w3tidemd
#endif
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    LOGICAL, INTENT(IN)     :: FLFRST
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER                 :: ISEA, IX, IY
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: D0, DN, DD, DT0N, DT0T, RD, CABS, CDIR
#ifdef W3_CRT2
    REAL                    :: RD2, CI2
#endif
#ifdef W3_TIDE
    INTEGER          :: J,K
    INTEGER(KIND=4)  :: TIDE_KD0, INT24, INTDYS       ! "Gregorian day constant"
    REAL             :: WCURTIDEX, WCURTIDEY, TIDE_ARGX, TIDE_ARGY
    REAL(KIND=8)     :: d1,h,tide_hour,hh,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau
    REAL             :: FX(44),UX(44),VX(44)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3UCUR')
#endif
    !
    ! 1.  Prepare auxiliary arrays
    !
    IF ( flfrst ) THEN
      DO isea=1, nsea
#ifdef W3_SMC
        !!Li  For sea-point SMC grid current, the 1-D current is stored on
        !!Li  2-D CX0(NSEA, 1) variable.
        IF( fswnd ) THEN
          ix = isea
          iy = 1
        ELSE
#endif
          ix        = mapsf(isea,1)
          iy        = mapsf(isea,2)
#ifdef W3_SMC
        ENDIF
#endif
 
        ca0(isea) = sqrt( cx0(ix,iy)**2 + cy0(ix,iy)**2 )
        cai(isea) = sqrt( cxn(ix,iy)**2 + cyn(ix,iy)**2 )
        IF ( ca0(isea) .GT. 1.e-7) THEN
          d0     = mod( tpi+atan2(cy0(ix,iy),cx0(ix,iy)) , tpi )
        ELSE
          d0     = 0
        END IF
        IF ( cai(isea) .GT. 1.e-7) THEN
          dn     = mod( tpi+atan2(cyn(ix,iy),cxn(ix,iy)) , tpi )
        ELSE
          dn     = d0
        END IF
        IF ( ca0(isea) .GT. 1.e-7) THEN
          cd0(isea) = d0
        ELSE
          cd0(isea) = dn
        END IF
        dd     = dn - cd0(isea)
        IF (abs(dd).GT.pi) dd = dd - tpi*sign(1.,dd)
        cdi(isea) = dd
        cai(isea) = cai(isea) - ca0(isea)
      END DO
    END IF
    !
    ! 2.  Calculate interpolation factor
    !
    dt0n    = dsec21( tc0, tcn )
    dt0t    = dsec21( tc0, time )
    !
#ifdef W3_CRT0
    rd     = 0.
#endif
#ifdef W3_CRT1
    rd     = dt0t / max( 1.e-7 , dt0n )
#endif
#ifdef W3_CRT2
    rd     = dt0t / max( 1.e-7 , dt0n )
    rd2    = 1. - rd
#endif
#ifdef W3_OASOCM
    rd     = 1.
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) dt0n, dt0t, rd
#endif
 
#ifdef W3_TIDE
    IF (flcurtide) THEN
      !          WRITE(6,*) 'TIME CUR:',TIME, '##',TC0, '##',TCN
      tide_hour = time2hours(time)
      !
      !*  THE ASTRONOMICAL ARGUMENTS ARE CALCULATED BY LINEAR APPROXIMATION
      !*  AT THE MID POINT OF THE ANALYSIS PERIOD.
      d1=tide_hour/24.d0
      tide_kd0= 2415020
      d1=d1-dfloat(tide_kd0)-0.5d0
      call astr(d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp)
      int24=24
      intdys=int((tide_hour+0.00001)/int24)
      hh=tide_hour-dfloat(intdys*int24)
      tau=hh/24.d0+h-s
    END IF
    !
    !  ONLY THE FRACTIONAL PART OF A SOLAR DAY NEED BE RETAINED FOR COMPU-
    !  TING THE LUNAR TIME TAU.
    !
#endif
 
    !
    ! 3.  Actual currents for all grid points
    !
    DO isea=1, nsea
#ifdef W3_TIDE
      IF (flcurtide) THEN  ! could move IF test outside of ISEA loop ...
        ! VUF should only be updated in latitude changes significantly ...
        ix        = mapsf(isea,1)
        iy        = mapsf(isea,2)
        CALL setvuf_fast(h,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau,real(ygrd(iy,ix)),fx,ux,vx)
        wcurtidex = cxtide(ix,iy,1,1)
        wcurtidey = cytide(ix,iy,1,1)
 
        DO j=2,tide_mf
          tide_argx=(vx(j)+ux(j))*twpi-cxtide(ix,iy,j,2)*dera
          tide_argy=(vx(j)+ux(j))*twpi-cytide(ix,iy,j,2)*dera
          wcurtidex = wcurtidex+fx(j)*cxtide(ix,iy,j,1)*cos(tide_argx)
          wcurtidey = wcurtidey+fx(j)*cytide(ix,iy,j,1)*cos(tide_argy)
        END DO
 
#endif
 
#ifdef W3_TIDET
        !Verification
        IF (isea.EQ.1) THEN
 
          tide_ampc(1:ntide,1)=cxtide(ix,iy,1:ntide,1)
          tide_phg(1:ntide,1 )=cxtide(ix,iy,1:ntide,2)
          tide_ampc(1:ntide,2)=cytide(ix,iy,1:ntide,1)
          tide_phg(1:ntide,2) =cytide(ix,iy,1:ntide,2)
 
          WRITE(993,'(A,F20.2,13F8.3)') 'TEST ISEA 0:',    &
               d1,h,s,tau,pp,s,p,enp,dh,dpp,ds,dp,dnp,real(ygrd(iy,ix))
 
          DO j=1,tide_mf
            WRITE(993,'(A,4I9,F12.0,3F8.3,I4,X,A)') 'TEST ISEA 1:',ix,j,time,tide_hour,    &
                 fx(j),ux(j),vx(j),tide_index2(j),tidecon_allnames(tide_index2(j))
          END DO
          DO k=1,2
            DO j=1,tide_mf
              WRITE(993,'(A,5I9,F12.0,5F8.3)') 'TEST ISEA 2:',ix,k,j,time,tide_hour,    &
                   fx(j),ux(j),vx(j),tide_ampc(j,k),tide_phg(j,k)
            END DO
          END DO
 
          WRITE(993,'(A,2F8.4,A,2F8.4)') '#:',cx0(ix,iy),cy0(ix,iy),'##',wcurtidex,wcurtidey
          CLOSE(993)
        END IF
        ! End of verification
#endif
#ifdef W3_TIDE
        cx(isea) = wcurtidex
        cy(isea) = wcurtidey
      ELSE
#endif
 
        cabs    = ca0(isea) + rd * cai(isea)
#ifdef W3_CRT2
        ci2      = sqrt( rd2 *      ca0(isea)**2 +             &
             rd  *(ca0(isea)+cai(isea))**2 )
        cabs    = cabs * min( 1.25 , ci2/max(1.e-7,cabs) )
#endif
        cdir    = cd0(isea) + rd * cdi(isea)
 
#ifdef W3_SMC
        !Li   Rotate curreent direction by ANGARC for Arctic part cells.  JGLi23Mar2016
        IF( arctc .AND. (isea .GT. nglo) ) THEN
          dn = cdir + angarc( isea - nglo )*dera
          cdir = mod( tpi + dn, tpi )
        ENDIF
#endif
 
        cx(isea) = cabs * cos(cdir)
        cy(isea) = cabs * sin(cdir)
#ifdef W3_TIDE
        !        IF (ISEA.EQ.1)  WRITE(6,'(A,4F8.4,A,4F8.4)') 'CUR#:',RD,CA0(ISEA),CAI(ISEA),CABS,'##', &
        !                                      CX(ISEA), CY(ISEA),WCURTIDEX, WCURTIDEY
      END IF
#endif
      !
    END DO
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3UCUR : DT0N, DT0T, RD :',2f8.1,f6.3)
#endif
    !/
    !/ End of W3UCUR ----------------------------------------------------- /
    !/
  END SUBROUTINE w3ucur
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE w3uwnd ( FLFRST, VGX, VGY )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         27-May-2014 |
    !/                  +-----------------------------------+
    !/
    !/    03-Dec-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    20-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    13-Nov-2002 : Add stress vector.                  ( version 3.00 )
    !/    26-Dec-2002 : Moving grid wind correction.        ( version 3.02 )
    !/    15-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    04-Jul-2006 : Consolidate stress arrays.          ( version 3.09 )
    !/    16-Sep-2013 : Rotating wind for Arctic part.      ( version 4.11 )
    !/    27-May-2014 : Adding OMPG parallelizations dir.   ( version 5.02 )
    !/    27-Aug-2015 : Rename DT0,DTT by DT0T,DT0N         ( version 5.10 )
    !/    26-Mar-2018 : Sea-point only wind for SMC grid.   ( version 6.07 )
    !/
    !  1. Purpose :
    !
    !     Interpolate wind fields to the given time.
    !
    !  2. Method :
    !
    !     Linear interpolation of wind speed and direction, with a simple
    !     correction to obtain quasi-conservation of energy.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       FLFRST  Log.  I   Flag for first pass through routine.
    !       VGX/Y   Real  I   Grid velocity                      (!/MGW)
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - Only winds over sea points are considered.
    !     - Time ranges checked in W3WAVE.
    !
    !  8. Structure :
    !
    !     --------------------------------------
    !      1.  Prepare auxiliary arrays.
    !      2.  Calculate interpolation factors
    !      3.  Get actual winds
    !      4.  Correct for currents
    !      5.  Convert to stresses
    !      6.  Stability correction
    !     --------------------------------------
    !
    !  9. Switches :
    !
    !     !/OMPG   OpenMP compiler directives.
    !
    !     !/WNT0   No wind interpolation.
    !     !/WNT1   Linear wind interpolation.
    !     !/WNT2   Energy conservation in wind interpolation.
    !
    !     !/RWND   Use wind speeds relative to currents.
    !     !/MGW    Moving grid wind correction.
    !
    !     !/STAB2  Calculate effective wind speed factor for stability
    !              to be used with !/ST2.
    !
    !     !/S      Enable subroutine tracing.
    !     !/T      Test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsf
#ifdef W3_WCOR
    USE w3gdatmd, ONLY:   wwcor
#endif
#ifdef W3_RWND
    USE w3gdatmd, ONLY:   rwindc
#endif
#ifdef W3_ST2
    USE w3gdatmd, ONLY: zwind, ofstab, ffng, ffps, ccng, ccps, shstab
#endif
#ifdef W3_SMC
    USE w3gdatmd, ONLY: narc, nglo, angarc, arctc, fswnd
#endif
    USE w3wdatmd, ONLY: time, asf
    USE w3adatmd, ONLY: dw, cx, cy, ua, ud, u10, u10d, as,          &
         ua0, uai, ud0, udi, as0, asi
    USE w3idatmd, ONLY: tw0, wx0, wy0, dt0, twn, wxn, wyn, dtn, flcur
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: VGX, VGY
    LOGICAL, INTENT(IN)     :: FLFRST
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER                 :: ISEA, IX, IY
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: D0, DN, DD, DT0N, DT0T, RD, UI2,      &
         UXR, UYR
#ifdef W3_WNT2
    REAL                    :: RD2
#endif
#ifdef W3_STAB2
    REAL                    :: STAB0, STAB, THARG1, THARG2, COR1, COR2
#endif
    REAL                    :: UDARC
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3UWND')
#endif
    !
    ! 1.  Prepare auxiliary arrays
    !
    IF ( flfrst ) THEN
      DO isea=1, nsea
#ifdef W3_SMC
        !!Li  For sea-point only SMC grid wind 1-D wind is stored on
        !!Li  2-D WX0(NSEA, 1) variable.
        IF( fswnd ) THEN
          ix = isea
          iy = 1
        ELSE
#endif
          ix        = mapsf(isea,1)
          iy        = mapsf(isea,2)
#ifdef W3_SMC
        ENDIF
#endif
 
        ua0(isea) = sqrt( wx0(ix,iy)**2 + wy0(ix,iy)**2 )
        uai(isea) = sqrt( wxn(ix,iy)**2 + wyn(ix,iy)**2 )
        IF ( ua0(isea) .GT. 1.e-7) THEN
          d0     = mod( tpi+atan2(wy0(ix,iy),wx0(ix,iy)) , tpi )
        ELSE
          d0     = 0
        END IF
        IF ( uai(isea) .GT. 1.e-7) THEN
          dn     = mod( tpi+atan2(wyn(ix,iy),wxn(ix,iy)) , tpi )
        ELSE
          dn     = d0
        END IF
        IF ( ua0(isea) .GT. 1.e-7) THEN
          ud0(isea) = d0
        ELSE
          ud0(isea) = dn
        END IF
        dd     = dn - ud0(isea)
        IF (abs(dd).GT.pi) dd = dd - tpi*sign(1.,dd)
        udi(isea) = dd
        uai(isea) = uai(isea) - ua0(isea)
        as0(isea) = dt0(ix,iy)
        asi(isea) = dtn(ix,iy) - dt0(ix,iy)
      END DO
    END IF
    !
    ! 2.  Calculate interpolation factor
    !
    dt0n    = dsec21( tw0, twn )
    dt0t   = dsec21( tw0, time )
    !
#ifdef W3_WNT0
    rd     = 0.
#endif
#ifdef W3_WNT1
    rd     = dt0t / max( 1.e-7 , dt0n )
#endif
#ifdef W3_WNT2
    rd     = dt0t / max( 1.e-7 , dt0n )
    rd2    =  1. - rd
#endif
#ifdef W3_OASACM
    rd     = 1.
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) dt0n, dt0t, rd
#endif
    !
    ! 3.  Actual wind for all grid points
    !
#ifdef W3_OMPG
    !$OMP PARALLEL DO PRIVATE (ISEA,UI2,UXR,UYR,UDARC)
#endif
    !
    DO isea=1, nsea
      !
      ua(isea) = ua0(isea) + rd * uai(isea)
#ifdef W3_WNT2
      ui2      = sqrt( rd2 *      ua0(isea)**2 +             &
           rd  *(ua0(isea)+uai(isea))**2 )
      ua(isea) = ua(isea) * min(1.25,ui2/max(1.e-7,ua(isea)))
#endif
      ud(isea) = ud0(isea) + rd * udi(isea)
#ifdef W3_MGW
      uxr        = ua(isea)*cos(ud(isea)) + vgx
      uyr        = ua(isea)*sin(ud(isea)) + vgy
      ua(isea) = max( 0.001 , sqrt(uxr**2+uyr**2) )
      ud(isea) = mod( tpi+atan2(uyr,uxr) , tpi )
#endif
#ifdef W3_SMC
      !Li   Rotate wind direction by ANGARC for Arctic part cells.
      IF( arctc .AND. (isea .GT. nglo) ) THEN
        udarc = ud(isea) + angarc( isea - nglo )*dera
        ud(isea) = mod( tpi + udarc, tpi )
      ENDIF
#endif
      !
      as(isea) = as0(isea) + rd * asi(isea)
      !        IF (UA(ISEA).NE.UA(ISEA)) WRITE(6,*) 'BUG WIND:',ISEA,UA(ISEA),MAPSF(ISEA,1), MAPSF(ISEA,2),UA0(ISEA),RD,UAI(ISEA)
      !        IF (UD(ISEA).NE.UD(ISEA)) WRITE(6,*) 'BUG WIN2:',ISEA,UD(ISEA),MAPSF(ISEA,1), MAPSF(ISEA,2)
      !
    END DO
    !
#ifdef W3_OMPG
    !$OMP END PARALLEL DO
#endif
    !
    ! 3.b  Bias correction ( !/WCOR )
#ifdef W3_WCOR
    WHERE ( ua .GE. wwcor(1) ) ua = ua+(ua-wwcor(1))*wwcor(2)
#endif
 
    !
    ! 4.  Correct for currents and grid motion
    !
#ifdef W3_RWND
    IF ( flcur ) THEN
#endif
      !
#ifdef W3_RWND
      DO isea=1, nsea
        uxr        = ua(isea)*cos(ud(isea)) - rwindc*cx(isea)
        uyr        = ua(isea)*sin(ud(isea)) - rwindc*cy(isea)
        u10(isea) = max( 0.001 , sqrt(uxr**2+uyr**2) )
        u10d(isea) = mod( tpi+atan2(uyr,uxr) , tpi )
      END DO
#endif
      !
#ifdef W3_RWND
    ELSE
#endif
      !
#ifdef W3_OMPG
      !$OMP PARALLEL DO PRIVATE (ISEA)
#endif
      !
      DO isea=1, nsea
        u10(isea) = max( ua(isea) , 0.001 )
        u10d(isea) = ud(isea)
      END DO
      !
#ifdef W3_OMPG
      !$OMP END PARALLEL DO
#endif
      !
#ifdef W3_RWND
    END IF
#endif
    !
    ! 5.  Stability correction ( !/STAB2 )
    !     Original settings :
    !
    !     SHSTAB =    1.4
    !     OFSTAB =   -0.01
    !     CCNG   =   -0.1
    !     CCPS   =    0.1
    !     FFNG   = -150.
    !     FFPS   =  150.
    !
#ifdef W3_STAB2
    stab0  = zwind * grav / 273.
#endif
    !
#ifdef W3_STAB2
    DO isea=1, nsea
      stab   = stab0 * as(isea) / max(5.,u10(isea))**2
      stab   = max( -1. , min( 1. , stab ) )
#endif
      !
#ifdef W3_STAB2
      tharg1 = max( 0. , ffng*(stab-ofstab))
      tharg2 = max( 0. , ffps*(stab-ofstab))
      cor1   = ccng * tanh(tharg1)
      cor2   = ccps * tanh(tharg2)
#endif
      !
#ifdef W3_STAB2
      asf(isea) = sqrt( (1.+cor1+cor2)/shstab )
      u10(isea) = u10(isea) / asf(isea)
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3UWND : DT0N, DT0T, RD :',2f8.1,f6.3)
#endif
    !/
    !/ End of W3UWND ----------------------------------------------------- /
    !/
  END SUBROUTINE w3uwnd
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3utau ( FLFRST )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           J. M. Castillo          |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         22-Mar-2021 |
    !/                  +-----------------------------------+
    !/
    !/    22-Mar-2021 : First implementation                ( version 7.13 )
    !/
    !  1. Purpose :
    !
    !     Interpolate atmosphere momentum fields to the given time.
    !
    !  2. Method :
    !
    !     Linear interpolation of momentum module and direction, with a simple
    !     correction to obtain quasi-conservation of energy.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       FLFRST  Log.  I   Flag for first pass through routine.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - Only momentum over sea points is considered.
    !     - Time ranges checked in W3WAVE.
    !
    !  8. Structure :
    !
    !     --------------------------------------
    !      1.  Prepare auxiliary arrays.
    !      2.  Calculate interpolation factors
    !      3.  Get actual momentum
    !     --------------------------------------
    !
    !  9. Switches :
    !
    !     !/OMPG   OpenMP compiler directives.
    !
    !     !/WNT0   No momentum interpolation.
    !     !/WNT1   Linear momentum interpolation.
    !     !/WNT2   Energy conservation in momentum interpolation.
    !
    !     !/S      Enable subroutine tracing.
    !     !/T      Test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nsea, mapsf
#ifdef W3_SMC
    USE w3gdatmd, ONLY: narc, nglo, angarc
    USE w3gdatmd, ONLY: fswnd, arctc
#endif
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: taua, tauadir, ma0, mai, md0, mdi
    USE w3idatmd, ONLY: tu0, ux0, uy0, tun, uxn, uyn
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    LOGICAL, INTENT(IN)     :: FLFRST
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER                 :: ISEA, IX, IY
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: D0, DN, DD, DT0N, DT0T, RD, MI2,      &
         MXR, MYR
#ifdef W3_WNT2
    REAL                    :: RD2
#endif
    REAL                    :: MDARC
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3UTAU')
#endif
    !
    ! 1.  Prepare auxiliary arrays
    !
    IF ( flfrst ) THEN
      DO isea=1, nsea
#ifdef W3_SMC
        !!Li  For sea-point only SMC grid momentum 1-D momentum is stored on
        !!Li  2-D UX0(NSEA, 1) variable.
        IF( fswnd ) THEN
          ix = isea
          iy = 1
        ELSE
#endif
          ix        = mapsf(isea,1)
          iy        = mapsf(isea,2)
#ifdef W3_SMC
        ENDIF
#endif
 
        ma0(isea) = sqrt( ux0(ix,iy)**2 + uy0(ix,iy)**2 )
        mai(isea) = sqrt( uxn(ix,iy)**2 + uyn(ix,iy)**2 )
        IF ( ma0(isea) .GT. 1.e-7) THEN
          d0     = mod( tpi+atan2(uy0(ix,iy),ux0(ix,iy)) , tpi )
        ELSE
          d0     = 0
        END IF
        IF ( mai(isea) .GT. 1.e-7) THEN
          dn     = mod( tpi+atan2(uyn(ix,iy),uxn(ix,iy)) , tpi )
        ELSE
          dn     = d0
        END IF
        IF ( ma0(isea) .GT. 1.e-7) THEN
          md0(isea) = d0
        ELSE
          md0(isea) = dn
        END IF
        dd     = dn - md0(isea)
        IF (abs(dd).GT.pi) dd = dd - tpi*sign(1.,dd)
        mdi(isea) = dd
        mai(isea) = mai(isea) - ma0(isea)
      END DO
    END IF
    !
    ! 2.  Calculate interpolation factor
    !
    dt0n    = dsec21( tu0, tun )
    dt0t   = dsec21( tu0, time )
    !
#ifdef W3_WNT0
    rd     = 0.
#endif
#ifdef W3_WNT1
    rd     = dt0t / max( 1.e-7 , dt0n )
#endif
#ifdef W3_WNT2
    rd     = dt0t / max( 1.e-7 , dt0n )
    rd2    =  1. - rd
#endif
#ifdef W3_OASACM
    rd     = 1.
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) dt0n, dt0t, rd
#endif
    !
    ! 3.  Actual momentum for all grid points
    !
#ifdef W3_OMPG
    !$OMP PARALLEL DO PRIVATE (ISEA,MI2,MXR,MYR,MDARC)
#endif
    !
    DO isea=1, nsea
      !
      taua(isea) = ma0(isea) + rd * mai(isea)
#ifdef W3_WNT2
      mi2      = sqrt( rd2 *      ma0(isea)**2 +             &
           rd  *(ma0(isea)+mai(isea))**2 )
      taua(isea) = taua(isea) * min(1.25,mi2/max(1.e-7,taua(isea)))
#endif
      tauadir(isea) = md0(isea) + rd * mdi(isea)
#ifdef W3_SMC
      !Li   Rotate momentum direction by ANGARC for Arctic part cells.
      IF( arctc .AND. (isea .GT. nglo) ) THEN
        mdarc = tauadir(isea) + angarc( isea - nglo )*dera
        tauadir(isea) = mod( tpi + mdarc, tpi )
      ENDIF
#endif
      !
    END DO
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3UTAU : DT0N, DT0T, RD :',2f8.1,f6.3)
#endif
    !/
    !/ End of W3UTAU ----------------------------------------------------- /
    !/
  END SUBROUTINE w3utau
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3uini ( A )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         06-Jun-2018 |
    !/                  +-----------------------------------+
    !/
    !/    19-Oct-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    20-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    24-Jan-2001 : Flat grid version.                  ( version 2.06 )
    !/    18-May-2001 : Fix CG declaration.                 ( version 2.11 )
    !/    15-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
    !/    06-Jun-2018 : use W3PARALL and INIT_GET_ISEA      ( version 6.04 )
    !/
    !  1. Purpose :
    !
    !     Initialize the wave field with fetch-limited spectra before the
    !     actual calculation start. (Named as an update routine due to
    !     placement in code.)
    !
    !  2. Method :
    !
    !     Fetch-limited JONSWAP spectra with a cosine^2 directional
    !     distribution and a mean direction taken from the wind.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.   O   Action density spectra.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - Wind speeds filtered by U10MIN and U10MAX (DATA statements)
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/SHRD  Switch for shared / distributed memory architecture.
    !       !/DIST  Id.
    !
    !       !/S     Enable subroutine tracing.
    !       !/T     General test output.
    !       !/T1    Parameters at grid points.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nx, ny, nsea, nseal, mapsf,                 &
         nk, nth, th, sig, dth, dsip, ungtype,       &
         rlgtype, clgtype, gtype, flagll,            &
         hpfac, hqfac
    USE w3adatmd, ONLY: u10, u10d, cg
    USE w3parall, only : init_get_jsea_isproc, init_get_isea
    USE w3parall, only : get_jsea_ibelong
#ifdef W3_T
    USE w3arrymd, ONLY : prtblk
#endif
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(OUT)       :: A(NTH,NK,0:NSEAL)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local variables
    !/
    INTEGER                 :: IX, IY, ISEA, JSEA, IK, ITH, ISPROC
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T
    INTEGER                 :: IX0, IXN, MAPOUT(NX,NY)
    INTEGER                 :: NXP = 60
#endif
    REAL                    :: ALFA(NSEAL), FP(NSEAL), YLN(NSEAL),  &
         AA, BB, CC
    REAL                    :: XGR, U10C, U10DIR, XSTAR, FSTAR,     &
         GAMMA, FR, D1(NTH), D1INT, F1, F2
    REAL                    :: ETOT, E1I
    REAL                    :: U10MIN =  1.
    REAL                    :: U10MAX = 20.
#ifdef W3_T
    REAL                    :: HSIG(NX,NY)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3UINI')
#endif
    !
    !
    ! Pre-process JONSWAP data for all grid points ----------------------- *
    !
#ifdef W3_T1
    WRITE (ndst,9010)
#endif
    !
    !  this is not clear what is going on betwen w3init and this ...
    a(:,:,:)=0
    DO jsea=1, nseal
      CALL init_get_isea(isea, jsea)
      IF (gtype.EQ.ungtype) THEN
        xgr=1.  ! to be fixed later
      ELSE
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        xgr    = 0.5 * sqrt(hpfac(iy,ix)**2+hqfac(iy,ix)**2)
      END IF
      IF ( flagll ) THEN
        xgr    = xgr * radius * dera
      END IF
      !
      u10c   = max( min(u10(isea),u10max) , u10min )
      !
      xstar  = grav * xgr / u10c**2
      fstar  = 3.5 / xstar**(0.33)
      gamma  = max( 1. , 7.0 / xstar**(0.143) )
      !
      alfa(jsea) = 0.076 / xstar**(0.22)
      fp(jsea) = fstar * grav / u10c
      yln(jsea) = log( gamma )
      !
#ifdef W3_T1
      WRITE (ndst,9011) isea, u10c, xstar,                      &
           alfa(jsea), fp(jsea), gamma
#endif
      !
    END DO
    !
    ! 1-D spectrum at location ITH = NTH --------------------------------- *
    !
    DO ik=1, nk
      fr     = sig(ik) * tpiinv
      DO jsea=1, nseal
        !
        !/ ----- INLINED EJ5P (REDUCED) -------------------------------------- /
        !
        aa     = alfa(jsea) * 0.06175/fr**5
        bb     = max( -50. , -1.25*(fp(jsea)/fr)**4 )
        cc     = max( -50. , -0.5*((fr-fp(jsea))/(0.07*fp(jsea)))**2 )
        a(nth,ik,jsea)                                              &
             = aa * exp(bb + exp(cc) * yln(jsea))
        !
        !/ ----- INLINED EJ5P (END) ------------------------------------------ /
        !
      END DO
    END DO
    !
    ! Apply directional distribution ------------------------------------- *
    !
    DO jsea=1, nseal
      CALL init_get_isea(isea, jsea)
      u10dir = u10d(isea)
      d1int  = 0.
      !
      DO ith=1, nth
        d1(ith) = ( max( 0. , cos(th(ith)-u10dir) ) )**2
        d1int   = d1int + d1(ith)
      END DO
      !
      d1int  = d1int * dth
      f1     = tpiinv / d1int
      !
      DO ik=1, nk
        f2     = f1 * a(nth,ik,jsea) * cg(ik,isea) / sig(ik)
        DO ith=1, nth
          a(ith,ik,jsea) = f2 * d1(ith)
        END DO
      END DO
      !
    END DO
    !
    ! Test output -------------------------------------------------------- *
    !
#ifdef W3_T
    hsig   = 0.
    mapout = 0
#endif
    !
#ifdef W3_T
    DO isea=iaproc, nsea, naproc
      jsea   = 1 + (isea-1)/naproc
      etot   = 0.
      DO ik=1, nk
        e1i    = 0.
        DO ith=1, nth
          e1i    = e1i + a(ith,ik,jsea)
        END DO
        etot   = etot + e1i * dsip(ik) * sig(ik) / cg(ik,isea)
      END DO
      ix            = mapsf(isea,1)
      iy            = mapsf(isea,2)
      hsig(ix,iy) = 4. * sqrt( etot * dth )
      mapout(ix,iy) = 1
    END DO
#endif
    !
#ifdef W3_T
    ix0    = 1
    DO
      ixn    = min( nx , ix0+nxp-1 )
      CALL prtblk (ndst, nx, ny, nx, hsig, mapout, 0, 0.,        &
           ix0, ixn, 1, 1, ny, 1, 'Hs', 'm')
      IF ( ixn .EQ. nx ) EXIT
      ix0    = ix0 + nxp
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3UINI : XGR = ',e10.3)
#endif
    !
#ifdef W3_T1
9010 FORMAT (' TEST W3UINI :  ISEA, U10C, XSTAR, ALPHA, FP, GAMMA')
9011 FORMAT ('   ',i6,f8.2,f10.1,2f6.3,f6.2)
#endif
    !/
    !/ End of W3UINI ----------------------------------------------------- /
    !/
  END SUBROUTINE w3uini
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3ubpt
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         06-Jun-2018 |
    !/                  +-----------------------------------+
    !/
    !/    19-Oct-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    20-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    15-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    07-Sep-2005 : Moving update to end of time step.  ( version 3.08 )
    !/    17-Aug-2010 : Add initialization ABPI0-N(:,0).  ( version 3.14.5 )
    !/    12-Jun-2012 : Add /RTD option or rotated grid option.
    !/                  (Jian-Guo Li)                       ( version 4.06 )
    !/    06-Jun-2018 : Add DEBUGIOBC/SETUP/DEBUGW3ULEV     ( version 6.04 )
    !/    13-Jun-2019 : Rotation only if POLAT<90 (C.Hansen)( version 7.11 )
    !/
    !  1. Purpose :
    !
    !     Update spectra at the active boundary points.
    !
    !  2. Method :
    !
    !     Spectra are read and interpolated in space and time from the
    !     data read by W3IOBC.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !       STRACE, DSEC21
    !                Service routines.
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - The data arrays contain sigma spectra to assure conservation
    !       when changing grids.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S     Enable subroutine tracing.
    !     !/T0    Test output of wave heights.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nspec, mapwn, sig2, dden
#ifdef W3_RTD
    !!   Use rotation angle and action conversion sub.  JGLi12Jun2012
    USE w3gdatmd, ONLY: nk, nth, nspec, angld, polat
    USE w3servmd, ONLY: w3acturn
#endif
    USE w3adatmd, ONLY: cg
    USE w3odatmd, ONLY: nbi, abpi0, abpin, isbpi, ipbpi, rdbpi,     &
         bbpi0, bbpin
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER                 :: IBI, ISP, ISEA
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T0
    REAL                    :: HS1, HS2
#endif
#ifdef W3_RTD
    !!    Declare a temporary spectr variable.  JGLi12Jun2012
    REAL :: Spectr(NSPEC), AnglBP
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3UBPT')
#endif
    !
    ! 1.  Process BBPI0 -------------------------------------------------- *
    ! 1.a First intialization
 
 
    !
    IF ( bbpi0(1,0) .EQ. -1. ) THEN
      !
      bbpi0(:,0) = 0.
      bbpin(:,0) = 0.
      abpi0(:,0) = 0.
      abpin(:,0) = 0.
      !
      DO ibi=1, nbi
        isea   = isbpi(ibi)
        DO isp=1, nspec
          bbpi0(isp,ibi) = cg(mapwn(isp),isea) / sig2(isp) *      &
               ( rdbpi(ibi,1) * abpi0(isp,ipbpi(ibi,1))   &
               + rdbpi(ibi,2) * abpi0(isp,ipbpi(ibi,2))   &
               + rdbpi(ibi,3) * abpi0(isp,ipbpi(ibi,3))   &
               + rdbpi(ibi,4) * abpi0(isp,ipbpi(ibi,4)) )
        END DO
      END DO
      !
      ! 1.b Shift BBPIN
      !
    ELSE
      bbpi0 = bbpin
    END IF
    !
    ! 2.  Process BBPIN -------------------------------------------------- *
    !
    DO ibi=1, nbi
      isea   = isbpi(ibi)
      DO isp=1, nspec
        bbpin(isp,ibi) = cg(mapwn(isp),isea) / sig2(isp) *          &
             ( rdbpi(ibi,1) * abpin(isp,ipbpi(ibi,1))       &
             + rdbpi(ibi,2) * abpin(isp,ipbpi(ibi,2))       &
             + rdbpi(ibi,3) * abpin(isp,ipbpi(ibi,3))       &
             + rdbpi(ibi,4) * abpin(isp,ipbpi(ibi,4)) )
      END DO
      !
#ifdef W3_RTD
      !!  Rotate the spectra if model is on rotated grid.  JGLi12Jun2012
      !!  PoLat == 90. if the grid is standard lat/lon (C. Hansen 20190613)
      IF ( polat < 90. ) THEN
        spectr = bbpin(:,ibi)
        anglbp = angld(isea)
        CALL  w3acturn( nth, nk, anglbp, spectr )
        bbpin(:,ibi) = spectr
      END IF
 
#endif
      !
    END DO
 
    ! 3.  Wave height test output ---------------------------------------- *
    !
#ifdef W3_T0
    WRITE (ndst,9000)
    DO ibi=1, nbi
      hs1    = 0.
      hs2    = 0.
      DO isp=1, nspec
        hs1    = hs1 + bbpi0(isp,ibi) * dden(mapwn(isp)) /       &
             cg(mapwn(isp),isbpi(ibi))
        hs2    = hs2 + bbpin(isp,ibi) * dden(mapwn(isp)) /       &
             cg(mapwn(isp),isbpi(ibi))
      END DO
      hs1    = 4. * sqrt( hs1 )
      hs2    = 4. * sqrt( hs2 )
      WRITE (ndst,9001) ibi, isbpi(ibi), hs1, hs2
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T0
9000 FORMAT ( ' TEST W3UBPT : WAVE HEIGHTS BBPI0/N (NO TAIL)')
9001 FORMAT ( '         ',2i8,2x,2f8.2)
#endif
 
    !/
    !/ End of W3UBPT ----------------------------------------------------- /
    !/
  END SUBROUTINE w3ubpt
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3uic1( FLFRST )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           C. Sevigny              |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         27-Aug-2015 |
    !/                  +-----------------------------------+
    !/
    !/    27-Aug-2015 : Creation                            ( version 5.10 )
    !/
    !  1. Purpose :
    !
    !     Update ice thickness in the wave model.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !      FLFRST   L.    I     Spectra in 1-D or 2-D representation
    !                           (points to same address).
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for shared / distributed memory architecture.
    !     !/DIST  Id.
    !
    !     !/S  Enable subroutine tracing.
    !     !/T  Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nsea, nsea, mapsf, iicehmin, iicehfac
    USE w3wdatmd, ONLY: time, tic1, iceh
    USE w3idatmd, ONLY: ti1, icep1, flic1
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    LOGICAL, INTENT(IN)     :: FLFRST
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local variables
    !/
    INTEGER                 :: IX, IY, ISEA
    !/
    !/
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Update times
    !
#ifdef W3_T
    WRITE (ndst,9010) time, tic1, ti1
#endif
    tic1(1) = ti1(1)
    tic1(2) = ti1(2)
 
    ! 2.  Main loop over sea points -------------------------------------- *
 
    DO isea=1, nsea
      !
      ix        = mapsf(isea,1)
      iy        = mapsf(isea,2)
      iceh(isea) = max(iicehmin,iicehfac*icep1(ix,iy))
    END DO
    !
    RETURN
#ifdef W3_T
9010 FORMAT ( ' TEST W3UIC1 : TIME     :',i9.8,i7.6/              &
         '               OLD TICE :',i9.8,i7.6/              &
         '               NEW TICE :',i9.8,i7.6)
#endif
    !/
    !/ End of W3UIC1 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3uic1
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3uic5( FLFRST )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |  C. Sevigny  & F. Ardhuin         |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jan-2016 |
    !/                  +-----------------------------------+
    !/
    !/    27-Aug-2015 : Creation                            ( version 5.08 )
    !/    13-Jan-2016 : Changed initial value of ICEDMAX    ( version 5.08 )
    !/
    !  1. Purpose :
    !
    !     Update ice floe mean and max diameters in the wave model.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !      FLFRST   L.    I     Spectra in 1-D or 2-D representation
    !                           (points to same address).
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for shared / distributed memory architecture.
    !     !/DIST  Id.
    !
    !     !/S  Enable subroutine tracing.
    !     !/T  Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3idatmd, ONLY: ti5, icep5
    USE w3gdatmd, ONLY: nsea, mapsf
    USE w3wdatmd, ONLY: time, tic5, ice, iceh, icef, icedmax
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    LOGICAL, INTENT(IN)     :: FLFRST
    !/
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local variables
    !/
    INTEGER                 :: IX, IY, ISEA
    LOGICAL                 :: FLFLOE
    !/
    !/
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Update times
    !
#ifdef W3_T
    WRITE (ndst,9010) time, tic5, ti5
#endif
    tic5(1) = ti5(1)
    tic5(2) = ti5(2)
 
    ! 2.  Main loop over sea points -------------------------------------- *
 
    DO isea=1, nsea
      !
      ix        = mapsf(isea,1)
      iy        = mapsf(isea,2)
      flfloe = ice(isea) .EQ. 0 .OR. iceh(isea) .EQ. 0
      IF ( flfloe) THEN
        icef(isea) = 0.0
        icedmax(isea) = 1000.0
      ELSE
        icef(isea) = icep5(ix,iy)
        icedmax(isea) = icep5(ix,iy)
      END IF
    END DO
    !
    RETURN
#ifdef W3_T
9010 FORMAT ( ' TEST W3UIC5 : TIME     :',i9.8,i7.6/              &
         '               OLD TICE :',i9.8,i7.6/              &
         '               NEW TICE :',i9.8,i7.6)
#endif
 
    !/
    !/
    !/ End of W3UIC5 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3uic5
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3uice ( VA )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         28-Mar-2014 |
    !/                  +-----------------------------------+
    !/
    !/    19-Oct-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    20-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    11-Jan-2002 : Sub-grid ice.                       ( version 2.15 )
    !/    15-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    28-Jun-2005 : Adding MAPST2.                      ( version 3.07 )
    !/                  Taking out initilization.
    !/    11-May-2007 : Adding NTPROC/NAPROC separation.    ( version 3.11 )
    !/    15-May-2010 : Adding second field for icebergs    ( version 3.14 )
    !/    13-Mar-2012 : Add initialization of UST on re-    ( version 4.07 )
    !/                  activation of grid point.
    !/    06-Jun-2012 : Porting bugfixes from 3.14 to 4.07  ( version 4.07 )
    !/    28-Mar-2014 : Adapting to ICx source terms        ( version 4.18 )
    !/
    !  1. Purpose :
    !
    !     Update ice map in the wave model.
    !
    !  2. Method :
    !
    !     Points with an ice concentration larger than FICEN are removed
    !     from the sea map in the wave model. Such points are identified
    !     by negative numbers is the grid status map MAPSTA. For ice
    !     points spectra are set to zero. Points from wich ice disappears
    !     are initialized with a "small" JONSWAP spectrum, based on the
    !     frequency SIG(NK-1) and the local wind direction.
    !
    !     In the case of icebergs, the iceberg attenuation coefficient is
    !     added to the subgrid obstruction map.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !      VA       R.A.  I/O   Spectra in 1-D or 2-D representation
    !                           (points to same address).
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for shared / distributed memory architecture.
    !     !/DIST  Id.
    !
    !     !/S  Enable subroutine tracing.
    !     !/T  Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsf, mapsta, mapst2, &
         nspec, ficen
    USE w3wdatmd, ONLY: time, tice, ice, berg, ust
    USE w3adatmd, ONLY: nsealm, charn
#if defined W3_ST3 || defined(W3_ST4)
    USE w3gdatmd, ONLY: aalpha
#endif
    USE w3idatmd, ONLY: tin, icei, bergi
    USE w3parall, ONLY: init_get_jsea_isproc, init_get_isea
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(INOUT)     :: VA(NSPEC,0:NSEALM)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER                 :: ISEA, JSEA, IX, IY
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER                 :: MAPICE(NY,NX), ISPROC
    LOGICAL                 :: LOCAL
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3UICE')
#endif
    !
    local   = iaproc .LE. naproc
    !
#ifdef W3_T
    WRITE (ndst,9000) ficen
    IF ( .NOT. local ) WRITE (ndst,9001)
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Update times
    !
#ifdef W3_T
    WRITE (ndst,9010) time, tice, tin
#endif
    tice(1) = tin(1)
    tice(2) = tin(2)
    !
    ! 1.b Process maps
    !
#ifdef W3_IC0
    mapice = mod(mapst2,2)
    mapst2 = mapst2 - mapice
#endif
    !
    ! 2.  Main loop over sea points -------------------------------------- *
    !
    DO isea=1, nsea
      !
      ! 2.a Get grid counters
      !
      ix        = mapsf(isea,1)
      iy        = mapsf(isea,2)
      ice(isea) = icei(ix,iy)
      berg(isea)= bergi(ix,iy)
      !
      ! 2.b Sea point to be de-activated..
      !
#ifdef W3_IC0
      IF ( icei(ix,iy).GE.ficen .AND. mapice(iy,ix).EQ.0 ) THEN
        mapsta(iy,ix) = - abs(mapsta(iy,ix))
        mapice(iy,ix) = 1
        CALL init_get_jsea_isproc(isea, jsea, isproc)
        IF (local .AND. (iaproc .eq. isproc)) THEN
#ifdef W3_T
          WRITE (ndst,9021) isea, ix, iy, mapsta(iy,ix),     &
               icei(ix,iy), 'ICE (NEW)'
#endif
          va(:,jsea) = 0.
#if defined W3_ST3 || defined(W3_ST4)
          charn(jsea) = aalpha
#else
          charn(jsea) = 0.
#endif
#ifdef W3_T
        ELSE
          WRITE (ndst,9021) isea, ix, iy, mapsta(iy,ix),     &
               icei(ix,iy), 'ICE (NEW X)'
#endif
        END IF
 
#ifdef W3_T
      ELSE IF ( icei(ix,iy).GE.ficen ) THEN
        WRITE (ndst,9021) isea, ix, iy, mapsta(iy,ix),         &
             icei(ix,iy), 'ICE'
#endif
      END IF
      !
      ! 2.b Ice point to be re-activated.
      !
      IF ( icei(ix,iy).LT.ficen .AND. mapice(iy,ix).EQ.1 ) THEN
 
        mapice(iy,ix) = 0
        ust(isea)     = 0.05
 
        IF ( mapst2(iy,ix) .EQ. 0 ) THEN
          mapsta(iy,ix) = abs(mapsta(iy,ix))
 
          CALL init_get_jsea_isproc(isea, jsea, isproc)
          IF ( local .AND. (iaproc .eq. isproc) ) THEN
#ifdef W3_T
            WRITE (ndst,9021) isea, ix, iy, mapsta(iy,ix), &
                 icei(ix,iy), 'SEA (NEW)'
#endif
            va(:,jsea) = 0.
#if defined W3_ST3 || defined(W3_ST4)
            charn(jsea) = aalpha
#else
            charn(jsea) = 0.
#endif
#ifdef W3_T
          ELSE
            WRITE (ndst,9021) isea, ix, iy, mapsta(iy,ix), &
                 icei(ix,iy), 'SEA (NEW X)'
#endif
          END IF
 
#ifdef W3_T
        ELSE
          WRITE (ndst,9021) isea, ix, iy, mapsta(iy,ix),     &
               icei(ix,iy), 'DIS'
#endif
        END IF
 
#ifdef W3_T
      ELSE IF ( icei(ix,iy).LT.ficen ) THEN
        WRITE (ndst,9021) isea, ix, iy, mapsta(iy,ix),     &
             icei(ix,iy), 'SEA'
#endif
 
      END IF
#endif
 
    END DO
    !
    ! 3.  Update MAPST2 -------------------------------------------------- *
    !
#ifdef W3_IC0
    mapst2 = mapst2 + mapice
#endif
    !
    RETURN
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3UICE : FICEN    :',f9.3)
9001 FORMAT ( ' TEST W3UICE : NO LOCAL SPECTRA')
9010 FORMAT ( ' TEST W3UICE : TIME     :',i9.8,i7.6/              &
         '               OLD TICE :',i9.8,i7.6/              &
         '               NEW TICE :',i9.8,i7.6)
9020 FORMAT ( ' TEST W3UICE : ISEA, IX, IY, MAP, ICE, STATUS :')
9021 FORMAT ( '           ',i8,3i4,f6.2,2x,a)
#endif
    !/
    !/ End of W3UICE ----------------------------------------------------- /
    !/
  END SUBROUTINE w3uice
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3ulev ( A, VA )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         26-Sep-2012 |
    !/                  +-----------------------------------+
    !/
    !/    15-Jan-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    21-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    30-Apr-2002 : Water level fixes.                  ( version 2.20 )
    !/    15-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    15-Jul-2005 : Adding drying out of points.        ( version 3.07 )
    !/    11-May-2007 : Adding NTPROC/NAPROC separation.    ( version 3.11 )
    !/    23-Aug-2011 : Bug fix for UG grids : new boundary ( version 4.04 )
    !/    13-Mar-2012 : Add initialization of UST on re-    ( version 4.07 )
    !/                  activation of grid point.
    !/    06-Jun-2012 : Porting bugfixes from 3.14 to 4.07  ( version 4.07 )
    !/    26-Sep-2012 : Adding update from tidal analysis   ( version 4.08 )
    !/
    !  1. Purpose :
    !
    !     Update the water level.
    !
    !  2. Method :
    !
    !     The wavenumber grid is modified without modyfying the spectrum
    !     (conservative linear interpolation to new grid).
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !      (V)A     R.A.  I/O  2-D and 1-D represetation of the spectra.
    !     ----------------------------------------------------------------
    !
    !     Local variables
    !     ----------------------------------------------------------------
    !       KDMAX   Real   Deep water cut-off for kd.
    !       WNO     R.A.   Old wavenumbers.
    !       CGO     R.A.   Old group velocities.
    !       OWN     R.A.   Old wavenumber band width.
    !       DWN     R.A.   New wavenumber band width.
    !       TA      R.A.   Auxiliary spectrum.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - The grid is updated only if KDmin > KDMAX.
    !     - The grid is updated for inactive points too.
    !     - The local wavenumber bandwidth is DSIGMA/CG.
    !     - The local spectrum is updated only if the grid is updated,
    !       the grid point is not disabled (MAPST2) and if the change of
    !       the lowest wavenumber exceeds RDKMIN times the band width.
    !     - No spectral initialization for newly wet points.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S     Enable subroutine tracing.
    !     !/T     Basic test output.
    !     !/T2    Output of minimum relative depth per grid point.
    !     !/T3    Spectra before and after
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nx, ny, nsea, nseal, mapsf, mapsta, mapst2, &
         zb, dmin, nk, nth, nspec, sig, dsip,        &
         mapwn, mapth, fachfa, gtype, ungtype, w3setref
    USE w3wdatmd, ONLY: time, tlev, wlv, ust
    USE w3adatmd, ONLY: cg, wn, dw, hs
    USE w3idatmd, ONLY: tln, wlev
    USE w3servmd, ONLY: extcde
    USE w3dispmd, ONLY: wavnu1
    USE w3timemd
    USE w3parall, only : init_get_jsea_isproc, init_get_isea
    USE w3parall, only : get_jsea_ibelong
#ifdef W3_PDLIB
    USE w3dispmd, ONLY: wavnu3
    USE pdlib_w3profsmd, ONLY : set_iobdp_pdlib
#endif
#ifdef W3_TIDE
    USE w3gdatmd, ONLY: ygrd
    USE w3idatmd, ONLY: fllevtide, wltide, ntide
    USE w3tidemd
#endif
#ifdef W3_SETUP
    USE w3wdatmd, ONLY: zeta_setup
    USE w3gdatmd, ONLY : do_change_wlv
#endif
 
 
#ifdef W3_T3
    USE w3arrymd, ONLY: prt2ds
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(INOUT)     :: A(NTH,NK,0:NSEAL), VA(NSPEC,0:NSEAL)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER                 :: ISEA, JSEA, IX, IY, IK, I1, I2,      &
         ISPEC, IK0, ITH
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER                 :: MAPDRY(NY,NX), ISPROC
    REAL                    :: DWO(NSEA), KDCHCK, WNO(0:NK+1),      &
         CGO(0:NK+1), DEPTH,                  &
         RDK, RD1, RD2, TA(NTH,NK),           &
         OWN(NK), DWN(NK)
    REAL                    :: KDMAX = 4., rdkmin = 0.05
    REAL                    :: WLVeff
#ifdef W3_T3
    REAL                    :: OUT(NK,NTH)
#endif
    LOGICAL                 :: LOCAL
    INTEGER                 :: IBELONG
    !
#ifdef W3_TIDE
    INTEGER          :: J
    INTEGER(KIND=4)  :: TIDE_KD0, INT24, INTDYS       ! "Gregorian day constant"
    REAL             :: WLEVTIDE, TIDE_ARG, WLEVTIDE2(1)
    REAL(KIND=8)     :: d1,h,tide_hour,hh,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau
    REAL             :: FX(44),UX(44),VX(44)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3ULEV')
#endif
    !
    local   = iaproc .LE. naproc
    !
#ifdef W3_T
    WRITE (ndst,9000) kdmax, rdkmin
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Check NK
    !
    IF ( nk .LT. 2 ) THEN
      IF ( iaproc .EQ. naperr ) WRITE (ndse,1000)
      CALL extcde ( 1 )
    END IF
    !
    ! 1.b Update times
    !
#ifdef W3_T
    WRITE (ndst,9010) time, tlev
#endif
    tlev = tln
#ifdef W3_T
    WRITE (ndst,9011) tlev
#endif
    !
    ! 1.c Extract dry point map, and residual MAPST2
    !
    mapdry = mod(mapst2/2,2)
    mapst2 = mapst2 - 2*mapdry
    !
    ! 1.d Update water levels and save old
    !
#ifdef W3_TIDE
    IF (fllevtide) THEN
      !          WRITE(6,*) 'TIME:',TIME
      tide_hour = time2hours(time)
      !
      !*  THE ASTRONOMICAL ARGUMENTS ARE CALCULATED BY LINEAR APPROXIMATION
      !*  AT THE MID POINT OF THE ANALYSIS PERIOD.
      d1=tide_hour/24.d0
      tide_kd0= 2415020
      d1=d1-dfloat(tide_kd0)-0.5d0
      call astr(d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp)
      int24=24
      intdys=int((tide_hour+0.00001)/int24)
      hh=tide_hour-dfloat(intdys*int24)
      tau=hh/24.d0+h-s
    END IF
    !
    !  ONLY THE FRACTIONAL PART OF A SOLAR DAY NEED BE RETAINED FOR COMPU-
    !  TING THE LUNAR TIME TAU.
    !
#endif
    DO isea=1, nsea
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      dwo(isea) = dw(isea)
      !
#ifdef W3_TIDE
      IF (fllevtide) THEN
        ! VUF should be updated only if latitude changes significantly ...
        CALL setvuf_fast(h,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau,real(ygrd(iy,ix)),fx,ux,vx)
        wlevtide = wltide(ix,iy,1,1)
        !Verification
        !          IF (ISEA.EQ.1) THEN
 
        tide_ampc(1:ntide,1)=wltide(ix,iy,1:ntide,1)
        tide_phg(1:ntide,1)=wltide(ix,iy,1:ntide,2)
        !
        !           WRITE(991,'(A,F20.2,13F8.3)') 'TEST ISEA 0:',    &
        !                       d1,H,S,TAU,pp,s,p,enp,dh,dpp,ds,dp,dnp,YGRD(IY,IX)
        j=1
        !           WRITE(991,'(A,4I9,F12.0,3F8.3,I4,X,A)') 'TEST ISEA 1:',IX,J,TIME,TIDE_HOUR,    &
        !                       FX(J),UX(J),VX(J),TIDE_INDEX2(J),TIDECON_ALLNAMES(TIDE_INDEX2(J))
        DO j=2,tide_mf
          tide_arg=(vx(j)+ux(j))*twpi-wltide(ix,iy,j,2)*dera
          wlevtide =wlevtide+fx(j)*wltide(ix,iy,j,1)*cos(tide_arg)
          !           WRITE(991,'(A,4I9,F12.0,3F8.3,I4,X,A)') 'TEST ISEA 1:',IX,J,TIME,TIDE_HOUR,    &
          !                       FX(J),UX(J),VX(J),TIDE_INDEX2(J),TIDECON_ALLNAMES(TIDE_INDEX2(J))
        END DO
        DO j=1,tide_mf
          !           WRITE(991,'(A,4I9,F12.0,5F8.3)') 'TEST ISEA 2:',IX,J,TIME,TIDE_HOUR,    &
          !                       FX(J),UX(J),VX(J),TIDE_AMPC(J,1),TIDE_PHG(J,1)
        END DO
        !         WRITE(991,'(A,3F7.3)') '#:',WLEV(IX,IY),WLEVTIDE,WLEV(IX,IY)-WLEVTIDE
#endif
 
#ifdef W3_TIDE
        !         CLOSE(991)
        !         END IF
        ! End of verification
        wlv(isea) = wlevtide
      ELSE
#endif
        !
        wlv(isea) = wlev(ix,iy)
        wlveff    = wlv(isea)
 
#ifdef W3_SETUP
       IF (do_change_wlv) THEN
         wlveff    = wlveff + zeta_setup(isea)
         wlv(isea) = wlveff
       END IF
#endif
#ifdef W3_TIDE
      ENDIF
#endif
      dw(isea) = max( 0. , wlveff-zb(isea) )
 
    END DO ! NSEA
 
    !
    ! 2.  Loop over all sea points --------------------------------------- *
    !
#ifdef W3_T2
    WRITE (ndst,9020)
#endif
    !
    DO isea=1, nsea
      !
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      !
      ! 2.a Check if deep water
      !
      kdchck = wn(1,isea) * min( dwo(isea) , dw(isea) )
      IF ( kdchck .LT. kdmax ) THEN
        !
        ! 2.b Update grid and save old grid
        !
        depth  = max( dmin, dw(isea) )
        !
        DO ik=0, nk+1
          wno(ik) = wn(ik,isea)
          cgo(ik) = cg(ik,isea)
          !
          !   Calculate wavenumbers and group velocities.
#ifdef W3_PDLIB
              CALL wavnu3(sig(ik),depth,wn(ik,isea),cg(ik,isea))
#else
              CALL wavnu1(sig(ik),depth,wn(ik,isea),cg(ik,isea))
#endif
        END DO
        !
        DO ik=1, nk
          own(ik) = dsip(ik) / cgo(ik)
          dwn(ik) = dsip(ik) / cg(ik,isea)
        END DO
        !
        ! 2.c Process dry points
        !
        IF ( wlv(isea)-zb(isea) .LE.0. ) THEN
          IF ( mapdry(iy,ix) .EQ. 0 ) THEN
            CALL get_jsea_ibelong(isea, jsea, ibelong)
            IF ( local .AND. (ibelong .eq. 1) ) THEN
              va(:,jsea) = 0.
            END IF
            mapdry(iy,ix) = 1
            mapsta(iy,ix) = -abs(mapsta(iy,ix))
#ifdef W3_T2
            WRITE (ndst,9021) isea, wlv(isea)-zb(isea),   &
                 0., 0., '  (NEW DRY)'
          ELSE
            WRITE (ndst,9021) isea, wlv(isea)-zb(isea),   &
                 0., 0., '  (DRY)'
#endif
          ENDIF
          cycle
        END IF
        !
        ! 2.d Process new wet point
        !
        IF (wlv(isea)-zb(isea).GT.0. .AND. mapdry(iy,ix).EQ.1) THEN
          mapdry(iy,ix) = 0
          !
          ! Resets the spectrum to zero
          !
          CALL get_jsea_ibelong(isea, jsea, ibelong)
          IF ( local .AND. (ibelong .eq. 1) ) THEN
            va(:,jsea) = 0.
          END IF
          !
          ust(isea)     = 0.05
          IF ( mapst2(iy,ix) .EQ. 0 ) THEN
            mapsta(iy,ix) = abs(mapsta(iy,ix))
#ifdef W3_T2
            WRITE (ndst,9021) isea, wlv(isea)-zb(isea),   &
                 0., 0., '  (NEW WET)'
          ELSE
            WRITE (ndst,9021) isea, wlv(isea)-zb(isea),   &
                 0., 0., '  (NEW WET INACTIVE)'
#endif
          END IF
          cycle
        END IF
        !
        ! 2.e Check if ice on grid point, or if grid changes negligible
        !
        rdk    = abs(wno(1)-wn(1,isea)) / dwn(1)
        !
#ifdef W3_T2
        IF ( mapsta(iy,ix) .LT. 0 ) THEN
          WRITE (ndst,9021)                                 &
               isea, dw(isea), kdchck, rdk, '  (INACTIVE)'
        ELSE IF ( rdk .LT. rdkmin ) THEN
          WRITE (ndst,9021)                                 &
               isea, dw(isea), kdchck, rdk, '  (NEGL)'
        ELSE
          WRITE (ndst,9021)                                 &
               isea, dw(isea), kdchck, rdk, ' '
        END IF
#endif
        !
        IF ( rdk.LT.rdkmin .OR. mapsta(iy,ix).LT.0 ) cycle
        CALL get_jsea_ibelong(isea, jsea, ibelong)
        IF ( ibelong .eq. 0) cycle
        !
        IF ( .NOT. local ) cycle
        !
        ! 2.d Save discrete actions and clean spectrum
        !
        DO ik=1, nk
          DO ith=1, nth
#ifdef W3_T3
            out(ik,ith) = a(ith,ik,jsea) * sig(ik) / cgo(ik)
#endif
            ta(ith,ik) = a(ith,ik,jsea) * own(ik)
          END DO
        END DO
        !
        va(:,jsea) = 0.
        !
#ifdef W3_T3
        CALL prt2ds ( ndst, nk, nk, nth, out, sig, ' ',    &
             tpi, 0., 1.e-5, 'F(f,th)', 'm2s', 'Before' )
#endif
        !
        ! 2.e Redistribute discrete action density
        !
        IF ( wno(1) .LT. wn(1,isea) ) THEN
          ik0    = 1
          i1     = 0
          i2     = 1
220       CONTINUE
          ik0    = ik0 + 1
          IF ( ik0 .GT. nk+1 ) GOTO 251
          IF ( wno(ik0) .GE. wn(1,isea) ) THEN
            ik0    = ik0 - 1
          ELSE
            GOTO 220
          END IF
        ELSE
          ik0    = 1
          i1     = 1
          i2     = 2
        END IF
        !
        DO ik=ik0, nk
          !
230       CONTINUE
          IF ( wno(ik) .GT. wn(i2,isea) ) THEN
            i1     = i1 + 1
            IF ( i1 .GT. nk ) GOTO 250
            i2     = i1 + 1
            GOTO 230
          END IF
          !
          IF ( i1 .EQ. 0 ) THEN
            rd1    = ( wn(1,isea) - wno(ik) ) / dwn(1)
            rd2    = 1. - rd1
          ELSE
            rd1    = ( wn(i2,isea) - wno(ik) ) /                &
                 ( wn(i2,isea) - wn(i1,isea) )
            rd2    = 1. - rd1
          END IF
          !
          IF ( i1 .GE. 1 ) THEN
            DO ith=1, nth
              a(ith,i1,jsea) = a(ith,i1,jsea) + rd1*ta(ith,ik)
            END DO
          END IF
          !
          IF ( i2 .LE. nk ) THEN
            DO ith=1, nth
              a(ith,i2,jsea) = a(ith,i2,jsea) + rd2*ta(ith,ik)
            END DO
          END IF
          !
250       CONTINUE
        END DO
251     CONTINUE
        !
        ! 2.f Convert discrete action densities to spectrum
        !
        DO ispec=1, nspec
          va(ispec,jsea) = va(ispec,jsea) / dwn(mapwn(ispec))
        END DO
        !
        ! 2.f Add tail if necessary
        !
        IF ( i2.LE.nk .AND. rd2.LE.0.95 ) THEN
          DO ik=max(i2,2), nk
            DO ith=1, nth
              a(ith,ik,jsea) = fachfa * a(ith,ik-1,jsea)
            END DO
          END DO
        END IF
        !
#ifdef W3_T3
        DO ispec=1, nspec
          ik          = mapwn(ispec)
          ith         = mapth(ispec)
          out(ik,ith) = a(ith,ik,jsea) * sig(ik) / cg(ik,isea)
        END DO
#endif
        !
#ifdef W3_T3
        CALL prt2ds ( ndst, nk, nk, nth, out, sig, ' ',    &
             tpi, 0., 1.e-5, 'F(f,th)', 'm2s', 'After' )
#endif
        !
#ifdef W3_T2
      ELSE
        WRITE (ndst,9021) isea, kdchck, '  (DEEP)'
#endif
      END IF
      !
    END DO ! NSEA
    !
    ! 3.  Reconstruct new MAPST2 ----------------------------------------- *
    !
    mapst2 = mapst2 + 2*mapdry
    !
    ! 4. Re-generates the boundary data ---------------------------------- *
    !
    IF (gtype.EQ.ungtype) THEN
#ifdef W3_PDLIB
      CALL set_iobdp_pdlib
#endif
#ifdef W3_REF1
    ELSE
      CALL w3setref
#endif
    ENDIF
    !
    RETURN
    !
    ! Formats
    !
1000 FORMAT (/' *** ERROR W3ULEV *** '/                              &
         '     THIS ROUTINE REQUIRES NK > 1 '/)
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3ULEV : KDMAX    :',f6.1/                   &
         '               RDKMIN   :',f8.3)
#endif
    !
#ifdef W3_T
9010 FORMAT ( ' TEST W3ULEV : TIME     :',i9.8,i7.6/              &
         '               OLD TLEV :',i9.8,i7.6)
9011 FORMAT ( '               NEW TLEV :',i9.8,i7.6)
#endif
    !
#ifdef W3_T2
9020 FORMAT ( ' TEST W3ULEV : LOOP OVER ALL POINTS:',            &
         ' ISEA, DW, KDMIN, RDK : ')
9021 FORMAT ( '             ',i6,f8.2,f6.2,f7.3,a)
#endif
    !/
    !/ End of W3ULEV ----------------------------------------------------- /
    !/
  END SUBROUTINE w3ulev
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3urho ( FLFRST )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           J. M. Castillo          |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Aug-2021 |
    !/                  +-----------------------------------+
    !/
    !/    22-Mar-2021 : First implementation                ( version 7.13 )
    !/    13-Aug-2021 : Enable time interpolation           ( version 7.14 )
    !/
    !  1. Purpose :
    !
    !     Interpolate air density field to the given time.
    !
    !  2. Method :
    !
    !     Linear interpolation.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       FLFRST  Log.  I   Flag for first pass through routine.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - Only air density over sea points is considered.
    !     - Time ranges checked in W3WAVE.
    !
    !  8. Structure :
    !
    !     --------------------------------------
    !      1.  Prepare auxiliary arrays
    !      2.  Calculate interpolation factors
    !      3.  Get actual air density
    !     --------------------------------------
    !
    !  9. Switches :
    !
    !     !/OMPG   OpenMP compiler directives
    !
    !     !/WNT0   No air density interpolation.
    !     !/WNT1   Linear air density interpolation.
    !     !/WNT2   Linear air density interpolation (and energy conservation for momentum).
    !
    !     !/S  Enable subroutine tracing.
    !     !/T  Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nsea, mapsf
#ifdef W3_SMC
    USE w3gdatmd, ONLY: fswnd
#endif
    USE w3wdatmd, ONLY: time, trho, rhoair
    USE w3idatmd, ONLY: tr0, trn, rh0, rhn
    USE w3adatmd, ONLY: ra0, rai
    USE w3odatmd, ONLY: iaproc, naproc
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    LOGICAL, INTENT(IN)     :: FLFRST
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER                 :: ISEA, IX, IY
#ifdef W3_S
    INTEGER, SAVE            :: IENT = 0
#endif
    REAL                    :: DT0N, DT0T, RD
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3URHO')
#endif
    !
    ! 1.  Prepare auxiliary arrays
    !
    IF ( flfrst ) THEN
      DO isea=1, nsea
#ifdef W3_SMC
        !!Li  For sea-point only SMC grid air density is stored on
        !!Li  2-D RH0(NSEA, 1) variable.
        IF( fswnd ) THEN
          ix = isea
          iy = 1
        ELSE
#endif
          ix        = mapsf(isea,1)
          iy        = mapsf(isea,2)
#ifdef W3_SMC
        ENDIF
#endif
 
        ra0(isea) = rh0(ix,iy)
        rai(isea) = rhn(ix,iy) - rh0(ix,iy)
      END DO
    END IF
    !
    ! 2.  Calculate interpolation factor
    !
    dt0n   = dsec21( tr0, trn )
    dt0t   = dsec21( tr0, time )
    !
#ifdef W3_WNT0
    rd     = 0.
#endif
#ifdef W3_WNT1
    rd     = dt0t / max( 1.e-7 , dt0n )
#endif
#ifdef W3_WNT2
    rd     = dt0t / max( 1.e-7 , dt0n )
#endif
#ifdef W3_OASACM
    rd     = 1.
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) dt0n, dt0t, rd
#endif
    !
    ! 3.  Actual momentum for all grid points
    !
#ifdef W3_OMPG
    !$OMP PARALLEL DO PRIVATE (ISEA,RA0,RAI)
#endif
    !
    DO isea=1, nsea
      !
      rhoair(isea) = ra0(isea) + rd * rai(isea)
      !
    END DO
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3URHO : DT0N, DT0T, RD :',2f8.1,f6.3)
#endif
    !/
    !/ End of W3URHO ----------------------------------------------------- /
    !/
  END SUBROUTINE w3urho
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3utrn ( TRNX, TRNY )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         30-Oct-2009 |
    !/                  +-----------------------------------+
    !/
    !/    02-Apr-2001 : Origination.                        ( version 2.10 )
    !/    11-Jan-2002 : Sub-grid ice.                       ( version 2.15 )
    !/    30-Apr-2002 : Change to ICE on storage grid.      ( version 2.20 )
    !/    15-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    11-Jan-2007 : Clean-up for boundary points.       ( version 3.10 )
    !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/
    !  1. Purpose :
    !
    !     Update cell boundary transparencies for general use in propagation
    !     routines.
    !
    !  2. Method :
    !
    !     Two arrays are generated with the size (NY*NX,-1:1). The value
    !     at (IXY,-1) indicates the transparency to be used if the lower
    !     or left boundary is an inflow boundary. (IXY,1) is used if the
    !     upper or right boundary is an inflow boundary. (IXY,0) is used
    !     for all other cases (by definition full transparency).
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       TRNX/Y  R.A.   I   Transparencies from model defintion file.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !       !/T      Basic test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsta, mapsf,                &
         trflag, fice0, ficen, ficel,                &
         rlgtype, clgtype, gtype, flagll,            &
         hpfac, hqfac, ffacberg
    USE w3wdatmd, ONLY: ice, berg
    USE w3adatmd, ONLY: atrnx, atrny
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: TRNX(NY*NX), TRNY(NY*NX)
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER                 :: ISEA, IX, IY, IXY, IXN, IXP, IYN, IYP
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T
    INTEGER                 :: ILEV, NLEV
#endif
 
    REAL                    :: TRIX(NY*NX), TRIY(NY*NX), DX, DY,    &
         LICE0, LICEN
#ifdef W3_T
    REAL                    :: LEVS(0:10)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3UTRN')
#endif
#ifdef W3_T
    WRITE (ndst,9000) trflag
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    !
    atrnx = 1.
    atrny = 1.
#ifdef W3_T
    WRITE (ndst,9001) 'INITIALIZING ATRNX/Y'
#endif
    !
    ! 2.  Filling arrays from TRNX/Y for obstructions -------------------- *
    ! 2.a TRFLAG = 0, no action needed
    IF ( trflag .EQ. 0 ) THEN
#ifdef W3_T
      WRITE (ndst,9001) 'NO FURTHER ACTION REQUIRED'
#endif
      RETURN
      !
      ! 2.b TRFLAG = 1,3: TRNX/Y defined at boundaries
      !
    ELSE IF ( trflag.EQ.1 .OR. trflag.EQ.3 .OR. trflag.EQ.5 ) THEN
#ifdef W3_T
      WRITE (ndst,9001) 'DATA APPLIED AT CELL BOUNDARIES'
      levs   = 0.
#endif
      !
      DO isea=1, nsea
        !
        ix            = mapsf(isea,1)
        iy            = mapsf(isea,2)
        ixy           = mapsf(isea,3)
        IF ( ix .EQ. 1 ) THEN
          atrnx(ixy,-1) = trnx(iy+(nx-1)*ny)
          atrnx(ixy, 1) = trnx(ixy)
        ELSE IF ( ix .EQ. nx ) THEN
          atrnx(ixy,-1) = trnx(ixy-ny)
          atrnx(ixy, 1) = trnx(iy)
        ELSE
          atrnx(ixy,-1) = trnx(ixy-ny)
          atrnx(ixy, 1) = trnx(ixy)
        END IF
        atrny(ixy,-1) = trny(ixy-1)
        atrny(ixy, 1) = trny(ixy)
        !
#ifdef W3_T
        ilev          = nint(10.*min(trnx(ixy),trny(ixy)))
        levs(ilev)    = levs(ilev) + 1.
#endif
        !
      END DO
      !
      ! 2.c TRFLAG = 2,4: TRNX/Y defined at cell centers
      !
    ELSE
#ifdef W3_T
      WRITE (ndst,9001) 'DATA APPLIED AT CELL CENTERS'
      levs   = 0.
#endif
      !
      DO isea=1, nsea
        !
        ix            = mapsf(isea,1)
        iy            = mapsf(isea,2)
        ixy           = mapsf(isea,3)
        !
        IF ( ix .EQ. 1 ) THEN
          ixn    = iy + (nx-1)*ny
          ixp    = iy +  ix   *ny
        ELSE IF ( ix .EQ. nx ) THEN
          ixn    = iy + (ix-2)*ny
          ixp    = iy
        ELSE
          ixn    = iy + (ix-2)*ny
          ixp    = iy +  ix   *ny
        END IF
        !
        IF ( iy .EQ. 1 ) THEN
          iyn    = ixy
          iyp    = ixy + 1
        ELSE IF ( iy .EQ. ny ) THEN
          iyn    = ixy - 1
          iyp    = ixy
        ELSE
          iyn    = ixy - 1
          iyp    = ixy + 1
        END IF
        !
        ! factors 0.5 in first term and 2. in second term cancel
        !
        atrnx(ixy,-1) = (1.+trnx(ixy)) * trnx(ixn)/(1.+trnx(ixn))
        atrnx(ixy, 1) = (1.+trnx(ixy)) * trnx(ixp)/(1.+trnx(ixp))
        atrny(ixy,-1) = (1.+trny(ixy)) * trny(iyn)/(1.+trny(iyn))
        atrny(ixy, 1) = (1.+trny(ixy)) * trny(iyp)/(1.+trny(iyp))
        !
        IF ( mapsta(iy,ix) .EQ. 2 ) THEN
          IF ( ix .EQ. 1  ) THEN
            atrnx(ixy,-1) = 1.
          ELSE IF ( mapsta( iy ,ix-1) .LE. 0 ) THEN
            atrnx(ixy,-1) = 1.
          END IF
          IF ( ix .EQ. nx ) THEN
            atrnx(ixy, 1) = 1.
          ELSE IF ( mapsta( iy ,ix+1) .LE. 0 ) THEN
            atrnx(ixy, 1) = 1.
          END IF
          IF ( iy .EQ. 1  ) THEN
            atrny(ixy,-1) = 1.
          ELSE IF ( mapsta(iy-1, ix ) .LE. 0 ) THEN
            atrny(ixy,-1) = 1.
          END IF
          IF ( iy .EQ. ny ) THEN
            atrny(ixy, 1) = 1.
          ELSE IF ( mapsta(iy+1, ix ) .LE. 0 ) THEN
            atrny(ixy, 1) = 1.
          END IF
        END IF
        !
#ifdef W3_T
        ilev          = nint(10.*min(trnx(ixy),trny(ixy)))
        levs(ilev)    = levs(ilev) + 1.
#endif
        !
      END DO
    END IF
    !
#ifdef W3_T
    WRITE(ndst,9010) 'ISLANDS'
    nlev   = 0
    DO ilev=0, 10
      WRITE (ndst,9011) ilev, levs(ilev)/real(nsea)
      nlev = nlev + nint(levs(ilev))
    END DO
#endif
    !
    ! 3.  Adding ice to obstructions ------------------------------------- *
    ! 3.a TRFLAG < 3, no action needed
    !
    IF ( trflag.LT.3 .OR. ficen-fice0.LT.1.e-6 ) THEN
#ifdef W3_T
      WRITE (ndst,9001) 'NO ICE ACTION REQUIRED'
#endif
      RETURN
      !
      ! 3.b TRFLAG = 3,4: Calculate ice transparencies
      !
    ELSE
#ifdef W3_T
      WRITE (ndst,9001) 'CALCULATE ICE TRANSPARENCIES'
      levs   = 0.
#endif
      trix   = 1.
      triy   = 1.
      !
      DO isea=1, nsea
        !
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
        !
        dx     = hpfac(iy,ix)
        dy     = hqfac(iy,ix)
        IF ( flagll ) THEN
          dx     = dx * radius * dera
          dy     = dy * radius * dera
        END IF
 
        !
#ifdef W3_IC0
        IF (ice(isea).GT.0) THEN
          IF (ficel.GT.0.) THEN
            trix(ixy) = exp(-ice(isea)*dx/ficel)
            triy(ixy) = exp(-ice(isea)*dy/ficel)
          ELSE
#endif
            ! Otherwise: original Tolman expression (Tolman 2003)
#ifdef W3_IC0
            lice0  = fice0*dx
            licen  = ficen*dx
            trix(ixy) = ( licen - ice(isea)*dx ) / ( licen - lice0 )
#endif
 
            ! begin temporary notes
            !                TRIX = (   LICEN    - ICE(ISEA)*DX ) / (   LICEN -      LICE0 )
            !    thus, it is TRIX=  ( (FICEN*DX) - ICE(ISEA)*DX ) / ( (FICEN*DX) - (FICE0*DX) )
            !    thus, it is TRIX=  (   FICEN -    ICE(ISEA) )   /  (  FICEN     -   FICE0 )
            !    in other words, the variables DX DY are not used
            !                    and the variables LICE0 LICEN are not necessary.
            ! end temporary notes
 
#ifdef W3_IC0
            lice0  = fice0*dy
            licen  = ficen*dy
            triy(ixy) = ( licen - ice(isea)*dy ) / ( licen - lice0 )
          END IF
#endif
          !
#ifdef W3_IC0
          trix(ixy) = max( 0. , min( 1. , trix(ixy) ) )
          triy(ixy) = max( 0. , min( 1. , triy(ixy) ) )
        END IF
#endif
        !
        !  Adding iceberg attenuation
        !
        IF (berg(isea).GT.0) THEN
          trix(ixy) = trix(ixy)*exp(-berg(isea)*ffacberg  *dx*0.0001)
          triy(ixy) = triy(ixy)*exp(-berg(isea)*ffacberg  *dy*0.0001)
        END IF
        !
#ifdef W3_T
        ilev          = nint(10.*min(trix(ixy),triy(ixy)))
        levs(ilev)    = levs(ilev) + 1.
#endif
        !
      END DO
      !
#ifdef W3_T
      WRITE(ndst,9010) 'ICE'
      nlev   = 0
      DO ilev=0, 10
        WRITE (ndst,9011) ilev, levs(ilev)/real(nsea)
        nlev = nlev + nint(levs(ilev))
      END DO
#endif
      !
      ! 3.c Combine transparencies, ice always defined at cell center !
      !
      DO isea=1, nsea
        !
        ix            = mapsf(isea,1)
        iy            = mapsf(isea,2)
        ixy           = mapsf(isea,3)
        !
        IF ( ix .EQ. 1 ) THEN
          ixn    = iy + (nx-1)*ny
          ixp    = iy +  ix   *ny
        ELSE IF ( ix .EQ. nx ) THEN
          ixn    = iy + (ix-2)*ny
          ixp    = iy
        ELSE
          ixn    = iy + (ix-2)*ny
          ixp    = iy +  ix   *ny
        END IF
        !
        IF ( iy .EQ. 1 ) THEN
          iyn    = ixy
          iyp    = ixy + 1
        ELSE IF ( iy .EQ. ny ) THEN
          iyn    = ixy - 1
          iyp    = ixy
        ELSE
          iyn    = ixy - 1
          iyp    = ixy + 1
        END IF
        !
        atrnx(ixy,-1) = atrnx(ixy,-1)                             &
             * (1.+trix(ixy)) * trix(ixn)/(1.+trix(ixn))
        atrnx(ixy, 1) = atrnx(ixy, 1)                             &
             * (1.+trix(ixy)) * trix(ixp)/(1.+trix(ixp))
        atrny(ixy,-1) = atrny(ixy,-1)                             &
             * (1.+triy(ixy)) * triy(iyn)/(1.+triy(iyn))
        atrny(ixy, 1) = atrny(ixy, 1)                             &
             * (1.+triy(ixy)) * triy(iyp)/(1.+triy(iyp))
        !
      END DO
      !
    END IF
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3UTRN : TRFLAG = ',i3)
9001 FORMAT ( ' TEST W3UTRN : ',a)
9010 FORMAT ( ' TEST W3UTRN : OBSTRICTION LEVELS FOR ',a,' :')
9011 FORMAT ( '             ',i4,f8.5)
#endif
    !/
    !/ End of W3UTRN ----------------------------------------------------- /
    !/
  END SUBROUTINE w3utrn
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3dzxy( ZZ, ZUNIT, DZZDX, DZZDY )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |         W. E. Rogers, NRL         |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         06-Dec-2010 |
    !/                  +-----------------------------------+
    !/
    !/    30-Oct-2009 : Origination.                        ( version 3.14 )
    !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
    !/                  specify index closure for a grid.   ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/
    !  1. Purpose :
    !
    !     Calculate derivatives of a field.
    !
    !  2. Method :
    !
    !     Derivatives are calculated in m/m from the longitude/latitude
    !     grid, central in space for iternal points, one-sided for
    !     coastal points.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ZZ      R.A.  I   Field to calculate derivatives of.
    !       ZUNIT   R.A.  I   Units of ZZ (used for test output).
    !       DZZDX   R.A.  O   Derivative in X-direction (W-E).
    !       DZZDY   R.A.  O   Derivative in Y-direction (S-N).
    !       IXP: IX plus 1 (with branch cut incorporated)
    !       IYP, IXM, IYM: ditto
    !       IXPS: value to use for IXP if IXPS is not masked.
    !             (use IX if masked)
    !       IYPS, IXMS, IYMS : ditto
    !       IXTRPL : in case of needing IY+1 for IY=NY, IX needs to be
    !                modified (tripole grid only)
    !       IXTRPLS : value to use for IXTRPL if IXTRPLS is not masked
    !             (use IX if masked)
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     This routine replaces the functionality of W3DDXY and W3DCXY.
    !     NB: subroutine "W3CGDM" has a similar purpose.
    !     Output arrays are always initialized to zero.
    !
    !  8. Structure :
    !
    !     ----------------------------------------
    !      1.  Preparations
    !        a Initialize arrays
    !        b Set constants
    !      2.  Derivatives in X-direction (W-E).
    !      3.  Derivatives in Y-direction (S-N).
    !     ----------------------------------------
    !
    !  9. Switches :
    !
    !       !/S   Enable subroutine tracing.
    !       !/T   Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsta, mapfs, mapfs, &
         dpdx, dpdy, dqdx, dqdy, flagll, iclose,          &
         iclose_none, iclose_smpl, iclose_trpl
    USE w3odatmd, ONLY: ndse, iaproc, naperr, naproc
    USE w3servmd, ONLY: extcde
#ifdef W3_T
    USE w3arrymd, ONLY : prtblk
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    REAL, INTENT(IN)        :: ZZ(NSEA)
    CHARACTER, INTENT(IN)   :: ZUNIT*(*)
    REAL, INTENT(OUT)       :: DZZDX(NY,NX), DZZDY(NY,NX)
    INTEGER                 :: ISEA, IX, IY, IXP, IXM, IYP, IYM
#ifdef W3_T
    INTEGER                 :: ISX, ISY, MAPOUT(NX,NY)
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T
    INTEGER, SAVE           :: NXS = 49
#endif
    REAL                    :: DFAC , STX, STY
    INTEGER                 :: IXPS,IYPS,IXMS,IYMS,IXTRPL,IXTRPLS
    INTEGER                 :: IXSTART,IXEND
#ifdef W3_T
    REAL                    :: XOUT(NX,NY)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3DZXY')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
 
    ! 1.a Initialize arrays
    !
    dzzdx = 0.
    dzzdy = 0.
    !
    ! 1.b Set constants
    !
 
    IF ( flagll ) THEN
      dfac   = 1. / ( dera * radius )
    ELSE
      dfac   = 1.
    END IF
 
    !
    ! 2.  Derivatives in X-direction (W-E) and Y-direction (S-N) ----- *
    !
 
    ! 2a. All points previously done in 2a,2b,2c of v4.18 done in 2a now:
    IF ( iclose.EQ.iclose_none ) THEN
      ixstart=2
      ixend=nx-1
    ELSE
      ixstart=1
      ixend=nx
    ENDIF
 
    DO iy=2, ny-1
      DO ix=ixstart,ixend
        IF ( mapsta(iy,ix) .NE. 0 ) THEN
          stx    = 0.5
          IF (ix.EQ.nx)THEN
            ixps=1
          ELSE
            ixps=ix+1
          ENDIF
 
          IF (mapsta(iy,ixps).EQ.0) THEN
            ixp    = ix
            stx    = 1.0
          ELSE
            ixp    = ixps
          END IF
 
          IF(ix.EQ.1)THEN
            ixms=nx
          ELSE
            ixms=ix-1
          ENDIF
 
          IF (mapsta(iy,ixms).EQ.0) THEN
            ixm    = ix
            stx    = 1.0
          ELSE
            ixm    = ixms
          END IF
          sty    = 0.5
          iyps=iy+1
          IF (mapsta(iyps,ix).EQ.0) THEN
            iyp    = iy
            sty    = 1.0
          ELSE
            iyp    = iyps
          END IF
          iyms=iy-1
          IF (mapsta(iyms,ix).EQ.0) THEN
            iym    = iy
            sty    = 1.0
          ELSE
            iym    = iyms
          END IF
          dzzdx(iy,ix) = (zz(mapfs(iy ,ixp))-zz(mapfs(iy ,ixm))) * stx * dpdx(iy,ix) &
               + (zz(mapfs(iyp,ix ))-zz(mapfs(iym,ix ))) * sty * dqdx(iy,ix)
          dzzdy(iy,ix) = (zz(mapfs(iy ,ixp))-zz(mapfs(iy ,ixm))) * stx * dpdy(iy,ix) &
               + (zz(mapfs(iyp,ix ))-zz(mapfs(iym,ix ))) * sty * dqdy(iy,ix)
          dzzdx(iy,ix) = dzzdx(iy,ix) * dfac
          dzzdy(iy,ix) = dzzdy(iy,ix) * dfac
        END IF
      END DO
    END DO
 
    ! 2b. column IY=NY for tripole case
    ! This is more complex, since for these two points: (IYP,IX) (IYM,IX),
    !     not only is the first index different (IYP.NE.IYM), but also the
    !     second index is different (IX.NE.IX)!
    IF ( iclose.EQ.iclose_trpl ) THEN
 
      iy=ny
      DO ix=1, nx
        IF ( mapsta(iy,ix) .NE. 0 ) THEN
 
          stx    = 0.5
 
          IF (ix.EQ.nx)THEN
            ixps=1
          ELSE
            ixps=ix+1
          ENDIF
          IF (mapsta(iy,ixps).EQ.0) THEN
            ixp    = ix
            stx    = 1.0
          ELSE
            ixp    = ixps
          END IF
 
          IF(ix.EQ.1)THEN
            ixms=nx
          ELSE
            ixms=ix-1
          ENDIF
          IF (mapsta(iy,ixms).EQ.0) THEN
            ixm    = ix
            stx    = 1.0
          ELSE
            ixm    = ixms
          END IF
 
          sty    = 0.5
 
          !..............next point: j+1: tripole: j==>j+1==>j and i==>ni-i+1
          !..............i.e. target point is MAPFS(IY,(NX-IX+1))
          ixtrpls=nx-ix+1
          IF (mapsta(iy,ixtrpls).EQ.0) THEN
            ixtrpl = ix
            sty    = 1.0
          ELSE
            ixtrpl=ixtrpls
          END IF
 
          iyms=iy-1
          IF (mapsta(iyms,ix).EQ.0) THEN
            iym    = iy
            sty    = 1.0
          ELSE
            iym    = iyms
          END IF
 
          ! tripole grid: (IYP,IX) is replaced with (IY,IXTRPL)
          dzzdx(iy,ix) = (zz(mapfs(iy ,ixp))-zz(mapfs(iy ,ixm))) * stx * dpdx(iy,ix) &
               + (zz(mapfs(iy,ixtrpl))-zz(mapfs(iym,ix ))) * sty * dqdx(iy,ix)
          dzzdy(iy,ix) = (zz(mapfs(iy ,ixp))-zz(mapfs(iy ,ixm))) * stx * dpdy(iy,ix) &
               + (zz(mapfs(iy,ixtrpl))-zz(mapfs(iym,ix ))) * sty * dqdy(iy,ix)
          dzzdx(iy,ix) = dzzdx(iy,ix) * dfac
          dzzdy(iy,ix) = dzzdy(iy,ix) * dfac
        END IF
      END DO
 
    END IF ! IF ( ICLOSE.EQ.ICLOSE_TRPL ) THEN
 
    !
    ! 3.  Test output of fields ------------------------------------------ *
    !
#ifdef W3_T
    WRITE (ndst,9010)
    isx   = 1 + nx/nxs
    isy   = 1 + ny/nxs
    DO iy=1, ny
      DO ix=1, nx
        mapout(ix,iy) = mapsta(iy,ix)
        IF ( mapfs(iy,ix) .NE. 0 )                               &
             xout(ix,iy) = zz(mapfs(iy,ix))
      END DO
    END DO
    CALL prtblk (ndst, nx, ny, nx, xout, mapout, 0, 0.,          &
         1, nx, isx, 1, ny, isy, 'ZZ', zunit)
    DO iy=1, ny
      DO ix=1, nx
        xout(ix,iy) = dzzdx(iy,ix)
      END DO
    END DO
    CALL prtblk (ndst, nx, ny, nx, xout, mapout, 0, 0.,          &
         1, nx, isx, 1, ny, isy, 'DZZDX',trim(zunit)//'/m')
    DO iy=1, ny
      DO ix=1, nx
        xout(ix,iy) = dzzdy(iy,ix)
      END DO
    END DO
    CALL prtblk (ndst, nx, ny, nx, xout, mapout, 0, 0.,          &
         1, nx, isx, 1, ny, isy, 'DZZDY',trim(zunit)//'/m')
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3DZXY : DX0I, DY0I : ',2e12.5)
9010 FORMAT (' TEST W3DZXY : FIELDS ')
#endif
    !/
    !/ End of W3DZXY ----------------------------------------------------- /
    !/
  END SUBROUTINE w3dzxy
  !/ ------------------------------------------------------------------- /
  !/ End of module W3UPDTMD -------------------------------------------- /
  !/
END MODULE w3updtmd