#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 !/