 
#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 :
!
!       !/C90   Cray FORTRAN 90 compiler directives.
!       !/NEC   NEC SXF90 compiler directives.
!
!       !/OMPH  Ading OMP directves for hybrid paralellization.
!
!       !/S     Enable subroutine tracing.
!       !/Tn    Enable test output.
!
!  7. Source code :
!
!/ ------------------------------------------------------------------- /
!/
      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
      REAL                    :: CFL, QB, DQ, DQNZ, QCN, QBN, QBR, CFAC
      REAL                    :: FLA(1-MY:MY*MX)
!/
!/ ------------------------------------------------------------------- /
!/
!
! 1.  Initialize aux. array FLA and closure ------------------------- *
!
      FLA = 0.
!
      IF ( CLOSE ) THEN
          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 )
!
      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 ) )
!
        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
!
        END DO
!
! 3.  Fluxes for points with boundary above ------------------------- *
!     ( 1st order without limiter )
!
      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)
        END DO
!
! 4.  Fluxes for points with boundary below ------------------------- *
!     ( 1st order without limiter )
!
      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)
        END DO
!
! 5.  Global closure ----------------------------------------------- *
!
      IF ( CLOSE ) THEN
          DO IY=1, NY
            FLA (IY+IAD00) = FLA (IY+IADN0)
            END DO
        END IF
!
! 6.  Propagation -------------------------------------------------- *
!
      DO IP=1, NACT
        IXY    = MAPACT(IP)
        Q(IXY) = MAX ( 0. , Q(IXY) + FLA(IXY-INC) - FLA(IXY) )
        END DO
!
      RETURN
!
! Formats
!
!/
!/ 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
      REAL                    :: CFL, VEL, QB, DQ, DQNZ, QCN, QBN,    &
                                 QBR, CFAC, FLA(1-MY:MY*MX)
!/
!/ ------------------------------------------------------------------- /
!/
!
! 1.  Initialize aux. array FLA and closure ------------------------- *
!
      FLA = 0.
!
      IF ( CLOSE ) THEN
          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 )
!
      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 ) )
!
        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
!
        END DO
!
! 3.  Fluxes for points with boundary above ------------------------- *
!     ( 1st order without limiter )
!
      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)
        END DO
!
! 4.  Fluxes for points with boundary below ------------------------- *
!     ( 1st order without limiter )
!
      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)
        END DO
!
! 5.  Global closure ----------------------------------------------- *
!
      IF ( CLOSE ) THEN
          DO IY=1, NY
            FLA (IY+IAD00) = FLA (IY+IADN0)
            END DO
        END IF
!
! 6.  Propagation -------------------------------------------------- *
!
      DO IP=1, NACT
        IXY    = MAPACT(IP)
        Q(IXY) = MAX ( 0. , Q(IXY) + DT/DX1(IXY) *                    &
                              (FLA(IXY-INC)-FLA(IXY)) )
        END DO
!
      RETURN
!
! Formats
!
!/
!/ 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
      REAL                    :: CFL, QB, DQ, DQNZ, QCN, QBN, QBR, CFAC
      REAL                    :: FLA(1-MY:MY*MX)
!/
!/ ------------------------------------------------------------------- /
!/
!
! 1.  Initialize aux. array FLA and closure ------------------------- *
!
      FLA = 0.
!
      IF ( CLOSE ) THEN
          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) ! 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
!
        END IF
!
! 2.  Fluxes for central points ------------------------------------- *
!     ( 3rd order + limiter )
!
      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 ) )
!
        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
!
        END DO
!
! 3.  Fluxes for points with boundary above ------------------------- *
!     ( 1st order without limiter )
!
!!!/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)
        END DO
!!!
!!!/OMPH/!$OMP END PARALLEL DO
!
! 4.  Fluxes for points with boundary below ------------------------- *
!     ( 1st order without limiter )
!
!!!/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)
        END DO
!
!!!/OMPH/!$OMP END PARALLEL DO
!
! 5.  Global closure ----------------------------------------------- *
!
      IF ( CLOSE ) THEN
          DO IY=1, NY
            FLA (IY+IAD00) = FLA (IY+IADN0)
            END DO
        END IF
!
! 6.  Propagation -------------------------------------------------- *
!
      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
!
        Q(IXY) = MAX ( 0. , Q(IXY) + TRANS(IXY,JN) * FLA(IXY-INC)     &
                                   - TRANS(IXY,JP) * FLA(IXY) )
 
        END DO
!
      RETURN
!
! Formats
!
!/
!/ End of W3QCK3 ----------------------------------------------------- /
!/
      END SUBROUTINE W3QCK3
!/
!/ End of module W3UQCKMD -------------------------------------------- /
!/
      END MODULE W3UQCKMD
