w3pro1md.ftn

#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      MODULE W3PRO1MD
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         05-Jun-2018 |
!/                  +-----------------------------------+
!/
!/    04-Feb-2000 : Origination                         ( version 2.00 )
!/    28-Mar-2001 : Partial time step bug fix (proper   ( version 2.10 )
!/                  ingest of boundaries).
!/    02-Apr-2001 : Sub-grid obstructions.              ( version 2.10 )
!/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
!/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/    07-Sep-2005 : Improved XY boundary conditions.    ( version 3.08 )
!/    10-Jan-2007 : Clean-up FACVX/Y compute.           ( version 3.10 )
!/    05-Mar-2008 : Added NEC sxf90 compiler directives
!/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
!/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
!/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
!/                  (W. E. Rogers & T. J. Campbell, NRL)
!/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
!/                  specify index closure for a grid.   ( version 3.14 )
!/                  (T. J. Campbell, NRL)
!/    29-May-2014 : Adding OMPH switch.                 ( version 5.02 )
!/    08-May-2014 : Implement tripolar grid for first order propagation
!/                  scheme                              ( version 5.03 )
!/                  (W. E. Rogers, NRL)
!/    05-Jun-2018 : Add DEBUG                           ( version 6.04 )
!/
!/    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 routines for first order propagation scheme in single
!     module.
!
!  2. Variables and types :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  3. Subroutines and functions :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!      W3MAP1    Subr. Public   Set up auxiliary maps.
!      W3XYP1    Subr. Public   First order spatial propagation.
!      W3KTP1    Subr. Public   First order spectral propagation.
!     ----------------------------------------------------------------
!
!  4. Subroutines and functions used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      DSEC21    Func. W3TIMEMD Time difference.
!      STRACE    Subr. W3SERVMD Subroutine tracing.
!     ----------------------------------------------------------------
!
!  5. Remarks :
!
!  6. Switches :
!
!       !/C90   Cray FORTRAN 90 compiler directives.
!       !/NEC   NEC SXF90 compiler directives.
!
!       !/S     Enable subroutine tracing.
!       !/Tn    Enable test output.
!
!  7. Source code :
!
!/ ------------------------------------------------------------------- /
      CONTAINS
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3MAP1 ( MAPSTA )
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         06-Dec-2010 |
!/                  +-----------------------------------+
!/
!/    19-Dec-1996 : Final FORTRAN 77                    ( version 1.18 )
!/    14-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
!/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/    10-Jan-2007 : Clean-up FACVX/Y compute.           ( version 3.10 )
!/    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 :
!
!     Generate 'map' arrays for the first order upstream scheme.
!
!  2. Method :
!
!     See section 3.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       MAPSTA  I.A.   I   Status map.
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!     See module documentation.
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ---------------------------------------------------------------- 
!      W3WAVE    Subr. W3WAVEMD Wave model routine.
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!  7. Remarks :
!
!  8. Structure :
!
!     ------------------------------------------------------
!      1.   Initialize arrays.
!      2.   Fill arrays.
!      3.   Invert arrays.
!     ------------------------------------------------------
!
!  9. Switches :
!
!     !/S   Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3GDATMD, ONLY: NTH, NSPEC, NX, NY, ICLOSE,                 &
                          ICLOSE_NONE, ICLOSE_SMPL, ICLOSE_TRPL
      USE W3ADATMD, ONLY: IS0, IS2, FACVX, FACVY
      USE W3ODATMD, ONLY: NDSE, IAPROC, NAPERR
      USE W3SERVMD, ONLY: EXTCDE
!/S      USE W3SERVMD, ONLY: STRACE
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: MAPSTA(NY*NX)
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: IX, IY, IXY, ISP, IXNEXT
!/S      INTEGER, SAVE           :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'W3MAP1')
!

! 1.  Initialize x-y arrays ------------------------------------------ *
!
      FACVX = 0.
      FACVY = 0.
!
! 2.  Fill x-y arrays ------------------------------------------------ *
!
!.....FACVY
      DO IX=1, NX
        DO IY=1, NY-1
          IXY    = IY +(IX-1)*NY
          IF ( MAPSTA( IXY ) .NE. 0 ) FACVY(IXY) = FACVY(IXY) + 1.
!.........next point : j+1 : increment IXY by 1
          IF ( MAPSTA(IXY+1) .NE. 0 ) FACVY(IXY) = FACVY(IXY) + 1.
          END DO
        END DO
!
!.....FACVY for IY=NY
      IF ( ICLOSE.EQ.ICLOSE_TRPL ) THEN
         IY=NY
         DO IX=1, NX
            IXY    = IY +(IX-1)*NY
            IF ( MAPSTA( IXY ) .NE. 0 ) FACVY(IXY) = FACVY(IXY) + 1.
!...........next point: j+1: tripole: j==>j+1==>j and i==>ni-i+1
            IXNEXT=NX-IX+1
            IXY    = IY +(IXNEXT-1)*NY
            IF ( MAPSTA( IXY ) .NE. 0 ) FACVY(IXY) = FACVY(IXY) + 1.
         END DO
!BGR: Adding the following lines to compute FACVX over all
!      IX for IY=NY (this allows along-seam propagation).
!      Located here since already inside "TRPL" if-block.
!{
         DO IX=1, NX-1
            IXY    = IY +(IX-1)*NY
            IF ( MAPSTA( IXY ) .NE. 0 ) FACVX(IXY) = FACVX(IXY) + 1.
            IF ( MAPSTA(IXY+NY) .NE. 0 ) FACVX(IXY) = FACVX(IXY) + 1.
         END DO
!}
      END IF
!
!.....FACVX
      DO IX=1, NX-1
        DO IY=2, NY-1
          IXY    = IY +(IX-1)*NY
          IF ( MAPSTA( IXY  ) .NE. 0 ) FACVX(IXY) = FACVX(IXY) + 1.
!.........next point : i+1 : increment IXY by NY
          IF ( MAPSTA(IXY+NY) .NE. 0 ) FACVX(IXY) = FACVX(IXY) + 1.
          END DO
        END DO
!
!.....FACVX for IX=NX
      IF ( ICLOSE.NE.ICLOSE_NONE ) THEN
          DO IY=2, NY-1
            IXY    = IY +(NX-1)*NY
            IF ( MAPSTA(IXY) .NE. 0 ) FACVX(IXY) = FACVX(IXY) + 1.
!...........next point : i+1 : increment IXY by NY
!...........IXY+NY=IY+(IX-1)*NY+NY = IY+IX*NY = IY+NX*NY ==> wrap to IY
            IF ( MAPSTA(IY ) .NE. 0 ) FACVX(IXY) = FACVX(IXY) + 1.
            END DO
        END IF
!
! 3.  Invert x-y arrays ---------------------------------------------- *
!
      DO IXY=1, NX*NY
        IF ( FACVX(IXY) .NE. 0. ) FACVX(IXY) = 1. / FACVX(IXY)
        IF ( FACVY(IXY) .NE. 0. ) FACVY(IXY) = 1. / FACVY(IXY)
        END DO
!
! 4.  Fill theta arrays ---------------------------------------------- *
!
      DO ISP=1, NSPEC
        IS2  (ISP) = ISP + 1
        IS0  (ISP) = ISP - 1
        END DO
!
      DO ISP=NTH, NSPEC, NTH
        IS2(ISP) = IS2(ISP) - NTH
        END DO
!
      DO ISP=1, NSPEC, NTH
        IS0(ISP) = IS0(ISP) + NTH
        END DO
!
      RETURN
!/
!/ End of W3MAP1 ----------------------------------------------------- /
!/
      END SUBROUTINE W3MAP1
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3XYP1 ( ISP, DTG, MAPSTA, FIELD, VGX, VGY )
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         29-May-2014 |
!/                  +-----------------------------------+
!/
!/    07-Jul-1998 : Final FORTRAN 77                    ( version 1.18 )
!/    14-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
!/    28-Mar-2001 : Partial time step bug fix.          ( version 2.10 )
!/    02-Apr-2001 : Sub-grid obstructions.              ( version 2.10 )
!/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
!/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/    07-Sep-2005 : Improved XY boundary conditions.    ( version 3.08 )
!/    05-Mar-2008 : Added NEC sxf90 compiler directives
!/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
!/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
!/                  (W. E. Rogers & T. J. Campbell, NRL)
!/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
!/                  specify index closure for a grid.   ( version 3.14 )
!/                  (T. J. Campbell, NRL)
!/    29-May-2014 : Adding OMPH switch.                 ( version 5.02 )
!/
!  1. Purpose :
!
!     Propagation in physical space for a given spectral component.
!
!  2. Method :
!
!     First order scheme with flux formulation.
!     Curvilinear grid implementation: Fluxes are computed in index space
!       and then transformed back into physical space.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       ISP     Int.   I   Number of spectral bin (IK-1)*NTH+ITH
!       DTG     Real   I   Total time step.
!       MAPSTA  I.A.   I   Grid point status map.
!       FIELD   R.A.  I/O  Wave action spectral densities on full
!                          grid.
!       VGX/Y   Real   I   Speed of grid.
!     ----------------------------------------------------------------
!
!     Local variables.
!     ----------------------------------------------------------------
!       NTLOC   Int.  Number of local steps.
!       DTLOC   Real  Local propagation time step.
!       VCX     R.A.  Propagation velocities in index space.
!       VCY     R.A.
!       CXTOT   R.A.  Propagation velocities in physical space.
!       CYTOT   R.A.
!       VFLX    R.A.  Discrete fluxes between grid points in index space.
!       VFLY    R.A.
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!     See module documentation.
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ---------------------------------------------------------------- 
!      W3WAVE    Subr. W3WAVEMD Wave model routine.
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!       None.
!
!  7. Remarks :
!
!     - The local work arrays are initialized on the first entry to
!       the routine.
!     - Curvilinear grid implementation. Variables FACX, FACY, CCOS, CSIN,
!       CCURX, CCURY are not needed and have been removed.  FACX is accounted
!       for as approriate in this subroutine.  FACX is also accounted for in
!       the case of .NOT.FLCX.  Since FACX is removed, there is now a check for
!       .NOT.FLCX in this subroutine.  In CFL calcs dx and dy are omitted,
!       since dx=dy=1 in index space.  Curvilinear grid derivatives
!       (DPDY, DQDX, etc.) and metric (GSQRT) are brought in via W3GDATMD.
!     - Standard VCB calculation for Y is:
!               VCB       = FACVY(IXY) * ( VCY2D(IY,IX) + VCY2D(IY+1,IX) )
!       This is to calculate the flux VCY(IY+0.5). For the tripole grid,
!       we cannot do it this way, since the sign of VCY flips as we jump
!       over the seam. If we were to do it this way, VCY(IY) and VCY(IY+1)
!       are two numbers of similar magnitude and opposite sign, so the 
!       average of the two gives something close to zero, so energy does 
!       not leave via VCY(IY+0.5). One alternative is:
!               VCB       = VCY2D(IY,IX) 
!       Another alternative is :
!               VCB       = FACVY(IXY) * ( VCY2D(IY,IX) - VCY2D(IY+1,IX) )
!       Both appear to give correct results for ww3_tp2.13. We use the 
!       second alternative.
!
!  8. Structure :
!
!     ---------------------------------------
!       1.  Preparations
!         a Set constants
!         b Initialize arrays
!       2.  Calculate local discrete fluxes
!       3.  Calculate propagation fluxes
!       4.  Propagate
!       5.  Update boundary conditions
!     ---------------------------------------
!
!  9. Switches :
!
!     !/NEC Enable NEC SXF90 compiler directives.
!
!     !/S   Enable subroutine tracing.
!
!     !/OMPH  Hybrid OpenMP directives.
!
!     !/T   Enable general test output.
!     !/T1  Test output local fluxes (V)FX-YL.
!     !/T2  Test output propagation fluxes (V)FLX-Y.
!     !/T3  Test output propagation.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS
!
      USE W3TIMEMD, ONLY: DSEC21
!
      USE W3GDATMD, ONLY: NK, NTH, SIG, ECOS, ESIN, NX, NY, NSEA,     &
                          MAPSF, DTCFL, ICLOSE, CLATS, FLCX, FLCY,    &
                          ICLOSE_NONE, ICLOSE_SMPL, ICLOSE_TRPL,      &
                          FLAGLL, DPDX, DPDY, DQDX, DQDY, GSQRT
      USE W3WDATMD, ONLY: TIME
      USE W3ADATMD, ONLY: CG, CX, CY, ATRNX, ATRNY, FACVX, FACVY
      USE W3IDATMD, ONLY: FLCUR
      USE W3ODATMD, ONLY: NDST, FLBPI, NBI, TBPI0, TBPIN, ISBPI,      &
                          BBPI0, BBPIN, NDSE, IAPROC, NAPERR
      USE W3SERVMD, ONLY: EXTCDE
!/S      USE W3SERVMD, ONLY: STRACE
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: ISP, MAPSTA(NY*NX)
      REAL, INTENT(IN)        :: DTG, VGX, VGY
      REAL, INTENT(INOUT)     :: FIELD(1-NY:NY*(NX+2))
!/
!/ ------------------------------------------------------------------ /
!/ Local parameters
!/
      INTEGER                 :: IK, ITH, NTLOC, ITLOC, ISEA, IXY,    &
                                 IY0, IX, IY, JXN, JXP, JYN, JYP,     &
                                 IBI, NYMAX
!/T3      INTEGER                 ::  IXF, IYF
!/S      INTEGER, SAVE           :: IENT = 0
      REAL                    :: CG0, CGL, CGA, CC, CGN
      REAL                    :: DTLOC,DTRAD, VCB
      REAL                    :: RD1, RD2
      REAL                    :: CP, CQ
!/T3      REAL                    :: AOLD
!/
!/ Automatic work arrays
!/
      REAL                    :: CXTOT2D(NY,NX)
      REAL                    :: CYTOT2D(NY,NX)
      REAL                    :: FLD2D(NY+1,NX+1)
      REAL                    :: VCX2D(NY,NX+1)
      REAL                    :: VCY2D(NY+1,NX)
      REAL                    :: VFLX2D(1:NY,0:NX)
      REAL                    :: VFLY2D(NY,NX)

!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'W3XYP1')
!
! 1.  Preparations --------------------------------------------------- *

! 1.a Set constants
!
      ITH    = 1 + MOD(ISP-1,NTH)
      IK     = 1 + (ISP-1)/NTH
!
      CG0    = 0.575 * GRAV / SIG(1)
      CGL    = 0.575 * GRAV / SIG(IK)
!
      IF ( FLCUR ) THEN
          CGA    = SQRT(MAXVAL((CGL*ECOS(ITH)+CX(1:NSEA))**2          &
                              +(CGL*ESIN(ITH)+CY(1:NSEA))**2))
          CC     = SQRT(MAXVAL(CX(1:NSEA)**2+CY(1:NSEA)**2))
!/MGP          CGA    = SQRT(MAXVAL((CGL*ECOS(ITH)+CX(1:NSEA)-VGX)**2      &
!/MGP                              +(CGL*ESIN(ITH)+CY(1:NSEA)-VGY)**2))
!/MGP          CC     = SQRT(MAXVAL((CX(1:NSEA)-VGX)**2+(CY(1:NSEA)-VGY)**2))
        ELSE
          CGA    = CGL
!/MGP          CGA    = SQRT((CGL*ECOS(ITH)-VGX)**2+(CGL*ESIN(ITH)-VGY)**2)
          CC     = 0.
        END IF
!
      CGN    = 0.9999 * MAX ( CGA, CC, 0.001*CG0 )
!
      NTLOC  = 1 + INT(DTG/(DTCFL*CG0/CGN))
      DTLOC  = DTG / REAL(NTLOC)
      DTRAD  = DTLOC
      IF ( FLAGLL ) DTRAD=DTRAD/(DERA*RADIUS)

!
!/T      WRITE (NDST,9000) NTLOC
!/T      WRITE (NDST,9001) ISP, ITH, IK
!
! ====================== Loop partial ================================ *
!
      DO ITLOC=1, NTLOC
!
! 1.b Initialize arrays
!
!/T1      WRITE (NDST,9010) ITLOC
!
        VCX2D = 0.
        VCY2D = 0.
        CXTOT2D  = 0.
        CYTOT2D  = 0.
        FLD2D  = 0.
        VFLX2D  = 0.
        VFLY2D  = 0.
!
! 2.  Calculate field and velocities --------------------------------- *
!
!     FIELD = A / CG * CLATS
!     VCX   = COS*CG / CLATS
!     VCY   = SIN*CG
!
!/T1      WRITE (NDST,9020)
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IXY, IX, IY)
!
        DO ISEA=1, NSEA
          IX     = MAPSF(ISEA,1)
          IY     = MAPSF(ISEA,2)
          IXY    = MAPSF(ISEA,3)

          FLD2D(IY,IX) = FIELD(IXY) / CG(IK,ISEA) * CLATS(ISEA)

          CXTOT2D(IY,IX) = ECOS(ITH) * CG(IK,ISEA) / CLATS(ISEA)
          CYTOT2D(IY,IX) = ESIN(ITH) * CG(IK,ISEA)
!/MGP          CXTOT2D(IY,IX) = CXTOT2D(IY,IX) - VGX/CLATS(ISEA)
!/MGP          CYTOT2D(IY,IX) = CYTOT2D(IY,IX) - VGY

!/T1          WRITE (NDST,9021) ISEA, IXY, FLD2D(IY,IX), &
!/T1                            CXTOT2D(IY,IX), CYTOT2D(IY,IX)
        END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
        IF ( FLCUR ) THEN
            DO ISEA=1, NSEA
              IX     = MAPSF(ISEA,1)
              IY     = MAPSF(ISEA,2)

              CXTOT2D(IY,IX) = CXTOT2D(IY,IX) + CX(ISEA)/CLATS(ISEA)
              CYTOT2D(IY,IX) = CYTOT2D(IY,IX) + CY(ISEA)

            END DO
        END IF
          
        IF ( FLCX ) THEN
           DO ISEA=1, NSEA
              IX     = MAPSF(ISEA,1)
              IY     = MAPSF(ISEA,2)
              CP=CXTOT2D(IY,IX)*DPDX(IY,IX)+CYTOT2D(IY,IX)*DPDY(IY,IX)
              VCX2D(IY,IX) = CP*DTRAD
           END DO
        ELSE
           VCX2D=0.0
        ENDIF

        IF ( FLCY ) THEN
           DO ISEA=1, NSEA
              IX     = MAPSF(ISEA,1)
              IY     = MAPSF(ISEA,2)
              CQ=CXTOT2D(IY,IX)*DQDX(IY,IX)+CYTOT2D(IY,IX)*DQDY(IY,IX)
              VCY2D(IY,IX) = CQ*DTRAD
           END DO
        ELSE
           VCY2D=0.0
        ENDIF

! Transform FIELD to index space, i.e. straightened space
! Bugfix: This is now done *before* adding the ghost row, so that ghost
!   row will be in index space (bug applied only to global, irregular
!   grids, so it did not apply to any test case that existed w/v4.18)
        FLD2D(1:NY,1:NX)=FLD2D(1:NY,1:NX)*GSQRT(1:NY,1:NX)

!
! Deal with longitude closure by duplicating one row *to the right*
!   in FIELD/FLD2D, VCX
        IF ( ICLOSE.NE.ICLOSE_NONE ) THEN
!/T1            WRITE (NDST,9024)
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
            DO IY=1, NY
              FLD2D(IY,NX+1)=FLD2D(IY,1)
              VCX2D(IY,NX+1)=VCX2D(IY,1)
!/T1              WRITE (NDST,9025) IY, FLD2D(IY,NX+1), VCX2D(IY,NX+1)
            END DO
        END IF

! Deal with tripole closure by duplicating one row *at the top*
!   in FIELD/FLD2D, VCY
        IF ( ICLOSE.EQ.ICLOSE_TRPL ) THEN
            DO IX=1,NX
!...........next point: j+1: tripole: j==>j+1==>j and i==>ni-i+1
              FLD2D(NY+1,IX)=FLD2D(NY,NX-IX+1)
              VCY2D(NY+1,IX)=VCY2D(NY,NX-IX+1)
            END DO
        END IF

!
! 3.  Calculate propagation fluxes ----------------------------------- *
!
        NYMAX=NY-1
        IF ( ICLOSE.EQ.ICLOSE_TRPL ) NYMAX=NY
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (IX, IY, IXY)
!
        DO IX=1, NX
           DO IY=1, NYMAX
              IXY    = IY +(IX-1)*NY
              VCB       = FACVX(IXY) * ( VCX2D(IY,IX) + VCX2D(IY,IX+1) )
              VFLX2D(IY,IX) = MAX ( VCB , 0. ) * FLD2D(IY,IX)          &
                   + MIN ( VCB , 0. ) * FLD2D(IY,IX+1)
           END DO
        END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
! Deal with longitude closure by duplicating one row *to the left*
!    in VFLX. Note that a similar action is not take for tripole grid,
!    since tripole seam is only: IY=NY communicating with other points
!    at IY=NY, not a case of IY=NY communicating with IY=1
        IF ( ICLOSE.NE.ICLOSE_NONE ) THEN
!/T2            WRITE (NDST,9032)
            DO IY=1, NY
              VFLX2D(IY,0) = VFLX2D(IY,NX)
!/T2              WRITE (NDST,9033) IY, VFLX2D(IY,0)
            END DO
        END IF
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (IX, IY, IXY, VCB)
!
        DO IX=1, NX
           DO IY=1, NY-1
              IXY    = IY +(IX-1)*NY
              VCB       = FACVY(IXY) * ( VCY2D(IY,IX) + VCY2D(IY+1,IX) )
              VFLY2D(IY,IX) = MAX ( VCB , 0. ) * FLD2D(IY,IX)          &
                            + MIN ( VCB , 0. ) * FLD2D(IY+1,IX)
           END DO
        END DO
!
!/OMPH/!$OMP END PARALLEL DO
!

! For tripole grid, include IY=NY in calculation. VCB is handled
!    differently. See notes in Section "7. Remarks" above.
        IF ( ICLOSE.EQ.ICLOSE_TRPL ) THEN
           IY=NY
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (IXY, VCB, IX, IY)
!
           DO IX=1, NX
              IXY    = IY +(IX-1)*NY
              VCB       = FACVY(IXY) * ( VCY2D(IY,IX) - VCY2D(IY+1,IX) )
              VFLY2D(IY,IX) = MAX ( VCB , 0. ) * FLD2D(IY,IX)          &
                            + MIN ( VCB , 0. ) * FLD2D(IY+1,IX)
           END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
        END IF

! 4.  Propagate ------------------------------------------------------ *
!
!/T3        WRITE (NDST,9040)
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IXY, JXN, JXP, JYN, JYP, IX, IY)
!
        DO ISEA=1, NSEA
!
          IX    = MAPSF(ISEA,1)
          IY    = MAPSF(ISEA,2)
          IXY   = MAPSF(ISEA,3)

!/T3          AOLD   = FLD2D(IY,IX) * CG(IK,ISEA) / CLATS(ISEA)
!
          IF (MAPSTA(IXY).EQ.1) THEN
!
             IF ( VFLX2D(IY,IX-1) .GT. 0. ) THEN
                JXN   = -1
             ELSE
                JXN   =  0
             END IF
             IF ( VFLX2D(IY,IX) .LT. 0. ) THEN
                JXP   =  1
             ELSE
                JXP   =  0
             END IF
             IF ( VFLY2D(IY-1,IX) .GT. 0. ) THEN
                JYN   = -1
             ELSE
                JYN   =  0
             END IF
             IF ( VFLY2D(IY,IX) .LT. 0. ) THEN
                JYP   =  1
             ELSE
                JYP   =  0
             END IF
!
             FLD2D(IY,IX) =   FLD2D(IY,IX)                            &
                                 + ATRNX(IXY,JXN) * VFLX2D(IY,IX-1)   &
                                 - ATRNX(IXY,JXP) * VFLX2D(IY,IX)     &
                                 + ATRNY(IXY,JYN) * VFLY2D(IY-1,IX)   &
                                 - ATRNY(IXY,JYP) * VFLY2D(IY,IX)   

!/T3              WRITE (NDST,9041) ISEA, IXY, IXY-NY, IXY-1,         &
!/T3                 VFLX2D(IY,IX-1), VFLX2D(IY,IX), VFLY2D(IY-1,IX), &
!/T3                 VFLY2D(IY,IX) , CG(IK,ISEA)/CLATS(ISEA),AOLD,    &
!/T3                 FLD2D(IY,IX)
!
!
!/T3              WRITE (NDST,9042) ISEA, MAPSTA(IXY), AOLD,FLD2D(IY,IX)
!
          END IF !  IF (MAPSTA(IXY).EQ.1) THEN

          FLD2D(IY,IX) = CG(IK,ISEA) / CLATS(ISEA) * FLD2D(IY,IX)

        END DO !  DO ISEA=1, NSEA
!
!/OMPH/!$OMP END PARALLEL DO
!

! Transform FIELD back to physical space, i.e. may be curvilinear
        FLD2D(1:NY,1:NX)=FLD2D(1:NY,1:NX)/GSQRT(1:NY,1:NX)
!
! 5.  Update boundary conditions ------------------------------------- *
!
        IF ( FLBPI ) THEN
           RD1    = DSEC21 ( TBPI0, TIME ) - DTG *                   &
                REAL(NTLOC-ITLOC)/REAL(NTLOC)
           RD2    = DSEC21 ( TBPI0, TBPIN )
           IF ( RD2 .GT. 0.001 ) THEN
              RD2    = MIN(1.,MAX(0.,RD1/RD2))
              RD1    = 1. - RD2
           ELSE
              RD1    = 0.
              RD2    = 1.
           END IF
           DO IBI=1, NBI
              IX    = MAPSF(ISBPI(IBI),1)
              IY    = MAPSF(ISBPI(IBI),2)
              FLD2D(IY,IX) = RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI)
           END DO
        END IF
!
! 6.  Put back in 1d shape ------------------------------------------- *
!
        DO ISEA=1, NSEA
           IX     = MAPSF(ISEA,1)
           IY     = MAPSF(ISEA,2)
           IXY    = MAPSF(ISEA,3)
           FIELD(IXY) = FLD2D(IY,IX)
        END DO
!
! ... End of partial time step loop
!
      END DO !   DO ITLOC=1, NTLOC
!
      RETURN
!
! Formats
!
!/T 9000 FORMAT (' TEST W3XYP1 : NTLOC :',I4)
!/T 9001 FORMAT (' TEST W3XYP1 : ISP, ITH, IK :',I8,2I4)
!
!/T1 9010 FORMAT (' TEST W3XYP1 : INIT. VFX-YL, ITLOC =',I3)
!
!/T1 9020 FORMAT (' TEST W3XYP1 : ISEA, IXY, FIELD, VCX, VCY')
!/T1 9021 FORMAT ('           ',2I8,3E12.4)
!/T1 9024 FORMAT (' TEST W3XYP1 : GLOBAL CLOSURE: IY, FIELD, VCX ')
!/T1 9025 FORMAT ('               ',I4,2E12.4)
!
!/T2 9032 FORMAT (' TEST W3XYP1 : CLOSE. : IY, VFLX')
!/T2 9033 FORMAT ('            ',I4,E12.4)
!
!/T3 9040 FORMAT (' TEST W3XYP1 : PROPAGATION '/                      &
!/T3              '      ISEA, IXY(3), , FLX(2), FLY(2), FAC, A(2)')
!/T3 9041 FORMAT (2X,4I5,1X,4E10.3,1X,E10.3,1X,2E10.3)
!/T3 9042 FORMAT (2X,I5,'( MAP = ',I2,' )',56X,2E10.3)
!/
!/ End of W3XYP1 ----------------------------------------------------- /
!/
      END SUBROUTINE W3XYP1
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3KTP1 ( ISEA, FACTH, FACK, CTHG0, CG, WN, DEPTH,    &
                          DDDX, DDDY, CX, CY, DCXDX, DCXDY, DCYDX,    &
                          DCYDY, DCDX, DCDY, VA )
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         20-Dec-2004 |
!/                  +-----------------------------------+
!/
!/    29-Aug-1997 : Final FORTRAN 77                    ( version 1.18 )
!/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
!/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
!/
!  1. Purpose :
!
!     Propagation in spectral space.
!
!  2. Method :
!
!     First order scheme.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       ISEA    Int.   I   Number of sea point.
!       FACTH/K Real   I   Factor in propagation velocity.
!       CTHG0   Real   I   Factor in great circle refracftion term.
!       CG      R.A.   I   Local group velocities.
!       WN      R.A.   I   Local wavenumbers.
!       DEPTH   Real   I   Depth.
!       DDDx    Real   I   Depth gradients.
!       CX/Y    Real   I   Current components.
!       DCxDx   Real   I   Current gradients.
!       DCDX-Y  Real   I   Phase speed gradients.
!       VA      R.A.  I/O  Spectrum.
!     ----------------------------------------------------------------
!
!     Local variables.
!     ----------------------------------------------------------------
!       DSDD    R.A.  Partial derivative of sigma for depth.
!       FRK, FRG, FKC
!               R.A.  Partial velocity terms.
!       DWNI    R.A.  Inverse band width.
!       CTH-WN  R.A.  Propagation velocities of local fluxes.
!       FLTH-WN R.A.  Propagation fluxes.
!       AA      R.A.  Extracted spectrum
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!     See module documentation.
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ---------------------------------------------------------------- 
!      W3WAVE    Subr. W3WAVEMD Wave model routine.
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!       None.
!
!  8. Structure :
!
!     -----------------------------------------------------------------
!       1.  Preparations
!         a Calculate DSDD
!         b Extract spectrum
!       2.  Refraction velocities
!         a Filter level depth reffraction.
!         b Depth refratcion velocity.
!         c Current refraction velocity.
!       3.  Wavenumber shift velocities
!         a Prepare directional arrays
!         b Calcuate velocity.
!       4.  Refraction
!         a Discrete fluxes.
!         b Propagation fluxes.
!         c Refraction.
!       5.  Wavenumber shifts.
!         a Discrete fluxes.
!         b Propagation fluxes.
!         c Refraction.
!     -----------------------------------------------------------------
!
!  9. Switches :
!
!     C/S   Enable subroutine tracing.
!     C/T   Enable general test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS
      USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, DSIP, ECOS, ESIN, ES2, &
                          ESC, EC2, FACHFA, MAPWN, FLCTH, FLCK, CTMAX
      USE W3ADATMD, ONLY: IS0, IS2
      USE W3IDATMD, ONLY: FLCUR
      USE W3ODATMD, ONLY: NDST
!/DEBUG      USE W3ODATMD, only : IAPROC
!/S      USE W3SERVMD, ONLY: STRACE
!/
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: ISEA
      REAL, INTENT(IN)        :: FACTH, FACK, CTHG0, CG(0:NK+1),      &
                                 WN(0:NK+1), DEPTH, DDDX, DDDY,       &
                                 CX, CY, DCXDX, DCXDY, DCYDX, DCYDY
      REAL, INTENT(IN)        :: DCDX(0:NK+1), DCDY(0:NK+1)
      REAL, INTENT(INOUT)     :: VA(NSPEC)
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: ITH, IK, ISP, ITH0
!/S      INTEGER, SAVE           :: IENT = 0
      REAL                    :: FDDMAX, FDG, DCYX, DCXXYY, DCXY,     &
                                 DCXX, DCXYYX, DCYY, FKD, FKD0, CTHB, &
                                 CWNB
      REAL                    :: VCTH(NSPEC), VCWN(1-NTH:NSPEC+NTH),  &
                                 VAA(1-NTH:NSPEC+NTH), VFLTH(NSPEC),  &
                                 VFLWN(1-NTH:NSPEC), DSDD(0:NK+1),    &
                                 FRK(NK), FRG(NK), FKC(NTH), DWNI(NK)
!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'W3KTP1')
!
!/T      WRITE (NDST,9000) FLCTH, FLCK, FACTH, FACK, CTMAX
!/T      WRITE (NDST,9001) ISEA, DEPTH, CX, CY,                       &
!/T             DDDX, DDDY, DCXDX, DCXDY, DCYDX, DCYDY
!
! 1.  Preparations --------------------------------------------------- *
! 1.a Array with partial derivative of sigma versus depth
!
      DO IK=0, NK+1
        IF ( DEPTH*WN(IK) .LT. 5. ) THEN
            DSDD(IK) = MAX ( 0. ,                                     &
                CG(IK)*WN(IK)-0.5*SIG(IK) ) / DEPTH
          ELSE
            DSDD(IK) = 0.
          END IF
        END DO
!
!/T      WRITE (NDST,9010)
!/T      DO IK=1, NK+1
!/T        WRITE (NDST,9011) IK, TPI/SIG(IK), TPI/WN(IK),             &
!/T                              CG(IK), DSDD(IK)
!/T        END DO
!
! 1.b Extract spectrum
!
      DO ISP=1, NSPEC
        VAA(ISP) = VA(ISP)
        END DO
!
! 2.  Refraction velocities ------------------------------------------ *
!
      IF ( FLCTH ) THEN
!
! 2.a Set slope filter for depth refraction
!
          FDDMAX = 0.
          FDG    = FACTH * CTHG0
!
          DO ITH=1, NTH
            FDDMAX = MAX ( FDDMAX , ABS (                             &
               ESIN(ITH)*DDDX - ECOS(ITH)*DDDY ) )
            END DO
!
          DO IK=1, NK
            FRK(IK) = FACTH * DSDD(IK) / WN(IK)
            FRK(IK) = FRK(IK) / MAX ( 1. , FRK(IK)*FDDMAX/CTMAX )
            FRG(IK) = FDG * CG(IK)
            END DO
!
! 2.b Depth refraction and great-circle propagation
!
          DO ISP=1, NSPEC
            VCTH(ISP) = FRG(MAPWN(ISP)) * ECOS(ISP)                   &
              + FRK(MAPWN(ISP)) * ( ESIN(ISP)*DDDX - ECOS(ISP)*DDDY )
            END DO
!/DEBUG      WRITE(740+IAPROC,*) 'pro1 FACTH=', FACTH
!/DEBUG      WRITE(740+IAPROC,*) 'pro1 CTHG0=', CTHG0
!/DEBUG      WRITE(740+IAPROC,*) 'pro1 FDG=', FDG
!/DEBUG      WRITE(740+IAPROC,*) 'pro1 FDDMAX=', FDDMAX
!/DEBUG      WRITE(740+IAPROC,*) 'pro1 sum(FRK)=', sum(FRK)
!/DEBUG      WRITE(740+IAPROC,*) 'pro1 sum(FRG)=', sum(FRG)
!/DEBUG      WRITE(740+IAPROC,*) 'pro1 sum(DSDD)=', sum(DSDD)
!/DEBUG      WRITE(740+IAPROC,*) 'ISEA=', ISEA, ' sum(VCTH)=', sum(VCTH)
!/DEBUG      FLUSH(740+IAPROC)

!
!/REFRX! 2.c @C/@x refraction and great-circle propagation
!/REFRX          VCTH = 0.
!/REFRX          FRK  = 0.
!/REFRX          FDDMAX = 0.
!
!/REFRX          DO ISP=1, NSPEC
!/REFRX            FDDMAX = MAX ( FDDMAX , ABS (                      &
!/REFRX               ESIN(ISP)*DCDX(MAPWN(ISP)) - ECOS(ISP)*DCDY(MAPWN(ISP)) ) )
!/REFRX            END DO
!
!/REFRX          DO IK=1, NK
!/REFRX            FRK(IK) = FACTH * CG(IK) * WN(IK) / SIG(IK)
!/REFRX            FRK(IK) = FRK(IK) / MAX ( 1. , FRK(IK)*FDDMAX/CTMAX )
!/REFRX            FRG(IK) = FDG * CG(IK)
!/REFRX          END DO
!/REFRX          DO ISP=1, NSPEC
!/REFRX            VCTH(ISP) = FRG(MAPWN(ISP)) * ECOS(ISP)            &
!/REFRX              + FRK(MAPWN(ISP)) * ( ESIN(ISP)*DCDX(MAPWN(ISP)) &
!/REFRX              - ECOS(ISP)*DCDY(MAPWN(ISP)) )
!/REFRX          END DO
!
! 2.d Current refraction
!
          IF ( FLCUR ) THEN
!
              DCYX   = FACTH *   DCYDX
              DCXXYY = FACTH * ( DCXDX - DCYDY )
              DCXY   = FACTH *   DCXDY
!
              DO ISP=1, NSPEC
                VCTH(ISP) = VCTH(ISP) + ES2(ISP)*DCYX                 &
                  + ESC(ISP)*DCXXYY - EC2(ISP)*DCXY
                END DO
!
            END IF
!
        END IF
!
! 3.  Wavenumber shift velocities ------------------------------------ *
!
      IF ( FLCK ) THEN
!
          DCXX   =  - FACK *   DCXDX
          DCXYYX =  - FACK * ( DCXDY + DCYDX )
          DCYY   =  - FACK *   DCYDY
          FKD    =    FACK * ( CX*DDDX + CY*DDDY )
!
          DO ITH=1, NTH
            FKC(ITH) = EC2(ITH)*DCXX +                                &
                       ESC(ITH)*DCXYYX + ES2(ITH)*DCYY
            END DO
!
          ISP    = -NTH
          DO IK=0, NK+1
            FKD0   = FKD / CG(IK) * DSDD(IK)
            DO ITH=1, NTH
              ISP    = ISP + 1
              VCWN(ISP) = FKD0 + WN(IK)*FKC(ITH)
              END DO
            END DO
!
          ITH0   = NSPEC - NTH
          DO ITH=1, NTH
            VAA(ITH+NSPEC) = FACHFA * VAA(ITH+ITH0)
            VAA(ITH- NTH ) = 0.
            END DO
!
          DO IK=1, NK
            DWNI(IK) = CG(IK) / DSIP(IK)
            END DO
!
        END IF
!
! 4.  Refraction ----------------------------------------------------- *
!
      IF ( FLCTH ) THEN
!
! 4.a Boundary velocity and fluxes
!
          DO ISP=1, NSPEC
            CTHB       = 0.5 * ( VCTH(ISP) + VCTH(IS2(ISP)) )
            VFLTH(ISP) = MAX ( CTHB , 0. ) * VAA(ISP)                 &
                       + MIN ( CTHB , 0. ) * VAA(IS2(ISP))
            END DO
!
! 4.b Propagation
!
          DO ISP=1, NSPEC
            VA(ISP) = VA(ISP) + VFLTH(IS0(ISP)) - VFLTH(ISP )
            END DO
!
        END IF
!
! 5.  Wavenumber shifts ---------------------------------------------- *
!
      IF ( FLCK ) THEN
!
! 5.a Boundary velocity and fluxes
!
          DO ISP=1-NTH, NSPEC
            CWNB       = 0.5 * ( VCWN(ISP) + VCWN(ISP+NTH) )
            VFLWN(ISP) = MAX ( CWNB , 0. ) * VAA(  ISP  )             &
                       + MIN ( CWNB , 0. ) * VAA(ISP+NTH)
            END DO
!
! 5.c Propagation
!
          DO ISP=1, NSPEC
            VA(ISP) = VA(ISP) + DWNI(MAPWN(ISP)) *                    &
                           ( VFLWN(ISP-NTH) - VFLWN(ISP) )
            END DO
!
        END IF
!
      RETURN
!
! Formats
!
!/T 9000 FORMAT (' TEST W3KTP1 : FLCTH-K, FACTH-K, CTMAX  :',         &
!/T                2L2,2E10.3,F7.3)
!/T 9001 FORMAT (' TEST W3KTP1 : LOCAL DATA :',I7,F7.1,2F6.2,1X,      &
!/T                6E10.3)
!/T 9010 FORMAT (' TEST W3KTP1 : IK, T, L, CG, DSDD : ')
!/T 9011 FORMAT ('              ',I3,F7.2,F7.1,F7.2,E11.3)
!/
!/ End of W3KTP1 ----------------------------------------------------- /
!/
      END SUBROUTINE W3KTP1
!/
!/ End of module W3PRO1MD -------------------------------------------- /
!/
      END MODULE W3PRO1MD