WW3/w3pro1md_8F90_source.html

#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 :

  !

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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    INTEGER, INTENT(IN)     :: MAPSTA(NY*NX)

    !/

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

    !/ Local parameters

    !/

    INTEGER                 :: IX, IY, IXY, ISP, IXNEXT

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3MAP1')

#endif

    !


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

    !

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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    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

#ifdef W3_T3

    INTEGER                 ::  IXF, IYF

#endif

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    REAL                    :: CG0, CGL, CGA, CC, CGN

    REAL                    :: DTLOC,DTRAD, VCB

    REAL                    :: RD1, RD2

    REAL                    :: CP, CQ

#ifdef W3_T3

    REAL                    :: AOLD

#endif

    !/

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


    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3XYP1')

#endif

    !

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

#ifdef W3_MGP

      cga    = sqrt(maxval((cgl*ecos(ith)+cx(1:nsea)-vgx)**2      &

           +(cgl*esin(ith)+cy(1:nsea)-vgy)**2))

      cc     = sqrt(maxval((cx(1:nsea)-vgx)**2+(cy(1:nsea)-vgy)**2))

#endif

    ELSE

      cga    = cgl

#ifdef W3_MGP

      cga    = sqrt((cgl*ecos(ith)-vgx)**2+(cgl*esin(ith)-vgy)**2)

#endif

      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)


    !

#ifdef W3_T

    WRITE (ndst,9000) ntloc

    WRITE (ndst,9001) isp, ith, ik

#endif

    !

    ! ====================== Loop partial ================================ *

    !

    DO itloc=1, ntloc

      !

      ! 1.b Initialize arrays

      !

#ifdef W3_T1

      WRITE (ndst,9010) itloc

#endif

      !

      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

      !

#ifdef W3_T1

      WRITE (ndst,9020)

#endif

      !

#ifdef W3_OMPH

      !$OMP PARALLEL DO PRIVATE (ISEA, IXY, IX, IY)

#endif

      !

      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)

#ifdef W3_MGP

        cxtot2d(iy,ix) = cxtot2d(iy,ix) - vgx/clats(isea)

        cytot2d(iy,ix) = cytot2d(iy,ix) - vgy

#endif


#ifdef W3_T1

        WRITE (ndst,9021) isea, ixy, fld2d(iy,ix), &

             cxtot2d(iy,ix), cytot2d(iy,ix)

#endif

      END DO

      !

#ifdef W3_OMPH

      !$OMP END PARALLEL DO

#endif

      !

      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

#ifdef W3_T1

        WRITE (ndst,9024)

#endif

        DO iy=1, ny

          fld2d(iy,nx+1)=fld2d(iy,1)

          vcx2d(iy,nx+1)=vcx2d(iy,1)

#ifdef W3_T1

          WRITE (ndst,9025) iy, fld2d(iy,nx+1), vcx2d(iy,nx+1)

#endif

        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

      !

#ifdef W3_OMPH

      !$OMP PARALLEL DO PRIVATE (IX, IY, IXY, VCB)

#endif

      !

      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

      !

#ifdef W3_OMPH

      !$OMP END PARALLEL DO

#endif

      !

      ! 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

#ifdef W3_T2

        WRITE (ndst,9032)

#endif

        DO iy=1, ny

          vflx2d(iy,0) = vflx2d(iy,nx)

#ifdef W3_T2

          WRITE (ndst,9033) iy, vflx2d(iy,0)

#endif

        END DO

      END IF

      !

#ifdef W3_OMPH

      !$OMP PARALLEL DO PRIVATE (IX, IY, IXY, VCB)

#endif

      !

      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

      !

#ifdef W3_OMPH

      !$OMP END PARALLEL DO

#endif

      !


      ! 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

        !

#ifdef W3_OMPH

        !$OMP PARALLEL DO PRIVATE (IXY, VCB, IX)

#endif

        !

        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

        !

#ifdef W3_OMPH

        !$OMP END PARALLEL DO

#endif

        !

      END IF


      ! 4.  Propagate ------------------------------------------------------ *

      !

#ifdef W3_T3

      WRITE (ndst,9040)

#endif

      !

#ifdef W3_OMPH

      !$OMP PARALLEL DO PRIVATE (ISEA, IXY, JXN, JXP, JYN, JYP, IX, IY)

#endif

      !

      DO isea=1, nsea

        !

        ix    = mapsf(isea,1)

        iy    = mapsf(isea,2)

        ixy   = mapsf(isea,3)


#ifdef W3_T3

        aold   = fld2d(iy,ix) * cg(ik,isea) / clats(isea)

#endif

        !

        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)


#ifdef W3_T3

          WRITE (ndst,9041) isea, ixy, ixy-ny, ixy-1,         &

               vflx2d(iy,ix-1), vflx2d(iy,ix), vfly2d(iy-1,ix), &

               vfly2d(iy,ix) , cg(ik,isea)/clats(isea),aold,    &

               fld2d(iy,ix)

#endif

          !

          !

#ifdef W3_T3

          WRITE (ndst,9042) isea, mapsta(ixy), aold,fld2d(iy,ix)

#endif

          !

        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

      !

#ifdef W3_OMPH

      !$OMP END PARALLEL DO

#endif

      !


      ! 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

    !

#ifdef W3_T

9000 FORMAT (' TEST W3XYP1 : NTLOC :',i4)

9001 FORMAT (' TEST W3XYP1 : ISP, ITH, IK :',i8,2i4)

#endif

    !

#ifdef W3_T1

9010 FORMAT (' TEST W3XYP1 : INIT. VFX-YL, ITLOC =',i3)

9020 FORMAT (' TEST W3XYP1 : ISEA, IXY, FIELD, VCX, VCY')

9021 FORMAT ('           ',2i8,3e12.4)

9024 FORMAT (' TEST W3XYP1 : GLOBAL CLOSURE: IY, FIELD, VCX ')

9025 FORMAT ('               ',i4,2e12.4)

#endif

    !

#ifdef W3_T2

9032 FORMAT (' TEST W3XYP1 : CLOSE. : IY, VFLX')

9033 FORMAT ('            ',i4,e12.4)

#endif

    !

#ifdef W3_T3

9040 FORMAT (' TEST W3XYP1 : PROPAGATION '/                      &

         '      ISEA, IXY(3), , FLX(2), FLY(2), FAC, A(2)')

9041 FORMAT (2x,4i5,1x,4e10.3,1x,e10.3,1x,2e10.3)

9042 FORMAT (2x,i5,'( MAP = ',i2,' )',56x,2e10.3)

#endif

    !/

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

#ifdef W3_DEBUG

    USE w3odatmd, only : iaproc

#endif

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    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

    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)

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

    CALL strace (ient, 'W3KTP1')

#endif

    !/

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

    !/

#ifdef W3_T

    WRITE (ndst,9000) flcth, flck, facth, fack, ctmax

    WRITE (ndst,9001) isea, depth, cx, cy,                       &

         dddx, dddy, dcxdx, dcxdy, dcydx, dcydy

#endif

    !

    ! 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

    !

#ifdef W3_T

    WRITE (ndst,9010)

    DO ik=1, nk+1

      WRITE (ndst,9011) ik, tpi/sig(ik), tpi/wn(ik),             &

           cg(ik), dsdd(ik)

    END DO

#endif

    !

    ! 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

      !

#ifdef W3_REFRX

      ! 2.c @C/@x refraction and great-circle propagation

      vcth = 0.

      frk  = 0.

      fddmax = 0.

      !

      DO isp=1, nspec

        fddmax = max( fddmax , abs(                      &

             esin(isp)*dcdx(mapwn(isp)) - ecos(isp)*dcdy(mapwn(isp)) ) )

      END DO

      !

      DO ik=1, nk

        frk(ik) = facth * cg(ik) * wn(ik) / sig(ik)

        frk(ik) = frk(ik) / max( 1. , frk(ik)*fddmax/ctmax )

        frg(ik) = fdg * cg(ik)

      END DO

      DO isp=1, nspec

        vcth(isp) = frg(mapwn(isp)) * ecos(isp)            &

             + frk(mapwn(isp)) * ( esin(isp)*dcdx(mapwn(isp)) &

             - ecos(isp)*dcdy(mapwn(isp)) )

      END DO

#endif

      !

      ! 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

    !

#ifdef W3_T

9000 FORMAT (' TEST W3KTP1 : FLCTH-K, FACTH-K, CTMAX  :',         &

         2l2,2e10.3,f7.3)

9001 FORMAT (' TEST W3KTP1 : LOCAL DATA :',i7,f7.1,2f6.2,1x,      &

         6e10.3)

9010 FORMAT (' TEST W3KTP1 : IK, T, L, CG, DSDD : ')

9011 FORMAT ('              ',i3,f7.2,f7.1,f7.2,e11.3)

#endif

    !/

    !/ End of W3KTP1 ----------------------------------------------------- /

    !/

  END SUBROUTINE w3ktp1

  !/

  !/ End of module W3PRO1MD -------------------------------------------- /

  !/

END MODULE w3pro1md