w3uqckmd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3uqckmd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         27-May-2014 |
  !/                  +-----------------------------------+
  !/
  !/    08-Feb-2001 : Origination of module. Routines     ( version 2.08 )
  !/                  taken out of w3pro2md.ftn
  !/    13-Nov-2001 : Version with obstacles added.       ( version 2.14 )
  !/    16-Oct-2002 : Fix par list W3QCK3.                ( version 3.00 )
  !/    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 : Fixed a couple of doc lines.        ( version 3.14 )
  !/                  (T. J. Campbell, NRL)
  !/    27-May-2014 : Added OMPH switches in W3QCK3.      ( 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 :
  !
  !     Portable ULTIMATE QUICKEST schemes.
  !
  !  2. Variables and types :
  !
  !     None.
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3QCK1    Subr. Public   Original ULTIMATE QUICKEST scheme.
  !      W3QCK2    Subr. Public   UQ scheme for irregular grid.
  !      W3QCK3    Subr. Public   Original ULTIMATE QUICKEST with obst.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !     - STRACE and !/S irrelevant for running code. The module is
  !       therefore fully portable to any other model.
  !
  !  6. Switches :
  !
  !       !/OMPH  Ading OMP directves for hybrid paralellization.
  !
  !       !/S     Enable subroutine tracing.
  !       !/Tn    Enable test output.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
#ifdef W3_S
  USE w3servmd, ONLY: strace
#endif
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3qck1 (MX, MY, NX, NY, CFLL, Q, CLOSE, INC,         &
       MAPACT, NACT, MAPBOU, NB0, NB1, NB2,         &
       NDSE, NDST )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         30-Oct-2009 |
    !/                  +-----------------------------------+
    !/
    !/    11-Mar-1997 : Final FORTRAN 77                    ( version 1.18 )
    !/    15-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    15-Feb-2001 : Unit numbers added to par list.     ( version 2.08 )
    !/    05-Mar-2008 : Added NEC sxf90 compiler directives.
    !/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
    !/    30-Oct-2009 : Fixed "Called by" doc line.         ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/
    !  1. Purpose :
    !
    !     Preform one-dimensional propagation in a two-dimensional space
    !     with irregular boundaries and regular grid.
    !
    !  2. Method :
    !
    !     ULTIMATE QUICKEST scheme (see manual).
    !
    !     Note that the check on monotonous behavior of QCN is performed
    !     using weights CFAC, to avoid the need for IF statements.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       MX,MY   Int.   I   Field dimensions, if grid is 'closed' or
    !                          circular, MX is the closed dimension.
    !       NX,NY   Int.   I   Part of field actually used.
    !       CFLL    R.A.   I   Local Courant numbers.         (MY,  MX+1)
    !       Q       R.A.  I/O  Propagated quantity.           (MY,0:MX+2)
    !       CLOSE   Log.   I   Flag for closed 'X' dimension'
    !       INC     Int.   I   Increment in 1-D array corresponding to
    !                          increment in 2-D space.
    !       MAPACT  I.A.   I   List of active grid points.
    !       NACT    Int.   I   Size of MAPACT.
    !       MAPBOU  I.A.   I   Map with boundary information (see W3MAP2).
    !       NBn     Int.   I   Counters in MAPBOU.
    !       NDSE    Int.   I   Error output unit number.
    !       NDST    Int.   I   Test output unit number.
    !     ----------------------------------------------------------------
    !           - CFLL amd Q need only bee filled in the (MY,MX) range,
    !             extension is used internally for closure.
    !           - CFLL and Q are defined as 1-D arrays internally.
    !
    !  4. Subroutines used :
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3KTP2   Propagation in spectral space
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - This routine can be used independently from WAVEWATCH III.
    !
    !  8. Structure :
    !
    !     ------------------------------------------------------
    !       1. Initialize aux. array FLA.
    !       2. Fluxes for central points (3rd order + limiter).
    !       3. Fluxes boundary point above (1st order).
    !       4. Fluxes boundary point below (1st order).
    !       5. Closure of 'X' if required
    !       6. Propagate.
    !     ------------------------------------------------------
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable test output.
    !     !/T0  Test output input/output fields.
    !     !/T1  Test output fluxes.
    !     !/T2  Test output integration.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: MX, MY, NX, NY, INC, MAPACT(MY*MX),  &
         NACT, MAPBOU(MY*MX), NB0, NB1, NB2,  &
         NDSE, NDST
    REAL, INTENT(INOUT)     :: CFLL(MY*(MX+1)), Q(1-MY:MY*(MX+2))
    LOGICAL, INTENT(IN)     :: CLOSE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IXY, IP, IXYC, IXYU, IXYD, IY, IX,   &
         IAD00, IAD02, IADN0, IADN1, IADN2
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T1
    INTEGER                 :: IX2, IY2
#endif
    REAL                    :: CFL, QB, DQ, DQNZ, QCN, QBN, QBR, CFAC
    REAL                    :: FLA(1-MY:MY*MX)
#ifdef W3_T0
    REAL                    :: QMAX
#endif
#ifdef W3_T1
    REAL                    :: QBO, QN
#endif
#ifdef W3_T2
    REAL                    :: QOLD
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3QCK1')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) mx, my, nx, ny, CLOSE, inc, nb0, nb1, nb2
#endif
    !
#ifdef W3_T0
    qmax   = 0.
    DO iy=1, ny
      DO ix=1, nx
        qmax   = max( qmax , q(iy+(ix-1)*my) )
      END DO
    END DO
    qmax   = max( 0.01*qmax , 1.e-10 )
#endif
    !
#ifdef W3_T0
    WRITE (ndst,9001) 'CFLL'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(100.*cfll(iy+(ix-1)*my)),ix=1,nx)
    END DO
    WRITE (ndst,9001) 'Q'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(q(iy+(ix-1)*my)/qmax),ix=1,nx)
    END DO
    WRITE (ndst,9001) 'MAPACT'
    WRITE (ndst,9003) (mapact(ixy),ixy=1,nact)
#endif
    !
    ! 1.  Initialize aux. array FLA and closure ------------------------- *
    !
    fla = 0.
    !
    IF ( CLOSE ) THEN
#ifdef W3_T
      WRITE (ndst,9005)
#endif
      iad00  = -my
      iad02  =  my
      iadn0  = iad00 + my*nx
      iadn1  =         my*nx
      iadn2  = iad02 + my*nx
      DO iy=1, ny
        q(iy+iad00) = q(iy+iadn0)
        q(iy+iadn1) = q(   iy   )
        q(iy+iadn2) = q(iy+iad02)
        cfll(iy+iadn1) = cfll(   iy   )
      END DO
    END IF
    !
    ! 2.  Fluxes for central points ------------------------------------- *
    !     ( 3rd order + limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9010)
    WRITE (ndst,9011) nb0, 'CENTRAL'
#endif
    !
    DO ip=1, nb0
      !
      ixy    = mapbou(ip)
      cfl    = 0.5 * ( cfll(ixy) + cfll(ixy+inc) )
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,cfl) ) )
      qb     = 0.5 * ( (1.-cfl)*q(ixy+inc) + (1.+cfl)*q(ixy) )      &
           - (1.-cfl**2)/6. * (q(ixyc-inc)-2.*q(ixyc)+q(ixyc+inc))
      !
      ixyu   = ixyc - inc * int( sign(1.1,cfl) )
      ixyd   = 2*ixyc - ixyu
      dq     = q(ixyd) - q(ixyu)
      dqnz   = sign( max(1.e-15,abs(dq)) , dq )
      qcn    = ( q(ixyc) - q(ixyu) ) / dqnz
      qcn    = min( 1.1, max( -0.1 , qcn ) )
      !
#ifdef W3_T1
      qbo    = qb
#endif
      qbn    = max( (qb-q(ixyu))/dqnz , qcn )
      qbn    = min( qbn , 1. , qcn/max(1.e-10,abs(cfl)) )
      qbr    = q(ixyu) + qbn*dq
      cfac   = real( int( 2. * abs(qcn-0.5) ) )
      qb     = (1.-cfac)*qbr + cfac*q(ixyc)
      !
      fla(ixy) = cfl * qb
      !
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( qb, qbo, q(ixy-inc), q(   ixy   ),         &
           q(ixy+inc), q(ixy+2*inc) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9012) ip, ix, iy, ix2, iy2,               &
             cfl, cfll(ixy), cfll(ixy+inc),      &
             qbo*qn, qb*qn, q(ixy-inc)*qn, q(   ixy   )*qn, &
             q(ixy+inc)*qn, q(ixy+2*inc)*qn
      END IF
#endif
      !
    END DO
    !
    ! 3.  Fluxes for points with boundary above ------------------------- *
    !     ( 1st order without limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9011) nb1-nb0, 'BOUNDARY ABOVE'
#endif
    !
    DO ip=nb0+1, nb1
      ixy    = mapbou(ip)
      cfl    = cfll(ixy)
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,cfl) ) )
      fla(ixy) = cfl * q(ixyc)
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( q(ixy+inc), q(ixy) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9013) ip, ix, iy, ix2, iy2, cfl,          &
             cfll(ixy), q(ixyc)*qn, q(ixy)*qn, q(ixy+inc)*qn
      END IF
#endif
    END DO
    !
    ! 4.  Fluxes for points with boundary below ------------------------- *
    !     ( 1st order without limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9011) nb2-nb1, 'BOUNDARY BELOW'
#endif
    !
    DO ip=nb1+1, nb2
      ixy    = mapbou(ip)
      cfl    = cfll(ixy+inc)
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,cfl) ) )
      fla(ixy) = cfl * q(ixyc)
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( q(ixy+inc), q(ixy) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9014) ip, ix, iy, ix2, iy2, cfl,         &
             cfll(ixy+inc), q(ixyc)*qn, q(ixy)*qn, q(ixy+inc)*qn
      END IF
#endif
    END DO
    !
    ! 5.  Global closure ----------------------------------------------- *
    !
    IF ( CLOSE ) THEN
#ifdef W3_T
      WRITE (ndst,9015)
#endif
      DO iy=1, ny
        fla(iy+iad00) = fla(iy+iadn0)
      END DO
    END IF
    !
    ! 6.  Propagation -------------------------------------------------- *
    !
#ifdef W3_T2
    WRITE (ndst,9020)
#endif
    DO ip=1, nact
      ixy    = mapact(ip)
#ifdef W3_T2
      qold   = q(ixy)
#endif
      q(ixy) = max( 0. , q(ixy) + fla(ixy-inc) - fla(ixy) )
#ifdef W3_T2
      IF ( qold + q(ixy) .GT. 1.e-10 )                          &
           WRITE (ndst,9021) ip, ixy, qold, q(ixy),             &
           fla(ixy-inc), fla(ixy)
#endif
    END DO
    !
#ifdef W3_T0
    WRITE (ndst,9001) 'Q'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(q(iy+(ix-1)*my)/qmax),ix=1,nx)
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3QCK1 : ARRAY DIMENSIONS  :',2i6/           &
         '               USED              :',2i6/           &
         '               CLOSE, INC        :',l6,i6/         &
         '               NB0, NB1, NB2     :',3i6)
#endif
#ifdef W3_T0
9001 FORMAT ( ' TEST W3QCK1 : DUMP ARRAY ',a,' :')
9002 FORMAT ( 1x,43i3)
9003 FORMAT ( 1x,21i6)
#endif
#ifdef W3_T
9005 FORMAT (' TEST W3QCK1 : GLOBAL CLOSURE (1)')
#endif
    !
#ifdef W3_T1
9010 FORMAT (' TEST W3QCK1 : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
         ' Q (b,b,i-1,i,i+1,i+2)')
9011 FORMAT (' TEST W3QCK1 :',i6,' POINTS OF TYPE ',a)
9012 FORMAT (10x,i6,4i4,1x,3f6.2,1x,f7.2,f6.2,1x,4f6.2)
9013 FORMAT (10x,i6,4i4,1x,f6.2,f6.2,'  --- ',1x,f7.2,1x,'  --- ',&
         2f6.2,'  --- ')
9014 FORMAT (10x,i6,4i4,1x,f6.2,'  --- ',f6.2,1x,f7.2,1x,'  --- ',&
         2f6.2,'  --- ')
#endif
#ifdef W3_T
9015 FORMAT (' TEST W3QCK1 : GLOBAL CLOSURE (2)')
#endif
    !
#ifdef W3_T2
9020 FORMAT (' TEST W3QCK1 : IP, IXY, 2Q, 2FL')
9021 FORMAT ('            ',2i6,2(1x,2e11.3))
#endif
    !/
    !/ End of W3QCK1 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3qck1
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3qck2 (MX, MY, NX, NY, VELO, DT, DX1, DX2, Q, CLOSE,&
       INC,  MAPACT, NACT, MAPBOU, NB0, NB1, NB2,    &
       NDSE, NDST )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         30-Oct-2009 |
    !/                  +-----------------------------------+
    !/
    !/    07-Sep-1997 : Final FORTRAN 77                    ( version 1.18 )
    !/    16-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    14-Feb-2001 : Unit numbers added to par list.     ( version 2.08 )
    !/    05-Mar-2008 : Added NEC sxf90 compiler directives.
    !/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
    !/    30-Oct-2009 : Fixed "Called by" doc line.         ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/
    !  1. Purpose :
    !
    !     Like W3QCK1 with variable grid spacing.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       MX,MY   Int.   I   Field dimensions, if grid is 'closed' or
    !                          circular, MX is the closed dimension.
    !       NX,NY   Int.   I   Part of field actually used.
    !       VELO    R.A.   I   Local velocities.              (MY,  MX+1)
    !       DT      Real   I   Time step.
    !       DX1     R.A.  I/O  Band width at points.          (MY,  MX+1)
    !       DX2     R.A.  I/O  Band width between points.     (MY,0:MX+1)
    !                          (local counter and counter+INC)
    !       Q       R.A.  I/O  Propagated quantity.           (MY,0:MX+2)
    !       CLOSE   Log.   I   Flag for closed 'X' dimension'
    !       INC     Int.   I   Increment in 1-D array corresponding to
    !                          increment in 2-D space.
    !       MAPACT  I.A.   I   List of active grid points.
    !       NACT    Int.   I   Size of MAPACT.
    !       MAPBOU  I.A.   I   Map with boundary information (see W3MAP2).
    !       NBn     Int.   I   Counters in MAPBOU.
    !       NDSE    Int.   I   Error output unit number.
    !       NDST    Int.   I   Test output unit number.
    !     ----------------------------------------------------------------
    !           - VELO amd Q need only bee filled in the (MY,MX) range,
    !             extension is used internally for closure.
    !           - VELO and Q are defined as 1-D arrays internally.
    !
    !  4. Subroutines used :
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3KTP2   Propagation in spectral space
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - This routine can be used independently from WAVEWATCH III.
    !
    !  8. Structure :
    !
    !     ------------------------------------------------------
    !       1. Initialize aux. array FLA.
    !       2. Fluxes for central points (3rd order + limiter).
    !       3. Fluxes boundary point above (1st order).
    !       4. Fluxes boundary point below (1st order).
    !       5. Closure of 'X' if required
    !       6. Propagate.
    !     ------------------------------------------------------
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable test output.
    !     !/T0  Test output input/output fields.
    !     !/T1  Test output fluxes.
    !     !/T2  Test output integration.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: MX, MY, NX, NY, INC, MAPACT(MY*MX),  &
         NACT, MAPBOU(MY*MX), NB0, NB1, NB2,  &
         NDSE, NDST
    REAL, INTENT(IN)        :: DT
    REAL, INTENT(INOUT)     :: VELO(MY*(MX+1)), DX1(MY*(MX+1)),     &
         DX2(1-MY:MY*(MX+1)), Q(1-MY:MY*(MX+2))
    LOGICAL, INTENT(IN)     :: CLOSE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IXY, IP, IXYC, IXYU, IXYD, IY, IX,   &
         IAD00, IAD02, IADN0, IADN1, IADN2
#ifdef W3_S
    INTEGER, SAVE           :: IENT
#endif
#ifdef W3_T1
    INTEGER                 :: IX2, IY2
#endif
    REAL                    :: CFL, VEL, QB, DQ, DQNZ, QCN, QBN,    &
         QBR, CFAC, FLA(1-MY:MY*MX)
#ifdef W3_T0
    REAL                    :: QMAX
#endif
#ifdef W3_T1
    REAL                    :: QBO, QN, XCFL
#endif
#ifdef W3_T2
    REAL                    :: QOLD
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3QCK2')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) mx, my, nx, ny, dt, CLOSE, inc, nb0, nb1, nb2
#endif
    !
#ifdef W3_T0
    qmax   = 0.
    DO iy=1, ny
      DO ix=1, nx
        qmax   = max( qmax , q(iy+(ix-1)*my) )
      END DO
    END DO
    qmax   = max( 0.01*qmax , 1.e-10 )
#endif
    !
#ifdef W3_T0
    WRITE (ndst,9001) 'VELO'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(100.*velo(iy+(ix-1)*my)           &
           *dt/dx1(iy+(ix-1)*my)),ix=1,nx)
    END DO
    WRITE (ndst,9001) 'Q'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(q(iy+(ix-1)*my)/qmax),ix=1,nx)
    END DO
    WRITE (ndst,9001) 'MAPACT'
    WRITE (ndst,9003) (mapact(ixy),ixy=1,nact)
#endif
    !
    ! 1.  Initialize aux. array FLA and closure ------------------------- *
    !
    fla = 0.
    !
    IF ( CLOSE ) THEN
#ifdef W3_T
      WRITE (ndst,9005)
#endif
      iad00  = -my
      iad02  =  my
      iadn0  = iad00 + my*nx
      iadn1  =         my*nx
      iadn2  = iad02 + my*nx
      DO iy=1, ny
        q(iy+iad00) = q(iy+iadn0)
        q(iy+iadn1) = q(   iy   )
        q(iy+iadn2) = q(iy+iad02)
        velo(iy+iadn1) = velo(   iy   )
        dx1(iy+iadn1) = dx1(   iy   )
        dx2(iy+iad00) = dx1(iy+iadn0)
        dx2(iy+iadn1) = dx1(   iy   )
      END DO
    END IF
    !
    ! 2.  Fluxes for central points ------------------------------------- *
    !     ( 3rd order + limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9010)
    WRITE (ndst,9011) nb0, 'CENTRAL'
#endif
    !
    DO ip=1, nb0
      !
      ixy    = mapbou(ip)
      vel    = 0.5 * ( velo(ixy) + velo(ixy+inc) )
      cfl    = dt *  vel / dx2(ixy)
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,cfl) ) )
      qb     = 0.5 * ( (1.-cfl)*q(ixy+inc) + (1.+cfl)*q(ixy) )      &
           - dx2(ixy)**2 / dx1(ixyc) * (1.-cfl**2) / 6.      &
           * ( (q(ixyc+inc)-q(ixyc))/dx2(ixyc)           &
           - (q(ixyc)-q(ixyc-inc))/dx2(ixyc-inc) )
      !
      ixyu   = ixyc - inc * int( sign(1.1,cfl) )
      ixyd   = 2*ixyc - ixyu
      dq     = q(ixyd) - q(ixyu)
      dqnz   = sign( max(1.e-15,abs(dq)) , dq )
      qcn    = ( q(ixyc) - q(ixyu) ) / dqnz
      qcn    = min( 1.1, max( -0.1 , qcn ) )
      !
#ifdef W3_T1
      qbo    = qb
#endif
      qbn    = max( (qb-q(ixyu))/dqnz , qcn )
      qbn    = min( qbn , 1. , qcn/max(1.e-10,abs(cfl)) )
      qbr    = q(ixyu) + qbn*dq
      cfac   = real( int( 2. * abs(qcn-0.5) ) )
      qb     = (1.-cfac)*qbr + cfac*q(ixyc)
      !
      fla(ixy) = vel * qb
      !
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( qb, qbo, q(ixy-inc), q(   ixy   ),         &
           q(ixy+inc), q(ixy+2*inc) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9012) ip, ix, iy, ix2, iy2,               &
             cfl, dt*velo(ixy)/dx1(ixy),                 &
             dt*velo(ixy+inc)/dx1(ixy+inc),         &
             qbo*qn, qb*qn, q(ixy-inc)*qn, q(   ixy   )*qn,  &
             q(ixy+inc)*qn, q(ixy+2*inc)*qn
      END IF
#endif
      !
    END DO
    !
    ! 3.  Fluxes for points with boundary above ------------------------- *
    !     ( 1st order without limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9011) nb1-nb0, 'BOUNDARY ABOVE'
#endif
    !
    DO ip=nb0+1, nb1
      ixy    = mapbou(ip)
      vel    = velo(ixy)
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,vel) ) )
      fla(ixy) = vel * q(ixyc)
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( q(ixy+inc), q(ixy) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9013) ip, ix, iy, ix2, iy2, xcfl,             &
             dt*velo(ixy)/dx2(ixy),                  &
             q(ixyc)*qn, q(ixy)*qn, q(ixy+inc)*qn
      END IF
#endif
    END DO
    !
    ! 4.  Fluxes for points with boundary below ------------------------- *
    !     ( 1st order without limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9011) nb2-nb1, 'BOUNDARY BELOW'
#endif
    !
    DO ip=nb1+1, nb2
      ixy    = mapbou(ip)
      vel    = velo(ixy+inc)
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,vel) ) )
      fla(ixy) = vel * q(ixyc)
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( q(ixy+inc), q(ixy) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9014) ip, ix, iy, ix2, iy2, xcfl,         &
             dt*velo(ixy+inc)/dx2(ixy),          &
             q(ixyc)*qn, q(ixy)*qn, q(ixy+inc)*qn
      END IF
#endif
    END DO
    !
    ! 5.  Global closure ----------------------------------------------- *
    !
    IF ( CLOSE ) THEN
#ifdef W3_T
      WRITE (ndst,9015)
#endif
      DO iy=1, ny
        fla(iy+iad00) = fla(iy+iadn0)
      END DO
    END IF
    !
    ! 6.  Propagation -------------------------------------------------- *
    !
#ifdef W3_T2
    WRITE (ndst,9020)
#endif
    DO ip=1, nact
      ixy    = mapact(ip)
#ifdef W3_T2
      qold   = q(ixy)
#endif
      q(ixy) = max( 0. , q(ixy) + dt/dx1(ixy) *                    &
           (fla(ixy-inc)-fla(ixy)) )
#ifdef W3_T2
      IF ( qold + q(ixy) .GT. 1.e-10 )                          &
           WRITE (ndst,9021) ip, ixy, qold, q(ixy),             &
           dt*fla(ixy-inc)/dx1(ixy),          &
           dt*fla(ixy)/dx1(ixy)
#endif
    END DO
    !
#ifdef W3_T0
    WRITE (ndst,9001) 'Q'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(q(iy+(ix-1)*my)/qmax),ix=1,nx)
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3QCK2 : ARRAY DIMENSIONS  :',2i6/           &
         '               USED              :',2i6/           &
         '               TIME STEP         :',f8.1/          &
         '               CLOSE, INC        :',l6,i6/         &
         '               NB0, NB1, NB2     :',3i6)
#endif
#ifdef W3_T0
9001 FORMAT ( ' TEST W3QCK2 : DUMP ARRAY ',a,' :')
9002 FORMAT ( 1x,43i3)
9003 FORMAT ( 1x,21i6)
#endif
#ifdef W3_T
9005 FORMAT (' TEST W3QCK2 : GLOBAL CLOSURE (1)')
#endif
    !
#ifdef W3_T1
9010 FORMAT (' TEST W3QCK2 : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
         ' Q (b,b,i-1,i,i+1,i+2)')
9011 FORMAT (' TEST W3QCK2 :',i6,' POINTS OF TYPE ',a)
9012 FORMAT (10x,i6,4i4,1x,3f6.2,1x,f7.2,f6.2,1x,4f6.2)
9013 FORMAT (10x,i6,4i4,1x,f6.2,f6.2,'  --- ',1x,f7.2,1x,'  --- ',&
         2f6.2,'  --- ')
9014 FORMAT (10x,i6,4i4,1x,f6.2,'  --- ',f6.2,1x,f7.2,1x,'  --- ',&
         2f6.2,'  --- ')
#endif
#ifdef W3_T
9015 FORMAT (' TEST W3QCK2 : GLOBAL CLOSURE (2)')
#endif
    !
#ifdef W3_T2
9020 FORMAT (' TEST W3QCK2 : IP, IXY, 2Q, 2FL')
9021 FORMAT ('            ',2i6,2(1x,2e11.3))
#endif
    !/
    !/ End of W3QCK2 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3qck2
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3qck3 (MX, MY, NX, NY, CFLL, TRANS, Q, CLOSE,       &
       INC, MAPACT, NACT, MAPBOU, NB0, NB1, NB2,    &
       NDSE, NDST )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         27-May-2014 |
    !/                  +-----------------------------------+
    !/
    !/    13_nov-2001 : Origination.                        ( version 2.14 )
    !/    16-Oct-2002 : Fix INTENT for TRANS.               ( version 3.00 )
    !/    05-Mar-2008 : Added NEC sxf90 compiler directives.
    !/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
    !/    27-May-2014 : Added OMPH switches in W3QCK3.      ( version 5.02 )
    !/
    !  1. Purpose :
    !
    !     Like W3QCK1 with cell transparencies added.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       MX,MY   Int.   I   Field dimensions, if grid is 'closed' or
    !                          circular, MX is the closed dimension.
    !       NX,NY   Int.   I   Part of field actually used.
    !       CFLL    R.A.   I   Local Courant numbers.         (MY,  MX+1)
    !       Q       R.A.  I/O  Propagated quantity.           (MY,0:MX+2)
    !       CLOSE   Log.   I   Flag for closed 'X' dimension'
    !       INC     Int.   I   Increment in 1-D array corresponding to
    !                          increment in 2-D space.
    !       MAPACT  I.A.   I   List of active grid points.
    !       NACT    Int.   I   Size of MAPACT.
    !       MAPBOU  I.A.   I   Map with boundary information (see W3MAP2).
    !       NBn     Int.   I   Counters in MAPBOU.
    !       NDSE    Int.   I   Error output unit number.
    !       NDST    Int.   I   Test output unit number.
    !     ----------------------------------------------------------------
    !           - CFLL amd Q need only bee filled in the (MY,MX) range,
    !             extension is used internally for closure.
    !           - CFLL and Q are defined as 1-D arrays internally.
    !
    !  4. Subroutines used :
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3XYP2   Propagation in physical space
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - This routine can be used independently from WAVEWATCH III.
    !
    !  8. Structure :
    !
    !     ------------------------------------------------------
    !       1. Initialize aux. array FLA.
    !       2. Fluxes for central points (3rd order + limiter).
    !       3. Fluxes boundary point above (1st order).
    !       4. Fluxes boundary point below (1st order).
    !       5. Closure of 'X' if required
    !       6. Propagate.
    !     ------------------------------------------------------
    !
    !  9. Switches :
    !
    !     !/OMPH  Ading OMP directves for hybrid paralellization.
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable test output.
    !     !/T0  Test output input/output fields.
    !     !/T1  Test output fluxes.
    !     !/T2  Test output integration.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: MX, MY, NX, NY, INC, MAPACT(MY*MX),  &
         NACT, MAPBOU(MY*MX), NB0, NB1, NB2,  &
         NDSE, NDST
    REAL, INTENT(IN)        :: TRANS(MY*MX,-1:1)
    REAL, INTENT(INOUT)     :: CFLL(MY*(MX+1)), Q(1-MY:MY*(MX+2))
    LOGICAL, INTENT(IN)     :: CLOSE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IXY, IP, IXYC, IXYU, IXYD, IY, IX,   &
         IAD00, IAD02, IADN0, IADN1, IADN2,   &
         JN, JP
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T1
    INTEGER                 :: IX2, IY2
#endif
    REAL                    :: CFL, QB, DQ, DQNZ, QCN, QBN, QBR, CFAC
    REAL                    :: FLA(1-MY:MY*MX)
#ifdef W3_T0
    REAL                    :: QMAX
#endif
#ifdef W3_T1
    REAL                    :: QBO, QN
#endif
#ifdef W3_T2
    REAL                    :: QOLD
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3QCK3')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) mx, my, nx, ny, CLOSE, inc, nb0, nb1, nb2
#endif
    !
#ifdef W3_T0
    qmax   = 0.
    DO iy=1, ny
      DO ix=1, nx
        qmax   = max( qmax , q(iy+(ix-1)*my) )
      END DO
    END DO
    qmax   = max( 0.01*qmax , 1.e-10 )
#endif
    !
#ifdef W3_T0
    WRITE (ndst,9001) 'CFLL'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(100.*cfll(iy+(ix-1)*my)),ix=1,nx)
    END DO
    WRITE (ndst,9001) 'Q'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(q(iy+(ix-1)*my)/qmax),ix=1,nx)
    END DO
    WRITE (ndst,9001) 'MAPACT'
    WRITE (ndst,9003) (mapact(ixy),ixy=1,nact)
#endif
    !
    ! 1.  Initialize aux. array FLA and closure ------------------------- *
    !
    fla = 0.
    !
    IF ( CLOSE ) THEN
#ifdef W3_T
      WRITE (ndst,9005)
#endif
      iad00  = -my
      iad02  =  my
      iadn0  = iad00 + my*nx
      iadn1  =         my*nx
      iadn2  = iad02 + my*nx
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (IY)
#endif
      !
      DO iy=1, ny
        q(iy+iad00) = q(iy+iadn0) ! 1 ghost column to left
        q(iy+iadn1) = q(   iy   ) ! 1st ghost column to right
        q(iy+iadn2) = q(iy+iad02) ! 2nd ghost column to right
        cfll(iy+iadn1) = cfll(   iy   ) ! as for Q above, 1st to rt
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
    END IF
    !
    ! 2.  Fluxes for central points ------------------------------------- *
    !     ( 3rd order + limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9010)
    WRITE (ndst,9011) nb0, 'CENTRAL'
#endif
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (IP, IXY, CFL, IXYC, QB, IXYU, IXYD, &
#ifdef W3_T1
    !$OMP                      QBO, QN, IX, IY, IX2, IY2,              &
#endif
    !$OMP&                     DQ, DQNZ, QCN, QBN, QBR, CFAC )
#endif
    !
    DO ip=1, nb0
      !
      ixy    = mapbou(ip)
      cfl    = 0.5 * ( cfll(ixy) + cfll(ixy+inc) )
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,cfl) ) )
      qb     = 0.5 * ( (1.-cfl)*q(ixy+inc) + (1.+cfl)*q(ixy) )      &
           - (1.-cfl**2)/6. * (q(ixyc-inc)-2.*q(ixyc)+q(ixyc+inc))
      !
      ixyu   = ixyc - inc * int( sign(1.1,cfl) )
      ixyd   = 2*ixyc - ixyu
      dq     = q(ixyd) - q(ixyu)
      dqnz   = sign( max(1.e-15,abs(dq)) , dq )
      qcn    = ( q(ixyc) - q(ixyu) ) / dqnz
      qcn    = min( 1.1, max( -0.1 , qcn ) )
      !
#ifdef W3_T1
      qbo    = qb
#endif
      qbn    = max( (qb-q(ixyu))/dqnz , qcn )
      qbn    = min( qbn , 1. , qcn/max(1.e-10,abs(cfl)) )
      qbr    = q(ixyu) + qbn*dq
      cfac   = real( int( 2. * abs(qcn-0.5) ) )
      qb     = (1.-cfac)*qbr + cfac*q(ixyc)
      !
      fla(ixy) = cfl * qb
      !
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( qb, qbo, q(ixy-inc), q(   ixy   ),         &
           q(ixy+inc), q(ixy+2*inc) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9012) ip, ix, iy, ix2, iy2,               &
             cfl, cfll(ixy), cfll(ixy+inc),      &
             qbo*qn, qb*qn, q(ixy-inc)*qn, q(   ixy   )*qn,  &
             q(ixy+inc)*qn, q(ixy+2*inc)*qn
      END IF
#endif
      !
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    ! 3.  Fluxes for points with boundary above ------------------------- *
    !     ( 1st order without limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9011) nb1-nb0, 'BOUNDARY ABOVE'
#endif
    !
!!!/OMPH/!$OMP PARALLEL DO PRIVATE (IP, IXY, CFL, IXYC)
!!!
    DO ip=nb0+1, nb1
      ixy    = mapbou(ip)
      cfl    = cfll(ixy)
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,cfl) ) )
      fla(ixy) = cfl * q(ixyc)
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( q(ixy+inc), q(ixy) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9013) ip, ix, iy, ix2, iy2, cfl,          &
             cfll(ixy), q(ixyc)*qn, q(ixy)*qn, q(ixy+inc)*qn
      END IF
#endif
    END DO
!!!
!!!/OMPH/!$OMP END PARALLEL DO
    !
    ! 4.  Fluxes for points with boundary below ------------------------- *
    !     ( 1st order without limiter )
    !
#ifdef W3_T1
    WRITE (ndst,9011) nb2-nb1, 'BOUNDARY BELOW'
#endif
    !
!!!/OMPH/!$OMP PARALLEL DO PRIVATE (IP, IXY, CFL, IXYC)
!!!
    DO ip=nb1+1, nb2
      ixy    = mapbou(ip)
      cfl    = cfll(ixy+inc)
      ixyc   = ixy - inc * int( min( 0. , sign(1.1,cfl) ) )
      fla(ixy) = cfl * q(ixyc)
#ifdef W3_T1
      iy     = mod( ixy , my )
      ix     = 1 + ixy/my
      iy2    = mod( ixy+inc , my )
      ix2    = 1 + (ixy+inc)/my
      qn     = max( q(ixy+inc), q(ixy) )
      IF ( qn .GT. 1.e-10 ) THEN
        qn     = 1. /qn
        WRITE (ndst,9014) ip, ix, iy, ix2, iy2, cfl, cfll(ixy+inc), &
             q(ixyc)*qn, q(ixy)*qn, q(ixy+inc)*qn
      END IF
#endif
    END DO
    !
!!!/OMPH/!$OMP END PARALLEL DO
    !
    ! 5.  Global closure ----------------------------------------------- *
    !
    IF ( CLOSE ) THEN
#ifdef W3_T
      WRITE (ndst,9015)
#endif
      DO iy=1, ny
        fla(iy+iad00) = fla(iy+iadn0)
      END DO
    END IF
    !
    ! 6.  Propagation -------------------------------------------------- *
    !
#ifdef W3_T2
    WRITE (ndst,9020)
#endif
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (IP, IXY, JN, JP )
#endif
    !
    DO ip=1, nact
      !
      ixy    = mapact(ip)
      IF ( fla(ixy-inc) .GT. 0. ) THEN
        jn    = -1
      ELSE
        jn    =  0
      END IF
      IF ( fla(ixy    ) .LT. 0. ) THEN
        jp    =  1
      ELSE
        jp    =  0
      END IF
      !
#ifdef W3_T2
      qold   = q(ixy)
#endif
      q(ixy) = max( 0. , q(ixy) + trans(ixy,jn) * fla(ixy-inc)     &
           - trans(ixy,jp) * fla(ixy) )
 
#ifdef W3_T2
      IF ( qold + q(ixy) .GT. 1.e-10 )                          &
           WRITE (ndst,9021) ip, ixy, qold, q(ixy),             &
           fla(ixy-inc), fla(ixy)
#endif
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
#ifdef W3_T0
    WRITE (ndst,9001) 'Q'
    DO iy=ny,1,-1
      WRITE (ndst,9002) (nint(q(iy+(ix-1)*my)/qmax),ix=1,nx)
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3QCK3 : ARRAY DIMENSIONS  :',2i6/           &
         '               USED              :',2i6/           &
         '               CLOSE, INC        :',l6,i6/         &
         '               NB0, NB1, NB2     :',3i6)
#endif
#ifdef W3_T0
9001 FORMAT ( ' TEST W3QCK3 : DUMP ARRAY ',a,' :')
9002 FORMAT ( 1x,43i3)
9003 FORMAT ( 1x,21i6)
#endif
#ifdef W3_T
9005 FORMAT (' TEST W3QCK3 : GLOBAL CLOSURE (1)')
#endif
    !
#ifdef W3_T1
9010 FORMAT (' TEST W3QCK3 : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
         ' Q (b,b,i-1,i,i+1,i+2)')
9011 FORMAT (' TEST W3QCK3 :',i6,' POINTS OF TYPE ',a)
9012 FORMAT (10x,i6,4i4,1x,3f6.2,1x,f7.2,f6.2,1x,4f6.2)
9013 FORMAT (10x,i6,4i4,1x,f6.2,f6.2,'  --- ',1x,f7.2,1x,'  --- ',&
         2f6.2,'  --- ')
9014 FORMAT (10x,i6,4i4,1x,f6.2,'  --- ',f6.2,1x,f7.2,1x,'  --- ',&
         2f6.2,'  --- ')
#endif
#ifdef W3_T
9015 FORMAT (' TEST W3QCK3 : GLOBAL CLOSURE (2)')
#endif
    !
#ifdef W3_T2
9020 FORMAT (' TEST W3QCK3 : IP, IXY, 2Q, 2FL')
9021 FORMAT ('            ',2i6,2(1x,2e11.3))
#endif
    !/
    !/ End of W3QCK3 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3qck3
  !/
  !/ End of module W3UQCKMD -------------------------------------------- /
  !/
END MODULE w3uqckmd