#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3PSMCMD !/ !/ +------------------------------------+ !/ | Spherical Multiple-Cell (SMC) grid | !/ | Adv, GCT, Rfr, Dif subroutines. | !/ | Jian-Guo Li | !/ | First created: 8 Nov 2010 | !/ | Last modified: 18 Apr 2018 | !/ +------------------------------------+ !/ !/ 08-Nov-2010 : Coding started by adapting w3pro2md.ftn. !/ 18-Nov-2010 : Refraction and GCT by rotation and k-shift. !/ 12-Apr-2011 : Restore x-advective flux for intermediate update. !/ 3-Jun-2011 : New refraction formulation using Cg only. !/ 8-Jun-2011 : Optimise classic refraction formulation. !/ 16-Jun-2011 : Add refraction limter to gradient direction. !/ 1-Jul-2011 : New refraction using Cg and gradient limiter. !/ 28-Jul-2011 : Finalise with old refraction scheme and gradient limiter. !/ 4-Nov-2011 : Separate x and y obstruction coefficients. !/ 5-Jan-2012 : Update to multi-resolution SMC grid with sub-time-steps. !/ 2-Feb-2012 : Separate single- and multi-resolution advection. !/ 6-Mar-2012 : Tidy up code and minor adjustments, CLATF. !/ 12-Mar-2012 : Remove net flux bug and optimise upstream code. !/ 16-Jan-2013 : Adapted for Version 4.08, removing FACX/Y. !/ 16-Sep-2013 : Add Arctic part for SMC grid in WW3 V4.11 !/ 3-Jan-2014 : Remove bug in SMCDHXY for AU/V as cell size. !/ 7-Jan-2014 : Remove bug in SMCGtCrfr for K definition. !/ 28-Jan-2014 : Move Arctic boundary condition update out. !/ 18-Aug-2015 : New gradient, average and 3rd order advection subs. !/ 3-Sep-2015 : UNO3 advection scheme by logical option FUNO3. !/ 14-Sep-2015 : Modify DHDX/Y for Arctic part refraction term. !/ 8-Aug-2017 : Update SMCGradn for 0 or 1 boundary conditions. !/ 9-Jan-2018 : Parallelization by adding OpenMP directives. !/ ! 1. Purpose : ! ! Bundles routines for SMC advection (UNO2) and diffusion schemes in ! single module, including GCT and refraction rotation schemes. ! ! 2. Variables and types : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! TRNMIN R.P. Private Minimum transparancy for local ! ---------------------------------------------------------------- ! ! 3. Subroutines and functions : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3PSMC Subr. Public Spatial propagation on SMC grid. ! W3KSMC Subr. Public Spectral modification by GCT and refraction. ! SMCxUNO2 Subr. Public Irregular grid mid-flux on U-faces by UNO2. ! SMCyUNO2 Subr. Public Irregular grid mid-flux on V-faces by UNO2. ! SMCxUNO2r Subr. Public Regular grid mid-flux on U-faces by UNO2. ! SMCyUNO2r Subr. Public Regular grid mid-flux on V-faces by UNO2. ! SMCkUNO2 Subr. Public Shift in k-space due to refraction by UNO2. ! SMCxUNO3 Subr. Public Irregular grid mid-flux on U-faces by UNO3. ! SMCyUNO3 Subr. Public Irregular grid mid-flux on V-faces by UNO3. ! SMCxUNO3r Subr. Public Regular grid 3rd order U-mid-flux by UNO3. ! SMCyUNO3r Subr. Public Regular grid 3rd order V-mid-flux by UNO3. ! SMCGtCrfr Subr. Public Refraction and GCT rotation in theta. ! SMCDHXY Subr. Public Evaluate depth gradient and refraction limiter. ! SMCGradn Subr. Public Evaluate local gradient for sea points. ! SMCAverg Subr. Public Numerical 1-2-1 average of sea point field. ! W3GATHSMC W3SCATSMC Gather and scatter spectral components. ! ---------------------------------------------------------------- ! ! 4. Subroutines and functions used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! W3ACTURN Subr. W3SERVMD Subroutine rotating action spectrum. ! ---------------------------------------------------------------- ! ! 5. Remarks : ! ! 6. Switches : ! ! !/MGP Correct for motion of grid. ! !/S Enable subroutine tracing. ! !/Tn Enable test output. ! !/ARC Enable Arctic part. ! ! 7. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ Use omp_lib for OpenMP functions if switched on. JGLi10Jan2018 !$ USE omp_lib !/ !/ Public variables !/ PUBLIC !/ !/ Private data !/ REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95 !/ CONTAINS !/ ------------------------------------------------------------------- / !/ SUBROUTINE W3PSMC ( ISP, DTG, VQ ) !/ !/ +------------------------------------+ !/ | Spherical Multiple-Cell (SMC) grid | !/ | Advection and diffusion sub. | !/ | First created: JG Li 8 Nov 2010 | !/ | Last modified: JG Li 18 Apr 2018 | !/ +------------------------------------+ !/ !/ 08-Nov-2010 : Origination. JGLi ( version 1.00 ) !/ 16-Dec-2010 : Check U/V CFL values. JGLi ( version 1.10 ) !/ 18-Mar-2011 : Check MPI communication. JGLi ( version 1.20 ) !/ 16-May-2011 : Tidy up diagnosis lines. JGLi ( version 1.30 ) !/ 4 Nov-2011 : Separate x and y obstruc. JGLi ( version 1.40 ) !/ 5 Jan-2012 : Multi-resolution SMC grid. JGLi ( version 1.50 ) !/ 2 Feb-2012 : Separate single multi adv. JGLi ( version 1.60 ) !/ 6 Mar-2012 : Minor adjustments of CLATF. JGLi ( version 1.70 ) !/ 12 Feb-2012 : Remove net flux bug. JGLi ( version 1.80 ) !/ 16 Sep-2013 : Add Arctic part. JGLi ( version 2.00 ) !/ 3 Sep-2015 : Gradient, UNO3 and Average. JGLi ( version 2.10 ) !/ 26 Feb-2016 : Update boundary spectra. JGLi ( version 2.20 ) !/ 23 Mar-2016 : Add current option. JGLi ( version 2.30 ) !/ 18 Apr-2018 : Refined sub-grid blocking. JGLi ( version 2.40 ) !/ ! 1. Purpose : ! ! Propagation in phyiscal space for a given spectral component. ! ! 2. Method : ! ! Unstructured SMC grid, point-oriented face and cell loops. ! UNO2 advection scheme and isotropic FTCS diffusion scheme. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ISP Int. I Number of spectral bin (IK-1)*NTH+ITH ! FACX/Y Real I Factor in propagation velocity (1 or 0 *DT/DX) ! DTG Real I Total time step. ! MAPSTA I.A. I Grid point status map. ! MAPFS I.A. I Storage map. ! VQ R.A. I/O Field to propagate. ! ---------------------------------------------------------------- ! ! Local variables. ! ---------------------------------------------------------------- ! NTLOC Int Number of local time steps. ! DTLOC Real Local propagation time step. ! CGD Real Deep water group velocity. ! DSSD, DNND Deep water diffusion coefficients. ! ULCFLX R.A. Local courant numbers in 'x' (norm. velocities) ! VLCFLY R.A. Id. in 'y'. ! CXTOT R.A. Propagation velocities in physical space. ! CYTOT R.A. ! DFRR Real Relative frequency increment. ! DX0I Real Inverted grid incremenent in meters (longitude, eq.). ! DY0I Real Inverted grid incremenent in meters (latitude). ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! See module documentation. ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3WAVE Subr. W3WAVEMD Wave model routine. ! --------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! --------------------------------------------- ! 1. Preparations ! a Set constants ! b Initialize arrays ! 2. Prepare arrays ! a Velocities and 'Q' ! b diffusion coefficients ! 3. Loop over sub-steps ! ---------------------------------------- ! a Propagate ! b Update boundary conditions ! c Diffusion correction ! ---------------------------------------- ! 4. Store Q field in spectra ! --------------------------------------------- ! ! 9. Switches : ! ! !/MGP Correct for motion of grid. ! ! !/TDYN Dynamic increase of DTME ! !/DSS0 Disable diffusion in propagation direction ! !/XW0 Propagation diffusion only. ! !/XW1 Growth diffusion only. ! ! !/S Enable subroutine tracing. ! ! !/T Enable general test output. ! !/T1 Dump of input field and fluxes. ! !/T2 Dump of output field. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS ! USE W3TIMEMD, ONLY: DSEC21 ! USE W3GDATMD, ONLY: NK, NTH, DTH, XFR, ESIN, ECOS, SIG, NX, NY, & NSEA, SX, SY, MAPSF, FUNO3, FVERG, & IJKCel, IJKUFc, IJKVFc, NCel, NUFc, NVFc, & NLvCel, NLvUFc, NLvVFc, NRLv, MRFct, & DTCFL, CLATS, DTME, CLATMN, CTRNX, CTRNY !/ARC USE W3GDATMD, ONLY: NGLO, ANGARC USE W3WDATMD, ONLY: TIME USE W3ADATMD, ONLY: CG, WN, U10, CX, CY, ATRNX, ATRNY, ITIME ! USE W3IDATMD, ONLY: FLCUR USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, TBPI0, TBPIN, & ISBPI, BBPI0, BBPIN ! !/S USE W3SERVMD, ONLY: STRACE !/ IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: ISP !/ REAL, INTENT(IN) :: FACX, FACY, DTG REAL, INTENT(IN) :: DTG REAL, INTENT(INOUT) :: VQ(NSEA) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: ITH, IK, NTLOC, ITLOC, ISEA, IXY, & IY, IY0, IP, IBI, LvR INTEGER :: i, j, k, L, M, N, LL, MM, NN, LMN, & iuf, juf, ivf, jvf, icl, jcl !/S INTEGER, SAVE :: IENT = 0 REAL :: CG0, CGA, CGN, CGX, CGY, FMR, RD1, & RD2, CXMIN, CXMAX, CYMIN, CYMAX, & CXC, CYC, DTLDX, DTLDY REAL :: DTLOC, CGCOS, CGSIN, FUTRN, FVTRN, & DFRR, DX0I, DY0I, CGD, DSSD, & DNND, DCELL, XWIND, TFAC, DSS, DNN !/ARC REAL :: PCArea, ARCTH LOGICAL :: YFIRST !/ !/ Automatic work arrays ! REAL, Dimension(-9:NCel) :: FCNt, AFCN, BCNt, UCFL, VCFL, CQ, & CQA, CXTOT, CYTOT REAL, Dimension( NUFc) :: FUMD, FUDIFX, ULCFLX REAL, Dimension( NVFc) :: FVMD, FVDIFY, VLCFLY !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3PSMC') ! ! 1. Preparations --------------------------------------------------- * !!Li Spectral bin direction and frequency indices ITH = 1 + MOD(ISP-1,NTH) IK = 1 + (ISP-1)/NTH !!Li Maximum group speed for 1st and the transported frequency bin !!Li A factor of 1.2 is added to account for the shallow water peak. CG0 = 0.6 * GRAV / SIG(1) CGA = 0.6 * GRAV / SIG(IK) !!Li Maximum group speed for given spectral bin. First bin speed is !!Li used to avoid zero speed component. ! CGX = ABS( CGA * ECOS(ITH) ) ! CGY = ABS( CGA * ESIN(ITH) ) CGX = CGA * ECOS(ITH) CGY = CGA * ESIN(ITH) !!Li Add maximum current components to maximum group components. IF ( FLCUR ) THEN CXMIN = MINVAL ( CX(1:NSEA) ) CXMAX = MAXVAL ( CX(1:NSEA) ) CYMIN = MINVAL ( CY(1:NSEA) ) CYMAX = MAXVAL ( CY(1:NSEA) ) IF ( ABS(CGX+CXMIN) .GT. ABS(CGX+CXMAX) ) THEN CGX = CGX + CXMIN ELSE CGX = CGX + CXMAX END IF IF ( ABS(CGY+CYMIN) .GT. ABS(CGY+CYMAX) ) THEN CGY = CGY + CYMIN ELSE CGY = CGY + CYMAX END IF CXC = MAX ( ABS(CXMIN) , ABS(CXMAX) ) CYC = MAX ( ABS(CYMIN) , ABS(CYMAX) ) ELSE CXC = 0. CYC = 0. END IF !!Li Base-cell grid lenth at Equator (size-4 on SMC625 grid). DX0I = 1.0/(SX * DERA * RADIUS) DY0I = 1.0/(SY * DERA * RADIUS) !!Li Miminum time step determined by Courant number < 0.8 !!Li Note, minimum x grid length is half the Equator value. !!Li Minimum time step should not be less than sub w3init requirement, !!Li where IAPPRO array is initialised for propagation parallization. CGN = 0.9999 * MAX( ABS(CGX), ABS(CGY), CXC, CYC, 0.001*CG0 ) DTLOC = DTCFL*CG0/CGN NTLOC = 1 + INT(DTG/DTLOC - 0.001) DTLOC = DTG / REAL(NTLOC) !!Li Group speed component common factors, FACX=DTG*DX0I !!Li FACX and FACY are evaluated here directly. JGLi16Jan2013 ! CGCOS = FACX * ECOS(ITH) / REAL(NTLOC) ! CGSIN = FACY * ESIN(ITH) / REAL(NTLOC) CGCOS = ECOS(ITH) CGSIN = ESIN(ITH) DTLDX = DTLOC * DX0I DTLDY = DTLOC * DY0I ! YFIRST = MOD(ITIME,2) .EQ. 0 ! !Li Homogenous diffusion Fourier number DNND and DSSD will be used. !Li They have to be divided by base-cell size for size-1 stability. !Li So they are equivalent to the Fourier number in size-1 cell at !Li the sub-time step DTLOC/MRFct. IF ( DTME .GT. 0. ) THEN DFRR = XFR - 1. CGD = 0.5 * GRAV / SIG(IK) DNN = ((DTH*CGD)**2)*DTME / 12. DNND = DNN*DTLOC*(DX0I*DX0I) DSSD = DNN*DTLOC*(DY0I*DY0I) ELSE DSSD = 0.0 DNND = 0.0 END IF ! ! 1.b Initialize arrays ! !/T WRITE (NDST,9010) ! ULCFLX = 0. VLCFLY = 0. !Li Pass spectral element VQ to CQ and define size-1 cell CFL !$OMP Parallel DO Private(ISEA) DO ISEA=1, NSEA !Li Transported variable is divided by CG as in WW3 (???) CQ(ISEA) = VQ(ISEA)/CG(IK,ISEA) !Li Resetting NaNQ VQ to zero if any. JGLi18Mar2013 IF( .NOT. (CQ(ISEA) .EQ. CQ(ISEA)) ) CQ(ISEA) = 0.0 END DO !$OMP END Parallel DO !Li Add current components if any to wave velocity. IF ( FLCUR ) THEN !$OMP Parallel DO Private(ISEA) DO ISEA=1, NSEA CXTOT(ISEA) = (CGCOS * CG(IK,ISEA) + CX(ISEA)) CYTOT(ISEA) = (CGSIN * CG(IK,ISEA) + CY(ISEA)) ENDDO !$OMP END Parallel DO ELSE !Li No current case use group speed only. !$OMP Parallel DO Private(ISEA) DO ISEA=1, NSEA CXTOT(ISEA) = CGCOS * CG(IK,ISEA) CYTOT(ISEA) = CGSIN * CG(IK,ISEA) END DO !$OMP END Parallel DO !Li End of IF( FLCUR ) block. ENDIF !/ARC !Li Arctic cell velocity components need to be rotated !/ARC !Li back to local east referenence system for propagation. !/ARC DO ISEA=NGLO+1, NSEA !/ARC ARCTH = ANGARC(ISEA-NGLO)*DERA !/ARC CXC = CXTOT(ISEA)*COS(ARCTH) + CYTOT(ISEA)*SIN(ARCTH) !/ARC CYC = CYTOT(ISEA)*COS(ARCTH) - CXTOT(ISEA)*SIN(ARCTH) !/ARC CXTOT(ISEA) = CXC !/ARC CYTOT(ISEA) = CYC !/ARC END DO !/ARC !Li Polar cell area factor for V-flux update !/ARC PCArea = DY0I/(MRFct*PI*DX0I*FLOAT(IJKCel(4,NSEA))) !/ARC !Li V-component is reset to zero for Polar cell as direction !/ARC !Li is undefined there. !/ARC CYTOT(NSEA) = 0.0 !/ARC !Li Convert velocity components into CFL factors. !$OMP Parallel DO Private(ISEA) DO ISEA=1, NSEA UCFL(ISEA) = DTLDX*CXTOT(ISEA)/CLATS(ISEA) VCFL(ISEA) = DTLDY*CYTOT(ISEA) ENDDO !$OMP END Parallel DO !Li Initialise boundary cell CQ and Velocity values. CQ(-9:0)=0.0 UCFL(-9:0)=0.0 VCFL(-9:0)=0.0 ! ! 3. Loop over frequency-dependent sub-steps -------------------------* ! DO ITLOC=1, NTLOC ! ! Initialise net flux arrays. FCNt(-9:NCel) = 0.0 AFCN(-9:NCel) = 0.0 BCNt(-9:NCel) = 0.0 ! ! Single-resolution SMC grid uses regular grid advection with ! partial blocking enabled when NRLv = 1 IF ( NRLv .EQ. 1 ) THEN IF( FUNO3 ) THEN ! Use 3rd order UNO3 scheme. JGLi20Aug2015 CALL SMCxUNO3r(1, NUFc, CQ, UCFL, ULCFLX, DNND, FUMD, FUDIFX) ELSE ! Call SMCxUNO2 to calculate MFx value CALL SMCxUNO2r(1, NUFc, CQ, UCFL, ULCFLX, DNND, FUMD, FUDIFX) ENDIF ! Store conservative flux in FCNt advective one in AFCN !$OMP Parallel DO Private(i, M, N, FUTRN) DO i=1, NUFc M=IJKUFc(5,i) N=IJKUFc(6,i) FUTRN = FUMD(i)*ULCFLX(i) - FUDIFX(i) !! Add sub-grid transparency for input flux update. JGLi16May2011 !! Transparency is also applied on diffusion flux. JGLi12Mar2012 !$OMP CRITICAL IF( (CTRNX(M)+CTRNX(N)) .GE. 1.96 ) THEN FCNt(M) = FCNt(M) - FUTRN FCNt(N) = FCNt(N) + FUTRN ELSE IF( ULCFLX(i) .GE. 0.0 ) THEN FCNt(M) = FCNt(M) - FUTRN*CTRNX(M) FCNt(N) = FCNt(N) + FUTRN*CTRNX(M)*CTRNX(N) ELSE FCNt(M) = FCNt(M) - FUTRN*CTRNX(N)*CTRNX(M) FCNt(N) = FCNt(N) + FUTRN*CTRNX(N) ENDIF ! Also divided by another cell length as UCFL is in basic unit. AFCN(M) = AFCN(M) - FUMD(i)*UCFL(M) + FUDIFX(i) AFCN(N) = AFCN(N) + FUMD(i)*UCFL(N) - FUDIFX(i) !$OMP END CRITICAL ENDDO !$OMP END Parallel DO ! Store conservative update in CQA and advective update in CQ ! The side length in MF value has to be cancelled with cell length ! Note ULCFLX has been divided by the cell size inside SMCxUNO2. !$OMP Parallel DO Private(n) DO n=1, NSEA CQA(n)=CQ(n) + FCNt(n)/FLOAT(IJKCel(3,n)) CQ (n)=CQ(n) + AFCN(n)/FLOAT(IJKCel(3,n)) ENDDO !$OMP END Parallel DO ! Call advection subs. IF( FUNO3 ) THEN ! Use 3rd order UNO3 scheme. JGLi20Aug2015 CALL SMCyUNO3r(1, NVFc, CQ, VCFL, VLCFLY, DSSD, FVMD, FVDIFY) ELSE ! Call SMCyUNO2 to calculate MFy value CALL SMCyUNO2r(1, NVFc, CQ, VCFL, VLCFLY, DSSD, FVMD, FVDIFY) ENDIF !$OMP Parallel DO Private(j, M, N, FVTRN) DO j=1, NVFc M=IJKVFc(5,j) N=IJKVFc(6,j) FVTRN = FVMD(j)*VLCFLY(j) - FVDIFY(j) !! Add sub-grid transparency for input flux update. JGLi16May2011 !! Transparency is also applied on diffusion flux. JGLi12Mar2012 !$OMP CRITICAL IF( (CTRNY(M)+CTRNY(N)) .GE. 1.96 ) THEN BCNt(M) = BCNt(M) - FVTRN BCNt(N) = BCNt(N) + FVTRN ELSE IF( VLCFLY(j) .GE. 0.0 ) THEN BCNt(M) = BCNt(M) - FVTRN*CTRNY(M) BCNt(N) = BCNt(N) + FVTRN*CTRNY(M)*CTRNY(N) ELSE BCNt(M) = BCNt(M) - FVTRN*CTRNY(N)*CTRNY(M) BCNt(N) = BCNt(N) + FVTRN*CTRNY(N) ENDIF !$OMP END CRITICAL ENDDO !$OMP END Parallel DO ! Store conservative update of CQA in CQ ! The v side length in MF value has to be cancelled with cell length !! One cosine factor is also needed to be divided for SMC grid !$OMP Parallel DO Private(n) DO n=1, NSEA CQ(n)=CQA(n) + BCNt(n)/( CLATS(n)*FLOAT(IJKCel(3,n)) ) ENDDO !$OMP END Parallel DO !/ARC !Li Polar cell needs a special area factor, one-level case. !/ARC CQ(NSEA) = CQA(NSEA) + BCNt(NSEA)*PCArea ! End of single-resolution advection and diffusion. ELSE ! Multi-resolution SMC grid uses irregular grid advection scheme ! without partial blocking when NRLv > 1 ! ! 3.a Multiresolution sub-steps depend on refined levels MRFct DO LMN = 1, MRFct ! ! 3.b Loop over all levels, starting from the finest level. DO LL= 1, NRLv ! Cell size of this level LvR=2**(LL - 1) FMR=FLOAT( LvR ) ! ! 3.c Calculate this level only if size is factor of LMN IF( MOD(LMN, LvR) .EQ. 0 ) THEN ! ! 3.d Select cell and face ranges icl=NLvCel(LL-1)+1 iuf=NLvUFc(LL-1)+1 ivf=NLvVFc(LL-1)+1 jcl=NLvCel(LL) juf=NLvUFc(LL) jvf=NLvVFc(LL) ! ! Use 3rd order UNO3 scheme. JGLi03Sep2015 IF( FUNO3 ) THEN CALL SMCxUNO3(iuf, juf, CQ, UCFL, ULCFLX, DNND, FUMD, FUDIFX, FMR) ELSE ! Call SMCxUNO2 to calculate finest level (size-1) MFx value CALL SMCxUNO2(iuf, juf, CQ, UCFL, ULCFLX, DNND, FUMD, FUDIFX, FMR) ENDIF ! Store fineset level conservative flux in FCNt advective one in AFCN !$OMP Parallel DO Private(i, L, M, FUTRN) DO i=iuf, juf L=IJKUFc(5,i) M=IJKUFc(6,i) FUTRN = FUMD(i)*ULCFLX(i) - FUDIFX(i) !$OMP CRITICAL !! Add sub-grid blocking for refined cells. JGLi18Apr2018 IF( (CTRNX(M)+CTRNX(L)) .GE. 1.96 ) THEN FCNt(L) = FCNt(L) - FUTRN FCNt(M) = FCNt(M) + FUTRN ELSE IF( ULCFLX(i) .GE. 0.0 ) THEN FCNt(L) = FCNt(L) - FUTRN*CTRNX(L) FCNt(M) = FCNt(M) + FUTRN*CTRNX(M)*CTRNX(L) ELSE FCNt(L) = FCNt(L) - FUTRN*CTRNX(L)*CTRNX(M) FCNt(M) = FCNt(M) + FUTRN*CTRNX(M) ENDIF AFCN(L) = AFCN(L) - FUMD(i)*UCFL(L)*FMR + FUDIFX(i) AFCN(M) = AFCN(M) + FUMD(i)*UCFL(M)*FMR - FUDIFX(i) !$OMP END CRITICAL ENDDO !$OMP END Parallel DO ! Store conservative update in CQA and advective update in CQ ! The side length in MF value has to be cancelled with cell y-length. ! Also divided by another cell x-size as UCFL is in size-1 unit. !$OMP Parallel DO Private(n) DO n=icl, jcl CQA(n)=CQ(n) + FCNt(n)/FLOAT( IJKCel(3, n)*IJKCel(4, n) ) CQ (n)=CQ(n) + AFCN(n)/FLOAT( IJKCel(3, n)*IJKCel(4, n) ) FCNt(n)=0.0 AFCN(n)=0.0 ENDDO !$OMP END Parallel DO ! ! Use 3rd order UNO3 scheme. JGLi03Sep2015 IF( FUNO3 ) THEN CALL SMCyUNO3(ivf, jvf, CQ, VCFL, VLCFLY, DSSD, FVMD, FVDIFY, FMR) ELSE ! Call SMCyUNO2 to calculate MFy value CALL SMCyUNO2(ivf, jvf, CQ, VCFL, VLCFLY, DSSD, FVMD, FVDIFY, FMR) ENDIF ! ! Store conservative flux in BCNt !$OMP Parallel DO Private(j, L, M, FVTRN) DO j=ivf, jvf L=IJKVFc(5,j) M=IJKVFc(6,j) FVTRN = FVMD(j)*VLCFLY(j) - FVDIFY(j) !$OMP CRITICAL !! Add sub-grid blocking for refined cells. JGLi18Apr2018 IF( (CTRNY(M)+CTRNY(L)) .GE. 1.96 ) THEN BCNt(L) = BCNt(L) - FVTRN BCNt(M) = BCNt(M) + FVTRN ELSE IF( VLCFLY(j) .GE. 0.0 ) THEN BCNt(L) = BCNt(L) - FVTRN*CTRNY(L) BCNt(M) = BCNt(M) + FVTRN*CTRNY(M)*CTRNY(L) ELSE BCNt(L) = BCNt(L) - FVTRN*CTRNY(L)*CTRNY(M) BCNt(M) = BCNt(M) + FVTRN*CTRNY(M) ENDIF !$OMP END CRITICAL ENDDO !$OMP END Parallel DO ! Store conservative update of CQA in CQ ! The v side length in MF value has to be cancelled with x-size. ! Also divided by cell y-size as VCFL is in size-1 unit. !! One cosine factor is also needed to be divided for SMC grid. !$OMP Parallel DO Private(n) DO n=icl, jcl CQ(n)=CQA(n) + BCNt(n)/( CLATS(n)* & & FLOAT( IJKCel(3, n)*IJKCel(4, n) ) ) BCNt(n)=0.0 ENDDO !$OMP END Parallel DO !/ARC !Li Polar cell needs a special area factor, multi-level case. !/ARC IF( jcl .EQ. NSEA ) THEN !/ARC CQ(NSEA) = CQA(NSEA) + BCNt(NSEA)*PCArea !/ARC ENDIF ! ! End of refine level if block MOD(LMN, LvR) .EQ. 0 ENDIF ! End of refine level loop LL=1, NRLv ENDDO !! !! END of multi-resolution sub-step loop LMN = 1, MRFct ENDDO ! End of multi-resolution advection ELSE block of NRLv > 1 ENDIF !! Update boundary spectra if any. JGLi26Feb2016 ! IF ( FLBPI ) THEN RD1 = DSEC21(TBPI0, TIME)-DTG*REAL(NTLOC-ITLOC)/REAL(NTLOC) RD2 = DSEC21(TBPI0, TBPIN) IF ( RD2 .GT. 0.001 ) THEN RD2 = MIN(1.,MAX(0.,RD1/RD2)) RD1 = 1. - RD2 ELSE RD1 = 0. RD2 = 1. END IF DO IBI=1, NBI ISEA = ISBPI(IBI) CQ(ISEA) = (RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI)) & /CG(IK,ISEA) END DO ENDIF ! !! End of ITLOC DO ENDDO ! Average with 1-2-1 scheme. JGLi20Aug2015 IF(FVERG) CALL SMCAverg(CQ) ! ! 4. Store results in VQ in proper format --------------------------- * ! !$OMP Parallel DO Private(ISEA) DO ISEA=1, NSEA VQ(ISEA) = MAX ( 0. , CQ(ISEA)*CG(IK,ISEA) ) END DO !$OMP END Parallel DO ! RETURN ! ! Formats ! !/T 9001 FORMAT (' TEST W3PSMC : ISP, ITH, IK, COS-SIN :',I8,2I4,2F7.3) !/T 9003 FORMAT (' TEST W3PSMC : NO DISPERSION CORRECTION ') ! !/T 9010 FORMAT (' TEST W3PSMC : INITIALIZE ARRAYS') ! !/T 9020 FORMAT (' TEST W3PSMC : CALCULATING LCFLX/Y AND DSS/NN (NSEA=', & !/T I6,')') !/T1 9021 FORMAT (1X,I6,2I5,E12.4,2f7.3) !/T 9022 FORMAT (' TEST W3PSMC : CORRECTING FOR CURRENT') ! !/T 9040 FORMAT (' TEST W3PSMC : FIELD AFTER PROP. (NSEA=',I6,')') !/T2 9041 FORMAT (1X,I6,2I5,E12.4) !/ !/ End of W3PSMC ----------------------------------------------------- / !/ END SUBROUTINE W3PSMC !/ !/ ------------------------------------------------------------------- / !/ SUBROUTINE W3KRTN ( ISEA, FACTH, FACK, CTHG0, CG, WN, DEPTH, & DDDX, DDDY, ALFLMT, CX, CY, DCXDX, DCXDY, & DCYDX, DCYDY, DCDX, DCDY, VA ) !/ !/ +------------------------------------+ !/ | Spherical Multiple-Cell (SMC) grid | !/ | Refraction and great-cirle turning | !/ | Jian-Guo Li | !/ | First created: 8 Nov 2010 | !/ | Last modified: 06-Jun-2018 | !/ +------------------------------------+ !/ !/ 08-Nov-2010 : Origination. ( version 1.00 ) !/ 10-Jun-2011 : New refraction formulation. ( version 1.10 ) !/ 16-Jun-2011 : Add refraction limiter to gradient. ( version 1.20 ) !/ 21-Jul-2011 : Old refraction formula + limiter. ( version 1.30 ) !/ 26-Jul-2011 : Tidy up refraction schemes. ( version 1.40 ) !/ 28-Jul-2011 : Finalise with old refraction. ( version 1.50 ) !/ 23-Mar-2016 : Add current option in refraction. ( version 2.30 ) !/ 06-Jun-2018 : Add DEBUGDCXDX ( version 6.04 ) !/ !/ ! 1. Purpose : ! ! Refraction and great-circle turning by spectral rotation. ! ! 2. Method : ! ! Linear interpolation equivalent to 1st order upstream scheme ! but without restriction on rotation angle. However, refraction ! is limited towards the depth gradient direction (< 90 degree). ! Refraction induced spectral shift in the k-space will remain ! to be advected using the UNO2 scheme. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ISEA Int. I Number of sea point. ! FACTH/K Real I Factor in propagation velocity. ! CTHG0 Real I Factor in great circle refracftion term. ! MAPxx2 I.A. I Propagation and storage maps. ! CG R.A. I Local group velocities. ! WN R.A. I Local wavenumbers. ! DEPTH R.A. I Depth. ! DDDx Real I Depth gradients. ! CX/Y Real I Current components. ! DCxDx Real I Current gradients. ! DCDx Real I Phase speed gradients. ! VA R.A. I/O Spectrum. ! ---------------------------------------------------------------- ! ! Local variables. ! ---------------------------------------------------------------- ! DPH2K R.A. 2*Depth*Wave_number_K ! SNH2K R.A. SINH(2*Depth*Wave_number_K) ! FDD, FDU, FGC, FCD, FCU ! R.A. Directionally varying part of depth, current and ! great-circle refraction terms and of consit. ! of Ck term. ! CFLT-K R.A. Propagation velocities of local fluxes. ! DB R.A. Wavenumber band widths at cell centers. ! DM R.A. Wavenumber band widths between cell centers and ! next cell center. ! Q R.A. Extracted spectrum ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! SMCGtCrfr Refraction and GCT rotation in theta. ! SMCkUNO2 Refraction shift in k-space by UNO2. ! STRACE Service routine. ! ! 5. Called by : ! ! W3WAVE Wave model routine. ! ! 6. Error messages : ! ! None. ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Preparations ! a Initialize arrays ! b Set constants and counters ! 2. Point preparations ! a Calculate SNH2K ! b Extract spectrum ! 3. Refraction velocities ! a Filter level depth reffraction. ! b Depth refratcion velocity. ! c Current refraction velocity. ! 4. Wavenumber shift velocities ! a Prepare directional arrays ! b Calcuate velocity. ! 5. Propagate. ! 6. Store results. ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable general test output. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, DSIP, ECOS, ESIN, & EC2, ESC, ES2, FLCTH, FLCK, CTMAX, DTH USE W3ADATMD, ONLY: ITIME USE W3IDATMD, ONLY: FLCUR USE W3ODATMD, ONLY: NDSE, NDST !/S USE W3SERVMD, ONLY: STRACE !/ IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: ISEA !/S INTEGER, SAVE :: IENT = 0 REAL, INTENT(IN) :: FACTH, FACK, CTHG0, CG(0:NK+1), & WN(0:NK+1), DEPTH, DDDX, DDDY, & ALFLMT(NTH), CX, CY, DCXDX, DCXDY, & DCYDX, DCYDY, DCDX(0:NK+1), DCDY(0:NK+1) REAL, INTENT(INOUT) :: VA(NSPEC) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: ITH, IK, ISP REAL :: FGC, FKD, FKS, FRK(NK), FRG(NK), DDNorm(NTH), & FKC(NTH), VQ(NSPEC), VCFLT(NSPEC), DEPTH30, & DB(0:NK+1), DM(-1:NK+1), CFLK(NTH,0:NK), & !Li For new refraction scheme using Cg. JGLi26Jul2011 ! DPH2K(0:NK+1), SNH2K(0:NK+1) !Li For old refraction scheme using phase speed. JGLi26Jul2011 SIGSNH(0:NK+1) !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3KRTN') ! ! 1. Preparation for point ------------------------------------------ * ! Array with partial derivative of sigma versus depth !Li Use of minimum depth 30 m for refraction factor. JGLi12Feb2014 DEPTH30=MAX(30.0, DEPTH) DO IK=0, NK+1 !Li For old refraction scheme using phase speed. JGLi8Jun2011 ! DPH2K(IK) = 2.0*WN(IK)*DEPTH ! SNH2K(IK) = SINH( DPH2K(IK) ) !Li For new refraction scheme using Cg. JGLi3Jun2011 !! SIGSNH(IK) = SIG(IK)/SINH(2.0*WN(IK)*DEPTH) !!AC Replacing SIGSINH with a delimiter to prevent the SINH value from !!AC becoming significantly large. Right now set to a max around 1E21 ! SIGSNH(IK) = SIG(IK)/SINH(MIN(2.0*WN(IK)*DEPTH,50.0)) !Li Refraction factor uses minimum depth of 30 m. JGLi12Feb2014 SIGSNH(IK) = SIG(IK)/SINH(MIN(2.0*WN(IK)*DEPTH30,50.0)) END DO ! ! 2. Extract spectrum without mapping ! VQ = VA ! 3. Refraction velocities ------------------------------------------ * ! ! IF ( FLCTH ) THEN ! ! 3.a Set slope filter for depth refraction ! !Li Lift theta-refraction limit and use new formulation. 25 Nov 2010 !Li FACTH = DTG / DTH / REAL(NTLOC), CTHG0 = - TAN(DERA*Y) / RADIUS FGC = FACTH * CTHG0 ! DO IK=1, NK !Li New refraction formulation using Cg only. JGLi3Jun2011 ! FRK(IK) = FACTH*2.0*SIG(IK)*(1.-DEPTH*SIG(IK)*SIG(IK)/GRAV) & ! & /(DPH2K(IK)+SNH2K(IK)) !Li Old refraction formulation using phase speed. JGLi8Jun2011 FRK(IK) = FACTH * SIGSNH(IK) !Li FRG(IK) = FGC * CG(IK) END DO ! !Li Current induced refraction stored in FKC. JGLi30Mar2016 IF ( FLCUR ) THEN DO ITH=1, NTH !Li Put a CTMAX limit on current theta rotation. JGLi02Mar2017 ! FKC(ITH) = FACTH*( DCYDX*ES2(ITH) - DCXDY*EC2(ITH) + & FGC = FACTH*( DCYDX*ES2(ITH) - DCXDY*EC2(ITH) + & (DCXDX - DCYDY)*ESC(ITH) ) FKC(ITH) = SIGN( MIN(ABS(FGC), CTMAX), FGC ) END DO ELSE FKC(:)=0.0 END IF ! ! 3.b Depth refraction and great-circle turning. ! DO ITH=1, NTH DDNorm(ITH)=ESIN(ITH)*DDDX-ECOS(ITH)*DDDY DO IK=1, NK ISP = (IK-1)*NTH + ITH !Li Apply depth gradient limited refraction, current and GCT term VCFLT(ISP)=FRG(IK)*ECOS(ITH) + FKC(ITH) + & SIGN( MIN(ABS(FRK(IK)*DDNorm(ITH)), ALFLMT(ITH)), & !Li For new refraction scheme using Cg. JGLi26Jul2011 ! FRK(IK)*DDNorm(ITH) ) !Li For old refraction scheme using phase speed. JGLi26Jul2011 DDNorm(ITH) ) END DO END DO END IF ! ! 4. Wavenumber shift velocities due to current refraction ---------- * ! IF ( FLCK ) THEN ! ! 4.a Directionally dependent part ! DO ITH=1, NTH !Li Depth induced refraction is suspended as it is absorbed in !Li the fixed frequency bin used for wave spectrum. JGLi30Mar2016 ! FKC(ITH) = ( ECOS(ITH)*DDDX + ESIN(ITH)*DDDY ) FKC(ITH) = -DCXDX*EC2(ITH) -DCYDY*ES2(ITH) & -(DCXDY + DCYDX)*ESC(ITH) END DO FKD = CX*DDDX + CY*DDDY ! ! 4.b Band widths ! !Li Cell and side indices for k-dimension are arranged as ! Cell: | -1 | 0 | 1 | 2 | ... | NK | NK+1 | NK+2 | ! Side: -1 0 1 2 ... NK NK+1 !Li DSIP = SIG(K+1) - SIG(K), radian frequency increment ! DO IK=0, NK DB(IK) = DSIP(IK) / CG(IK) DM(IK) = WN(IK+1) - WN(IK) END DO DB(NK+1) = DSIP(NK+1) / CG(NK+1) DM(NK+1) = DM(NK) DM( -1) = DM( 0) ! 4.c Courant number of k-velocity without dividing by dk !!Li FACK = DTG / REAL(NTLOC) ! DO IK=0, NK !Li For new refraction scheme using Cg. JGLi3Jun2011 ! FKS = - FACK*WN(IK)*SIG(IK)/SNH2K(IK) !Li Old refraction formulation using phase speed. JGLi8Jun2011 ! FKS = - FACK*WN(IK)*SIGSNH(IK) !Li Current induced k-shift. JGLi30Mar2016 FKS = MAX( 0.0, CG(IK)*WN(IK)-0.5*SIG(IK) )*FKD / & ( DEPTH30*CG(IK) ) DO ITH=1, NTH CFLK(ITH,IK) = FACK*( FKS + FKC(ITH)*WN(IK) ) END DO END DO !Li No CFL limiter is required here as it is applied in SMCkUNO2. ! END IF ! ! 5. Propagate ------------------------------------------------------ * ! IF ( MOD(ITIME,2) .EQ. 0 ) THEN IF ( FLCK ) THEN !!Li Refraction caused k-space shift. CALL SMCkUNO2(CFLK, VQ, DB, DM) END IF IF ( FLCTH ) THEN !!Li GCT and refraction by rotation in theta direction. CALL SMCGtCrfr(VCFLT, VQ) END IF ELSE IF ( FLCTH ) THEN !!Li GCT and refraction by rotation in theta direction. CALL SMCGtCrfr(VCFLT, VQ) END IF IF ( FLCK ) THEN !!Li Refraction caused k-space shift. CALL SMCkUNO2(CFLK, VQ, DB, DM) END IF END IF ! ! 6. Store reults --------------------------------------------------- * ! VA = VQ ! RETURN ! !/ End of W3KRTN ----------------------------------------------------- / !/ END SUBROUTINE W3KRTN ! Subroutine that calculate mid-flux values for x dimension SUBROUTINE SMCxUNO2(NUA, NUB, CF, UC, UFLX, AKDif, FU, FX, FTS) !!Li CALL SMCxUNO2(iuf, juf, CQ, UCFL, ULCFLX, DNND, FUMD, FUDIFX, FMR) USE CONSTANTS USE W3GDATMD, ONLY: NCel, MRFct, NUFc, IJKCel, IJKUFc, CLATS USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER, INTENT( IN):: NUA, NUB REAL, INTENT( IN):: CF(-9:NCel), UC(-9:NCel), AKDif, FTS REAL, INTENT(Out):: UFLX(NUFc), FU(NUFc), FX(NUFc) ! INTEGER :: i, j, k, L, M, N, ij REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, CNST8, CNST9 ! Two layer of boundary cells are added to each boundary cell face ! with all boundary cell values CF(-9:0)=0.0. ! Diffusion Fourier no. at sub-time-step, proportional to face size, ! which is also equal to the sub-time-step factor FTS. CNST0=AKDif*FTS*FTS ! Uniform diffusion coefficient for all sizes. JGLi24Feb2012 ! CNST0=AKDif*MRFct*FTS !$OMP Parallel Default(Shared), Private(i, ij, K, L, M, N), & !$OMP& Private(CNST,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST8,CNST9) ! Notice an extra side length L is multiplied to mid-flux to give correct ! proportion of flux into the cells. This length will be removed by the ! cell length when the tracer concentration is updated. !$OMP DO DO i=NUA, NUB ! Select Upstream, Central and Downstream cells K=IJKUFc(4,i) L=IJKUFc(5,i) M=IJKUFc(6,i) N=IJKUFc(7,i) ! Face bounding cell lengths and central gradient CNST2=FLOAT( IJKCel(3,L) ) CNST3=FLOAT( IJKCel(3,M) ) CNST5=(CF(M)-CF(L))/( CNST2 + CNST3 ) ! Courant number in local size-1 cell, arithmetic average. CNST6=0.5*( UC(L)+UC(M) )*FTS UFLX(i) = CNST6 ! Multi-resolution SMC grid requires flux multiplied by face factor. CNST8 = FLOAT( IJKUFc(3,i) ) ! Diffusion factor in local size-1 cell, plus the cosine factors. ! 2.0 factor to cancel that in gradient CNST5. JGLi08Mar2012 ! The maximum cell number is used to avoid the boundary cell number ! in selection of the cosine factor. ij= MAX(L, M) CNST9 = 2.0/( CLATS( ij )*CLATS( ij ) ) ! For positive velocity case IF(CNST6 >= 0.0) THEN ! Use central cell velocity for boundary flux. JGLi06Apr2011 IF( M .LE. 0) UFLX(i) = UC(L)*FTS ! Upstream cell length and gradient, depending on UFLX sign. CNST1=FLOAT( IJKCel(3,K) ) CNST4=(CF(L)-CF(K))/( CNST2 + CNST1 ) ! Use minimum gradient all region. CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ! Mid-flux value inside central cell FU(i)=(CF(L) + CNST*(CNST2 - UFLX(i)))*CNST8 ! For negative velocity case ELSE ! Use central cell velocity for boundary flux. JGLi06Apr2011 IF( L .LE. 0) UFLX(i) = UC(M)*FTS ! Upstream cell length and gradient, depending on UFLX sign. CNST1=FLOAT( IJKCel(3,N) ) CNST4=(CF(N)-CF(M))/( CNST1 + CNST3 ) ! Use minimum gradient outside monotonic region. CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ! Mid-flux value inside central cell M FU(i)=(CF(M) - CNST*(CNST3+UFLX(i)))*CNST8 ENDIF ! Diffusion flux by face gradient x DT FX(i)=CNST0*CNST5*CNST8*CNST9 END DO !$OMP END DO !$OMP END Parallel ! 999 PRINT*, ' Sub SMCxUNO2 ended.' RETURN END SUBROUTINE SMCxUNO2 ! Subroutine that calculate mid-flux values for x dimension SUBROUTINE SMCyUNO2(NVA, NVB, CF, VC, VFLY, AKDif, FV, FY, FTS) !!Li CALL SMCyUNO2(ivf, jvf, CQ, VCFL, VLCFLY, DNND, FVMD, FVDIFY, FMR) USE CONSTANTS USE W3GDATMD, ONLY: NCel, MRFct, NVFc, IJKCel, IJKVFc, CLATF USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER, INTENT( IN):: NVA, NVB REAL, INTENT( IN):: CF(-9:NCel), VC(-9:NCel), AKDif, FTS REAL, INTENT(Out):: VFLY(NVFc), FV(NVFc), FY(NVFc) INTEGER :: i, j, k, L, M, N, ij REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, CNST8 ! Notice an extra side length L is multiplied to mid-flux to give correct ! proportion of flux into the cells. This length will be removed by the ! cell length when the tracer concentration is updated. ! Diffusion Fourier no. at sub-time-step, proportional to face size, ! which is also equal to the sub-time-step factor FTS. ! CNST0=AKDif*FTS*FTS ! 2.0 factor to cancel that in gradient CNST5. JGLi08Mar2012 CNST0=AKDif*FTS*FTS*2.0 ! Uniform diffusion coefficient for all sizes. JGLi24Feb2012 ! CNST0=AKDif*MRFct*FTS !$OMP Parallel Default(Shared), Private(j, K, L, M, N), & !$OMP& Private(CNST,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST8) !$OMP DO DO j=NVA, NVB ! Select Upstream, Central and Downstream cells K=IJKVFc(4,j) L=IJKVFc(5,j) M=IJKVFc(6,j) N=IJKVFc(7,j) ! Face bounding cell lengths and gradient CNST2=FLOAT( IJKCel(4,L) ) CNST3=FLOAT( IJKCel(4,M) ) CNST5=(CF(M)-CF(L))/( CNST2 + CNST3 ) ! Courant number in local size-1 cell unit ! Multiply by multi-resolution time step factor FTS CNST6=0.5*( VC(L)+VC(M) )*FTS VFLY(j) = CNST6 ! Face size integer and cosine factor. ! CLATF is defined on V-face for SMC grid. JGLi28Feb2012 CNST8=CLATF(j)*FLOAT( IJKVFc(3,j) ) ! For positive velocity case IF(CNST6 >= 0.0) THEN ! Boundary cell y-size is set equal to central cell y-size ! as y-boundary cell sizes are not proportional to refined ! inner cells but constant of the base cell y-size, and ! Use central cell speed for face speed. JGLi06Apr2011 IF( M .LE. 0 ) THEN VFLY(j) = VC(L)*FTS CNST3 = CNST2 ENDIF ! Upstream cell size and irregular grid gradient, depending on VFLY. CNST1=FLOAT( IJKCel(4,K) ) CNST4=(CF(L)-CF(K))/( CNST2 + CNST1 ) ! Use minimum gradient outside monotonic region CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ! Mid-flux value multiplied by face width and cosine factor FV(j)=( CF(L) + CNST*(CNST2 - VFLY(j)) )*CNST8 ! For negative velocity case ELSE ! Set boundary cell y-size equal to central cell y-size and ! Use central cell speed for flux face speed. JGLi06Apr2011 IF( L .LE. 0 ) THEN VFLY(j) = VC(M)*FTS CNST2 = CNST3 ENDIF ! Upstream cell size and gradient, depending on VFLY sign. ! Side gradients for central cell includs 0.5 factor. CNST1=FLOAT( IJKCel(4,N) ) CNST4=(CF(N)-CF(M))/( CNST1 + CNST3 ) ! Use minimum gradient outside monotonic region CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ! Mid-flux value multiplied by face width and cosine factor FV(j)=( CF(M) - CNST*(CNST3 + VFLY(j)) )*CNST8 ENDIF ! Diffusion flux by face gradient x DT x face_width x cos(lat) ! Multiply by multi-resolution time step factor FTS FY(j)=CNST0*CNST5*CNST8 END DO !$OMP END DO !$OMP END Parallel ! 999 PRINT*, ' Sub SMCyUNO2 ended.' RETURN END SUBROUTINE SMCyUNO2 ! Subroutine that calculate mid-flux values for x dimension SUBROUTINE SMCxUNO2r(NUA, NUB, CF, UC, UFLX, AKDif, FU, FX) !!Li CALL SMCxUNO2r(1, NUFc, CQ, UCFL, ULCFLX, DNND, FUMD, FUDIFX) USE CONSTANTS USE W3GDATMD, ONLY: NSEA, NY, NCel, NUFc, IJKCel, IJKUFc, CLATS USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER, INTENT( IN):: NUA, NUB REAL, INTENT( IN):: CF(-9:NCel), UC(-9:NCel), AKDif REAL, INTENT(Out):: UFLX(NUFc), FU(NUFc), FX(NUFc) ! INTEGER :: i, j, k, L, M, N, ij REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6 ! Two layer of boundary cells are added to each boundary cell face ! with all boundary cell values CF(-9:0)=0.0. ! Notice an extra side length L is multiplied to mid-flux to give correct ! proportion of flux into the cells. This length will be removed by the ! cell length when the tracer concentration is updated. !$OMP Parallel Default(Shared), Private(i, ij, K, L, M, N), & !$OMP& Private(CNST,CNST0,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6) !$OMP DO DO i=NUA, NUB ! Select Upstream, Central and Downstream cells K=IJKUFc(4,i) L=IJKUFc(5,i) M=IJKUFc(6,i) N=IJKUFc(7,i) ! Face bounding cell lengths and gradient CNST2=FLOAT( IJKCel(3,L) ) CNST3=FLOAT( IJKCel(3,M) ) CNST5=(CF(M)-CF(L)) ! Averaged Courant number for base-level cell face CNST6= 0.5*( UC(L)+UC(M) ) UFLX(i) = CNST6 ! Diffusion Fourier number in local cell size ! To avoid boundary cell number, use maximum of L and M. ij= MAX(L, M) CNST0 = 2.0/( CLATS(ij)*CLATS(ij) ) ! For positive velocity case IF(CNST6 >= 0.0) THEN ! Use central cell velocity for boundary flux. JGLi06Apr2011 IF( M .LE. 0) UFLX(i) = UC(L) ! Side gradient for upstream cell as regular grid. CNST4=(CF(L)-CF(K)) ! Use minimum gradient all region with 0.5 factor CNST=Sign(0.5, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ! Mid-flux value inside central cell FU(i)=(CF(L) + CNST*(1.0-UFLX(i)/CNST2)) ! For negative velocity case ELSE ! Use central cell velocity for boundary flux. JGLi06Apr2011 IF( L .LE. 0) UFLX(i) = UC(M) ! Side gradient for upstream cell, depneding on UFLX sign. CNST4=(CF(N)-CF(M)) ! Use minimum gradient outside monotonic region, include 0.5 factor CNST=Sign(0.5, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ! Mid-flux value inside central cell M FU(i)=(CF(M) - CNST*(1.0+UFLX(i)/CNST3)) ENDIF ! Diffusion flux by face gradient x DT FX(i)=AKDif*CNST0*CNST5/(CNST2 + CNST3) END DO !$OMP END DO !$OMP END Parallel ! 999 PRINT*, ' Sub SMCxUNO2r ended.' RETURN END SUBROUTINE SMCxUNO2r ! Subroutine that calculate mid-flux values for y dimension SUBROUTINE SMCyUNO2r(NVA, NVB, CF, VC, VFLY, AKDif, FV, FY) !!Li CALL SMCyUNO2r(1, NVFc, CQ, VCFL, VLCFLY, DNND, FVMD, FVDIFY) USE CONSTANTS USE W3GDATMD, ONLY: NSEA, NY, NCel, NVFc, IJKCel, IJKVFc, CLATF USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER, INTENT( IN):: NVA, NVB REAL, INTENT( IN):: CF(-9:NCel), VC(-9:NCel), AKDif REAL, INTENT(Out):: VFLY(NVFc), FV(NVFc), FY(NVFc) INTEGER :: i, j, k, L, M, N, ij REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, CNST8 ! Notice an extra side length L is multiplied to mid-flux to give correct ! proportion of flux into the cells. This length will be removed by the ! cell length when the tracer concentration is updated. !$OMP Parallel Default(Shared), Private(j, K, L, M, N), & !$OMP& Private(CNST,CNST4,CNST5,CNST6,CNST8) !$OMP DO DO j=NVA, NVB ! Select Upstream, Central and Downstream cells K=IJKVFc(4,j) L=IJKVFc(5,j) M=IJKVFc(6,j) N=IJKVFc(7,j) ! Central face gradient. CNST5=(CF(M)-CF(L)) ! Courant number in basic cell unit as dy is constant CNST6=0.5*( VC(L)+VC(M) ) VFLY(j) = CNST6 ! Face size integer and cosine factor ! CLATF is defined on V-face for SMC grid. JGLi28Feb2012 CNST8=CLatF(j)*FLOAT( IJKVFc(3,j) ) ! For positive velocity case IF(CNST6 >= 0.0) THEN ! Use central cell speed for flux face speed. JGLi06Apr2011 IF( M .LE. 0 ) VFLY(j) = VC(L) ! Upstream face gradient, depending on VFLY sign. CNST4=(CF(L)-CF(K)) ! Use minimum gradient, including 0.5 factor and central sign. CNST=Sign(0.5, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ! Mid-flux value multiplied by face width and cosine factor FV(j)=( CF(L) + CNST*(1.0-VFLY(j)) )*CNST8 ! For negative velocity case ELSE ! Use central cell speed for flux face speed. JGLi06Apr2011 IF( L .LE. 0 ) VFLY(j) = VC(M) ! Side gradients for upstream face, depending on VFLY sign. CNST4=(CF(N)-CF(M)) ! Use minimum gradient, including 0.5 factor and central sign. CNST=Sign(0.5, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ! Mid-flux value multiplied by face width and cosine factor FV(j)=( CF(M) - CNST*(1.0+VFLY(j)) )*CNST8 ENDIF ! Diffusion flux by face gradient x DT x face_width x cos(lat) FY(j)=AKDif*CNST5*CNST8 END DO !$OMP END DO !$OMP END Parallel ! 999 PRINT*, ' Sub SMCyUNO2r ended.' RETURN END SUBROUTINE SMCyUNO2r ! Subroutine that calculate mid-flux values for x dimension with UNO3 scheme SUBROUTINE SMCxUNO3(NUA, NUB, CF, UC, UFLX, AKDif, FU, FX, FTS) !!Li CALL SMCxUNO3(iuf, juf, CQ, UCFL, ULCFLX, DNND, FUMD, FUDIFX, FMR) USE CONSTANTS USE W3GDATMD, ONLY: NCel, MRFct, NUFc, IJKCel, IJKUFc, CLATS USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER, INTENT( IN):: NUA, NUB REAL, INTENT( IN):: CF(-9:NCel), UC(-9:NCel), AKDif, FTS REAL, INTENT(Out):: UFLX(NUFc), FU(NUFc), FX(NUFc) ! INTEGER :: i, j, k, L, M, N, ij REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, & & CNST7, CNST8, CNST9 ! Two layer of boundary cells are added to each boundary cell face ! with all boundary cell values CF(-9:0)=0.0. ! Diffusion Fourier no. at sub-time-step, proportional to face size, ! which is also equal to the sub-time-step factor FTS. ! CNST0=AKDif*FTS*FTS ! 2.0 factor to cancel that in gradient CNST5. JGLi03Sep2015 CNST0=AKDif*FTS*FTS*2.0 ! Notice an extra side length L is multiplied to mid-flux to give correct ! proportion of flux into the cells. This length will be removed by the ! cell length when the tracer concentration is updated. !$OMP Parallel Default(Shared), Private(i, ij, K, L, M, N), & !$OMP& Private(CNST,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST7,CNST8,CNST9) !$OMP DO DO i=NUA, NUB ! Select Upstream, Central and Downstream cells K=IJKUFc(4,i) L=IJKUFc(5,i) M=IJKUFc(6,i) N=IJKUFc(7,i) ! Face bounding cell lengths and central gradient CNST2=FLOAT( IJKCel(3,L) ) CNST3=FLOAT( IJKCel(3,M) ) CNST5=(CF(M)-CF(L))/( CNST2 + CNST3 ) ! Courant number in local size-1 cell, arithmetic average. CNST6=0.5*( UC(L)+UC(M) )*FTS UFLX(i) = CNST6 ! Multi-resolution SMC grid requires flux multiplied by face factor. CNST8 = FLOAT( IJKUFc(3,i) ) ! Diffusion factor in local size-1 cell, plus the cosine factors. ! 2.0 factor to cancel that in gradient CNST5. JGLi08Mar2012 ! The maximum cell number is used to avoid the boundary cell number ! in selection of the cosine factor. ij= MAX(L, M) ! For positive velocity case IF(CNST6 >= 0.0) THEN ! Use central cell velocity for boundary flux. JGLi06Apr2011 IF( M .LE. 0) UFLX(i) = UC(L)*FTS ! Upstream cell length and gradient, depending on UFLX sign. CNST1=FLOAT( IJKCel(3,K) ) CNST4=(CF(L)-CF(K))/( CNST2 + CNST1 ) ! Second order gradient CNST7 = CNST5 - CNST4 CNST9 = 2.0/( CNST3+CNST2+CNST2+CNST1 ) ! Use 3rd order scheme IF( Abs(CNST7) .LT. 0.6*CNST9*Abs(CF(M)-CF(K)) ) THEN CNST= CNST5 - ( CNST3+UFLX(i) )*CNST7*CNST9/1.5 ! Use doubled UNO2 scheme ELSE IF( DBLE(CNST4)*DBLE(CNST5) .GT. 0.d0 ) THEN CNST=Sign(2.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ELSE ! Use minimum gradient UNO2 scheme CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ENDIF ! Mid-flux value inside central cell FU(i)=(CF(L) + CNST*(CNST2 - UFLX(i)))*CNST8 ! For negative velocity case ELSE ! Use central cell velocity for boundary flux. JGLi06Apr2011 IF( L .LE. 0) UFLX(i) = UC(M)*FTS ! Upstream cell length and gradient, depending on UFLX sign. CNST1=FLOAT( IJKCel(3,N) ) CNST4=(CF(N)-CF(M))/( CNST1 + CNST3 ) ! Second order gradient CNST7 = CNST4 - CNST5 CNST9 = 2.0/( CNST2+CNST3+CNST3+CNST1 ) ! Use 3rd order scheme IF( Abs(CNST7) .LT. 0.6*CNST9*Abs(CF(N)-CF(L)) ) THEN CNST= CNST5 + ( CNST2-UFLX(i) )*CNST7*CNST9/1.5 ! Use doubled UNO2 scheme ELSE IF( DBLE(CNST4)*DBLE(CNST5) .GT. 0.d0 ) THEN CNST=Sign(2.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ELSE ! Use minimum gradient UNO2 scheme. CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ENDIF ! Mid-flux value inside central cell M FU(i)=(CF(M) - CNST*(CNST3+UFLX(i)))*CNST8 ENDIF ! Diffusion flux by face gradient x DT FX(i)=CNST0*CNST5*CNST8/( CLATS( ij )*CLATS( ij ) ) END DO !$OMP END DO !$OMP END Parallel ! 999 PRINT*, ' Sub SMCxUNO3 ended.' RETURN END SUBROUTINE SMCxUNO3 ! Subroutine that calculate mid-flux values for y dimension with UNO3 scheme SUBROUTINE SMCyUNO3(NVA, NVB, CF, VC, VFLY, AKDif, FV, FY, FTS) !!Li CALL SMCyUNO3(ivf, jvf, CQ, VCFL, VLCFLY, DNND, FVMD, FVDIFY, FMR) USE CONSTANTS USE W3GDATMD, ONLY: NCel, MRFct, NVFc, IJKCel, IJKVFc, CLATF USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER, INTENT( IN):: NVA, NVB REAL, INTENT( IN):: CF(-9:NCel), VC(-9:NCel), AKDif, FTS REAL, INTENT(Out):: VFLY(NVFc), FV(NVFc), FY(NVFc) INTEGER :: i, j, k, L, M, N, ij REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, & & CNST7, CNST8, CNST9 ! Notice an extra side length L is multiplied to mid-flux to give correct ! proportion of flux into the cells. This length will be removed by the ! cell length when the tracer concentration is updated. ! Diffusion Fourier no. at sub-time-step, proportional to face size, ! which is also equal to the sub-time-step factor FTS. ! CNST0=AKDif*FTS*FTS ! 2.0 factor to cancel that in gradient CNST5. JGLi08Mar2012 CNST0=AKDif*FTS*FTS*2.0 ! Uniform diffusion coefficient for all sizes. JGLi24Feb2012 ! CNST0=AKDif*MRFct*FTS !$OMP Parallel Default(Shared), Private(j, K, L, M, N), & !$OMP& Private(CNST,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST7,CNST8,CNST9) !$OMP DO DO j=NVA, NVB ! Select Upstream, Central and Downstream cells K=IJKVFc(4,j) L=IJKVFc(5,j) M=IJKVFc(6,j) N=IJKVFc(7,j) ! Face bounding cell lengths and gradient CNST2=FLOAT( IJKCel(4,L) ) CNST3=FLOAT( IJKCel(4,M) ) CNST5=(CF(M)-CF(L))/( CNST2 + CNST3 ) ! Courant number in local size-1 cell unit ! Multiply by multi-resolution time step factor FTS CNST6=0.5*( VC(L)+VC(M) )*FTS VFLY(j) = CNST6 ! Face size integer and cosine factor. ! CLATF is defined on V-face for SMC grid. JGLi28Feb2012 CNST8=CLATF(j)*FLOAT( IJKVFc(3,j) ) ! For positive velocity case IF(CNST6 >= 0.0) THEN ! Boundary cell y-size is set equal to central cell y-size ! as y-boundary cell sizes are not proportional to refined ! inner cells but constant of the base cell y-size, and ! Use central cell speed for face speed. JGLi06Apr2011 IF( M .LE. 0 ) THEN VFLY(j) = VC(L)*FTS CNST3 = CNST2 ENDIF ! Upstream cell size and irregular grid gradient, depending on VFLY. CNST1=FLOAT( IJKCel(4,K) ) CNST4=(CF(L)-CF(K))/( CNST2 + CNST1 ) ! Second order gradient CNST7 = CNST5 - CNST4 CNST9 = 2.0/( CNST3+CNST2+CNST2+CNST1 ) ! Use 3rd order scheme IF( Abs(CNST7) .LT. 0.6*CNST9*Abs(CF(M)-CF(K)) ) THEN CNST= CNST5 - ( CNST3+VFLY(j) )*CNST7*CNST9/1.5 ! Use doubled UNO2 scheme ELSE IF( DBLE(CNST4)*DBLE(CNST5) .GT. 0.d0 ) THEN CNST=Sign(2.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ELSE ! Use minimum gradient outside monotonic region CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ENDIF ! Mid-flux value multiplied by face width and cosine factor FV(j)=( CF(L) + CNST*(CNST2 - VFLY(j)) )*CNST8 ! For negative velocity case ELSE ! Set boundary cell y-size equal to central cell y-size and ! Use central cell speed for flux face speed. JGLi06Apr2011 IF( L .LE. 0 ) THEN VFLY(j) = VC(M)*FTS CNST2 = CNST3 ENDIF ! Upstream cell size and gradient, depending on VFLY sign. ! Side gradients for central cell includs 0.5 factor. CNST1=FLOAT( IJKCel(4,N) ) CNST4=(CF(N)-CF(M))/( CNST1 + CNST3 ) ! Second order gradient CNST7 = CNST4 - CNST5 CNST9 = 2.0/( CNST2+CNST3+CNST3+CNST1 ) ! Use 3rd order scheme IF( Abs(CNST7) .LT. 0.6*CNST9*Abs(CF(N)-CF(L)) ) THEN CNST= CNST5 + ( CNST2-VFLY(j) )*CNST7*CNST9/1.5 ! Use doubled UNO2 scheme ELSE IF( DBLE(CNST4)*DBLE(CNST5) .GT. 0.d0 ) THEN CNST=Sign(2.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ELSE ! Use minimum gradient outside monotonic region CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ENDIF ! Mid-flux value multiplied by face width and cosine factor FV(j)=( CF(M) - CNST*(CNST3 + VFLY(j)) )*CNST8 ENDIF ! Diffusion flux by face gradient x DT x face_width x cos(lat) ! Multiply by multi-resolution time step factor FTS FY(j)=CNST0*CNST5*CNST8 END DO !$OMP END DO !$OMP END Parallel ! 999 PRINT*, ' Sub SMCyUNO3 ended.' RETURN END SUBROUTINE SMCyUNO3 ! Subroutine that calculate mid-flux values for x dimension with UNO3 SUBROUTINE SMCxUNO3r(NUA, NUB, CF, UC, UFLX, AKDif, FU, FX) !!Li CALL SMCxUNO3r(1, NUFc, CQ, UCFL, ULCFLX, DNND, FUMD, FUDIFX) USE CONSTANTS USE W3GDATMD, ONLY: NSEA, NY, NCel, NUFc, IJKCel, IJKUFc, CLATS USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER, INTENT( IN):: NUA, NUB REAL, INTENT( IN):: CF(-9:NCel), UC(-9:NCel), AKDif REAL, INTENT(Out):: UFLX(NUFc), FU(NUFc), FX(NUFc) ! INTEGER :: i, j, k, L, M, N, ij REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, & CNST7, CNST8, CNST9 ! Two layer of boundary cells are added to each boundary cell face ! with all boundary cell values CF(-9:0)=0.0. ! Notice an extra side length L is multiplied to mid-flux to give correct ! proportion of flux into the cells. This length will be removed by the ! cell length when the tracer concentration is updated. !$OMP Parallel Default(Shared), Private(i, ij, K, L, M, N), & !$OMP& Private(CNST,CNST0,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST7,CNST8,CNST9) !$OMP DO DO i=NUA, NUB ! Select Upstream, Central and Downstream cells K=IJKUFc(4,i) L=IJKUFc(5,i) M=IJKUFc(6,i) N=IJKUFc(7,i) ! Face bounding cell lengths and gradient CNST2=FLOAT( IJKCel(3,L) ) CNST3=FLOAT( IJKCel(3,M) ) CNST5=(CF(M)-CF(L)) ! Averaged Courant number for base-level cell face CNST6= 0.5*( UC(L)+UC(M) ) UFLX(i) = CNST6 ! Diffusion Fourier number in local cell size ! To avoid boundary cell number, use maximum of L and M. ij= MAX(L, M) CNST0 = 2.0/( CLATS(ij)*CLATS(ij) ) ! For positive velocity case IF(CNST6 >= 0.0) THEN ! Use central cell velocity for boundary flux. JGLi06Apr2011 IF( M .LE. 0) UFLX(i) = UC(L) ! Side gradient for upstream cell as regular grid. CNST4=(CF(L)-CF(K)) CNST8=(CF(M)-CF(K)) CNST9=(CNST5-CNST4) IF( Abs(CNST9) <= 0.6*Abs(CNST8) ) THEN ! Use 3rd order scheme in limited zone, note division by 2 grid sizes CNST=0.5*CNST5 - (1.0+UFLX(i)/CNST2)*CNST9/6.0 ELSE IF( DBLE(CNST4)*DBLE(CNST5) .GT. 0.d0 ) THEN ! Use doubled minimum gradient in rest of monotonic region CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ELSE ! Use minimum gradient all region with 0.5 factor CNST=Sign(0.5, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ENDIF ! Mid-flux value inside central cell FU(i)=(CF(L) + CNST*(1.0-UFLX(i)/CNST2)) ! For negative velocity case ELSE ! Use central cell velocity for boundary flux. JGLi06Apr2011 IF( L .LE. 0) UFLX(i) = UC(M) ! Side gradient for upstream cell, depneding on UFLX sign. CNST4=(CF(N)-CF(M)) CNST8=(CF(N)-CF(L)) CNST9=(CNST4-CNST5) IF( Abs(CNST9) <= 0.6*Abs(CNST8) ) THEN ! Use 3rd order scheme in limited zone, note division by 2 grid sizes CNST=0.5*CNST5 + (1.0-UFLX(i)/CNST3)*CNST9/6.0 ELSE IF( DBLE(CNST4)*DBLE(CNST5) .GT. 0.d0 ) THEN ! Use doubled minimum gradient in rest of monotonic region CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ELSE ! Use minimum gradient outside monotonic region, include 0.5 factor CNST=Sign(0.5, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ENDIF ! Mid-flux value inside central cell M FU(i)=(CF(M) - CNST*(1.0+UFLX(i)/CNST3)) ENDIF ! Diffusion flux by face gradient x DT FX(i)=AKDif*CNST0*CNST5/(CNST2 + CNST3) END DO !$OMP END DO !$OMP END Parallel ! 999 PRINT*, ' Sub SMCxUNO3r ended.' RETURN END SUBROUTINE SMCxUNO3r ! Subroutine that calculate mid-flux values for y dimension with UNO3 SUBROUTINE SMCyUNO3r(NVA, NVB, CF, VC, VFLY, AKDif, FV, FY) !!Li CALL SMCyUNO3r(1, NVFc, CQ, VCFL, VLCFLY, DNND, FVMD, FVDIFY) USE CONSTANTS USE W3GDATMD, ONLY: NSEA, NY, NCel, NVFc, IJKCel, IJKVFc, CLATF USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER, INTENT( IN):: NVA, NVB REAL, INTENT( IN):: CF(-9:NCel), VC(-9:NCel), AKDif REAL, INTENT(Out):: VFLY(NVFc), FV(NVFc), FY(NVFc) INTEGER :: i, j, k, L, M, N, ij REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, & & CNST7, CNST8, CNST9 ! Notice an extra side length L is multiplied to mid-flux to give correct ! proportion of flux into the cells. This length will be removed by the ! cell length when the tracer concentration is updated. !$OMP Parallel Default(Shared), Private(j, K, L, M, N), & !$OMP& Private(CNST,CNST4,CNST5,CNST6,CNST7,CNST8,CNST9) !$OMP DO DO j=NVA, NVB ! Select Upstream, Central and Downstream cells K=IJKVFc(4,j) L=IJKVFc(5,j) M=IJKVFc(6,j) N=IJKVFc(7,j) ! Central face gradient. CNST5=(CF(M)-CF(L)) ! Courant number in basic cell unit as dy is constant CNST6=0.5*( VC(L)+VC(M) ) VFLY(j) = CNST6 ! Face size integer and cosine factor ! CLATF is defined on V-face for SMC grid. JGLi28Feb2012 CNST7=CLatF(j)*FLOAT( IJKVFc(3,j) ) ! For positive velocity case IF(CNST6 >= 0.0) THEN ! Use central cell speed for flux face speed. JGLi06Apr2011 IF( M .LE. 0 ) VFLY(j) = VC(L) ! Upstream face gradient, depending on VFLY sign. CNST4=(CF(L)-CF(K)) ! Second gradient for 3rd scheme CNST8=(CF(M)-CF(K)) CNST9=(CNST5-CNST4) IF( Abs(CNST9) <= 0.6*Abs(CNST8) ) THEN ! Use 3rd order scheme in limited zone, note division by 2 grid sizes CNST=0.5*CNST5-(1.0+VFLY(j))*CNST9/6.0 ELSE IF( DBLE(CNST4)*DBLE(CNST5) .GT. 0.d0 ) THEN ! Use doubled minimum gradient in rest of monotonic region CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ELSE ! Use minimum gradient, including 0.5 factor and central sign. CNST=Sign(0.5, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ENDIF ! Mid-flux value multiplied by face width and cosine factor FV(j)=( CF(L) + CNST*(1.0-VFLY(j)) )*CNST7 ! For negative velocity case ELSE ! Use central cell speed for flux face speed. JGLi06Apr2011 IF( L .LE. 0 ) VFLY(j) = VC(M) ! Side gradients for upstream face, depending on VFLY sign. CNST4=(CF(N)-CF(M)) ! Second gradient for 3rd scheme CNST8=(CF(N)-CF(L)) CNST9=(CNST4-CNST5) IF( Abs(CNST9) <= 0.6*Abs(CNST8) ) THEN ! Use 3rd order scheme in limited zone, note division by 2 grid sizes CNST=0.5*CNST5+(1.0-VFLY(j))*CNST9/6.0 ELSE IF( DBLE(CNST4)*DBLE(CNST5) .GT. 0.d0 ) THEN ! Use doubled minimum gradient in rest of monotonic region CNST=Sign(1.0, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ELSE ! Use minimum gradient, including 0.5 factor and central sign. CNST=Sign(0.5, CNST5)*min( Abs(CNST4), Abs(CNST5) ) ENDIF ! Mid-flux value multiplied by face width and cosine factor FV(j)=( CF(M) - CNST*(1.0+VFLY(j)) )*CNST7 ENDIF ! Diffusion flux by face gradient x DT x face_width x cos(lat) FY(j)=AKDif*CNST5*CNST7 END DO !$OMP END DO !$OMP END Parallel ! 999 PRINT*, ' Sub SMCyUNO3r ended.' RETURN END SUBROUTINE SMCyUNO3r ! ! Subroutine that calculate cell centre gradient for any input variable. ! Nemerical average is applied to size-changing faces and the gradients ! are along the lat-lon local east-north directions. JGLi18Aug2015 ! Add optional zero-gradient bounday conditions. JGLi08Aug2017 ! SUBROUTINE SMCGradn(CVQ, GrdX, GrdY, L0r1) USE CONSTANTS USE W3GDATMD, ONLY: NSEA, NUFc, NVFc, MRFct, & & IJKCel, IJKUFc, IJKVFC, CLATS, SX, SY !/ARC USE W3GDATMD, ONLY: NGLO USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE !! New boundary conditions depending on user defined L0r1. !! L0r1 = 0 will set zero at land points while L0r1 > 0 invokes !! the zero-gradient boundary condition. JGLi08Aug2017 REAL, INTENT( IN):: CVQ(NSEA) REAL, INTENT(Out):: GrdX(NSEA), GrdY(NSEA) INTEGER, INTENT( IN):: L0r1 ! INTEGER :: I, J, K, L, M, N REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6 REAL :: DX0I, DY0I ! Use a few working arrays REAL, Dimension(-9:NSEA):: CVF, AUN, AVN ! Two layer of boundary cells are added to each boundary cell face ! with all boundary cell default values CVF(-9:0)= 0.0. CVF(-9:0) = 0.0 CVF(1:NSEA)=CVQ(1:NSEA) !! Initialize arrays AUN = 0. AVN = 0. GrdX = 0. GrdY = 0. !! Multi-resolution base-cell size defined by refined levels. !! So the MRFct converts the base cell SX, SY into size-1 cell lenth. !! Constant size-1 dy=DY0 and dx on Equator DX0, inverted. DX0I = MRFct/ ( SX * DERA * RADIUS ) DY0I = MRFct/ ( SY * DERA * RADIUS ) !$OMP Parallel Default(Shared), Private(i, j, K, L, M, N), & !$OMP& Private(CNST,CNST0,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6) !$OMP DO !! Calculate x-gradient by averaging U-face gradients. DO i=1, NUFc ! Select Upstream, Central and Downstream cells L=IJKUFc(5,i) M=IJKUFc(6,i) !! For zero-gradient boundary conditions, simply skip boundary faces. IF( L0r1 .EQ. 0 .OR. (L > 0 .AND. M > 0) ) THEN ! Multi-resolution SMC grid requires flux multiplied by face factor. CNST1=FLOAT( IJKUFc(3,i) ) ! Face bounding cell lengths and central gradient CNST2=FLOAT( IJKCel(3,L) ) CNST3=FLOAT( IJKCel(3,M) ) ! Side gradients over 2 cell lengths for central cell. ! Face size factor is also included for average. CNST5=CNST1*(CVF(M)-CVF(L))/(CNST2+CNST3) !$OMP CRITICAL ! Store side gradient in two neighbouring cells AUN(L) = AUN(L) + CNST5 AUN(M) = AUN(M) + CNST5 !$OMP END CRITICAL ENDIF END DO !$OMP END DO ! Assign averaged side-gradient to GrdX, plus latitude factor ! Note averaging over 2 times of cell y-width factor but AUN ! has already been divied by two cell lengths. !$OMP DO DO n=1, NSEA ! Cell y-size IJKCel(4,i) is used to cancel the face size-factor in AUN. ! Plus the actual physical length scale for size-1 cell. ! Note polar cell (if any) AUN = 0.0 as it has no U-face. GrdX(n)=DX0I*AUN(n)/( CLats(n)*IJKCel(4,n) ) ENDDO !$OMP END DO !$OMP DO !! Calculate y-gradient by averaging V-face gradients. DO j=1, NVFc ! Select Central and Downstream cells L=IJKVFc(5,j) M=IJKVFc(6,j) !! For zero-gradient boundary conditions, simply skip boundary faces. IF( L0r1 .EQ. 0 .OR. (L > 0 .AND. M > 0) ) THEN ! Face size is required for multi-resolution grid. CNST1=Real( IJKVFc(3,j) ) ! Cell y-length of UCD cells CNST2=Real( IJKCel(4,L) ) CNST3=Real( IJKCel(4,M) ) ! Side gradients over 2 cell lengths for central cell. ! Face size factor is also included for average. CNST6=CNST1*(CVF(M)-CVF(L))/(CNST2+CNST3) !$OMP CRITICAL ! Store side gradient in two neighbouring cells AVN(L) = AVN(L) + CNST6 AVN(M) = AVN(M) + CNST6 !$OMP END CRITICAL ENDIF END DO !$OMP END DO !$OMP DO ! Assign averaged side-gradient to GrdY. DO n=1, NSEA ! AV is divided by the cell x-size IJKCel(3,i) to cancel face ! size-factor, and physical y-distance of size-1 cell. GrdY(n)=DY0I*AVN(n)/Real( IJKCel(3,n) ) END DO !$OMP END DO !$OMP END Parallel !/ARC !!Li Polar cell (if any) y-gradient is set to zero. !/ARC IF( NSEA .GT. NGLO ) GrdY(NSEA) = 0.0 ! 999 PRINT*, ' Sub SMCGradn ended.' RETURN END SUBROUTINE SMCGradn ! ! Subroutine that average sea point values with a 1-2-1 scheme. ! SUBROUTINE SMCAverg(CVQ) USE CONSTANTS USE W3GDATMD, ONLY: NSEA, NUFc, NVFc, & & IJKCel, IJKUFc, IJKVFC !/ARC USE W3GDATMD, ONLY: NGLO USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE REAL, INTENT(INOUT):: CVQ(-9:NSEA) ! INTEGER :: I, J, K, L, M, N REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6 ! Use a few working arrays REAL, Dimension(-9:NSEA):: CVF, AUN, AVN ! Two layer of boundary cells are added to each boundary cell face ! with all boundary cell values stored in CVF(-9:0). CVF=CVQ !! Initialize arrays AUN = 0. AVN = 0. !/ARC !!Li Save polar cell value !/ARC CNST0 = CVQ(NSEA) !$OMP Parallel Default(Shared), Private(i, j, L, M, n), & !$OMP& Private(CNST3,CNST4,CNST5,CNST6) !$OMP DO !! Calculate x-gradient by averaging U-face gradients. DO i=1, NUFc ! Select Upstream, Central and Downstream cells L=IJKUFc(5,i) M=IJKUFc(6,i) ! Multi-resolution SMC grid requires flux multiplied by face factor. CNST5=Real( IJKUFc(3,i) )*(CVF(M)+CVF(L)) !$OMP CRITICAL ! Store side gradient in two neighbouring cells AUN(L) = AUN(L) + CNST5 AUN(M) = AUN(M) + CNST5 !$OMP END CRITICAL END DO !$OMP END DO !$OMP DO !! Calculate y-gradient by averaging V-face gradients. DO j=1, NVFc ! Select Central and Downstream cells L=IJKVFc(5,j) M=IJKVFc(6,j) ! Face size is required for multi-resolution grid. CNST6=Real( IJKVfc(3,j) )*(CVF(M)+CVF(L)) !$OMP CRITICAL ! Store side gradient in two neighbouring cells AVN(L) = AVN(L) + CNST6 AVN(M) = AVN(M) + CNST6 !$OMP END CRITICAL END DO !$OMP END DO !$OMP DO ! Assign averaged value back to CVQ. DO n=1, NSEA CNST3=0.125/Real( IJKCel(3,n) ) CNST4=0.125/Real( IJKCel(4,n) ) ! AUN is divided by the cell y-size IJKCel(4,n) and AVN by ! the cell x-size IJKCel(3,n) to cancel face size factors. CVQ(n)= AUN(n)*CNST4 + AVN(n)*CNST3 END DO !$OMP END DO !$OMP END Parallel !/ARC !!Li Polar cell (if any) keep original value. !/ARC IF( NSEA .GT. NGLO ) CVQ(NSEA) = CNST0 ! 999 PRINT*, ' Sub SMCAverg ended.' RETURN END SUBROUTINE SMCAverg !Li ! Subroutine that calculate great circle turning (GCT) and refraction. ! The refraction and GCT terms are equivalent to a simgle rotation by each ! element and does not need to be calculated as advection. A simple rotation ! scheme similar to the 1st order upstream scheme but without any restriction ! on the rotation angle or the CFL limit by an Eulerian advection scheme. ! Jian-Guo Li 12 Nov 2010 !Li SUBROUTINE SMCGtCrfr(CoRfr, SpeTHK) USE CONSTANTS USE W3GDATMD, ONLY: NK, NTH, DTH, CTMAX IMPLICIT NONE REAL, INTENT(IN) :: CoRfr(NTH, NK) REAL, INTENT(INOUT):: SpeTHK(NTH, NK) INTEGER :: I, J, K, L, M, N REAL, Dimension(NTH):: SpeGCT, Spectr REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6 ! Loop through NK spectral bins. DO n=1, NK !! Asign cell spectrum to temporary variable Spcetr Spectr=SpeTHK(1:NTH,n) SpeGCT=0.0 !! Loop through NTH directional bins for each cell spectrum DO j=1, NTH ! GCT + refraction Courant number for this dirctional bin CNST6=CoRfr(j,n) ! Work out integer number of bins to be skipped. ! If K is great than NTH, full circle turning is removed. CNST5=ABS( CNST6 ) K= MOD( INT(CNST5), NTH ) ! Partitioning faraction of the spectral component CNST1=CNST5 - FLOAT( INT(CNST5) ) CNST2=1.0 - CNST1 ! For positive turning case IF(CNST6 > 0.0) THEN ! Select the upstream and downstream bins to rotate in, wrap at end L=j+K M=j+K+1 IF( L .GT. NTH ) L = L - NTH IF( M .GT. NTH ) M = M - NTH !! Divid the j bin energy by fraction of CNST6 and store in SpeGCT SpeGCT(L)=SpeGCT(L)+Spectr(j)*CNST2 SpeGCT(M)=SpeGCT(M)+Spectr(j)*CNST1 ! For negative or no turning case ELSE ! Select the upstream and downstream bins to rotate in, wrap at end L=j-K M=j-K-1 IF( L .LT. 1 ) L = L + NTH IF( M .LT. 1 ) M = M + NTH !! Divid the bin energy by fraction of CNST6 and store in SpeGCT SpeGCT(L)=SpeGCT(L)+Spectr(j)*CNST2 SpeGCT(M)=SpeGCT(M)+Spectr(j)*CNST1 ENDIF !! End of directional loop j END DO !! Store GCT spectrum SpeTHK(1:NTH,n) = SpeGCT !! End of cell loop n END DO ! 999 PRINT*, ' Sub SMCGtCrfr ended.' RETURN END SUBROUTINE SMCGtCrfr !Li ! Subroutine that calculates refraction induced shift in k-space. ! The term is equivalent to advection on an irregular k-space grid. ! The UNO2 scheme on irregular grid is used for this term. ! Jian-Guo Li 15 Nov 2010 ! Fix bug on CFL limiter and add positive filter. JGLi28Jun2017 !Li SUBROUTINE SMCkUNO2(CoRfr, SpeTHK, DKC, DKS) !!Li CALL SMCkUNO2(CFLK, VQ, DB, DM) USE CONSTANTS USE W3GDATMD, ONLY: NK, NK2, NTH, DTH, XFR, CTMAX IMPLICIT NONE REAL, INTENT(IN) :: CoRfr(NTH, 0:NK), DKC(0:NK+1), DKS(-1:NK+1) REAL, INTENT(INOUT):: SpeTHK(NTH, NK) INTEGER :: I, J, K, L, M, N REAL, Dimension(-1:NK+2):: SpeRfr, Spectr, SpeFlx REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6 !Li Cell and side indices for k-dimension are arranged as ! Cell: | -1 | 0 | 1 | 2 | ... | NK | NK+1 | NK+2 | ! Side: -1 0 1 2 ... NK NK+1 ! The wave action in k-space is extended at the high-wavenumber (frequency) end ! by the (m+2)th negative power of frequency for boundary conditions. Outside ! low-wavenumber (frequncy) end, wave action is assumed to be zero. CNST=XFR**(-7) DO n=1, NTH !! Asign cell spectrum to temporary variable Spcetr Spectr(-1) =0.0 Spectr( 0) =0.0 Spectr(1:NK)=SpeTHK(n,1:NK) Spectr(NK+1)=Spectr(NK )*CNST Spectr(NK+2)=Spectr(NK+1)*CNST !! Calculate k-space gradient for NK+2 faces by UNO2 scheme SpeRfr(-1)= 0.0 DO j=-1, NK+1 SpeRfr(j)=(Spectr(j+1)-Spectr(j))/DKS(j) ENDDO !! Calculate k-space fluxes for NK+1 faces by UNO2 scheme DO j=0, NK ! Final refraction Courant number for this k-space face CNST6=CoRfr(n,j) !! Note CoRfr is CFL for k but without dividing dk. ! For positive case IF(CNST6 > 0.0) THEN CNST0 = MIN( CTMAX*DKC(j), CNST6 ) SpeFlx(j) = CNST0*( Spectr(j) + SIGN(0.5, SpeRfr(j))*(DKC(j)-CNST0) & *MIN( ABS(SpeRfr(j-1)), ABS(SpeRfr(j)) ) ) ! For negative or no turning case ELSE CNST0 = MIN( CTMAX*DKC(j+1), -CNST6 ) SpeFlx(j) = -CNST0*( Spectr(j+1) - SIGN(0.5, SpeRfr(j))*(DKC(j+1)-CNST0) & *MIN( ABS(SpeRfr(j+1)), ABS(SpeRfr(j)) ) ) ENDIF !! End of flux loop j END DO !! Update spectrum for the given direction DO j=1, NK ! Final refraction Courant number for this k-space face ! SpeTHK(n, j) = Spectr(j) + (SpeFlx(j-1) - SpeFlx(j))/DKC(j) ! Add positive filter in case negative values. JGLi27Jun2017 SpeTHK(n, j) = MAX( 0.0, Spectr(j)+(SpeFlx(j-1)-SpeFlx(j))/DKC(j) ) END DO !! End of directional loop n END DO ! 999 PRINT*, ' Sub SMCkUNO2 ended.' RETURN END SUBROUTINE SMCkUNO2 ! Subroutine that calculates water-depth gradient for refraction. ! For consistency with the lat-lon grid, full grid DDDX, DDDY are ! also assigned here. DHDX, DHDY are used for refraction at present. ! It has to be rotated to map-east system in the Arctic part. SUBROUTINE SMCDHXY USE CONSTANTS USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSTA, MAPFS, MRFct, IJKCel, & CLATS, NTH, DTH, ESIN, ECOS, Refran, DMIN !/ARC USE W3GDATMD, ONLY: NGLO, ANGARC USE W3ADATMD, ONLY: DW, DDDX, DDDY, DHDX, DHDY, DHLMT USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER :: I, J, K, L, M, N REAL :: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6 REAL, Dimension(NSEA) :: HCel, GrHx, GrHy ! REAL, Dimension(-9:NSEA) :: HCel !! Assign water depth to HCel from DW integer values. !! set half the minimum depth DMIN for negative cells. ! HCel(-9:0) = 0.5*DMIN HCel(1:NSEA)= DW(1:NSEA) !! Reset shallow water depth with minimum depth !$OMP Parallel DO Private(k) DO k=1, NSEA IF(DW(k) .LT. DMIN) HCel(k)=DMIN ENDDO !$OMP END Parallel DO !! Initialize full grid gradient arrays DDDX = 0. DDDY = 0. !! Use zero-gradient boundary condition or L0r1 > 0 L = 1 !! Calculate sea point water depth gradient CALL SMCGradn(HCel, GrHx, GrHy, L) !! Pass gradient values to DHDX, DHDY DHDX(1:NSEA) = GrHx DHDY(1:NSEA) = GrHy !! Apply limiter to depth-gradient and copy to full grid. !$OMP Parallel DO Private(i,j,k,m,n, CNST0, CNST1, CNST2) DO n=1,NSEA ! A limiter of gradient <= 0.1 is applied. IF( ABS( DHDX(n) ) .GT. 0.1) DHDX(n)=SIGN( 0.1, DHDX(n) ) IF( ABS( DHDY(n) ) .GT. 0.1) DHDY(n)=SIGN( 0.1, DHDY(n) ) !! Asign DHDX value to full grid variable DDDX i= IJKCel(1,n)/MRFct + 1 j= IJKCel(2,n)/MRFct + 1 k= MAX(1, IJKCel(3,n)/MRFct) m= MAX(1, IJKCel(4,n)/MRFct) DDDX(j:j+m-1,i:i+k-1) = DHDX(n) DDDY(j:j+m-1,i:i+k-1) = DHDY(n) !/ARC !Li Depth gradient in the Arctic part has to be rotated into !/ARC !Li the map-east system for calculation of refraction. !/ARC IF( n .GT. NGLO ) THEN !/ARC CNST0 = ANGARC(n - NGLO)*DERA !/ARC CNST1 = DHDX(n)*COS(CNST0) - DHDY(n)*SIN(CNST0) !/ARC CNST2 = DHDX(n)*SIN(CNST0) + DHDY(n)*COS(CNST0) !/ARC DHDX(n) = CNST1 !/ARC DHDY(n) = CNST2 !/ARC ENDIF END DO !$OMP END Parallel DO !! Calculate the depth gradient limiter for refraction. L = 0 !$OMP Parallel DO Private(i, n, CNST4, CNST6) DO n=1,NSEA !Li Work out magnitude of depth gradient CNST4 = 1.0001*SQRT(DHDX(n)*DHDX(n) + DHDY(n)*DHDY(n)) ! !Li Directional depedent depth gradient limiter. JGLi16Jun2011 IF ( CNST4 .GT. 1.0E-5 ) THEN !$OMP ATOMIC Update L = L + 1 !$OMP END ATOMIC DO i=1, NTH !Li Refraction is done only when depth gradient is non-zero. !Li Note ACOS returns value between [0, Pi), always positive. CNST6 = ACOS(-(DHDX(n)*ECOS(i)+DHDY(n)*ESIN(i))/CNST4 ) !Li User-defined refraction limiter added. JGLi09Jan2012 DHLMT(i,n)=MIN(Refran, 0.75*MIN(CNST6,ABS(PI-CNST6)))/DTH END DO !Li Output some values for inspection. JGLi22Jul2011 !/T IF( MOD(n, 1000) .EQ. 0 ) & !/T & WRITE(NDST,'(i8,18F5.1)' ) n, (DHLMT(i,n), i=1,18) ELSE DHLMT(:,n) = 0.0 ENDIF ENDDO !$OMP END Parallel DO !/T WRITE(NDST,*) ' No. Refraction points =', L !/T 999 PRINT*, ' Sub SMCDHXY ended.' RETURN END SUBROUTINE SMCDHXY ! Subroutine that calculates current velocity gradient for refraction. ! For consistency with the lat-lon grid, full grid DCXDXY, DCYDXY are ! assigned here. They are rotated to map-east system in the Arctic part. ! JGLi23Mar2016 ! SUBROUTINE SMCDCXY USE CONSTANTS USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSTA, MAPFS, MRFct, IJKCel !/ARC USE W3GDATMD, ONLY: NGLO, ANGARC USE W3ADATMD, ONLY: CX, CY, DCXDX, DCXDY, DCYDX, DCYDY USE W3ODATMD, ONLY: NDSE, NDST IMPLICIT NONE INTEGER :: I, J, K, L, M, N REAL :: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6 REAL, Dimension(NSEA) :: CXCY, GrHx, GrHy ! REAL, Dimension(-9:NSEA) :: CXCY !! Assign current CX speed to CXCY and set negative cells. ! CXCY(-9:0) = 0.0 !! Use zero-gradient boundary condition or L0r1 > 0 L = 1 CXCY(1:NSEA)= CX(1:NSEA) !! Initialize full grid gradient arrays !/DEBUGDCXDX WRITE(740+IAPROC,*) 'Before assigning DCXDX to ZERO' DCXDX = 0.0 DCXDY = 0.0 !! Calculate sea point water depth gradient CALL SMCGradn(CXCY, GrHx, GrHy, L) !! Apply limiter to CX-gradient and copy to full grid. !$OMP Parallel DO Private(i, j, k, m, n, CNST0, CNST1, CNST2) DO n=1,NSEA ! A limiter of gradient <= 0.01 is applied. IF( ABS( GrHx(n) ) .GT. 0.01) GrHx(n)=SIGN( 0.01, GrHx(n) ) IF( ABS( GrHy(n) ) .GT. 0.01) GrHy(n)=SIGN( 0.01, GrHy(n) ) !/ARC !Li Current gradient in the Arctic part has to be rotated into !/ARC !Li the map-east system for calculation of refraction. !/ARC IF( n .GT. NGLO ) THEN !/ARC CNST0 = ANGARC(n - NGLO)*DERA !/ARC CNST1 = GrHx(n)*COS(CNST0) - GrHy(n)*SIN(CNST0) !/ARC CNST2 = GrHx(n)*SIN(CNST0) + GrHy(n)*COS(CNST0) !/ARC GrHx(n) = CNST1 !/ARC GrHy(n) = CNST2 !/ARC ENDIF !! Asign CX gradients to full grid variable DCXDX/Y i= IJKCel(1,n)/MRFct + 1 j= IJKCel(2,n)/MRFct + 1 k= MAX(1, IJKCel(3,n)/MRFct) m= MAX(1, IJKCel(4,n)/MRFct) DCXDX(j:j+m-1,i:i+k-1) = GrHx(n) DCXDY(j:j+m-1,i:i+k-1) = GrHy(n) END DO !$OMP END Parallel DO !/DEBUGDCXDX WRITE(740+IAPROC,*) 'After non-trivial assination to DCXDX array' !! Assign current CY speed to CXCY and set negative cells. ! CXCY(-9:0) = 0.0 !! Use zero-gradient boundary condition or L0r1 > 0 L = 1 CXCY(1:NSEA)= CY(1:NSEA) !! Initialize full grid gradient arrays DCYDX = 0.0 DCYDY = 0.0 !! Calculate sea point water depth gradient CALL SMCGradn(CXCY, GrHx, GrHy, L) !! Apply limiter to CX-gradient and copy to full grid. !$OMP Parallel DO Private(i, j, k, m, n, CNST0, CNST1, CNST2) DO n=1,NSEA ! A limiter of gradient <= 0.1 is applied. IF( ABS( GrHx(n) ) .GT. 0.01) GrHx(n)=SIGN( 0.01, GrHx(n) ) IF( ABS( GrHy(n) ) .GT. 0.01) GrHy(n)=SIGN( 0.01, GrHy(n) ) !/ARC !Li Current gradient in the Arctic part has to be rotated into !/ARC !Li the map-east system for calculation of refraction. !/ARC IF( n .GT. NGLO ) THEN !/ARC CNST0 = ANGARC(n - NGLO)*DERA !/ARC CNST1 = GrHx(n)*COS(CNST0) - GrHy(n)*SIN(CNST0) !/ARC CNST2 = GrHx(n)*SIN(CNST0) + GrHy(n)*COS(CNST0) !/ARC GrHx(n) = CNST1 !/ARC GrHy(n) = CNST2 !/ARC ENDIF !! Asign CX gradients to full grid variable DCXDX/Y i= IJKCel(1,n)/MRFct + 1 j= IJKCel(2,n)/MRFct + 1 k= MAX(1, IJKCel(3,n)/MRFct) m= MAX(1, IJKCel(4,n)/MRFct) DCYDX(j:j+m-1,i:i+k-1) = GrHx(n) DCYDY(j:j+m-1,i:i+k-1) = GrHy(n) END DO !$OMP END Parallel DO !/T 999 PRINT*, ' Sub SMCDCXY ended.' RETURN END SUBROUTINE SMCDCXY !/ !/ ------------------------------------------------------------------- / !/ Two modified subs for SMC grid. JGLi 15Mar2011 !/ ------------------------------------------------------------------- / !/ SUBROUTINE W3GATHSMC ( ISPEC, FIELD ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH-III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 13-Jun-2006 | !/ +-----------------------------------+ !/ !/ 04-Jan-1999 : Distributed FORTRAN 77 version. ( version 1.18 ) !/ 13-Jan-2000 : Upgrade to FORTRAN 90 ( version 2.00 ) !/ Major changes to logistics. !/ 29-Dec-2004 : Multiple grid version. ( version 3.06 ) !/ 13-Jun-2006 : Split STORE in G/SSTORE ( version 3.09 ) !/ 9-Dec-2010 : Adapted for SMC grid propagtion. JGLi !/ ! 1. Purpose : ! ! Gather spectral bin information into a propagation field array. ! ! 2. Method : ! ! Direct copy or communication calls (MPP version). ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ISPEC Int. I Spectral bin considered. ! FIELD R.A. O Full field to be propagated. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ! MPI_STARTALL, MPI_WAITALL ! Subr. mpif.h MPI persistent comm. routines (!/MPI). ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3WAVE Subr. W3WAVEMD Actual wave model routine. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! - The field is extracted but not converted. ! - Array FIELD is not initialized. ! - MPI version requires posing of send and receive calls in ! W3WAVE to match local calls. ! - MPI version does not require an MPI_TESTALL call for the ! posted gather operation as MPI_WAITALL is mandatory to ! reset persistent communication for next time step. ! - MPI version allows only two new pre-fetch postings per ! call to minimize chances to be slowed down by gathers that ! are not yet needed, while maximizing the pre-loading ! during the early (low-frequency) calls to the routine ! where the amount of calculation needed for proagation is ! the largest. ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/SHRD Switch for message passing method. ! !/MPI Id. ! ! !/S Enable subroutine tracing. ! !/MPIT MPI test output. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/S USE W3SERVMD, ONLY: STRACE !/ USE W3GDATMD, ONLY: NSPEC, NX, NY, NSEA, NSEAL, NCel, MAPSF USE W3WDATMD, ONLY: A => VA !/MPI USE W3ADATMD, ONLY: MPIBUF, BSTAT, IBFLOC, ISPLOC, BISPL, & !/MPI NSPLOC, NRQSG2, IRQSG2, GSTORE !/MPI USE W3ODATMD, ONLY: NDST, IAPROC, NAPROC !/ IMPLICIT NONE ! !/MPI INCLUDE "mpif.h" !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: ISPEC REAL, INTENT(OUT) :: FIELD(NCel) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ !/SHRD INTEGER :: ISEA, IXY !/MPI INTEGER :: STATUS(MPI_STATUS_SIZE,NSPEC), & !/MPI IOFF, IERR_MPI, JSEA, ISEA, & !/MPI IXY, IS0, IB0, NPST, J !/S INTEGER, SAVE :: IENT !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3GATH') ! ! FIELD = 0. ! ! 1. Shared memory version ------------------------------------------ / ! !/SHRD DO ISEA=1, NSEA !/SHRD FIELD(ISEA) = A(ISPEC,ISEA) !/SHRD END DO ! !/SHRD RETURN ! ! 2. Distributed memory version ( MPI ) ----------------------------- / ! 2.a Update counters ! !/MPI ISPLOC = ISPLOC + 1 !/MPI IBFLOC = IBFLOC + 1 !/MPI IF ( IBFLOC .GT. MPIBUF ) IBFLOC = 1 ! ! 2.b Check status of present buffer ! 2.b.1 Scatter (send) still in progress, wait to end ! !/MPI IF ( BSTAT(IBFLOC) .EQ. 2 ) THEN !/MPI IOFF = 1 + (BISPL(IBFLOC)-1)*NRQSG2 !/MPI IF ( NRQSG2 .GT. 0 ) CALL & !/MPI MPI_WAITALL ( NRQSG2, IRQSG2(IOFF,2), & !/MPI STATUS, IERR_MPI ) !/MPI BSTAT(IBFLOC) = 0 !/MPI END IF ! ! 2.b.2 Gather (recv) not yet posted, post now ! !/MPI IF ( BSTAT(IBFLOC) .EQ. 0 ) THEN !/MPI BSTAT(IBFLOC) = 1 !/MPI BISPL(IBFLOC) = ISPLOC !/MPI IOFF = 1 + (ISPLOC-1)*NRQSG2 !/MPI IF ( NRQSG2 .GT. 0 ) CALL & !/MPI MPI_STARTALL ( NRQSG2, IRQSG2(IOFF,1), IERR_MPI ) !/MPI END IF ! ! 2.c Put local spectral densities in store ! !/MPI DO JSEA=1, NSEAL !/MPI GSTORE(IAPROC+(JSEA-1)*NAPROC,IBFLOC) = A(ISPEC,JSEA) !/MPI END DO ! ! 2.d Wait for remote spectral densities ! !/MPI IOFF = 1 + (BISPL(IBFLOC)-1)*NRQSG2 !/MPI IF ( NRQSG2 .GT. 0 ) CALL & !/MPI MPI_WAITALL ( NRQSG2, IRQSG2(IOFF,1), STATUS, IERR_MPI ) ! ! 2.e Convert storage array to field. ! !/MPI DO ISEA=1, NSEA !/MPI FIELD(ISEA) = GSTORE(ISEA,IBFLOC) !/MPI END DO ! ! 2.f Pre-fetch data in available buffers ! !/MPI IS0 = ISPLOC !/MPI IB0 = IBFLOC !/MPI NPST = 0 ! !/MPI DO J=1, MPIBUF-1 !/MPI IS0 = IS0 + 1 !/MPI IF ( IS0 .GT. NSPLOC ) EXIT !/MPI IB0 = 1 + MOD(IB0,MPIBUF) !/MPI IF ( BSTAT(IB0) .EQ. 0 ) THEN !/MPI BSTAT(IB0) = 1 !/MPI BISPL(IB0) = IS0 !/MPI IOFF = 1 + (IS0-1)*NRQSG2 !/MPI IF ( NRQSG2 .GT. 0 ) CALL & !/MPI MPI_STARTALL ( NRQSG2, IRQSG2(IOFF,1), IERR_MPI ) !/MPI NPST = NPST + 1 !/MPI END IF !/MPI IF ( NPST .GE. 2 ) EXIT !/MPI END DO ! !/MPI RETURN ! !/ End of W3GATHSMC ----------------------------------------------------- / !/ END SUBROUTINE W3GATHSMC ! !/ ------------------------------------------------------------------- / SUBROUTINE W3SCATSMC ( ISPEC, MAPSTA, FIELD ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH-III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 13-Jun-2006 | !/ +-----------------------------------+ !/ !/ 04-Jan-1999 : Distributed FORTRAN 77 version. ( version 1.18 ) !/ 13-Jan-2000 : Upgrade to FORTRAN 90 ( version 2.00 ) !/ Major changes to logistics. !/ 28-Dec-2004 : Multiple grid version. ( version 3.06 ) !/ 07-Sep-2005 : Updated boundary conditions. ( version 3.08 ) !/ 13-Jun-2006 : Split STORE in G/SSTORE ( version 3.09 ) !/ 9-Dec-2010 : Adapted for SMC grid propagtion. JGLi09Dec2010 !/ 16-Jan-2012 : Remove MAPSTA checking for SMC grid. JGLi16Jan2012 !/ !/ ! 1. Purpose : ! ! 'Scatter' data back to spectral storage after propagation. ! ! 2. Method : ! ! Direct copy or communication calls (MPP version). ! See also W3GATH. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ISPEC Int. I Spectral bin considered. ! MAPSTA I.A. I Status map for spatial grid. ! FIELD R.A. I SMC grid field to be propagated. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ! MPI_STARTALL, MPI_WAITALL, MPI_TESTALL ! Subr. mpif.h MPI persistent comm. routines (!/MPI). ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3WAVE Subr. W3WAVEMD Actual wave model routine. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! - The field is put back but not converted ! ! - MPI persistent communication calls initialize in W3MPII. ! - See W3GATH and W3MPII for additional comments on data ! buffering. ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/SHRD Switch for message passing method. ! !/MPI Id. ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / !/S USE W3SERVMD, ONLY: STRACE !/ USE W3GDATMD, ONLY: NSPEC, NX, NY, NSEA, NCel, NSEAL, MAPSF USE W3WDATMD, ONLY: A => VA !/MPI USE W3ADATMD, ONLY: MPIBUF, BSTAT, IBFLOC, ISPLOC, BISPL, & !/MPI NSPLOC, NRQSG2, IRQSG2, SSTORE USE W3ODATMD, ONLY: NDST !/MPI USE W3ODATMD, ONLY: IAPROC, NAPROC !/ IMPLICIT NONE ! !/MPI INCLUDE "mpif.h" !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: ISPEC, MAPSTA(NY*NX) REAL, INTENT(IN) :: FIELD(NCel) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ !/SHRD INTEGER :: ISEA, IXY !/MPI INTEGER :: ISEA, IXY, IOFF, IERR_MPI, J, & !/MPI STATUS(MPI_STATUS_SIZE,NSPEC), & !/MPI JSEA, IB0 !/S INTEGER, SAVE :: IENT !/MPI LOGICAL :: DONE !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3SCAT') ! ! 1. Shared memory version ------------------------------------------ * ! !/SHRD DO ISEA=1, NSEA !/SHRD IXY = MAPSF(ISEA,3) !/SHRD IF ( MAPSTA(IXY) .GE. 1 ) A(ISPEC,ISEA) = FIELD(ISEA) !/SHRD END DO ! !/SHRD RETURN ! ! 2. Distributed memory version ( MPI ) ----------------------------- * ! 2.a Initializations ! ! 2.b Convert full grid to sea grid, active points only ! !/MPI DO ISEA=1, NSEA !/MPI IXY = MAPSF(ISEA,3) !/MPI IF ( MAPSTA(IXY) .GE. 1 ) SSTORE(ISEA,IBFLOC) = FIELD(ISEA) !/MPI END DO ! ! 2.c Send spectral densities to appropriate remote ! !/MPI IOFF = 1 + (ISPLOC-1)*NRQSG2 !/MPI IF ( NRQSG2 .GT. 0 ) CALL & !/MPI MPI_STARTALL ( NRQSG2, IRQSG2(IOFF,2), IERR_MPI ) !/MPI BSTAT(IBFLOC) = 2 ! ! 2.d Save locally stored results ! !/MPI DO JSEA=1, NSEAL !/MPI !!Li ISEA = IAPROC+(JSEA-1)*NAPROC !/MPI ISEA = MIN( IAPROC+(JSEA-1)*NAPROC, NSEA ) !/MPI A(ISPEC,JSEA) = SSTORE(ISEA,IBFLOC) !/MPI END DO ! ! 2.e Check if any sends have finished ! !/MPI IB0 = IBFLOC ! !/MPI DO J=1, MPIBUF !/MPI IB0 = 1 + MOD(IB0,MPIBUF) !/MPI IF ( BSTAT(IB0) .EQ. 2 ) THEN !/MPI IOFF = 1 + (BISPL(IB0)-1)*NRQSG2 !/MPI IF ( NRQSG2 .GT. 0 ) THEN !/MPI CALL MPI_TESTALL ( NRQSG2, IRQSG2(IOFF,2), DONE, & !/MPI STATUS, IERR_MPI ) !/MPI ELSE !/MPI DONE = .TRUE. !/MPI END IF !/MPI IF ( DONE .AND. NRQSG2.GT.0 ) CALL & !/MPI MPI_WAITALL ( NRQSG2, IRQSG2(IOFF,2), & !/MPI STATUS, IERR_MPI ) !/MPI IF ( DONE ) THEN !/MPI BSTAT(IB0) = 0 !/MPI END IF !/MPI END IF !/MPI END DO ! ! 2.f Last component, finish message passing, reset buffer control ! !/MPI IF ( ISPLOC .EQ. NSPLOC ) THEN ! !/MPI DO IB0=1, MPIBUF !/MPI IF ( BSTAT(IB0) .EQ. 2 ) THEN !/MPI IOFF = 1 + (BISPL(IB0)-1)*NRQSG2 !/MPI IF ( NRQSG2 .GT. 0 ) CALL & !/MPI MPI_WAITALL ( NRQSG2, IRQSG2(IOFF,2), & !/MPI STATUS, IERR_MPI ) !/MPI BSTAT(IB0) = 0 !/MPI END IF !/MPI END DO ! !/MPI ISPLOC = 0 !/MPI IBFLOC = 0 ! !/MPI END IF ! !/MPI RETURN ! ! Formats ! !/ !/ End of W3SCATSMC ----------------------------------------------------- / !/ END SUBROUTINE W3SCATSMC !/ !/ End of two new subs for SMC grid. JGLi 15Mar2011 !/ !/ End of module W3PSMCMD -------------------------------------------- / !/ END MODULE W3PSMCMD