w3pro2md.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
 
MODULE w3pro2md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         29-May-2014 |
  !/                  +-----------------------------------+
  !/
  !/    04-Feb-2000 : Origination.                        ( version 2.00 )
  !/    24-Jan-2001 : Flat grid version                   ( version 2.06 )
  !/    08-Feb-2001 : UQ routines moved to own module     ( version 2.08 )
  !/    09-Feb-2001 : Clean up of parameter lists         ( version 2.08 )
  !/    14-Feb-2001 : Unit numbers in UQ routines         ( version 2.08 )
  !/    13-Nov-2001 : Sub-grid obstacles added.           ( version 2.14 )
  !/    26-Dec-2002 : Moving grid option.                 ( version 3.02 )
  !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    07-Sep-2005 : Improved XY boundary conditions.    ( version 3.08 )
  !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
  !/    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)
  !/    23-Dec-2010 : Fix HPFAC and HQFAC by including the COS(YGRD)
  !/                  factor with DXDP and DXDQ terms.    ( version 3.14 )
  !/                  (T. J. Campbell, NRL)
  !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
  !/                  and second order schemes.           ( version 4.12 )
  !/    29-May-2014 : Adding OMPH switch.                 ( version 5.02 )
  !/
  !/    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 third order porpagation scheme in single
  !     module.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      TRNMIN    R.P.  Private   Minimum transparancy for local
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3MAP2    Subr. Public   Set up auxiliary maps.
  !      W3XYP2    Subr. Public   Third order spatial propagation.
  !      W3KTP2    Subr. Public   Third order spectral propagation.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      W3QCK1    Subr. W3UQCKMD Regular grid UQ scheme.
  !      W3QCK2    Subr.   Id.    Irregular grid UQ scheme.
  !      W3QCK3    Subr.   Id.    Regular grid UQ scheme + obstructions.
  !      W3UNO2    Subr. W3UNO2MD UNO2 scheme for irregular grid.
  !      W3UNO2r   Subr.   Id.    UNO2 scheme reduced to regular grid.
  !      W3UNO2s   Subr.   Id.    UNO2 regular grid with subgrid
  !                               obstruction.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !       !/UQ    3rd order UQ propagation scheme.
  !       !/UNO   2nd order UNO propagation scheme.
  !
  !       !/MGP   Correct for motion of grid.
  !
  !       !/OMPH  Hybrid OpenMP directives.
  !
  !       !/TDYN, !/DSS0, !/XW0, !/XW1
  !               Diffusion options in W3XYP2
  !
  !       !/S     Enable subroutine tracing.
  !       !/Tn    Enable test output.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !/
  !/ Public variables
  !/
  PUBLIC
  !/
  !/ Private data
  !/
  REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3map2
    !/
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Nov-2005 |
    !/                  +-----------------------------------+
    !/
    !/    19-Dec-1996 : Final FORTRAN 77                    ( version 1.18 )
    !/    15-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    09-Feb-2001 : Clean up of parameter lists         ( version 2.08 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
    !/
    !  1. Purpose :
    !
    !     Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
    !
    !  2. Method :
    !
    !     MAPX2, MAPY2, MAPTH2 and MAPWN2 contain consecutive 1-D spatial
    !     grid counters (e.g., IXY = (IX-1)*MY + IY). The arrays are
    !     devided in three parts. For MAPX2, these ranges are :
    !
    !         1    - NMX0  Counters IXY for which grid point (IX,IY) and
    !                      (IX+1,IY) both are active grid points.
    !       NMX0+1 - NMX1  Id. only (IX,IY) active.
    !       NMX1+1 - NMX2  Id. only (IX+1,IY) active.
    !
    !     The array MAPY2 has a similar layout varying IY instead of IX.
    !     MAPXY contains similar information for the cross term in the
    !     diffusion correction (counter NMXY only).
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  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.   Map MAPX2
    !        a  Range 1 to NMX0
    !        b  Range NMX0+1 to NMX1
    !        c  Range NMX1+1 to NMX2
    !      2.   Map MAPY2
    !        a  Range 1 to NMY0
    !        b  Range NMY0+1 to NMY1
    !        c  Range NMY1+1 to NMY2
    !      3.   Map MAPAXY and MAPXY
    !      4.   Maps for intra-spectral propagation
    !        a  MAPTH2, MAPATK
    !        b  MAPWN2
    !     ------------------------------------------------------
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable test output.
    !
    ! 10. Source code :
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nk, nth, nspec, nx, ny, nsea, mapsta
    USE w3adatmd, ONLY: nmx0, nmx1, nmx2, nmy0, nmy1, nmy2, nact,   &
         nmxy, mapx2, mapy2, mapaxy, mapxy,          &
         mapth2, mapwn2
    USE w3odatmd, ONLY: ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IX, IY, IXY0, IX2, IY2, IX0, IY0,    &
         IK, ITH, ISP, ISP0, ISP2
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T
    INTEGER                 :: MAPTXY(NY,NX), I, IXY
    INTEGER                 :: MAPTST(NK+2,NTH)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3MAP2')
#endif
    !
    ! 1.  Map MAPX2 ------------------------------------------------------ *
    ! 1.a Range 1 to NMX0
    !
#ifdef W3_T
    maptxy = 0.
#endif
    !
    nmx0   = 0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      ix2    = 1 + mod(ix,nx)
      DO iy=2, ny-1
        IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy,ix2).EQ.1 ) THEN
          nmx0   = nmx0 + 1
          mapx2(nmx0) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 1
#endif
        END IF
      END DO
    END DO
    !
    ! 1.b Range NMX0+1 to NMX1
    !
    nmx1   = nmx0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      ix2    = 1 + mod(ix,nx)
      DO iy=2, ny-1
        IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy,ix2).NE.1 ) THEN
          nmx1   = nmx1 + 1
          mapx2(nmx1) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 2
#endif
        END IF
      END DO
    END DO
    !
    ! 1.c Range NMX1+1 to NMX2
    !
    nmx2   = nmx1
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      ix2    = 1 + mod(ix,nx)
      DO iy=2, ny-1
        IF ( mapsta(iy,ix).NE.1 .AND. mapsta(iy,ix2).EQ.1 ) THEN
          nmx2   = nmx2 + 1
          mapx2(nmx2) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 4
#endif
        END IF
      END DO
    END DO
    !
 
#ifdef W3_T
    WRITE (ndst,9000) 'MAPX2', nmx0, nmx1-nmx0,                  &
         nmx2-nmx1, nmx2
    DO iy=ny, 1, -1
      WRITE (ndst,9001) (maptxy(iy,ix),ix=1, nx)
    END DO
#endif
    !
    ! 2.  Map MAPY2 ------------------------------------------------------ *
    ! 2.a Range 1 to NMY0
    !
#ifdef W3_T
    maptxy = 0.
#endif
    !
    nmy0   = 0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      DO iy=1, ny-1
        iy2    = iy + 1
        IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy2,ix).EQ.1 ) THEN
          nmy0   = nmy0 + 1
          mapy2(nmy0) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 1
#endif
        END IF
      END DO
    END DO
    !
    ! 2.b Range NMY0+1 to NMY1
    !
    nmy1   = nmy0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      DO iy=1, ny-1
        iy2    = iy + 1
        IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy2,ix).NE.1 ) THEN
          nmy1   = nmy1 + 1
          mapy2(nmy1) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 2
#endif
        END IF
      END DO
    END DO
    !
    ! 2.c Range NMY1+1 to NMY2
    !
    nmy2   = nmy1
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      DO iy=1, ny-1
        iy2    = iy + 1
        IF ( mapsta(iy,ix).NE.1 .AND. mapsta(iy2,ix).EQ.1 ) THEN
          nmy2   = nmy2 + 1
          mapy2(nmy2) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 4
#endif
        END IF
      END DO
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9000) 'MAPY2', nmy0, nmy1-nmy0,                  &
         nmy2-nmy1, nmy2
    DO iy=ny, 1, -1
      WRITE (ndst,9001) (maptxy(iy,ix),ix=1, nx)
    END DO
#endif
    !
    ! 3.  Map MAPAXY and MAPXY ------------------------------------------- *
    !
    nact   = 0
    DO ix=1, nx
      iy0    = (ix-1)*ny
      DO iy=2, ny-1
        IF ( mapsta(iy,ix).EQ.1 ) THEN
          nact         = nact + 1
          mapaxy(nact) = iy0 + iy
        END IF
      END DO
    END DO
    !
    nmxy   = 0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      ix2    = ix+1
      IF ( ix .EQ. nx ) ix2 =  1
      DO iy=2, ny-1
        IF ( mapsta( iy ,ix ).GE.1 .AND.                            &
             mapsta( iy ,ix2).GE.1 .AND.                            &
             mapsta(iy+1,ix ).GE.1 .AND.                            &
             mapsta(iy+1,ix2).GE.1 ) THEN
          nmxy   = nmxy + 1
          mapxy(nmxy) = ixy0 + iy
        END IF
      END DO
    END DO
    !
    ! 4.  Maps for intra-spectral propagation ---------------------------- *
    !
    IF ( mapth2(1) .NE. 0 ) RETURN
    !
#ifdef W3_T
    maptst = 0
#endif
    !
    ! 4.a MAPTH2 and MAPBTK
    !
    DO ik=1, nk
      DO ith=1, nth
        isp    = ith + (ik-1)*nth
        isp2   = (ik+1) + (ith-1)*(nk+2)
        mapth2(isp) = isp2
#ifdef W3_T
        maptst(ik+1,ith) = maptst(ik+1,ith) + 1
#endif
      END DO
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9000) 'MAPTH2', isp, 0, 0, isp
    DO ik=nk+2, 1, -1
      WRITE (ndst,9001) (maptst(ik,ith),ith=1, nth)
    END DO
    maptst = 0
#endif
    !
    ! 4.b MAPWN2
    !
    isp0   = 0
    DO ik=1, nk-1
      DO ith=1, nth
        isp0   = isp0 + 1
        isp2   = (ik+1) + (ith-1)*(nk+2)
        mapwn2(isp0) = isp2
#ifdef W3_T
        maptst(ik+1,ith) = maptst(ik+1,ith) + 1
#endif
      END DO
    END DO
    !
    DO ith=1, nth
      isp0   = isp0 + 1
      isp2   = nk+1 + (ith-1)*(nk+2)
      mapwn2(isp0) = isp2
#ifdef W3_T
      maptst(nk+1,ith) = maptst(nk+1,ith) + 2
#endif
    END DO
    !
    DO ith=1, nth
      isp0   = isp0 + 1
      isp2   = 1 + (ith-1)*(nk+2)
      mapwn2(isp0) = isp2
#ifdef W3_T
      maptst(1,ith) = maptst(1,ith) + 4
#endif
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9000) 'MAPWN2', nspec-nth, nth, nth, nspec+nth
    DO ik=nk+2, 1, -1
      WRITE (ndst,9001) (maptst(ik,ith),ith=1, nth)
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (/' TEST W3MAP2 : TEST MAP FOR PROPAGATION'/          &
         '      MAP     : ',a/                               &
         '      CENTRAL : ',i6/                              &
         '      ABOVE   : ',i6/                              &
         '      BELOW   : ',i6/                              &
         '      TOTAL   : ',i6/)
9001 FORMAT (1x,130i1)
9010 FORMAT (' TEST W3MAP2 : COMPOSITE MAPS TH2, WN2 AND BTK')
9011 FORMAT (2x,60i2)
#endif
    !/
    !/ End of W3MAP2 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3map2
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3xyp2 ( ISP, DTG, MAPSTA, MAPFS, VQ, 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 )
    !/    15-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    24-Jan-2001 : Flat grid version                   ( version 2.06 )
    !/    09-Feb-2001 : Clean up of parameter lists         ( version 2.08 )
    !/    14-Feb-2001 : Unit numbers in UQ routines         ( version 2.08 )
    !/    13-Nov-2001 : Sub-grid obstructions.              ( version 2.14 )
    !/    26-Dec-2002 : Moving grid option,                 ( version 3.02 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    07-Sep-2005 : Improved boundary conditions.       ( version 3.08 )
    !/    09-Nov-2005 : Removing soft boundary option.      ( 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)
    !/    23-Dec-2010 : Fix HPFAC and HQFAC by including the COS(YGRD)
    !/                  factor with DXDP and DXDQ terms.    ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
    !/                  and second order schemes.           ( version 4.12 )
    !/    29-May-2014 : Adding OMPH switch.                 ( version 5.02 )
    !/
    !  1. Purpose :
    !
    !     Propagation in physical space for a given spectral component.
    !
    !  2. Method :
    !
    !     Third-order ULTIMATE QUICKEST scheme and diffusion correction
    !     for linear dispersion (see manual).
    !     Curvilinear grid implementation: Fluxes are computed in index space
    !       and then transformed back into physical space.  The diffusion term
    !       is handled in 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.
    !       MAPFS   I.A.   I   Storage map.
    !       VQ      R.A.  I/O  Field to propagate.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       NTLOC   Int   Number of local time steps.
    !       DTLOC   Real  Local propagation time step.
    !       CGD     Real  Deep water group velocity.
    !       DSSD, DNND    Deep water diffusion coefficients.
    !       VLCFLX  R.A.  Local courant numbers in index space (norm. velocities)
    !       VLCFLY  R.A.
    !       CXTOT   R.A.  Propagation velocities in physical space.
    !       CYTOT   R.A.
    !       CELLP   Real  Cell Reynolds/Peclet number used to calculate
    !                     diffusion coefficient for growing spectral
    !                     components.
    !       DFRR    Real  Relative frequency increment.
    !     ----------------------------------------------------------------
    !
    !  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 :
    !
    !     - Note that the ULTIMATE limiter does not guarantee non-zero
    !       energies.
    !     - The present scheme shows a strong distortion when propaga-
    !       ting a field under an angle with the grid in a truly 2-D
    !       fashion. Propagation is therefore split along the two
    !       axes.
    !     - Two boundary treatments are available. The first uses real
    !       boundaries in each space. In this case waves will not
    !       penetrate in narrow straights under an angle with the grid.
    !       This behavior is improved by using a 'soft' option, in
    !       which the 'X' or 'Y' sweep allows for energy to go onto
    !       the land. This improves the above behavior, but implies
    !       that X-Y connenctions are required in barriers for them
    !       to become inpenetrable.
    !     - If TDYN is set to zero, ALL diffusion is skipped. Set TDYN
    !       to a small positive number to have growth diffusion only.
    !     - 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.
    !     - Factors VFDIFX_FAC VFDIFY_FAC VFDIFC_FAC are introduced so that results
    !       match for test case tp2.3.  Use of these factors is optional and removal
    !       can significantly reduce size/cost of code. These variants are marked as
    !       CURV1 or CURV2. NCEP will make final decision re: which version to adopt.
    !       CURV1 is the shorter version and results do not match the original code
    !       for all test cases. CURV2 is the longer version and results do match the
    !       original. DETAILED EXPLANATION: Discrepancies occur at the boundaries.
    !       This is because, at the boundaries, the pre-curvilinear version zeros out
    !       some terms in the diffusion calculation. Since they are zeroed out,
    !       they aren't there to *cancel* certain other terms: these "other terms"
    !       affect the result, so they have to be retained in the long vesion (CURV2)
    !       to get an exact match. In the short version, the "canceling out" is
    !       performed prior to coding the scheme, so both the canceled and canceling
    !       terms are always omitted.
    !
    !  8. Structure :
    !
    !     ---------------------------------------------
    !       1.  Preparations
    !         a Set constants
    !         b Initialize arrays
    !       2.  Prepare arrays
    !         a Velocities and 'Q'
    !         b diffusion coefficients
    !       3.  Loop over sub-steps
    !       ----------------------------------------
    !         a Propagate
    !         b Update boundary conditions
    !         c Diffusion correction
    !       ----------------------------------------
    !       4.  Store Q field in spectra
    !     ---------------------------------------------
    !
    !  9. Switches :
    !
    !       !/MGP   Correct for motion of grid.
    !
    !       !/TDYN  Dynamic increase of DTME
    !       !/DSS0  Disable diffusion in propagation direction
    !       !/XW0   Propagation diffusion only.
    !       !/XW1   Growth diffusion only.
    !
    !       !/OMPH  Hybrid OpenMP directives.
    !
    !       !/S     Enable subroutine tracing.
    !
    !       !/T     Enable general test output.
    !       !/T1    Dump of input field and fluxes.
    !       !/T2    Dump of output field.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !
    USE w3timemd, ONLY: dsec21
    !
    USE w3gdatmd, ONLY: nk, nth, dth, xfr, esin, ecos, sig, nx, ny, &
         nsea, sx, sy, mapsf, iclose, flcx, flcy,    &
         iclose_none, iclose_smpl, iclose_trpl,      &
         dtcfl, clats, dtme, clatmn, flagll,         &
         hpfac, hqfac, dpdx, dpdy, dqdx, dqdy, gsqrt
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cg, wn, u10, cx, cy, atrnx, atrny, itime,   &
         nmx0, nmx1, nmx2, nmy0, nmy1, nmy2, nact,   &
         nmxy, mapx2, mapy2, mapaxy, mapxy
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin,       &
         isbpi, bbpi0, bbpin, iaproc, naperr
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_UQ
    USE w3uqckmd
#endif
#ifdef W3_UNO
    USE w3uno2md
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISP, MAPSTA(NY*NX), MAPFS(NY*NX)
    REAL, INTENT(IN)        :: DTG, VGX, VGY
    REAL, INTENT(INOUT)     :: VQ(1-NY:NY*(NX+2))
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK, NTLOC, ITLOC, ISEA, IXY,    &
         IX,IY, IY0, IP, IBI
    INTEGER                 :: TTEST(2),DTTST
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: CG0, CGA, CGN, CGX, CGY, CXC, CYC,  &
         CXMIN, CXMAX, CYMIN, CYMAX
    REAL                    :: DTLOC, DTRAD,               &
         DFRR, CELLP,  CGD, DSSD,    &
         DNND, DCELL, XWIND, TFAC, DSS, DNN
    REAL                    :: RD1, RD2
    REAL                    :: RFAC, DFAC, DVQ, QXX, QXY, QYY
    REAL                    :: CP, CQ
    LOGICAL                 :: YFIRST
    LOGICAL                 :: GLOBAL
    !/
    !/ Automatic work arrays
    !/
    REAL                    :: VLCFLX((NX+1)*NY), VLCFLY((NX+1)*NY),&
         VFDIFX(1-NY:NX*NY), VFDIFY(NX*NY),   &
         VFDIFC(1-NY:NX*NY), VDXX((NX+1)*NY), &
         VDYY(NX*NY), VDXY((NX+1)*NY)
 
    REAL                    :: CXTOT((NX+1)*NY), CYTOT(NX*NY)
    REAL                    :: VQ_OLD(1-NY:NY*(NX+2))
    !CURV2: BEGIN -----------------------------------------------------------------
    REAL                    :: VFDIFX_FAC(1-NY:NX*NY),              &
         VFDIFY_FAC(1-NY:NX*NY),              &
         VFDIFC_FAC(1-NY:NX*NY)
    !CURV2: END -------------------------------------------------------------------
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3XYP2')
#endif
    !
    !     IF ( MAXVAL(VQ) .EQ. 0. ) THEN
    !         IF ( NBI .EQ. 0 ) THEN
    !             RETURN
    !           ELSE
    !             IF ( MAXVAL(BBPI0(ISP,:)) .EQ. 0. .AND.            &
    !                  MAXVAL(BBPIN(ISP,:)) .EQ. 0. ) RETURN
    !           END IF
    !       END IF
    !
    ! 1.  Preparations --------------------------------------------------- *
 
    IF ( iclose .EQ. iclose_trpl ) THEN
      IF (iaproc .EQ. naperr) &
           WRITE(ndse,*)'SUBROUTINE W3XYP2 IS NOT YET ADAPTED FOR '//    &
           'TRIPOLE GRIDS. STOPPING NOW.'
      CALL extcde ( 1 )
    END IF
 
    ! 1.a Set constants
    !
    global = iclose.NE.iclose_none
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
    !
    cg0    = 0.575 * grav / sig(1)
    cga    = 0.575 * grav / sig(ik)
    cgx    = cga * ecos(ith)
    cgy    = cga * esin(ith)
#ifdef W3_MGP
    cgx    = cgx - vgx
    cgy    = cgy - vgy
#endif
    !
    IF ( flcur ) THEN
      cxmin  = minval( cx(1:nsea) )
      cxmax  = maxval( cx(1:nsea) )
      cymin  = minval( cy(1:nsea) )
      cymax  = maxval( cy(1:nsea) )
      IF ( abs(cgx+cxmin) .GT. abs(cgx+cxmax) ) THEN
        cgx    = cgx + cxmin
      ELSE
        cgx    = cgx + cxmax
      END IF
      IF ( abs(cgy+cymin) .GT. abs(cgy+cymax) ) THEN
        cgy    = cgy + cymin
      ELSE
        cgy    = cgy + cymax
      END IF
      cxc    = max( abs(cxmin) , abs(cxmax) )
      cyc    = max( abs(cymin) , abs(cymax) )
#ifdef W3_MGP
      cxc    = max( abs(cxmin-vgx) , abs(cxmax-vgx) )
      cyc    = max( abs(cymin-vgy) , abs(cymax-vgy) )
#endif
    ELSE
      cxc    = 0.
      cyc    = 0.
    END IF
    !
    cgn    = 0.9999 * max( abs(cgx) , abs(cgy) , cxc, cyc, 0.001*cg0 )
    !
    ntloc  = 1 + int(dtg/(dtcfl*cg0/cgn))
    dtloc  = dtg / real(ntloc)
    dtrad  = dtloc
    IF ( flagll ) dtrad=dtrad/(dera*radius)
    !
    IF ( flagll ) THEN
      rfac = dera * radius
      dfac = 1. / rfac**2
    ELSE
      rfac = 1.
      dfac = 1.
    END IF
    !
    ttest(1) = time(1)
    ttest(2) = 0
    dttst = dsec21(ttest,time)
    yfirst = mod(nint(dttst/dtg),2) .EQ. 0
    !
#ifdef W3_T
    WRITE (ndst,9000) yfirst
    WRITE (ndst,9001) isp, ith, ik, ecos(ith), esin(ith)
#endif
    !
#ifdef W3_TDYN
    IF ( isp .EQ. 1 ) dtme = dtme + dtg
#endif
    !
    IF ( dtme .NE. 0. ) THEN
      dfrr   = xfr - 1.
      cellp  = 10.
      cgd    = 0.5 * grav / sig(ik)
      dssd   = ( dfrr * cgd )**2 * dtme / 12.
      dnnd   = ( cgd * dth )**2 * dtme / 12.
#ifdef W3_T
      WRITE (ndst,9002) dfrr, cellp, dtme
    ELSE
      WRITE (ndst,9003)
#endif
    END IF
 
    !
    ! 1.b Initialize arrays
    !
#ifdef W3_T
    WRITE (ndst,9010)
#endif
    !
    vlcflx = 0.
    vlcfly = 0.
    vfdifx = 0.
    vfdify = 0.
    vfdifc = 0.
    vdxx   = 0.
    vdyy   = 0.
    vdxy   = 0.
    cxtot  = 0.
    cytot  = 0.
    !
    ! 2.  Calculate velocities and diffusion coefficients ---------------- *
    ! 2.a Velocities
    !
    !     Q     = ( A / CG * CLATS )
    !     LCFLX = ( COS*CG / CLATS ) * DT / DX
    !     LCFLY = (     SIN*CG )     * DT / DX
    !
#ifdef W3_T
    WRITE (ndst,9020) nsea
#endif
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IXY)
#endif
    !
    DO isea=1, nsea
      ixy         = mapsf(isea,3)
      vq(ixy) = vq(ixy) / cg(ik,isea) * clats(isea)
      cxtot(ixy) = ecos(ith) * cg(ik,isea) / clats(isea)
      cytot(ixy) = esin(ith) * cg(ik,isea)
#ifdef W3_MGP
      cxtot(ixy) = cxtot(ixy) - vgx/clats(isea)
      cytot(ixy) = cytot(ixy) - vgy
#endif
#ifdef W3_T1
      IF ( .NOT. flcur )                                        &
           WRITE (ndst,9021) isea, mapsf(isea,1), mapsf(isea,2), &
           vq(ixy), cxtot(ixy), cytot(ixy)
#endif
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    IF ( flcur ) THEN
#ifdef W3_T
      WRITE (ndst,9022)
#endif
      DO isea=1, nsea
        ixy         = mapsf(isea,3)
        cxtot(ixy) = cxtot(ixy) + cx(isea)/clats(isea)
        cytot(ixy) = cytot(ixy) + cy(isea)
#ifdef W3_T1
        WRITE (ndst,9021) isea, mapsf(isea,1), mapsf(isea,2), &
             vq(ixy), cxtot(ixy), cytot(ixy)
#endif
      END DO
    END IF
 
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, CP, CQ)
#endif
    !
    DO isea=1, nsea
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      ixy    = mapsf(isea,3)
      cp=cxtot(ixy)*dpdx(iy,ix)+cytot(ixy)*dpdy(iy,ix)
      cq=cxtot(ixy)*dqdx(iy,ix)+cytot(ixy)*dqdy(iy,ix)
      vlcflx(ixy) = cp*dtrad
      vlcfly(ixy) = cq*dtrad
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    ! 2.b Diffusion coefficients
    !
    IF ( dtme .NE. 0. ) THEN
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, &
      !$OMP&                     DCELL, XWIND, TFAC, DSS, DNN)
#endif
      !
      DO isea=1, nsea
        ix        = mapsf(isea,1)
        iy        = mapsf(isea,2)
        ixy       = mapsf(isea,3)
        IF ( min( atrnx(ixy,1) , atrnx(ixy,-1) ,                 &
             atrny(ixy,1) , atrny(ixy,-1) ) .GT. trnmin ) THEN
          dcell     = cgd * min( hpfac(iy,ix)*rfac, &
               hqfac(iy,ix)*rfac ) / cellp
          xwind     = 3.3 * u10(isea)*wn(ik,isea)/sig(ik) - 2.3
          xwind     = max( 0. , min( 1. , xwind ) )
#ifdef W3_XW0
          xwind     = 0.
#endif
#ifdef W3_XW1
          xwind     = 1.
#endif
          tfac      = min( 1. , (clats(isea)/clatmn)**2 )
          dss       = xwind * dcell + (1.-xwind) * dssd * tfac
#ifdef W3_DSS0
          dss       = 0.
#endif
          dnn       = xwind * dcell + (1.-xwind) * dnnd * tfac
 
          vdxx(ixy) = dtloc * (dss*ecos(ith)**2+dnn*esin(ith)**2)
          vdyy(ixy) = dtloc * (dss*esin(ith)**2+dnn*ecos(ith)**2) &
               / clats(isea)**2
          vdxy(ixy) = dtloc * (dss-dnn) * esin(ith)*ecos(ith)     &
               / clats(isea)
 
        END IF
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
    END IF
    !
    ! 3.  Loop over sub-steps -------------------------------------------- *
    !
    DO itloc=1, ntloc
      !
      ! 3.a Propagate fields
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY )
#endif
      !
      DO isea=1, nsea
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
        vq(ixy)= vq(ixy) * gsqrt(iy,ix)
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
      IF ( yfirst ) THEN
        !
#ifdef W3_UQ
        IF ( flcy ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcfly, atrny, vq,       &
             .false., 1, mapaxy, nact, mapy2, nmy0,   &
             nmy1, nmy2, ndse, ndst )
        IF ( flcx ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,       &
             global, ny, mapaxy, nact, mapx2, nmx0,   &
             nmx1, nmx2, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        IF ( flcy ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcfly, atrny, vq,      &
             .false., 1, mapaxy, nact, mapy2, nmy0,  &
             nmy1, nmy2, ndse, ndst )
        IF ( flcx ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,      &
             global, ny, mapaxy, nact, mapx2, nmx0,  &
             nmx1, nmx2, ndse, ndst )
#endif
        !
      ELSE
        !
#ifdef W3_UQ
        IF ( flcx ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,       &
             global, ny, mapaxy, nact, mapx2, nmx0,   &
             nmx1, nmx2, ndse, ndst )
        IF ( flcy ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcfly, atrny, vq,       &
             .false., 1, mapaxy, nact, mapy2, nmy0,   &
             nmy1, nmy2, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        IF ( flcx ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,      &
             global, ny, mapaxy, nact, mapx2, nmx0,  &
             nmx1, nmx2, ndse, ndst )
        IF ( flcy ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcfly, atrny, vq,      &
             .false., 1, mapaxy, nact, mapy2, nmy0,  &
             nmy1, nmy2, ndse, ndst )
#endif
        !
      END IF
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY )
#endif
      !
      DO isea=1, nsea
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
        vq(ixy)= vq(ixy) / gsqrt(iy,ix)
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
      ! 3.b Update boundaries
      !
      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
          isea    = isbpi(ibi)
          ixy     = mapsf(isbpi(ibi),3)
          vq(ixy) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) )   &
               / cg(ik,isea) * clats(isea)
        END DO
      END IF
      !
      ! 3.c Diffusion correction
      !
      IF ( dtme .NE. 0. ) THEN
 
        IF ( global ) THEN
          DO iy=1, ny
            vq(iy+nx*ny) = vq(iy)
          END DO
        END IF
 
        !CURV2: BEGIN -----------------------------------------------------------------
        vfdifx_fac=0.0
        DO ip=1, nmx0
          ixy         = mapx2(ip)
          vfdifx_fac(ixy) = 1.0
        END DO
        vfdify_fac=0.0
        DO ip=1, nmy0
          ixy         = mapy2(ip)
          vfdify_fac(ixy) = 1.0
        END DO
        vfdifc_fac=0.0
        DO ip=1, nmxy
          ixy         = mapxy(ip)
          vfdifc_fac(ixy) = 1.0
        END DO
        IF ( global ) THEN
          iy0     = (nx-1)*ny
          DO iy=1, ny
            vfdifx_fac(iy-ny) = vfdifx_fac(iy+iy0)
            vfdifc_fac(iy-ny) = vfdifc_fac(iy+iy0)
          END DO
        END IF
        !CURV2: END -------------------------------------------------------------------
 
        vq_old = vq
        !
#ifdef W3_OMPH
        !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, &
        !$OMP&                     QXX, QYY, QXY, DVQ )
#endif
        !
        DO ip=1, nact
          ixy    = mapaxy(ip)
          isea   = mapfs(ixy)
          ix     = mapsf(isea,1)
          iy     = mapsf(isea,2)
 
          !CURV1: DOES NOT GIVE EXACT MATCH TO EARLIER WW3 VERSION FOR TEST CASE TP2.3
          !CURV1: NEAR THE BOUNDARY, NOTE THAT THIS VERSION USES NON-ACTIVE GRID POINTS
          !CURV1: IN ITS CALCS. THIS IS NO PROBLEM, AS LONG AS THE NON-ACTIVE GRID POINTS
          !CURV1: EXIST IN THE ARRAY AND ARE VQ=0. ALSO NOTE THAT WITH THE SHORT VERSION
          !CURV1: OF THE CODE, VFDIF?_FAC VARIABLES CAN BE REMOVED AND 3-4 DO LOOPS ABOVE
          !CURV1: CAN BE REMOVED.
          !CURV1: BEGIN -----------------------------------------------------------------
          !                QXX = VQ_OLD(IXY+NY) - 2.0*VQ_OLD(IXY) + VQ_OLD(IXY-NY)
          !                QYY = VQ_OLD(IXY+1)  - 2.0*VQ_OLD(IXY) + VQ_OLD(IXY-1)
          !                QXY = VQ_OLD(IXY+NY+1) - VQ_OLD(IXY-NY+1) &
          !                    - VQ_OLD(IXY+NY-1) + VQ_OLD(IXY-NY-1)
          !CURV1: END -------------------------------------------------------------------
 
          !CURV2: DOES GIVE EXACT MATCH TO EARLIER WW3 VERSION FOR TEST CASE TP2.3. NOTE
          !CURV2: THAT IF VFDIFC_FAC VARIABLES ARE ALL UNITY, MANY TERMS CANCEL OUT.
          !CURV2: HOWEVER, VFDIFC_FAC IS ZERO WHEN A RELATED VQ POINT IS NOT AN ACTIVE
          !CURV2: GRID POINT
          !CURV2: BEGIN -----------------------------------------------------------------
          qxx = vfdifx_fac(ixy)   *vq_old(ixy+ny)     &
               - vfdifx_fac(ixy)   *vq_old(ixy)        &
               - vfdifx_fac(ixy-ny)*vq_old(ixy)        &
               + vfdifx_fac(ixy-ny)*vq_old(ixy-ny)
          qyy = vfdify_fac(ixy)   *vq_old(ixy+1)      &
               - vfdify_fac(ixy)   *vq_old(ixy)        &
               - vfdify_fac(ixy-1) *vq_old(ixy)        &
               + vfdify_fac(ixy-1) *vq_old(ixy-1)
          qxy = vfdifc_fac(ixy)     *vq_old(ixy)      &
               + vfdifc_fac(ixy-ny-1)*vq_old(ixy)      &
               - vfdifc_fac(ixy-1)   *vq_old(ixy)      &
               - vfdifc_fac(ixy-ny)  *vq_old(ixy)      &
               + vfdifc_fac(ixy-ny)  *vq_old(ixy+1)    &
               - vfdifc_fac(ixy)     *vq_old(ixy+1)    &
               + vfdifc_fac(ixy-1)   *vq_old(ixy-1)    &
               - vfdifc_fac(ixy-ny-1)*vq_old(ixy-1)    &
               + vfdifc_fac(ixy-1)   *vq_old(ixy+ny)   &
               - vfdifc_fac(ixy)     *vq_old(ixy+ny)   &
               + vfdifc_fac(ixy-ny)  *vq_old(ixy-ny)   &
               - vfdifc_fac(ixy-ny-1)*vq_old(ixy-ny)   &
               + vfdifc_fac(ixy)     *vq_old(ixy+ny+1) &
               - vfdifc_fac(ixy-1)   *vq_old(ixy+ny-1) &
               + vfdifc_fac(ixy-ny-1)*vq_old(ixy-ny-1) &
               - vfdifc_fac(ixy-ny)  *vq_old(ixy-ny+1)
          !CURV2: END -------------------------------------------------------------------
          !
          qxy = 0.25*qxy
          !
          dvq = vdxx(ixy)*(  dpdx(iy,ix)*dpdx(iy,ix)*qxx          &
               + 2.0*dpdx(iy,ix)*dqdx(iy,ix)*qxy      &
               + dqdx(iy,ix)*dqdx(iy,ix)*qyy )        &
               + vdyy(ixy)*(  dpdy(iy,ix)*dpdy(iy,ix)*qxx          &
               + 2.0*dpdy(iy,ix)*dqdy(iy,ix)*qxy      &
               + dqdy(iy,ix)*dqdy(iy,ix)*qyy)         &
               + 2.0*vdxy(ixy)*(  dpdx(iy,ix)*dpdy(iy,ix)*qxx      &
               + dqdx(iy,ix)*dqdy(iy,ix)*qyy      &
               +      (  dpdx(iy,ix)*dqdy(iy,ix)  &
               + dqdx(iy,ix)*dpdy(iy,ix) )*qxy )
          !
          vq(ixy) = vq_old(ixy) + dvq * dfac
          !
        END DO
        !
#ifdef W3_OMPH
        !$OMP END PARALLEL DO
#endif
        !
      END IF
      !
      yfirst = .NOT. yfirst
    END DO
    !
    ! 4.  Store results in VQ in proper format --------------------------- *
    !
#ifdef W3_T
    WRITE (ndst,9040) nsea
#endif
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IXY )
#endif
    !
    DO isea=1, nsea
      ixy    = mapsf(isea,3)
      IF ( mapsta(ixy) .GT. 0 ) THEN
#ifdef W3_T2
        WRITE (ndst,9041) isea, mapsf(isea,1), mapsf(isea,2), vq(ixy)
#endif
        vq(ixy) =  max( 0. , cg(ik,isea) / clats(isea) * vq(ixy) )
        !         ELSE
        !           VQ(IXY) =  0.
      END IF
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3XYP2 : YFIRST  :',l2)
9001 FORMAT (' TEST W3XYP2 : ISP, ITH, IK, COS-SIN :',i8,2i4,2f7.3)
9002 FORMAT (' TEST W3XYP2 : DFRR, CELLP, DTME     :',3e10.3)
9003 FORMAT (' TEST W3XYP2 : NO DISPERSION CORRECTION ')
9010 FORMAT (' TEST W3XYP2 : INITIALIZE ARRAYS')
9020 FORMAT (' TEST W3XYP2 : CALCULATING LCFLX/Y AND DSS/NN (NSEA=', &
         i6,')')
9022 FORMAT (' TEST W3XYP2 : CORRECTING FOR CURRENT')
9040 FORMAT (' TEST W3XYP2 : FIELD AFTER PROP. (NSEA=',i6,')')
#endif
#ifdef W3_T1
9021 FORMAT (1x,i6,2i5,e12.4,2f7.3)
#endif
#ifdef W3_T2
9041 FORMAT (1x,i6,2i5,e12.4)
#endif
    !/
    !/ End of W3XYP2 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3xyp2
  !/
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3ktp2 ( 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 :         01-Jul-2013 |
    !/                  +-----------------------------------+
    !/
    !/    14-Feb-2000 : Origination.                        ( version 2.08 )
    !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
    !/                  and second order schemes.           ( version 4.12 )
    !/
    !  1. Purpose :
    !
    !     Propagation in spectral space.
    !
    !  2. Method :
    !
    !     Third order QUICKEST scheme with ULTIMATE limiter.
    !
    !     As with the spatial propagation, the two spaces are considered
    !     independently, but the propagation is performed in a 2-D space.
    !     Compared to the propagation in physical space, the directions
    !     rerpesent a closed space and are therefore comparable to the
    !     longitudinal or 'X' propagation. The wavenumber space has to be
    !     extended to allow for boundary treatment. Using a simple first
    !     order boundary treatment at both sided, two points need to
    !     be added. This implies that the spectrum needs to be extended,
    !     shifted and rotated, as is performed using MAPTH2 as set
    !     in W3MAP3.
    !
    !  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.
    !       MAPxx2  I.A.   I   Propagation and storage maps.
    !       CG      R.A.   I   Local group velocities.
    !       WN      R.A.   I   Local wavenumbers.
    !       DEPTH   R.A.   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.
    !       FDD, FDU, FDG, FCD, FCU
    !               R.A.  Directionally varying part of depth, current and
    !                     great-circle refraction terms and of consit.
    !                     of Ck term.
    !       CFLT-K  R.A.  Propagation velocities of local fluxes.
    !       DB      R.A.  Wavenumber band widths at cell centers.
    !       DM      R.A.  Wavenumber band widths between cell centers and
    !                     next cell center.
    !       Q       R.A.  Extracted spectrum
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       W3QCK1   Actual propagation routine.
    !       W3QCK2   Actual propagation routine.
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3WAVE   Wave model routine.
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  8. Structure :
    !
    !     -----------------------------------------------------------------
    !       1.  Preparations
    !         a Initialize arrays
    !         b Set constants and counters
    !       2.  Point  preparations
    !         a Calculate DSDD
    !         b Extract spectrum
    !       3.  Refraction velocities
    !         a Filter level depth reffraction.
    !         b Depth refratcion velocity.
    !         c Current refraction velocity.
    !       4.  Wavenumber shift velocities
    !         a Prepare directional arrays
    !         b Calcuate velocity.
    !       5.  Propagate.
    !       6.  Store results.
    !     -----------------------------------------------------------------
    !
    !  9. Switches :
    !
    !       !/S     Enable subroutine tracing.
    !       !/T     Enable general test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nk2, nth, nspec, sig, dsip, ecos, esin, &
         ec2, esc, es2, fachfa, mapwn, flcth, flck,  &
         ctmax
    USE w3adatmd, ONLY: mapth2, mapwn2, itime
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndse, ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_UQ
    USE w3uqckmd
#endif
#ifdef W3_UNO
    USE w3uno2md
#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
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: FDDMAX, FDG, FKD, FKD0, DCYX,        &
         DCXXYY, DCXY, DCXX, DCXYYX, DCYY
    REAL                    :: DSDD(0:NK+1), FRK(NK), FRG(NK),      &
         FKC(NTH), VQ(-NK-1:NK2*(NTH+2)),     &
         DB(NK2,NTH+1), DM(NK2,0:NTH+1),      &
         VCFLT(NK2*(NTH+1)), CFLK(NK2,NTH)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3KTP2')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Initialize arrays
    !
    IF ( flck ) vq = 0.
    !
#ifdef W3_T
    WRITE (ndst,9000) flcth, flck, facth, fack, ctmax
    WRITE (ndst,9010) isea, depth, cx, cy, dddx, dddy,           &
         dcxdx, dcxdy, dcydx, dcydy
#endif
    !
    ! 2.  Preparation for point ------------------------------------------ *
    ! 2.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,9020)
    DO ik=1, nk+1
      WRITE (ndst,9021) ik, tpi/sig(ik), tpi/wn(ik),             &
           cg(ik), dsdd(ik)
    END DO
#endif
    !
    ! 2.b Extract spectrum
    !
    DO isp=1, nspec
      vq(mapth2(isp)) = va(isp)
    END DO
    !
    ! 3.  Refraction velocities ------------------------------------------ *
    !
    IF ( flcth ) THEN
      !
      ! 3.a Set slope filter for depth refraction
      !
      fddmax = 0.
      fdg    = facth * cthg0
      !
      DO ith=1, nth/2
        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
      !
      ! 3.b Depth refraction and great-circle propagation
      !
      DO isp=1, nspec
        vcflt(mapth2(isp)) = frg(mapwn(isp)) * ecos(isp)          &
             + frk(mapwn(isp)) * ( esin(isp)*dddx - ecos(isp)*dddy )
      END DO
      !
#ifdef W3_REFRX
      ! 3.c @C/@x refraction and great-circle propagation
      vcflt  = 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)
      END DO
      DO isp=1, nspec
        vcflt(mapth2(isp)) = frg(mapwn(isp)) * ecos(isp)   &
             + frk(mapwn(isp)) * ( esin(isp)*dcdx(mapwn(isp)) &
             - ecos(isp)*dcdy(mapwn(isp)) )
      END DO
#endif
      !
      ! 3.d Current refraction
      !
      IF ( flcur ) THEN
        !
        dcyx   = facth *   dcydx
        dcxxyy = facth * ( dcxdx - dcydy )
        dcxy   = facth *   dcxdy
        !
        DO isp=1, nspec
          vcflt(mapth2(isp)) = vcflt(mapth2(isp)) +             &
               es2(isp)*dcyx  + esc(isp)*dcxxyy - ec2(isp)*dcxy
        END DO
        !
      END IF
      !
    END IF
    !
    ! 4.  Wavenumber shift velocities ------------------------------------ *
    !     FACK is just the time step, which is accounted for in W3QCK2
    !
    IF ( flck ) THEN
      !
      ! 4.a Directionally dependent part
      !
      dcxx   =  -   dcxdx
      dcxyyx =  - ( dcxdy + dcydx )
      dcyy   =  -   dcydy
      fkd    =    ( cx*dddx + cy*dddy )
      !
      DO ith=1, nth
        fkc(ith) = ec2(ith)*dcxx +                                &
             esc(ith)*dcxyyx + es2(ith)*dcyy
      END DO
      !
      ! 4.b Velocities
      !
      DO ik=0, nk+1
        fkd0   = fkd / cg(ik) * dsdd(ik)
        DO ith=1, nth
          cflk(ik+1,ith) = fkd0 + wn(ik)*fkc(ith)
        END DO
      END DO
      !
      ! 4.c Band widths
      !
      DO ik=0, nk
        db(ik+1,1) = dsip(ik) / cg(ik)
        dm(ik+1,1) = wn(ik+1) - wn(ik)
      END DO
      db(nk+2,1) = dsip(nk+1) / cg(nk+1)
      dm(nk+2,1) = 0.
      !
      DO ith=2, nth
        DO ik=1, nk+2
          db(ik,ith) = db(ik,1)
          dm(ik,ith) = dm(ik,1)
        END DO
      END DO
      !
    END IF
    !
    ! 5.  Propagate ------------------------------------------------------ *
    !
    IF ( mod(itime,2) .EQ. 0 ) THEN
      IF ( flck ) THEN
        DO ith=1, nth
          vq(nk+2+(ith-1)*nk2) = fachfa * vq(nk+1+(ith-1)*nk2)
        END DO
        !
#ifdef W3_UQ
        CALL w3qck2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,   &
             vq, .false., 1, mapth2, nspec,            &
             mapwn2, nspec-nth, nspec, nspec+nth,      &
             ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,  &
             vq, .false., 1, mapth2, nspec,           &
             mapwn2, nspec-nth, nspec, nspec+nth,     &
             ndse, ndst )
#endif
      END IF
      IF ( flcth ) THEN
        !
#ifdef W3_UQ
        CALL w3qck1 ( nth, nk2, nth, nk2, vcflt, vq, .true.,    &
             nk2, mapth2, nspec, mapth2, nspec, nspec, &
             nspec, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2r( nth, nk2, nth, nk2, vcflt, vq, .true.,   &
             nk2, mapth2, nspec, mapth2, nspec, nspec,&
             nspec, ndse, ndst )
#endif
        !
      END IF
    ELSE
      IF ( flcth ) THEN
        !
#ifdef W3_UQ
        CALL w3qck1 ( nth, nk2, nth, nk2, vcflt, vq, .true.,    &
             nk2, mapth2, nspec, mapth2, nspec, nspec, &
             nspec, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2r( nth, nk2, nth, nk2, vcflt, vq, .true.,   &
             nk2, mapth2, nspec, mapth2, nspec, nspec,&
             nspec, ndse, ndst )
#endif
        !
      END IF
      IF ( flck )  THEN
        DO ith=1, nth
          vq(nk+2+(ith-1)*nk2) = fachfa * vq(nk+1+(ith-1)*nk2)
        END DO
        !
#ifdef W3_UQ
        CALL w3qck2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,   &
             vq, .false., 1, mapth2, nspec,            &
             mapwn2, nspec-nth, nspec, nspec+nth,      &
             ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,  &
             vq, .false., 1, mapth2, nspec,           &
             mapwn2, nspec-nth, nspec, nspec+nth,     &
             ndse, ndst )
#endif
        !
      END IF
    END IF
    !
    ! 6.  Store reults --------------------------------------------------- *
    !
    DO isp=1, nspec
      va(isp) = vq(mapth2(isp))
    END DO
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3KTP2 : FLCTH-K, FACTH-K, CTMAX  :',        &
         2l2,2e10.3,f7.3)
9010 FORMAT ( ' TEST W3KTP2 : LOCAL DATA :',i7,f7.1,2f6.2,1x,6e10.2)
9020 FORMAT ( ' TEST W3KTP2 : IK, T, L, CG, DSDD : ')
9021 FORMAT ( '              ',i3,f7.2,f7.1,f7.2,e11.3)
#endif
    !
#ifdef W3_T0
9040 FORMAT (/' TEST W3KTP2 : NORMALIZED ',a/)
9041 FORMAT (1x,60(1x,i2))
9042 FORMAT (1x,60i3)
#endif
    !/
    !/ End of W3KTP2 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3ktp2
  !/
  !/ End of module W3PRO2MD -------------------------------------------- /
  !/
END MODULE w3pro2md