#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3UQCKMD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 27-May-2014 | !/ +-----------------------------------+ !/ !/ 08-Feb-2001 : Origination of module. Routines ( version 2.08 ) !/ taken out of w3pro2md.ftn !/ 13-Nov-2001 : Version with obstacles added. ( version 2.14 ) !/ 16-Oct-2002 : Fix par list W3QCK3. ( version 3.00 ) !/ 05-Mar-2008 : Added NEC sxf90 compiler directives. !/ (Chris Bunney, UK Met Office) ( version 3.13 ) !/ 29-May-2009 : Preparing distribution version. ( version 3.14 ) !/ 30-Oct-2009 : Fixed a couple of doc lines. ( version 3.14 ) !/ (T. J. Campbell, NRL) !/ 27-May-2014 : Added OMPH switches in W3QCK3. ( version 5.02 ) !/ !/ Copyright 2009-2014 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ ! 1. Purpose : ! ! Portable ULTIMATE QUICKEST schemes. ! ! 2. Variables and types : ! ! None. ! ! 3. Subroutines and functions : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3QCK1 Subr. Public Original ULTIMATE QUICKEST scheme. ! W3QCK2 Subr. Public UQ scheme for irregular grid. ! W3QCK3 Subr. Public Original ULTIMATE QUICKEST with obst. ! ---------------------------------------------------------------- ! ! 4. Subroutines and functions used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Remarks : ! ! - STRACE and !/S irrelevant for running code. The module is ! therefore fully portable to any other model. ! ! 6. Switches : ! ! !/C90 Cray FORTRAN 90 compiler directives. ! !/NEC NEC SXF90 compiler directives. ! ! !/OMPH Ading OMP directves for hybrid paralellization. ! ! !/S Enable subroutine tracing. ! !/Tn Enable test output. ! ! 7. Source code : ! !/ ------------------------------------------------------------------- / !/S USE W3SERVMD, ONLY: STRACE !/ CONTAINS !/ ------------------------------------------------------------------- / SUBROUTINE W3QCK1 (MX, MY, NX, NY, CFLL, Q, CLOSE, INC, & MAPACT, NACT, MAPBOU, NB0, NB1, NB2, & NDSE, NDST ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 30-Oct-2009 | !/ +-----------------------------------+ !/ !/ 11-Mar-1997 : Final FORTRAN 77 ( version 1.18 ) !/ 15-Dec-1999 : Upgrade to FORTRAN 90 ( version 2.00 ) !/ 15-Feb-2001 : Unit numbers added to par list. ( version 2.08 ) !/ 05-Mar-2008 : Added NEC sxf90 compiler directives. !/ (Chris Bunney, UK Met Office) ( version 3.13 ) !/ 30-Oct-2009 : Fixed "Called by" doc line. ( version 3.14 ) !/ (T. J. Campbell, NRL) !/ ! 1. Purpose : ! ! Preform one-dimensional propagation in a two-dimensional space ! with irregular boundaries and regular grid. ! ! 2. Method : ! ! ULTIMATE QUICKEST scheme (see manual). ! ! Note that the check on monotonous behavior of QCN is performed ! using weights CFAC, to avoid the need for IF statements. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! MX,MY Int. I Field dimensions, if grid is 'closed' or ! circular, MX is the closed dimension. ! NX,NY Int. I Part of field actually used. ! CFLL R.A. I Local Courant numbers. (MY, MX+1) ! Q R.A. I/O Propagated quantity. (MY,0:MX+2) ! CLOSE Log. I Flag for closed 'X' dimension' ! INC Int. I Increment in 1-D array corresponding to ! increment in 2-D space. ! MAPACT I.A. I List of active grid points. ! NACT Int. I Size of MAPACT. ! MAPBOU I.A. I Map with boundary information (see W3MAP2). ! NBn Int. I Counters in MAPBOU. ! NDSE Int. I Error output unit number. ! NDST Int. I Test output unit number. ! ---------------------------------------------------------------- ! - CFLL amd Q need only bee filled in the (MY,MX) range, ! extension is used internally for closure. ! - CFLL and Q are defined as 1-D arrays internally. ! ! 4. Subroutines used : ! ! STRACE Service routine. ! ! 5. Called by : ! ! W3KTP2 Propagation in spectral space ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! - This routine can be used independently from WAVEWATCH III. ! ! 8. Structure : ! ! ------------------------------------------------------ ! 1. Initialize aux. array FLA. ! 2. Fluxes for central points (3rd order + limiter). ! 3. Fluxes boundary point above (1st order). ! 4. Fluxes boundary point below (1st order). ! 5. Closure of 'X' if required ! 6. Propagate. ! ------------------------------------------------------ ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable test output. ! !/T0 Test output input/output fields. ! !/T1 Test output fluxes. ! !/T2 Test output integration. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: MX, MY, NX, NY, INC, MAPACT(MY*MX), & NACT, MAPBOU(MY*MX), NB0, NB1, NB2, & NDSE, NDST REAL, INTENT(INOUT) :: CFLL(MY*(MX+1)), Q(1-MY:MY*(MX+2)) LOGICAL, INTENT(IN) :: CLOSE !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IXY, IP, IXYC, IXYU, IXYD, IY, IX, & IAD00, IAD02, IADN0, IADN1, IADN2 !/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, 'W3QCK1') ! !/T WRITE (NDST,9000) MX, MY, NX, NY, CLOSE, 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 ( CLOSE ) 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) ) ) QB = 0.5 * ( (1.-CFL)*Q(IXY+INC) + (1.+CFL)*Q(IXY) ) & - (1.-CFL**2)/6. * (Q(IXYC-INC)-2.*Q(IXYC)+Q(IXYC+INC)) ! IXYU = IXYC - INC * INT ( SIGN (1.1,CFL) ) IXYD = 2*IXYC - IXYU DQ = Q(IXYD) - Q(IXYU) DQNZ = SIGN ( MAX(1.E-15,ABS(DQ)) , DQ ) QCN = ( Q(IXYC) - Q(IXYU) ) / DQNZ QCN = MIN ( 1.1, MAX ( -0.1 , QCN ) ) ! !/T1 QBO = QB QBN = MAX ( (QB-Q(IXYU))/DQNZ , QCN ) QBN = MIN ( QBN , 1. , QCN/MAX(1.E-10,ABS(CFL)) ) QBR = Q(IXYU) + QBN*DQ CFAC = REAL ( INT( 2. * ABS(QCN-0.5) ) ) QB = (1.-CFAC)*QBR + CFAC*Q(IXYC) ! FLA(IXY) = CFL * QB ! !/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 ( CLOSE ) 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 W3QCK1 : ARRAY DIMENSIONS :',2I6/ & !/T ' USED :',2I6/ & !/T ' CLOSE, INC :',L6,I6/ & !/T ' NB0, NB1, NB2 :',3I6) !/T0 9001 FORMAT ( ' TEST W3QCK1 : DUMP ARRAY ',A,' :') !/T0 9002 FORMAT ( 1X,43I3) !/T0 9003 FORMAT ( 1X,21I6) !/T 9005 FORMAT (' TEST W3QCK1 : GLOBAL CLOSURE (1)') ! !/T1 9010 FORMAT (' TEST W3QCK1 : IP, 2x(IX,IY), CFL (b,i,i+1), ', & !/T1 ' Q (b,b,i-1,i,i+1,i+2)') !/T1 9011 FORMAT (' TEST W3QCK1 :',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 W3QCK1 : GLOBAL CLOSURE (2)') ! !/T2 9020 FORMAT (' TEST W3QCK1 : IP, IXY, 2Q, 2FL') !/T2 9021 FORMAT (' ',2I6,2(1X,2E11.3)) !/ !/ End of W3QCK1 ----------------------------------------------------- / !/ END SUBROUTINE W3QCK1 !/ ------------------------------------------------------------------- / SUBROUTINE W3QCK2 (MX, MY, NX, NY, VELO, DT, DX1, DX2, Q, CLOSE,& INC, MAPACT, NACT, MAPBOU, NB0, NB1, NB2, & NDSE, NDST ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 30-Oct-2009 | !/ +-----------------------------------+ !/ !/ 07-Sep-1997 : Final FORTRAN 77 ( version 1.18 ) !/ 16-Dec-1999 : Upgrade to FORTRAN 90 ( version 2.00 ) !/ 14-Feb-2001 : Unit numbers added to par list. ( version 2.08 ) !/ 05-Mar-2008 : Added NEC sxf90 compiler directives. !/ (Chris Bunney, UK Met Office) ( version 3.13 ) !/ 30-Oct-2009 : Fixed "Called by" doc line. ( version 3.14 ) !/ (T. J. Campbell, NRL) !/ ! 1. Purpose : ! ! Like W3QCK1 with variable grid spacing. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! MX,MY Int. I Field dimensions, if grid is 'closed' or ! circular, MX is the closed dimension. ! NX,NY Int. I Part of field actually used. ! VELO R.A. I Local velocities. (MY, MX+1) ! DT Real I Time step. ! DX1 R.A. I/O Band width at points. (MY, MX+1) ! DX2 R.A. I/O Band width between points. (MY,0:MX+1) ! (local counter and counter+INC) ! Q R.A. I/O Propagated quantity. (MY,0:MX+2) ! CLOSE Log. I Flag for closed 'X' dimension' ! INC Int. I Increment in 1-D array corresponding to ! increment in 2-D space. ! MAPACT I.A. I List of active grid points. ! NACT Int. I Size of MAPACT. ! MAPBOU I.A. I Map with boundary information (see W3MAP2). ! NBn Int. I Counters in MAPBOU. ! NDSE Int. I Error output unit number. ! NDST Int. I Test output unit number. ! ---------------------------------------------------------------- ! - VELO amd Q need only bee filled in the (MY,MX) range, ! extension is used internally for closure. ! - VELO and Q are defined as 1-D arrays internally. ! ! 4. Subroutines used : ! ! STRACE Service routine. ! ! 5. Called by : ! ! W3KTP2 Propagation in spectral space ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! - This routine can be used independently from WAVEWATCH III. ! ! 8. Structure : ! ! ------------------------------------------------------ ! 1. Initialize aux. array FLA. ! 2. Fluxes for central points (3rd order + limiter). ! 3. Fluxes boundary point above (1st order). ! 4. Fluxes boundary point below (1st order). ! 5. Closure of 'X' if required ! 6. Propagate. ! ------------------------------------------------------ ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable test output. ! !/T0 Test output input/output fields. ! !/T1 Test output fluxes. ! !/T2 Test output integration. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: MX, MY, NX, NY, INC, MAPACT(MY*MX), & NACT, MAPBOU(MY*MX), NB0, NB1, NB2, & NDSE, NDST REAL, INTENT(IN) :: DT REAL, INTENT(INOUT) :: VELO(MY*(MX+1)), DX1(MY*(MX+1)), & DX2(1-MY:MY*(MX+1)), Q(1-MY:MY*(MX+2)) LOGICAL, INTENT(IN) :: CLOSE !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IXY, IP, IXYC, IXYU, IXYD, IY, IX, & IAD00, IAD02, IADN0, IADN1, IADN2 !/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, 'W3QCK2') ! !/T WRITE (NDST,9000) MX, MY, NX, NY, DT, CLOSE, 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 ( CLOSE ) 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 ------------------------------------- * ! ( 3rd order + limiter ) ! !/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) ) CFL = DT * VEL / DX2(IXY) IXYC = IXY - INC * INT( MIN ( 0. , SIGN(1.1,CFL) ) ) QB = 0.5 * ( (1.-CFL)*Q(IXY+INC) + (1.+CFL)*Q(IXY) ) & - DX2(IXY)**2 / DX1(IXYC) * (1.-CFL**2) / 6. & * ( (Q(IXYC+INC)-Q(IXYC))/DX2(IXYC) & - (Q(IXYC)-Q(IXYC-INC))/DX2(IXYC-INC) ) ! IXYU = IXYC - INC * INT ( SIGN (1.1,CFL) ) IXYD = 2*IXYC - IXYU DQ = Q(IXYD) - Q(IXYU) DQNZ = SIGN ( MAX(1.E-15,ABS(DQ)) , DQ ) QCN = ( Q(IXYC) - Q(IXYU) ) / DQNZ QCN = MIN ( 1.1, MAX ( -0.1 , QCN ) ) ! !/T1 QBO = QB QBN = MAX ( (QB-Q(IXYU))/DQNZ , QCN ) QBN = MIN ( QBN , 1. , QCN/MAX(1.E-10,ABS(CFL)) ) QBR = Q(IXYU) + QBN*DQ CFAC = REAL ( INT( 2. * ABS(QCN-0.5) ) ) QB = (1.-CFAC)*QBR + CFAC*Q(IXYC) ! FLA(IXY) = VEL * QB ! !/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 * 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 * 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 ( CLOSE ) 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) + DT/DX1(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 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 W3QCK2 : ARRAY DIMENSIONS :',2I6/ & !/T ' USED :',2I6/ & !/T ' TIME STEP :',F8.1/ & !/T ' CLOSE, INC :',L6,I6/ & !/T ' NB0, NB1, NB2 :',3I6) !/T0 9001 FORMAT ( ' TEST W3QCK2 : DUMP ARRAY ',A,' :') !/T0 9002 FORMAT ( 1X,43I3) !/T0 9003 FORMAT ( 1X,21I6) !/T 9005 FORMAT (' TEST W3QCK2 : GLOBAL CLOSURE (1)') ! !/T1 9010 FORMAT (' TEST W3QCK2 : IP, 2x(IX,IY), CFL (b,i,i+1), ', & !/T1 ' Q (b,b,i-1,i,i+1,i+2)') !/T1 9011 FORMAT (' TEST W3QCK2 :',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 W3QCK2 : GLOBAL CLOSURE (2)') ! !/T2 9020 FORMAT (' TEST W3QCK2 : IP, IXY, 2Q, 2FL') !/T2 9021 FORMAT (' ',2I6,2(1X,2E11.3)) !/ !/ End of W3QCK2 ----------------------------------------------------- / !/ END SUBROUTINE W3QCK2 !/ ------------------------------------------------------------------- / SUBROUTINE W3QCK3 (MX, MY, NX, NY, CFLL, TRANS, Q, CLOSE, & INC, MAPACT, NACT, MAPBOU, NB0, NB1, NB2, & NDSE, NDST ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 27-May-2014 | !/ +-----------------------------------+ !/ !/ 13_nov-2001 : Origination. ( version 2.14 ) !/ 16-Oct-2002 : Fix INTENT for TRANS. ( version 3.00 ) !/ 05-Mar-2008 : Added NEC sxf90 compiler directives. !/ (Chris Bunney, UK Met Office) ( version 3.13 ) !/ 27-May-2014 : Added OMPH switches in W3QCK3. ( version 5.02 ) !/ ! 1. Purpose : ! ! Like W3QCK1 with cell transparencies added. ! ! 2. Method : ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! MX,MY Int. I Field dimensions, if grid is 'closed' or ! circular, MX is the closed dimension. ! NX,NY Int. I Part of field actually used. ! CFLL R.A. I Local Courant numbers. (MY, MX+1) ! Q R.A. I/O Propagated quantity. (MY,0:MX+2) ! CLOSE Log. I Flag for closed 'X' dimension' ! INC Int. I Increment in 1-D array corresponding to ! increment in 2-D space. ! MAPACT I.A. I List of active grid points. ! NACT Int. I Size of MAPACT. ! MAPBOU I.A. I Map with boundary information (see W3MAP2). ! NBn Int. I Counters in MAPBOU. ! NDSE Int. I Error output unit number. ! NDST Int. I Test output unit number. ! ---------------------------------------------------------------- ! - CFLL amd Q need only bee filled in the (MY,MX) range, ! extension is used internally for closure. ! - CFLL and Q are defined as 1-D arrays internally. ! ! 4. Subroutines used : ! ! STRACE Service routine. ! ! 5. Called by : ! ! W3XYP2 Propagation in physical space ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! - This routine can be used independently from WAVEWATCH III. ! ! 8. Structure : ! ! ------------------------------------------------------ ! 1. Initialize aux. array FLA. ! 2. Fluxes for central points (3rd order + limiter). ! 3. Fluxes boundary point above (1st order). ! 4. Fluxes boundary point below (1st order). ! 5. Closure of 'X' if required ! 6. Propagate. ! ------------------------------------------------------ ! ! 9. Switches : ! ! !/OMPH Ading OMP directves for hybrid paralellization. ! ! !/S Enable subroutine tracing. ! !/T Enable test output. ! !/T0 Test output input/output fields. ! !/T1 Test output fluxes. ! !/T2 Test output integration. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: MX, MY, NX, NY, INC, MAPACT(MY*MX), & NACT, MAPBOU(MY*MX), NB0, NB1, NB2, & NDSE, NDST REAL, INTENT(IN) :: TRANS(MY*MX,-1:1) REAL, INTENT(INOUT) :: CFLL(MY*(MX+1)), Q(1-MY:MY*(MX+2)) LOGICAL, INTENT(IN) :: CLOSE !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IXY, IP, IXYC, IXYU, IXYD, IY, IX, & IAD00, IAD02, IADN0, IADN1, IADN2, & JN, JP !/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, 'W3QCK3') ! !/T WRITE (NDST,9000) MX, MY, NX, NY, CLOSE, 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 ( CLOSE ) 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) ! 1 ghost column to left Q (IY+IADN1) = Q ( IY ) ! 1st ghost column to right Q (IY+IADN2) = Q (IY+IAD02) ! 2nd ghost column to right CFLL(IY+IADN1) = CFLL( IY ) ! as for Q above, 1st to rt END DO ! !/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, QB, IXYU, IXYD, & !/OMPH/!/T1!$OMP QBO, QN, IX, IY, IX2, IY2, & !/OMPH/!$OMP& DQ, DQNZ, QCN, QBN, QBR, CFAC ) ! DO IP=1, NB0 ! IXY = MAPBOU(IP) CFL = 0.5 * ( CFLL(IXY) + CFLL(IXY+INC) ) IXYC = IXY - INC * INT( MIN ( 0. , SIGN(1.1,CFL) ) ) QB = 0.5 * ( (1.-CFL)*Q(IXY+INC) + (1.+CFL)*Q(IXY) ) & - (1.-CFL**2)/6. * (Q(IXYC-INC)-2.*Q(IXYC)+Q(IXYC+INC)) ! IXYU = IXYC - INC * INT ( SIGN (1.1,CFL) ) IXYD = 2*IXYC - IXYU DQ = Q(IXYD) - Q(IXYU) DQNZ = SIGN ( MAX(1.E-15,ABS(DQ)) , DQ ) QCN = ( Q(IXYC) - Q(IXYU) ) / DQNZ QCN = MIN ( 1.1, MAX ( -0.1 , QCN ) ) ! !/T1 QBO = QB QBN = MAX ( (QB-Q(IXYU))/DQNZ , QCN ) QBN = MIN ( QBN , 1. , QCN/MAX(1.E-10,ABS(CFL)) ) QBR = Q(IXYU) + QBN*DQ CFAC = REAL ( INT( 2. * ABS(QCN-0.5) ) ) QB = (1.-CFAC)*QBR + CFAC*Q(IXYC) ! FLA(IXY) = CFL * QB ! !/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' ! !!!/OMPH/!$OMP PARALLEL DO PRIVATE (IP, IXY, CFL, IXYC) !!! DO IP=NB0+1, NB1 IXY = MAPBOU(IP) CFL = CFLL(IXY) IXYC = IXY - INC * INT( MIN ( 0. , SIGN(1.1,CFL) ) ) FLA(IXY) = CFL * Q(IXYC) !/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 !!! !!!/OMPH/!$OMP END PARALLEL DO ! ! 4. Fluxes for points with boundary below ------------------------- * ! ( 1st order without limiter ) ! !/T1 WRITE (NDST,9011) NB2-NB1, 'BOUNDARY BELOW' ! !!!/OMPH/!$OMP PARALLEL DO PRIVATE (IP, IXY, CFL, IXYC) !!! DO IP=NB1+1, NB2 IXY = MAPBOU(IP) CFL = CFLL(IXY+INC) IXYC = IXY - INC * INT( MIN ( 0. , SIGN(1.1,CFL) ) ) FLA(IXY) = CFL * Q(IXYC) !/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 ! !!!/OMPH/!$OMP END PARALLEL DO ! ! 5. Global closure ----------------------------------------------- * ! IF ( CLOSE ) 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 W3QCK3 : ARRAY DIMENSIONS :',2I6/ & !/T ' USED :',2I6/ & !/T ' CLOSE, INC :',L6,I6/ & !/T ' NB0, NB1, NB2 :',3I6) !/T0 9001 FORMAT ( ' TEST W3QCK3 : DUMP ARRAY ',A,' :') !/T0 9002 FORMAT ( 1X,43I3) !/T0 9003 FORMAT ( 1X,21I6) !/T 9005 FORMAT (' TEST W3QCK3 : GLOBAL CLOSURE (1)') ! !/T1 9010 FORMAT (' TEST W3QCK3 : IP, 2x(IX,IY), CFL (b,i,i+1), ', & !/T1 ' Q (b,b,i-1,i,i+1,i+2)') !/T1 9011 FORMAT (' TEST W3QCK3 :',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 W3QCK3 : GLOBAL CLOSURE (2)') ! !/T2 9020 FORMAT (' TEST W3QCK3 : IP, IXY, 2Q, 2FL') !/T2 9021 FORMAT (' ',2I6,2(1X,2E11.3)) !/ !/ End of W3QCK3 ----------------------------------------------------- / !/ END SUBROUTINE W3QCK3 !/ !/ End of module W3UQCKMD -------------------------------------------- / !/ END MODULE W3UQCKMD