w3uno2md.F90

Go to the documentation of this file.

 
 
#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 :
  !
  !       !/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 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
#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, 'W3UNO2')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) mx, my, nx, ny, dt, bclose, 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 ( bclose ) 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 ------------------------------------- *
    !     ( 2rd order UNO2 scheme )
    !
#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) )
      !  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) )
      !
#ifdef W3_T1
      qbo    = qb
#endif
      !
      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, 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 * dt * 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 * dt * 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 ( bclose ) 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
      ! 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.
#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 W3UNO2 : ARRAY DIMENSIONS  :',2i6/           &
         '               USED              :',2i6/           &
         '               TIME STEP         :',f8.1/          &
         '              BCLOSE, INC        :',l6,i6/         &
         '               NB0, NB1, NB2     :',3i6)
#endif
#ifdef W3_T0
9001 FORMAT ( ' TEST W3UNO2 : DUMP ARRAY ',a,' :')
9002 FORMAT ( 1x,43i3)
9003 FORMAT ( 1x,21i6)
#endif
#ifdef W3_T
9005 FORMAT (' TEST W3UNO2 : GLOBAL CLOSURE (1)')
#endif
    !
#ifdef W3_T1
9010 FORMAT (' TEST W3UNO2 : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
         ' Q (b,b,i-1,i,i+1,i+2)')
9011 FORMAT (' TEST W3UNO2 :',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 W3UNO2 : GLOBAL CLOSURE (2)')
#endif
    !
#ifdef W3_T2
9020 FORMAT (' TEST W3UNO2 : IP, IXY, 2Q, 2FL')
9021 FORMAT ('            ',2i6,2(1x,2e11.3))
#endif
  END SUBROUTINE w3uno2
  !/
  !/ End of 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
#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, 'W3UNO2r')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) mx, my, nx, ny, bclose, 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 ( bclose ) 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) ) )
      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)) )
#ifdef W3_T1
      qbo    = qb
#endif
      !
      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 ( bclose ) 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 W3UNO2r : ARRAY DIMENSIONS  :',2i6/           &
         '               USED              :',2i6/           &
         '              BCLOSE, INC        :',l6,i6/         &
         '               NB0, NB1, NB2     :',3i6)
#endif
#ifdef W3_T0
9001 FORMAT ( ' TEST W3UNO2r : DUMP ARRAY ',a,' :')
9002 FORMAT ( 1x,43i3)
9003 FORMAT ( 1x,21i6)
#endif
#ifdef W3_T
9005 FORMAT (' TEST W3UNO2r : GLOBAL CLOSURE (1)')
#endif
    !
#ifdef W3_T1
9010 FORMAT (' TEST W3UNO2r : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
         ' Q (b,b,i-1,i,i+1,i+2)')
9011 FORMAT (' TEST W3UNO2r :',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 W3UNO2r : GLOBAL CLOSURE (2)')
#endif
    !
#ifdef W3_T2
9020 FORMAT (' TEST W3UNO2r : IP, IXY, 2Q, 2FL')
9021 FORMAT ('            ',2i6,2(1x,2e11.3))
#endif
    !/
  END SUBROUTINE w3uno2r
  !/
  !/ End of 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
#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, 'W3UNO2s')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) mx, my, nx, ny, bclose, 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 ( bclose ) 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)
        q(iy+iadn1) = q(   iy   )
        q(iy+iadn2) = q(iy+iad02)
        cfll(iy+iadn1) = cfll(   iy   )
      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,                       &
#ifdef W3_T1
    !$OMP QBO, IX, IY, IY2, IX2, QN                   &
#endif
    !$OMP IXYC, IXYD, QB)
#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) ) )
      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)) )
      !
#ifdef W3_T1
      qbo    = qb
#endif
      !
      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
    !
    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 ( bclose ) 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  &
#ifdef W3_T2
    !$OMP       PRIVATE(QOLD), &
#endif
    !$OMP       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 W3UNO2s : ARRAY DIMENSIONS  :',2i6/           &
         '               USED              :',2i6/           &
         '               BCLOSE, INC       :',l6,i6/         &
         '               NB0, NB1, NB2     :',3i6)
#endif
#ifdef W3_T0
9001 FORMAT ( ' TEST W3UNO2s : DUMP ARRAY ',a,' :')
9002 FORMAT ( 1x,43i3)
9003 FORMAT ( 1x,21i6)
#endif
#ifdef W3_T
9005 FORMAT (' TEST W3UNO2s : GLOBAL CLOSURE (1)')
#endif
    !
#ifdef W3_T1
9010 FORMAT (' TEST W3UNO2s : IP, 2x(IX,IY), CFL (b,i,i+1), ',    &
         ' Q (b,b,i-1,i,i+1,i+2)')
9011 FORMAT (' TEST W3UNO2s :',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 W3UNO2s : GLOBAL CLOSURE (2)')
#endif
    !
#ifdef W3_T2
9020 FORMAT (' TEST W3UNO2s : IP, IXY, 2Q, 2FL')
9021 FORMAT ('            ',2i6,2(1x,2e11.3))
#endif
    !/
    !/ End of W3UNO2s ---------------------------------------------------- /
    !/
  END SUBROUTINE w3uno2s
  !/
  !/ End of module W3UNO2MD -------------------------------------------- /
  !/
END MODULE w3uno2md
!/