w3uno2md.ftn

#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      MODULE W3UNO2MD
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           MetOffice |
!/                  |           Jian-Guo Li             |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         01-Jul-2013 |
!/                  +-----------------------------------+
!/
!/    Adapted from   WAVEWATCH-III  W3UQCKMD 
!/            for    UNO2 advection scheme. 
!/
!/    18-Mar-2008 : Origination.                        ( version 3.14 )
!/    ..-...-...  : .....                               ( version 3.14 )
!/    19-Mar-2008 : last modified by Jian-Guo           ( version 3.14 )
!/    01-Jul-2013 : Put in NCEP branch (Tolman).        ( version 4.12 )
!/    08-Jan-2018 : Added OMPH switches in W3UNO2.      ( version 6.02 )
!/
!  1. Purpose :
!
!     Portable UNO2 scheme on irregular grid.
!
!  2. Variables and types :
!
!     None.
!
!  3. Subroutines and functions :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!      W3UNO2    Subr. Public   UNO2 scheme for irregular grid.
!      W3UNO2r   Subr. Public   UNO2 scheme reduced to regular grid.
!      W3UNO2s   Subr. Public   UNO2 regular grid with subgrid obstruction.
!     ----------------------------------------------------------------
!
!  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 :
!
!/ ------------------------------------------------------------------- /
!/S      USE W3SERVMD, ONLY: STRACE
!/
      CONTAINS
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3UNO2 (MX, MY, NX, NY, VELO, DT, DX1, DX2, Q,BCLOSE,&
                        INC,  MAPACT, NACT, MAPBOU, NB0, NB1, NB2,    &
                        NDSE, NDST )
!/
!
!     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)
!       BCLOSE  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 :
!
!       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 :
!
!     !/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)     :: BCLOSE
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: IXY, IP, IXYC, IXYU, IXYD, IY, IX,   &
                                 IAD00, IAD02, IADN0, IADN1, IADN2
!/S      INTEGER, SAVE           :: IENT
!/T1      INTEGER                 :: IX2, IY2
      REAL                    :: CFL, VEL, QB, DQ, DQNZ, QCN, QBN,    &
                                 QBR, CFAC, FLA(1-MY:MY*MX)
!/T0      REAL                    :: QMAX
!/T1      REAL                    :: QBO, QN, XCFL
!/T2      REAL                    :: QOLD
!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'W3UNO2')
!
!/T      WRITE (NDST,9000) MX, MY, NX, NY, DT, BCLOSE, INC, NB0, NB1, NB2
!
!/T0      QMAX   = 0.
!/T0      DO IY=1, NY
!/T0        DO IX=1, NX
!/T0          QMAX   = MAX ( QMAX , Q(IY+(IX-1)*MY) )
!/T0          END DO
!/T0        END DO
!/T0      QMAX   = MAX ( 0.01*QMAX , 1.E-10 )
!
!/T0      WRITE (NDST,9001) 'VELO'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(100.*VELO(IY+(IX-1)*MY)           &
!/T0                                  *DT/DX1(IY+(IX-1)*MY)),IX=1,NX)
!/T0        END DO
!/T0      WRITE (NDST,9001) 'Q'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(Q(IY+(IX-1)*MY)/QMAX),IX=1,NX)
!/T0        END DO
!/T0      WRITE (NDST,9001) 'MAPACT'
!/T0      WRITE (NDST,9003) (MAPACT(IXY),IXY=1,NACT)
!
! 1.  Initialize aux. array FLA and closure ------------------------- *
!
      FLA = 0.
!
      IF ( BCLOSE ) THEN
!/T          WRITE (NDST,9005)
          IAD00  = -MY
          IAD02  =  MY
          IADN0  = IAD00 + MY*NX
          IADN1  =         MY*NX
          IADN2  = IAD02 + MY*NX
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
          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 ------------------------------------- *
!     ( 2rd order UNO2 scheme )
!
!/T1      WRITE (NDST,9010)
!/T1      WRITE (NDST,9011) NB0, 'CENTRAL'
!
      DO IP=1, NB0
!
        IXY  = MAPBOU(IP)
        VEL  = 0.5 * ( VELO(IXY) + VELO(IXY+INC) )
!  Assuming velocity is at cell centre, so face velocity is an average.
        CFL  = DT *  VEL 
!  Courant number without gradient distance (between IXY and IXY+INC cells)
        IXYC = IXY - INC * INT( MIN ( 0. , SIGN(1.1,CFL) ) )
!  Central cell index, depending on flow direction.
!          IXY for positive CFL, IXY+INC for negative CFL
!  Upstream and downstream cell numbers
        IXYD = IXYC + INC * INT ( SIGN (1.1,CFL) )
!  Minimum gradient is derived from the two sides of the central cell
!
        QB   = Q(IXYC)+SIGN(0.5, Q(IXYD)-Q(IXYC))*(DX1(IXYC)-ABS(CFL)) &
                      *MIN(ABS(Q(IXYC+INC)-Q(IXYC))/DX2(IXYC),         &
                           ABS(Q(IXYC)-Q(IXYC-INC))/DX2(IXYC-INC) )
!
!/T1        QBO    = QB
!
        FLA(IXY) = CFL * QB
!
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( QB, QBO, Q(IXY-INC), Q(   IXY   ),         &
!/T1                                Q(IXY+INC), Q(IXY+2*INC) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9012) IP, IX, IY, IX2, IY2,               &
!/T1                      CFL, DT*VELO(IXY)/DX1(IXY),                 &
!/T1                           DT*VELO(IXY+INC)/DX1(IXY+INC),         &
!/T1                  QBO*QN, QB*QN, Q(IXY-INC)*QN, Q(   IXY   )*QN,  &
!/T1                                 Q(IXY+INC)*QN, Q(IXY+2*INC)*QN
!/T1          END IF
!
      END DO
!
! 3.  Fluxes for points with boundary above ------------------------- *
!     ( 1st order without limiter )
!
!/T1      WRITE (NDST,9011) NB1-NB0, 'BOUNDARY ABOVE'
!
      DO IP=NB0+1, NB1
        IXY    = MAPBOU(IP)
        VEL    = VELO(IXY)
        IXYC   = IXY - INC * INT( MIN ( 0. , SIGN(1.1,VEL) ) )
        FLA(IXY) = VEL * DT * Q(IXYC)
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( Q(IXY+INC), Q(IXY) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9013) IP, IX, IY, IX2, IY2, XCFL,             &
!/T1                              DT*VELO(IXY)/DX2(IXY),                  &
!/T1                              Q(IXYC)*QN, Q(IXY)*QN, Q(IXY+INC)*QN
!/T1          END IF
      END DO
!
! 4.  Fluxes for points with boundary below ------------------------- *
!     ( 1st order without limiter )
!
!/T1      WRITE (NDST,9011) NB2-NB1, 'BOUNDARY BELOW'
!
      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 * DT * Q(IXYC)
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( Q(IXY+INC), Q(IXY) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9014) IP, IX, IY, IX2, IY2, XCFL,         &
!/T1                              DT*VELO(IXY+INC)/DX2(IXY),          &
!/T1                              Q(IXYC)*QN, Q(IXY)*QN, Q(IXY+INC)*QN
!/T1          END IF
      END DO
!
! 5.  Global closure ----------------------------------------------- *
!
      IF ( BCLOSE ) THEN
!/T          WRITE (NDST,9015)
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
          DO IY=1, NY
            FLA (IY+IAD00) = FLA (IY+IADN0)
          END DO
      END IF
!
! 6.  Propagation -------------------------------------------------- *
!
!/T2      WRITE (NDST,9020)
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
      DO IP=1, NACT
        IXY    = MAPACT(IP)
!/T2        QOLD   = Q(IXY)
! Li    Update transported quantity with fluxes
        Q(IXY) = MAX( 0., Q(IXY)+( FLA(IXY-INC)-FLA(IXY) )/DX1(IXY) )
! Li    This positive filter is not necessary for UNO2 scheme but kept here.
!/T2        IF ( QOLD + Q(IXY) .GT. 1.E-10 )                          &
!/T2             WRITE (NDST,9021) IP, IXY, QOLD, Q(IXY),             &
!/T2                               DT*FLA(IXY-INC)/DX1(IXY),          &
!/T2                               DT*FLA(IXY)/DX1(IXY)
      END DO
!
!/T0      WRITE (NDST,9001) 'Q'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(Q(IY+(IX-1)*MY)/QMAX),IX=1,NX)
!/T0        END DO
!
      RETURN
!
! Formats
!
!/T 9000 FORMAT ( ' TEST W3UNO2 : ARRAY DIMENSIONS  :',2I6/           &
!/T               '               USED              :',2I6/           &
!/T               '               TIME STEP         :',F8.1/          &
!/T               '              BCLOSE, INC        :',L6,I6/         &
!/T               '               NB0, NB1, NB2     :',3I6)
!/T0 9001 FORMAT ( ' TEST W3UNO2 : DUMP ARRAY ',A,' :')
!/T0 9002 FORMAT ( 1X,43I3)
!/T0 9003 FORMAT ( 1X,21I6)
!/T 9005 FORMAT (' TEST W3UNO2 : GLOBAL CLOSURE (1)')
!
!/T1 9010 FORMAT (' TEST W3UNO2 : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
!/T1              ' Q (b,b,i-1,i,i+1,i+2)')
!/T1 9011 FORMAT (' TEST W3UNO2 :',I6,' POINTS OF TYPE ',A)
!/T1 9012 FORMAT (10X,I6,4I4,1X,3F6.2,1X,F7.2,F6.2,1X,4F6.2)
!/T1 9013 FORMAT (10X,I6,4I4,1X,F6.2,F6.2,'  --- ',1X,F7.2,1X,'  --- ',&
!/T1              2F6.2,'  --- ')
!/T1 9014 FORMAT (10X,I6,4I4,1X,F6.2,'  --- ',F6.2,1X,F7.2,1X,'  --- ',&
!/T1              2F6.2,'  --- ')
!/T 9015 FORMAT (' TEST W3UNO2 : GLOBAL CLOSURE (2)')
!
!/T2 9020 FORMAT (' TEST W3UNO2 : IP, IXY, 2Q, 2FL')
!/T2 9021 FORMAT ('            ',2I6,2(1X,2E11.3))
!/
!/ End of W3UNO2 ----------------------------------------------------- /
!/
      END SUBROUTINE W3UNO2
!/
!/
      SUBROUTINE W3UNO2r (MX, MY, NX, NY, CFLL, Q, BCLOSE, INC,       &
                         MAPACT, NACT, MAPBOU, NB0, NB1, NB2,         &
                         NDSE, NDST )
!/
!/    Adapted from W3QCK1 for UNO2 regular grid scheme.
!/            First created:    19 Mar 2008   Jian-Guo Li
!/            Last modified:     8 Jan 2018   Jian-Guo Li
!/
!  1. Purpose :
!
!     Preform one-dimensional propagation in a two-dimensional space
!     with irregular boundaries and regular grid.
!
!  2. Method :
!
!     UNO2 regular grid scheme
!
!  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)
!       BCLOSE  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 :
!
!     !/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)     :: BCLOSE
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: IXY, IP, IXYC, IXYU, IXYD, IY, IX,   &
                                 IAD00, IAD02, IADN0, IADN1, IADN2
!/S      INTEGER, SAVE           :: IENT = 0
!/T1      INTEGER                 :: IX2, IY2
      REAL                    :: CFL, QB, DQ, DQNZ, QCN, QBN, QBR, CFAC
      REAL                    :: FLA(1-MY:MY*MX)
!/T0      REAL                    :: QMAX
!/T1      REAL                    :: QBO, QN
!/T2      REAL                    :: QOLD
!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'W3UNO2r')
!
!/T      WRITE (NDST,9000) MX, MY, NX, NY, BCLOSE, INC, NB0, NB1, NB2
!
!/T0      QMAX   = 0.
!/T0      DO IY=1, NY
!/T0        DO IX=1, NX
!/T0          QMAX   = MAX ( QMAX , Q(IY+(IX-1)*MY) )
!/T0          END DO
!/T0        END DO
!/T0      QMAX   = MAX ( 0.01*QMAX , 1.E-10 )
!
!/T0      WRITE (NDST,9001) 'CFLL'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(100.*CFLL(IY+(IX-1)*MY)),IX=1,NX)
!/T0        END DO
!/T0      WRITE (NDST,9001) 'Q'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(Q(IY+(IX-1)*MY)/QMAX),IX=1,NX)
!/T0        END DO
!/T0      WRITE (NDST,9001) 'MAPACT'
!/T0      WRITE (NDST,9003) (MAPACT(IXY),IXY=1,NACT)
!
! 1.  Initialize aux. array FLA and closure ------------------------- *
!
      FLA = 0.
!
      IF ( BCLOSE ) THEN
!/T          WRITE (NDST,9005)
          IAD00  = -MY
          IAD02  =  MY
          IADN0  = IAD00 + MY*NX
          IADN1  =         MY*NX
          IADN2  = IAD02 + MY*NX
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
          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 )
!
!/T1      WRITE (NDST,9010)
!/T1      WRITE (NDST,9011) NB0, 'CENTRAL'
!
      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) ) )
        IXYD = IXYC + INC * INT( SIGN (1.1,CFL) )
        QB   = Q(IXYC)+SIGN(0.5, Q(IXYD)-Q(IXYC))*(1.0-ABS(CFL)) &
                      *MIN(ABS(Q(IXYC+INC)-Q(IXYC)),             &
                           ABS(Q(IXYC)-Q(IXYC-INC)) )
!/T1        QBO    = QB
!
        FLA(IXY) = CFL * QB
!
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( QB, QBO, Q(IXY-INC), Q(   IXY   ),         &
!/T1                                Q(IXY+INC), Q(IXY+2*INC) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9012) IP, IX, IY, IX2, IY2,               &
!/T1                              CFL, CFLL(IXY), CFLL(IXY+INC),      &
!/T1                   QBO*QN, QB*QN, Q(IXY-INC)*QN, Q(   IXY   )*QN, &
!/T1                                  Q(IXY+INC)*QN, Q(IXY+2*INC)*QN
!/T1          END IF
!
        END DO
!
! 3.  Fluxes for points with boundary above ------------------------- *
!     ( 1st order without limiter )
!
!/T1      WRITE (NDST,9011) NB1-NB0, 'BOUNDARY ABOVE'
!
      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)
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( Q(IXY+INC), Q(IXY) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9013) IP, IX, IY, IX2, IY2, CFL,          &
!/T1                   CFLL(IXY), Q(IXYC)*QN, Q(IXY)*QN, Q(IXY+INC)*QN
!/T1          END IF
        END DO
!
! 4.  Fluxes for points with boundary below ------------------------- *
!     ( 1st order without limiter )
!
!/T1      WRITE (NDST,9011) NB2-NB1, 'BOUNDARY BELOW'
!
      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)
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( Q(IXY+INC), Q(IXY) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9014) IP, IX, IY, IX2, IY2, CFL,         &
!/T1               CFLL(IXY+INC), Q(IXYC)*QN, Q(IXY)*QN, Q(IXY+INC)*QN
!/T1          END IF
        END DO
!
! 5.  Global closure ----------------------------------------------- *
!
      IF ( BCLOSE ) THEN
!/T          WRITE (NDST,9015)
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
          DO IY=1, NY
            FLA (IY+IAD00) = FLA (IY+IADN0)
            END DO
        END IF
!
! 6.  Propagation -------------------------------------------------- *
!
!/T2      WRITE (NDST,9020)
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
      DO IP=1, NACT
        IXY    = MAPACT(IP)
!/T2        QOLD   = Q(IXY)
        Q(IXY) = MAX ( 0. , Q(IXY) + FLA(IXY-INC) - FLA(IXY) )
!/T2        IF ( QOLD + Q(IXY) .GT. 1.E-10 )                          &
!/T2             WRITE (NDST,9021) IP, IXY, QOLD, Q(IXY),             &
!/T2                               FLA(IXY-INC), FLA(IXY)
        END DO
!
!/T0      WRITE (NDST,9001) 'Q'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(Q(IY+(IX-1)*MY)/QMAX),IX=1,NX)
!/T0        END DO
!
      RETURN
!
! Formats
!
!/T 9000 FORMAT ( ' TEST W3UNO2r : ARRAY DIMENSIONS  :',2I6/           &
!/T               '               USED              :',2I6/           &
!/T               '              BCLOSE, INC        :',L6,I6/         &
!/T               '               NB0, NB1, NB2     :',3I6)
!/T0 9001 FORMAT ( ' TEST W3UNO2r : DUMP ARRAY ',A,' :')
!/T0 9002 FORMAT ( 1X,43I3)
!/T0 9003 FORMAT ( 1X,21I6)
!/T 9005 FORMAT (' TEST W3UNO2r : GLOBAL CLOSURE (1)')
!
!/T1 9010 FORMAT (' TEST W3UNO2r : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
!/T1              ' Q (b,b,i-1,i,i+1,i+2)')
!/T1 9011 FORMAT (' TEST W3UNO2r :',I6,' POINTS OF TYPE ',A)
!/T1 9012 FORMAT (10X,I6,4I4,1X,3F6.2,1X,F7.2,F6.2,1X,4F6.2)
!/T1 9013 FORMAT (10X,I6,4I4,1X,F6.2,F6.2,'  --- ',1X,F7.2,1X,'  --- ',&
!/T1              2F6.2,'  --- ')
!/T1 9014 FORMAT (10X,I6,4I4,1X,F6.2,'  --- ',F6.2,1X,F7.2,1X,'  --- ',&
!/T1              2F6.2,'  --- ')
!/T 9015 FORMAT (' TEST W3UNO2r : GLOBAL CLOSURE (2)')
!
!/T2 9020 FORMAT (' TEST W3UNO2r : IP, IXY, 2Q, 2FL')
!/T2 9021 FORMAT ('            ',2I6,2(1X,2E11.3))
!/
!/ End of W3UNO2r ---------------------------------------------------- /
!/
      END SUBROUTINE W3UNO2r
!/
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3UNO2s (MX, MY, NX, NY, CFLL, TRANS, Q, BCLOSE,     &
                         INC, MAPACT, NACT, MAPBOU, NB0, NB1, NB2,    &
                         NDSE, NDST )
!/
!/
!/    Adapted from W3QCK3 for UNO2 regular grid scheme with
!/            subgrid obstruction.
!/            First created:    19 Mar 2008   Jian-Guo Li
!/            Last modified:     8 Jan 2018   Jian-Guo Li
!/
!  1. Purpose :
!
!     Like W3UNO2r 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)
!       BCLOSE  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)     :: BCLOSE
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: IXY, IP, IXYC, IXYU, IXYD, IY, IX,   &
                                 IAD00, IAD02, IADN0, IADN1, IADN2,   &
                                 JN, JP
!/S      INTEGER, SAVE           :: IENT = 0
!/T1      INTEGER                 :: IX2, IY2
      REAL                    :: CFL, QB, DQ, DQNZ, QCN, QBN, QBR, CFAC
      REAL                    :: FLA(1-MY:MY*MX)
!/T0      REAL                    :: QMAX
!/T1      REAL                    :: QBO, QN
!/T2      REAL                    :: QOLD
!/
!/ ------------------------------------------------------------------- /
!/
!/S      CALL STRACE (IENT, 'W3UNO2s')
!
!/T      WRITE (NDST,9000) MX, MY, NX, NY, BCLOSE, INC, NB0, NB1, NB2
!
!/T0      QMAX   = 0.
!/T0      DO IY=1, NY
!/T0        DO IX=1, NX
!/T0          QMAX   = MAX ( QMAX , Q(IY+(IX-1)*MY) )
!/T0          END DO
!/T0        END DO
!/T0      QMAX   = MAX ( 0.01*QMAX , 1.E-10 )
!
!/T0      WRITE (NDST,9001) 'CFLL'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(100.*CFLL(IY+(IX-1)*MY)),IX=1,NX)
!/T0        END DO
!/T0      WRITE (NDST,9001) 'Q'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(Q(IY+(IX-1)*MY)/QMAX),IX=1,NX)
!/T0        END DO
!/T0      WRITE (NDST,9001) 'MAPACT'
!/T0      WRITE (NDST,9003) (MAPACT(IXY),IXY=1,NACT)
!
! 1.  Initialize aux. array FLA and closure ------------------------- *
!
      FLA = 0.
!
      IF ( BCLOSE ) THEN
!/T          WRITE (NDST,9005)
          IAD00  = -MY
          IAD02  =  MY
          IADN0  = IAD00 + MY*NX
          IADN1  =         MY*NX
          IADN2  = IAD02 + MY*NX
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (IY)
!
          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
!
!/OMPH/!$OMP END PARALLEL DO
!
        END IF
!
! 2.  Fluxes for central points ------------------------------------- *
!     ( 3rd order + limiter )
!
!/T1      WRITE (NDST,9010)
!/T1      WRITE (NDST,9011) NB0, 'CENTRAL'
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (IP, IXY, CFL, IXYC, IXYD, QB)  
!
      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) ) )
        IXYD = IXYC + INC * INT( SIGN (1.1,CFL) )
        QB   = Q(IXYC)+SIGN(0.5, Q(IXYD)-Q(IXYC))*(1.0-ABS(CFL)) &
                      *MIN(ABS(Q(IXYC+INC)-Q(IXYC)),             &
                           ABS(Q(IXYC)-Q(IXYC-INC)) )
!
!/T1        QBO    = QB
!
        FLA(IXY) = CFL * QB
!
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( QB, QBO, Q(IXY-INC), Q(   IXY   ),         &
!/T1                                Q(IXY+INC), Q(IXY+2*INC) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9012) IP, IX, IY, IX2, IY2,               &
!/T1                              CFL, CFLL(IXY), CFLL(IXY+INC),      &
!/T1                  QBO*QN, QB*QN, Q(IXY-INC)*QN, Q(   IXY   )*QN,  &
!/T1                                 Q(IXY+INC)*QN, Q(IXY+2*INC)*QN
!/T1          END IF
!
        END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
! 3.  Fluxes for points with boundary above ------------------------- *
!     ( 1st order without limiter )
!
!/T1      WRITE (NDST,9011) NB1-NB0, 'BOUNDARY ABOVE'
!
      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)
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( Q(IXY+INC), Q(IXY) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9013) IP, IX, IY, IX2, IY2, CFL,          &
!/T1                   CFLL(IXY), Q(IXYC)*QN, Q(IXY)*QN, Q(IXY+INC)*QN
!/T1          END IF
        END DO
!
! 4.  Fluxes for points with boundary below ------------------------- *
!     ( 1st order without limiter )
!
!/T1      WRITE (NDST,9011) NB2-NB1, 'BOUNDARY BELOW'
!
      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)
!/T1        IY     = MOD ( IXY , MY )
!/T1        IX     = 1 + IXY/MY
!/T1        IY2    = MOD ( IXY+INC , MY )
!/T1        IX2    = 1 + (IXY+INC)/MY
!/T1        QN     = MAX ( Q(IXY+INC), Q(IXY) )
!/T1        IF ( QN .GT. 1.E-10 ) THEN
!/T1            QN     = 1. /QN
!/T1            WRITE (NDST,9014) IP, IX, IY, IX2, IY2, CFL, CFLL(IXY+INC), &
!/T1                              Q(IXYC)*QN, Q(IXY)*QN, Q(IXY+INC)*QN
!/T1          END IF
        END DO
!
! 5.  Global closure ----------------------------------------------- *
!
      IF ( BCLOSE ) THEN
!/T          WRITE (NDST,9015)
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
          DO IY=1, NY
            FLA (IY+IAD00) = FLA (IY+IADN0)
            END DO
        END IF
!
! 6.  Propagation -------------------------------------------------- *
!
!/T2      WRITE (NDST,9020)
!/C90/!DIR$ IVDEP
!/NEC/!CDIR NODEP
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (IP, IXY, JN, JP )
!
      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
!
!/T2        QOLD   = Q(IXY)
        Q(IXY) = MAX ( 0. , Q(IXY) + TRANS(IXY,JN) * FLA(IXY-INC)     &
                                   - TRANS(IXY,JP) * FLA(IXY) )

!/T2        IF ( QOLD + Q(IXY) .GT. 1.E-10 )                          &
!/T2             WRITE (NDST,9021) IP, IXY, QOLD, Q(IXY),             &
!/T2                               FLA(IXY-INC), FLA(IXY)
        END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
!/T0      WRITE (NDST,9001) 'Q'
!/T0      DO IY=NY,1,-1
!/T0        WRITE (NDST,9002) (NINT(Q(IY+(IX-1)*MY)/QMAX),IX=1,NX)
!/T0        END DO
!
      RETURN
!
! Formats
!
!/T 9000 FORMAT ( ' TEST W3UNO2s : ARRAY DIMENSIONS  :',2I6/           &
!/T               '               USED              :',2I6/           &
!/T               '               BCLOSE, INC       :',L6,I6/         &
!/T               '               NB0, NB1, NB2     :',3I6)
!/T0 9001 FORMAT ( ' TEST W3UNO2s : DUMP ARRAY ',A,' :')
!/T0 9002 FORMAT ( 1X,43I3)
!/T0 9003 FORMAT ( 1X,21I6)
!/T 9005 FORMAT (' TEST W3UNO2s : GLOBAL CLOSURE (1)')
!
!/T1 9010 FORMAT (' TEST W3UNO2s : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
!/T1              ' Q (b,b,i-1,i,i+1,i+2)')
!/T1 9011 FORMAT (' TEST W3UNO2s :',I6,' POINTS OF TYPE ',A)
!/T1 9012 FORMAT (10X,I6,4I4,1X,3F6.2,1X,F7.2,F6.2,1X,4F6.2)
!/T1 9013 FORMAT (10X,I6,4I4,1X,F6.2,F6.2,'  --- ',1X,F7.2,1X,'  --- ',&
!/T1              2F6.2,'  --- ')
!/T1 9014 FORMAT (10X,I6,4I4,1X,F6.2,'  --- ',F6.2,1X,F7.2,1X,'  --- ',&
!/T1              2F6.2,'  --- ')
!/T 9015 FORMAT (' TEST W3UNO2s : GLOBAL CLOSURE (2)')
!
!/T2 9020 FORMAT (' TEST W3UNO2s : IP, IXY, 2Q, 2FL')
!/T2 9021 FORMAT ('            ',2I6,2(1X,2E11.3))
!/
!/ End of W3UNO2s ---------------------------------------------------- /
!/
      END SUBROUTINE W3UNO2s
!/
!/ End of module W3UNO2MD -------------------------------------------- /
!/
      END MODULE W3UNO2MD
!/