! $Id: tpcore_mod.f,v 1.10 2007/11/05 16:16:25 bmy Exp $
MODULE TPCORE_MOD 1
!
!******************************************************************************
! Module TPCORE_MOD contains the TPCORE transport subroutine package by
! S-J Lin, version 7.1. (bmy, 7/16/01, 9/18/07)
!
! Module Routines:
! ============================================================================
! (1 ) TPCORE : TPCORE driver program
! (2 ) COSA : TPCORE intermediate subroutine
! (3 ) COSC : TPCORE intermediate subroutine
! (4 ) FCT3D : TPCORE intermediate subroutine
! (5 ) FILEW : TPCORE intermediate subroutine
! (6 ) FILNS : TPCORE intermediate subroutine
! (7 ) FXPPM : TPCORE intermediate subroutine
! (8 ) FYPPM : TPCORE intermediate subroutine
! (9 ) FZPPM : TPCORE intermediate subroutine
! (10) HILO : TPCORE intermediate subroutine
! (11) HILO3D : TPCORE intermediate subroutine
! (12) LMTPPM : TPCORE intermediate subroutine
! (13) QCKXYZ : TPCORE intermediate subroutine
! (14) XADV : TPCORE intermediate subroutine
! (15) XMIST : TPCORE intermediate subroutine
! (16) XTP : TPCORE intermediate subroutine
! (17) YMIST : TPCORE intermediate subroutine
! (18) YTP : TPCORE intermediate subroutine
! (19) PRESS_FIX : Wrapper for pressure-fixer subroutine DYN0
! (20) DYN0 : Implements pressure fix for mass fluxes in TPCORE
! (21) PFILTR : Applies pressure filter to ALFA and BETA mass fluxes
! (22) LOCFLT : Local pressure filter -- called from PFILTR
! (23) POLFLT : Polar pressure filter -- called from PFILTR
! (24) DIAG_FLUX : Computes TPCORE mass fluxes for ND24, ND25, ND26 diags
!
! GEOS-CHEM modules referenced by tagged_co_mod.f
! ============================================================================
! (1 ) diag_mod.f : Module containing GEOS-CHEM diagnostic arrays
! (2 ) dao_mod.f : Module containing DAO met field arrays
! (3 ) global_ch4_mod.f : Module containing routines to read 3-D CH4 field
! (4 ) grid_mod.f : Module containing horizontal grid information
! (5 ) pressure_mod.f : Module containing routines to compute P(I,J,L)
! (6 ) time_mod.f : Module containing routines to compute date & time
!
! NOTES:
! (1 ) The TPCORE subroutines have not been modified, except to replace
! obsolete parallel loop directives. It is more convenient to place
! all of the TPCORE subroutines into a single module, this reduces
! clutter. (bmy, 7/16/01)
! (2 ) All parallel loops are now specified with OpenMP directives,
! for cross-platform compatibility. (bmy, 7/16/01)
! (3 ) The routines in TPCORE_MOD have been validated against the previous
! version (Code_4.16). (bmy, 7/16/01)
! (4 ) Updated comments (bmy, 9/4/01)
! (5 ) Removed obsolete code from 7/12/01. Also implemented pressure-fix
! subroutines PRESS_FIX, DYN0, PFLITR, LOCFLT, POLFLT. (bmy, 10/9/01)
! (6 ) Now use PSC2 instead of PS in subroutine DYN0. Also delineate the
! first-time header text with horizontal lines. (bdf, bmy, 4/15/02)
! (7 ) Now zero XMASS_PF and YMASS_PF arrays on every call to TPCORE.
! This will avoid floating-point exceptions on the Alpha platform.
! (bmy, 4/18/02)
! (8 ) Now divide module header into MODULE PRIVATE, MODULE VARIABLES, and
! MODULE ROUTINES sections. Updated comments (bmy, 5/28/02)
! (9 ) Deleted obsolete code from 4/02. (bdf, bmy, 8/22/02)
! (10) Minor bug fix for ALPHA platform: delete extra comma in format
! statement 2 in routine TPCORE. Bug fix: now stop the run if NDT is
! too large. This makes sure we don't violate the Courant limit.
! (bmy, 11/22/02)
! (11) Also add output for the SUN/Sparc platform. Rename DEC_COMPAQ to
! COMPAQ. Also assume that all platforms other than CRAY use OPENMP
! parallelization commands (bmy, 3/23/03)
! (12) Now references "grid_mod.f" and "time_mod.f" (bmy, 3/24/03)
! (13) Now print output for IBM/AIX platform in "tpcore" (gcc, bmy, 6/27/03)
! (14) Remove obsolete code for CO-OH parameterization (bmy, 6/24/05)
! (15) Bug fix in DIAG_FLUX: now dimension FX, FX properly (bmy, 7/21/05)
! (16) Now print output for IFORT compiler in "tpcore" (bmy, 10/18/05)
! (17) Remove support for LINUX_IFC & LINUX_EFC compilers (bmy, 8/4/06)
! (18) Corrected mass flux diagnostics (phs, 9/18/07)
!******************************************************************************
!
!=================================================================
! MODULE PRIVATE DECLARATIONS -- keep certain internal variables
! and routines from being seen outside "tpcore_mod.f"
!=================================================================
! Make everything PRIVATE ...
PRIVATE
! ... except this routine
PUBLIC :: TPCORE
!=================================================================
! MODULE ROUTINES -- follow below the "CONTAINS" statement
!=================================================================
CONTAINS
!------------------------------------------------------------------------------
C ****6***0*********0*********0*********0*********0*********0**********72
subroutine tpcore(IGD,Q,PS1,PS2,U,V,W,NDT,IORD,JORD,KORD,NC,IM, 1,10
& JM,j1,NL,AP,BP,PT,AE,FILL,MFCT,Umax)
C****6***0*********0*********0*********0*********0*********0**********72
C TransPort module for Goddard Chemistry Transport Model (G-CTM), Goddard
C Earth Observing System General Circulation Model (GEOS-GCM), and Data
C Assimilation System (GEOS-DAS).
C Purpose: perform the transport of 3-D mixing ratio fields using
C externally specified winds on the hybrid Eta-coordinate.
C One call to tpcore updates the 3-D mixing ratio
C fields for one time step (NDT). [vertical mass flux is computed
C internally using a center differenced hydrostatic mass
C continuity equation].
C Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) schemes
C (Lin and Rood 1996, MWR) with a modified MFCT option (Zalesak 1979).
C Multitasking version: 7.1
C Last modified: Sept 2, 1999
C Changes from version 7.m: large-time-step bug in xtp fixed.
C Suggested compiler options:
C CRAY f77 compiler: cf77 -Zp -c -Wd'-dec' -Wf' -a stack -exm'
C CRAY f90 compiler: f90 -c -eZ -DCRAY -Dmultitask
C SGI Origin: f77 -c -DSGI -Dmultitask -r8 -64 -O3 -mips4 -mp
C loader: f77 -64 -mp
C
C Send comments/suggestions to
C
C S.-J. Lin
C Address:
C Code 910.3, NASA/GSFC, Greenbelt, MD 20771
C Phone: 301-614-6161
C E-mail: slin@dao.gsfc.nasa.gov
C
C The algorithm is based on the following papers:
C 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi-
C Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070.
C
C 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of
C the van Leer-type transport schemes and its applications to the moist-
C ure transport in a General Circulation Model. Mon. Wea. Rev., 122,
C 1575-1593.
C
C 3. Lin, S.-J., and R. B. Rood, 1997: Multidimensional flux form semi-
C Lagrangian transport schemes- MFCT option. To be submitted.
C ======
C INPUT:
C ======
C IGD: (horizontal) grid type on which winds are defined.
C IGD = 0 A-Grid [all variables defined at the same point from south
C pole (j=1) to north pole (j=JM) ]
C IGD = 1 GEOS-GCM C-Grid (Max Suarez's center difference dynamical core)
C [North]
C V(i,j)
C |
C |
C |
C [WEST] U(i-1,j)---Q(i,j)---U(i,j) [EAST]
C |
C |
C |
C V(i,j-1)
C [South]
C U(i, 1) is defined at South Pole.
C V(i, 1) is half grid north of the South Pole.
C V(i,JM-1) is half grid south of the North Pole.
C
C V must be defined at j=1 and j=JM-1 if IGD=1
C V at JM need not be defined.
C Q(IM,JM,NL,NC): mixing ratios at current time (t)
C NC: total # of constituents
C IM: first (E-W) dimension; # of Grid intervals in E-W is IM
C JM: 2nd (N-S) dimension; # of Grid intervals in N-S is JM-1
C NL: 3rd dimension (# of layers); vertical index increases from 1 at
C the model top to NL near the surface (see fig. below).
C It is assumed that NL > 5.
C
C PS1(IM,JM): surface pressure at current time (t)
C PS2(IM,JM): surface pressure at mid-time-level (t+NDT/2)
C PS2 is replaced by the predicted PS (at t+NDT) on output.
C Note: surface pressure can have any unit or can be multiplied by any
C const.
C
C The hybrid ETA-coordinate:
C
C pressure at layer edges are defined as follows:
C
C p(i,j,k) = AP(k)*PT + BP(k)*PS(i,j) (1)
C
C Where PT is a constant having the same unit as PS.
C AP and BP are unitless constants given at layer edges.
C In all cases BP(1) = 0., BP(NL+1) = 1.
C The pressure at the model top is PTOP = AP(1)*PT
C
C *********************
C For pure sigma system
C *********************
C AP(k) = 1 for all k, PT = PTOP,
C BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP, where Psfc
C is the true surface pressure.
C
C /////////////////////////////////
C / \ ------ Model top P=PTOP --------- AP(1), BP(1)
C |
C delp(1) | ........... Q(i,j,1) ............
C |
C W(k=1) \ / --------------------------------- AP(2), BP(2)
C
C
C
C W(k-1) / \ --------------------------------- AP(k), BP(k)
C |
C delp(K) | ........... Q(i,j,k) ............
C |
C W(k) \ / --------------------------------- AP(k+1), BP(k+1)
C
C
C
C / \ --------------------------------- AP(NL), BP(NL)
C |
C delp(NL) | ........... Q(i,j,NL) .........
C |
C W(NL)=0 \ / -----Earth's surface P=Psfc ------ AP(NL+1), BP(NL+1)
C //////////////////////////////////
C U(IM,JM,NL) & V(IM,JM,NL):winds (m/s) at mid-time-level (t+NDT/2)
C Note that on return U and V are destroyed.
C NDT (integer): time step in seconds (need not be constant during the course of
C the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5
C (Lat-Lon) resolution. Smaller values maybe needed if the model
C has a well-resolved stratosphere and Max(V) > 225 m/s
C
C J1 determines the size of the polar cap:
C South polar cap edge is located at -90 + (j1-1.5)*180/(JM-1) deg.
C North polar cap edge is located at 90 - (j1-1.5)*180/(JM-1) deg.
C There are currently only two choices (j1=2 or 3).
C IM must be an even integer if j1 = 2. Recommended value: J1=3.
C
C IORD, JORD, and KORD are integers controlling various options in E-W, N-S,
C and vertical transport, respectively.
C
C
C _ORD=
C 1: 1st order upstream scheme (too diffusive, not a real option; it
C can be used for debugging purposes; this is THE only known "linear"
C monotonic advection scheme.).
C 2: 2nd order van Leer (full monotonicity constraint;
C see Lin et al 1994, MWR)
C 3: monotonic PPM* (Collela & Woodward 1984)
C 4: semi-monotonic PPM (same as 3, but overshoots are allowed)
C 5: positive-definite PPM (constraint on the subgrid distribution is
C only strong enough to prevent generation of negative values;
C both overshoots & undershootes are possible).
C 6: un-constrained PPM (nearly diffusion free; faster but
C positivity of the subgrid distribution is not quaranteed. Use
C this option only when the fields and winds are very smooth or
C when MFCT=.true.)
C 7: Huynh/Van Leer/Lin full monotonicity constraint
C Only KORD can be set to 7 to enable the use of Huynh's 2nd monotonicity
C constraint for piece-wise parabolic distribution.
C
C *PPM: Piece-wise Parabolic Method
C
C Recommended values:
C IORD=JORD=3 for high horizontal resolution.
C KORD=6 or 7 if MFCT=.true.
C KORD=3 or 7 if MFCT=.false.
C
C The implicit numerical diffusion decreases as _ORD increases.
C DO not use option 4 or 5 for non-positive definite scalars
C (such as Ertel Potential Vorticity).
C
C If numerical diffusion is a problem (particularly at low horizontal
C resolution) then the following setup is recommended:
C IORD=JORD=KORD=6 and MFCT=.true.
C
C AE: Radius of the sphere (meters).
C Recommended value for the planet earth: 6.371E6
C
C FILL (logical): flag to do filling for negatives (see note below).
C MFCT (logical): flag to do a Zalesak-type Multidimensional Flux
C correction. It shouldn't be necessary to call the
C filling routine when MFCT is true.
C
C Umax: Estimate (upper limit) of the maximum U-wind speed (m/s).
C (225 m/s is a good value for troposphere model; 300 m/s otherwise)
C
C ======
C Output
C ======
C
C Q: the updated mixing ratios at t+NDT (original values are over-written)
C W(;;NL): large-scale vertical mass flux as diagnosed from the hydrostatic
C relationship. W will have the same unit as PS1 and PS2 (eg, mb).
C W must be divided by NDT to get the correct mass-flux unit.
C The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND
C is the pressure thickness in the "upwind" direction. For example,
C C(k) = W(k)/delp(k) if W(k) > 0;
C C(k) = W(k)/delp(k+1) if W(k) < 0.
C ( W > 0 is downward, ie, toward surface)
C PS2: predicted PS at t+NDT (original values are over-written)
C
C Memory usage:
C This code is optimized for speed. it requres 18 dynamically allocated
C 3D work arrays (IM,JM,NL) regardless of the value of NC.
C Older versions (version 4 or 4.5) use less memory if NC is small.
C =====
C NOTES:
C =====
C
C This forward-in-time upstream-biased transport scheme degenerates to
C the 2nd order center-in-time center-in-space mass continuity eqn.
C if Q = 1 (constant fields will remain constant). This degeneracy ensures
C that the computed vertical velocity to be identical to GEOS-1 GCM
C for on-line transport.
C
C A larger polar cap is used if j1=3 (recommended for C-Grid winds or when
C winds are noisy near poles).
C
C The user needs to change the parameter Jmax or Kmax if the resolution
C is greater than 0.25 deg in N-S or 500 layers in the vertical direction.
C (this TransPort Core is otherwise resolution independent and can be used
C as a library routine).
C PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd
C order accurate for non-uniform grid (vertical sigma coord.).
C Time step is limitted only by transport in the meridional direction.
C (the FFSL scheme is not implemented in the meridional direction).
C Since only 1-D limiters are applied, negative values could
C potentially be generated when large time step is used and when the
C initial fields contain discontinuities.
C This does not necessarily imply the integration is unstable.
C These negatives are typically very small. A filling algorithm is
C activated if the user set "fill" to be true.
C Alternatively, one can use the MFCT option to enforce monotonicity.
! Added to pass C-preprocessor switches (bmy, 3/9/01)
# include "define.h"
C ****6***0*********0*********0*********0*********0*********0**********72
PARAMETER (Jmax = 721, kmax = 200)
C ****6***0*********0*********0*********0*********0*********0**********72
C Input-Output arrays
REAL Q(IM,JM,NL,NC),PS1(IM,JM),PS2(IM,JM),W(IM,JM,NL),
& U(IM,JM,NL),V(IM,JM,NL),AP(NL+1),BP(NL+1)
LOGICAL ZCROSS, FILL, MFCT, deform
C Local dynamic arrays
REAL CRX(IM,JM,NL),CRY(IM,JM,NL),delp(IM,JM,NL),delp1(IM,JM,NL),
& xmass(IM,JM,NL),ymass(IM,JM,NL),delp2(IM,JM,NL),
& DG1(IM),DG2(IM,JM),DPI(IM,JM,NL),qlow(IM,JM,NL),
& WK(IM,JM,NL),PU(IM,JM,NL),DQ(IM,JM,NL),
& fx(IM+1,JM,NL),fy(IM,JM,NL),fz(IM,JM,NL+1),
& qz(IM,JM,NL),Qmax(IM,JM,NL),Qmin(IM,JM,NL)
! bey, 6/20/00. for mass-flux diagnostic
REAL fx1_tp(IM,JM,NL), fy1_tp(IM,JM,NL), fz1_tp(IM,JM,NL)
INTEGER JS(NL),JN(NL)
C Local static arrays
REAL DTDX(Jmax), DTDX5(Jmax), acosp(Jmax),cosp(Jmax),
& cose(Jmax), DAP(kmax), DBK(kmax)
DATA NDT0, NSTEP /0, 0/
DATA ZCROSS /.true./
C Saved internal variables:
SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML, DTDX,
& DTDX5, acosp, COSP, COSE, DAP,DBK
! New variables for TPCORE pressure fixer (bdf, bmy, 10/11/01)
REAL YMASS_PF(IM,JM,NL), XMASS_PF(IM,JM,NL), TEMP(IM,JM,NL)
LOGICAL PRESSURE_FIX
PRESSURE_FIX = .TRUE.
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
deform = .false.
JM1 = JM -1
IMH = IM/2
j2 = JM - j1 + 1
NSTEP = NSTEP + 1
C****6***0*********0*********0*********0*********0*********0**********72
C Initialization
C****6***0*********0*********0*********0*********0*********0**********72
! Moved further down to be done for each tracer (phs, 30/8/07)
!! ! For mass flux diagnostics (bey, 6/20/00)
!! fx1_tp(:,:,:) = 0d0
!! fy1_tp(:,:,:) = 0d0
!! fz1_tp(:,:,:) = 0d0
!!
!! ! Also need to initialize these arrays, so that the flux diagnostics
!! ! will be identical for single or multi processor (bmy, 9/29/00)
!! fx(:,:,:) = 0d0
!! fy(:,:,:) = 0d0
!! fz(:,:,:) = 0d0
!!
! Need to initialize these arrays in order to avoid
! floating-point exceptions on Alpha (lyj, bmy, 4/19/02)
YMASS_PF(:,:,:) = 0d0
XMASS_PF(:,:,:) = 0d0
if(NSTEP.eq.1) then
! Updated output (bmy, 3/13/03)
WRITE( 6, '(a)' ) REPEAT( '=', 79 )
WRITE( 6, '(a)' ) 'T P C O R E -- FFSL TransPort Core v. 7.1'
WRITE( 6, '(a)' )
WRITE( 6, '(a)' ) 'Originally written by S-J Lin'
WRITE( 6, '(a)' )
WRITE( 6, '(a)' )
& 'Modified for GEOS-CHEM by Isabelle Bey, Brendan Field, and'
WRITE( 6, '(a)' )
& 'Bob Yantosca, with the addition of flux diagnostics and the'
WRITE( 6, '(a)' ) 'DYN0 pressure fixer from M. Prather'
WRITE( 6, '(a)' )
WRITE( 6, '(a)' ) 'Last Modification Date: 8/22/02'
WRITE( 6, '(a)' )
#if ( multitask )
WRITE( 6, '(a)' ) 'TPCORE was compiled for multitasking'
#if defined( CRAY )
WRITE( 6, '(a)' ) 'for CRAY'
#elif defined( SGI_MIPS )
WRITE( 6, '(a)' ) 'for SGI Origin/Power Challenge machines'
#elif defined( COMPAQ )
WRITE( 6, '(a)' ) 'for COMPAQ/HP RISC Alpha machines'
#elif defined( LINUX_PGI )
WRITE( 6, '(a)' ) 'for Linux environment w/ PGI compiler'
#elif defined( LINUX_IFORT )
WRITE( 6, '(a)' ) 'for Linux environment w/ Intel IFORT compiler'
#elif defined( SPARC )
WRITE( 6, '(a)' ) 'for SUN/Sparc machines'
#elif defined( IBM_AIX )
WRITE( 6, '(a)' ) 'for IBM/AIX machines'
#endif
#endif
! Added output on the first time TPCORE is called (bmy, 10/11/01)
IF ( PRESSURE_FIX ) THEN
WRITE( 6, '(a)' )
WRITE( 6, '(a)' ) 'TPCORE PRESSURE FIXER is turned ON!'
ENDIF
if( MFCT ) then
WRITE( 6, '(a)' )
WRITE( 6, '(a)' ) 'MFCT option is on!'
endif
! Updated output (bmy, 4/15/02)
WRITE( 6, '(a)' )
WRITE( 6, 2 ) IM, JM, NL, j1
2 FORMAT( 'IM= ', i3,1x,'JM= ', i3,1x,'NL= ',i3,1x,'J1= ',i3 )
! Updated output (bmy, 4/15/02)
WRITE( 6, 3 ) NC, IORD, JORD, KORD, NDT
3 FORMAT( 'NC= ',i3,1x,'IORD=',i3,1x,'JORD=',i3,1x,
& 'KORD=',i3,1x,'NDT= ',i8)
if(NL.LT.6) then
write(6,*) 'stop in module tpcore'
write(6,*) 'NL must be >=6'
stop
endif
if(Jmax.lt.JM .or. Kmax.lt.NL) then
write(6,*) 'stop in module tpcore'
write(6,*) 'Jmax or Kmax is too small; see documentation'
stop
endif
DO 5 k=1,NL
DAP(k) = (AP(k+1) - AP(k))*PT
5 DBK(k) = BP(k+1) - BP(k)
PI = 4. * ATAN(1.)
DL = 2.*PI / float(IM)
DP = PI / float(JM1)
if(IGD.eq.0) then
C Compute analytic cosine at cell edges
call cosa
(cosp,cose,JM,PI,DP)
else
C Define cosine consistent with GEOS-GCM (using dycore2.0 or later)
call cosc(cosp,cose,JM,PI,DP)
endif
do 15 J=2,JM1
15 acosp(j) = 1./cosp(j)
C Inverse of the Scaled polar cap area.
agle = (float(j1)-1.5)*DP
RCAP = DP / ( float(IM)*(1.-COS(agle)) )
acosp(1) = RCAP
acosp(JM) = RCAP
ENDIF
if(NDT0 .ne. NDT) then
DT = NDT
NDT0 = NDT
CR1 = abs(Umax*DT)/(DL*AE)
MaxDT = DP*AE / abs(Umax) + 0.5
! Updated output (bmy, 4/15/02)
WRITE( 6, '(a)' )
WRITE(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT
! Bug fix: Now stop the run if NDT is too large. This will make
! sure that we don't violate the Courant limit. (bmy, 11/22/02)
if(MaxDT .lt. abs(NDT)) then
write(6,*) 'Warning!!! NDT maybe too large!'
STOP
endif
if(CR1.ge.0.95) then
JS0 = 0
JN0 = 0
IML = IM-2
ZTC = 0.
else
ZTC = acos(CR1) * (180./PI)
JS0 = float(JM1)*(90.-ZTC)/180. + 2
JS0 = max(JS0, J1+1)
IML = min(6*JS0/(J1-1)+2, 4*IM/5)
JN0 = JM-JS0+1
endif
! Updated output (bmy, 4/15/02)
WRITE( 6, '(''ZTC= '', f13.6)') ZTC
WRITE( 6, 21 ) JS0, JN0, IML
21 FORMAT( 'JS= ',i3,1x, 'JN= ',i3,1x,'IML= ',i3 )
do 22 J=2,JM1
DTDX(j) = DT / ( DL*AE*COSP(J) )
DTDX5(j) = 0.5*DTDX(j)
22 continue
DTDY = DT /(AE*DP)
DTDY5 = 0.5*DTDY
! Updated output (bmy, 4/15/02)
WRITE( 6, 23 ) J1, J2
23 FORMAT( 'J1= ',i3,1x, 'J2= ',i3 )
! Fancy output to stdout (bmy, 3/13/03)
WRITE( 6, '(a)' ) REPEAT( '=', 79 )
ENDIF ! END INITIALIZATION.
C****6***0*********0*********0*********0*********0*********0**********72
C Compute Courant number
C****6***0*********0*********0*********0*********0*********0**********72
if(IGD.eq.0) then
C Convert winds on A-Grid to Courant # on C-Grid.
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all shared(NL,im,jm1,jm,U,V,dtdx5,dtdy5,CRX,CRY)
CMIC$* private(i,j,k)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K )
#endif
#endif
do k=1,NL
do 46 j=2,JM1
do 46 i=2,IM
46 CRX(i,j,k) = dtdx5(j)*(U(i,j,k)+U(i-1,j,k))
C for i=1
do 48 j=2,JM1
48 CRX(1,j,k) = dtdx5(j)*(U(1,j,k)+U(IM,j,k))
do 49 j=2,JM
do 49 i=1,IM
49 CRY(i,j,k) = DTDY5*(V(i,j,k)+V(i,j-1,k))
enddo
else
C Convert winds on C-grid to Courant #
C Beware of the index shifting!! (GEOS-GCM)
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all shared(NL,im,jm1,jm,U,V,dtdx,dtdy,CRX,CRY)
CMIC$* private(i,j,k)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K )
#endif
#endif
DO 65 k=1,NL
do 50 j=2,JM1
do 50 i=2,IM
50 CRX(i,j,k) = dtdx(j)*U(i-1,j,k)
do 55 j=2,JM1
55 CRX(1,j,k) = dtdx(j)*U(IM,j,k)
do 60 j=2,JM
do 60 i=1,IM
60 CRY(i,j,k) = DTDY*V(i,j-1,k)
65 continue
endif
!=================================================================
! ***** T P C O R E P R E S S U R E F I X E R *****
!
! Run pressure fixer to fix mass conservation problem. Pressure
! fixer routines PRESS_FIX, DYN0, PFILTR, LOCFLT, and POLFLT
! change the mass fluxes so they become consistant with met field
! pressures. (bdf, bmy, 10/11/01)
!
! NOTE: The pressure fixer is not 100% perfect; tracer mass will
! increase on the order of 0.5%/yr. However, this is much
! better than w/o the pressure fixer, where the mass may
! increase by as much as 40%/yr. (bdf, bmy, 10/22/01)
!=================================================================
IF ( PRESSURE_FIX ) THEN
! Loop over vertical levels --
! added parallel loop #if statements (bmy, 10/11/01)
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all shared(NL,IM,JM,JM1,COSE,XMASS_PF,YMASS_PF,DELP2,CRX,CRY)
CMIC$* private(I,J,K,D5)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K, D5 )
#endif
#endif
DO K = 1, NL
! DELP = pressure thickness:
! the pseudo-density in a hydrostatic system.
DO J = 1, JM
DO I = 1, IM
DELP2(I,J,K) = DAP(K) + DBK(K)*PS2(I,J)
ENDDO
ENDDO
! calculate mass fluxes for pressure fixer.
! N-S component
DO J = J1, J2+1
D5 = 0.5 * COSE(J)
DO I = 1, IM
YMASS_PF(I,J,K) =
& CRY(I,J,K) * D5 * (DELP2(I,J,K)+DELP2(I,J-1,K))
ENDDO
ENDDO
! Enlarged polar cap.
IF(J1.NE.2) THEN
DO I=1,IM
YMASS_PF(I,1,K) = 0
YMASS_PF(I,JM1+1,K) = 0
ENDDO
ENDIF
! E-W component
DO J = J1, J2
DO I = 2, IM
PU(I,J,K) = 0.5 * (DELP2(I,J,K) + DELP2(I-1,J,K))
ENDDO
ENDDO
DO J = J1, J2
PU(1,J,K) = 0.5 * (DELP2(1,J,K) + DELP2(IM,J,K))
ENDDO
DO J = J1, J2
DO I = 1, IM
XMASS_PF(I,J,K) = PU(I,J,K) * CRX(I,J,K)
ENDDO
ENDDO
ENDDO
!==============================================================
! Call PRESS_FIX to apply the pressure fix to the mass fluxes
! XMASS_PF, YMASS_PF. PRESS_FIX will call routine DYN0, etc.
!==============================================================
CALL PRESS_FIX( XMASS_PF, YMASS_PF, NDT, ACOSP, J1 )
! Loop over vertical levels --
! added parallel loop #if statements (bmy, 10/11/01)
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all shared(NL,IM,JM,XMASS_PF,PU,YMASS_PF,DELP2,COSE,CRX,CRY)
CMIC$* private(I,J,K,D5)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K, D5 )
#endif
#endif
DO K = 1, NL
! Recreate the CRX variable with the new values
! of XMASS_PF, which has been adjusted by DYN0
DO J = J1, J2
DO I = 1, IM
CRX(I,J,K) = XMASS_PF(I,J,K) / PU(I,J,K)
ENDDO
ENDDO
! Recreate the CRY variable with the new values
! of YMASS_PF, which has been adjusted by DYN0
DO J = J1, J2+1
D5 = 0.5 * COSE(J)
DO I = 1, IM
CRY(I,J,K) = YMASS_PF(I,J,K) /
& ( D5 * ( DELP2(I,J,K) + DELP2(I,J-1,K) ) )
ENDDO
ENDDO
ENDDO
ENDIF
!=================================================================
! End of TPCORE PRESSURE FIXER -- continue as usual
!=================================================================
C****6***0*********0*********0*********0*********0*********0**********72
C Find JN and JS
C****6***0*********0*********0*********0*********0*********0**********72
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope shared(JS,JN,CRX,CRY,PS2,U,V,DPI,ymass,delp2,PU)
CMIC$* shared(xmass)
CMIC$* private(i,j,k,sum1,sum2,D5)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K, SUM1, SUM2, D5 )
#endif
#endif
do 1000 k=1,NL
JS(k) = j1
JN(k) = j2
do 111 j=JS0,j1+1,-1
do 111 i=1,IM
if(abs(CRX(i,j,k)) .GT. 1.) then
JS(k) = j
go to 112
endif
111 continue
112 continue
do 122 j=JN0,j2-1
do 122 i=1,IM
if(abs(CRX(i,j,k)) .GT. 1.) then
JN(k) = j
go to 133
endif
122 continue
133 continue
C****6***0*********0*********0*********0*********0*********0**********72
C ***** Compute horizontal mass fluxes *****
C****6***0*********0*********0*********0*********0*********0**********72
C delp = pressure thickness: the psudo-density in a hydrostatic system.
do 30 j=1,JM
do 30 i=1,IM
30 delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j)
C N-S componenet
do j=j1,j2+1
D5 = 0.5 * COSE(j)
do i=1,IM
ymass(i,j,k) = CRY(i,j,k)*D5*(delp2(i,j,k) + delp2(i,j-1,k))
enddo
enddo
DO 75 j=j1,j2
DO 75 i=1,IM
75 DPI(i,j,k) = (ymass(i,j,k)-ymass(i,j+1,k)) * acosp(j)
if(j1.ne.2) then ! Enlarged polar cap.
do 95 i=1,IM
DPI(i, 2,k) = 0.
95 DPI(i,JM1,k) = 0.
endif
C Poles
sum1 = ymass(IM,j1 ,k)
sum2 = ymass(IM,j2+1,k)
do 98 i=1,IM-1
sum1 = sum1 + ymass(i,j1 ,k)
98 sum2 = sum2 + ymass(i,j2+1,k)
sum1 = - sum1 * RCAP
sum2 = sum2 * RCAP
do 100 i=1,IM
DPI(i, 1,k) = sum1
100 DPI(i,JM,k) = sum2
C E-W component
do j=j1,j2
do i=2,IM
PU(i,j,k) = 0.5 * (delp2(i,j,k) + delp2(i-1,j,k))
enddo
enddo
do j=j1,j2
PU(1,j,k) = 0.5 * (delp2(1,j,k) + delp2(IM,j,k))
enddo
DO 110 j=j1,j2
DO 110 i=1,IM
110 xmass(i,j,k) = PU(i,j,k)*CRX(i,j,k)
DO 120 j=j1,j2
DO 120 i=1,IM-1
120 DPI(i,j,k) = DPI(i,j,k) + xmass(i,j,k) - xmass(i+1,j,k)
DO 130 j=j1,j2
130 DPI(IM,j,k) = DPI(IM,j,k) + xmass(IM,j,k) - xmass(1,j,k)
C****6***0*********0*********0*********0*********0*********0**********72
C Compute Courant number at cell center
C****6***0*********0*********0*********0*********0*********0**********72
DO 135 j=2,JM1
do 135 i=1,IM-1
if(CRX(i,j,k)*CRX(i+1,j,k) .gt. 0.) then
if(CRX(i,j,k) .gt. 0.) then
U(i,j,k) = CRX(i,j,k)
else
U(i,j,k) = CRX(i+1,j,k)
endif
else
U(i,j,k) = 0.
endif
135 continue
i=IM
DO 136 j=2,JM1
if(CRX(i,j,k)*CRX(1,j,k) .gt. 0.) then
if(CRX(i,j,k) .gt. 0.) then
U(i,j,k) = CRX(i,j,k)
else
U(i,j,k) = CRX(1,j,k)
endif
else
U(i,j,k) = 0.
endif
136 continue
do 138 j=2,JM1
do 138 i=1,IM
if(CRY(i,j,k)*CRY(i,j+1,k) .gt. 0.) then
if(CRY(i,j,k) .gt. 0.) then
V(i,j,k) = CRY(i,j,k)
else
V(i,j,k) = CRY(i,j+1,k)
endif
else
V(i,j,k) = 0.
endif
138 continue
do 139 i=1,IMH
V(i, 1,k) = 0.5*(CRY(i,2,k)-CRY(i+IMH,2,k))
V(i+IMH, 1,k) = -V(i,1,k)
V(i, JM,k) = 0.5*(CRY(i,JM,k)-CRY(i+IMH,JM1,k))
139 V(i+IMH,JM,k) = -V(i,JM,k)
1000 continue
C****6***0*********0*********0*********0*********0*********0**********72
C Compute vertical mass flux (same dimensional unit as PS)
C****6***0*********0*********0*********0*********0*********0**********72
C compute total column mass CONVERGENCE.
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope shared(im,jm,DPI,PS1,PS2,W,DBK)
CMIC$* shared(DPI,PS1,PS2,W,DBK)
CMIC$* private(i,j,k,DG1)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K, DG1 )
#endif
#endif
do 395 j=1,jm
do 320 i=1,IM
320 DG1(i) = DPI(i,j,1)
do 330 k=2,NL
do 330 i=1,IM
DG1(i) = DG1(i) + DPI(i,j,k)
330 continue
do 360 i=1,IM
C Compute PS2 (PS at n+1) using the hydrostatic assumption.
C Changes (increases) to surface pressure = total column mass convergence
PS2(i,j) = PS1(i,j) + DG1(i)
C compute vertical mass flux from mass conservation principle.
W(i,j,1) = DPI(i,j,1) - DBK(1)*DG1(i)
W(i,j,NL) = 0.
360 continue
do 370 k=2,NL-1
do 370 i=1,IM
W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*DG1(i)
370 continue
395 continue
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all
CMIC$* shared(deform,NL,im,jm,delp,delp1,delp2,DPI,DAP,DBK,PS1,PS2)
CMIC$* private(i,j,k)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K )
#endif
#endif
DO 390 k=1,NL
DO 380 j=1,JM
DO 380 i=1,IM
delp1(i,j,k) = DAP(k) + DBK(k)*PS1(i,j)
delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j)
380 delp (i,j,k) = delp1(i,j,k) + DPI(i,j,k)
C Check deformation of the flow fields
if(deform) then
DO 385 j=1,JM
DO 385 i=1,IM
if(delp(i,j,k) .le. 0.) then
c write(6,*) k,'Noisy wind fields -> delp* is negative!'
c write(6,*) ' *** Smooth the wind fields or reduce NDT'
stop
endif
385 continue
endif
390 continue
C****6***0*********0*********0*********0*********0*********0**********72
C Do transport one tracer at a time.
C****6***0*********0*********0*********0*********0*********0**********72
DO 5000 IC=1,NC
! Moved initialization to 0 here (30/8/07, phs)
! For mass flux diagnostics (bey, 6/20/00)
fx1_tp(:,:,:) = 0d0
fy1_tp(:,:,:) = 0d0
fz1_tp(:,:,:) = 0d0
! Also need to initialize these arrays, so that the flux diagnostics
! will be identical for single or multi processor (bmy, 9/29/00)
fx(:,:,:) = 0d0
fy(:,:,:) = 0d0
fz(:,:,:) = 0d0
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* shared(q,DQ,delp1,U,V,j1,j2,JS,JN,im,jm,IML,IC,IORD,JORD)
CMIC$* shared(CRX,CRY,PU,xmass,ymass,fx,fy,acosp,rcap,qz)
CMIC$* shared(fx1_tp, fy1_tp)
CMIC$* private(i,j,k,jt,wk,DG2)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K, JT, WK, DG2 )
#endif
#endif
do 2500 k=1,NL
if(j1.ne.2) then
DO 405 I=1,IM
q(I, 2,k,IC) = q(I, 1,k,IC)
405 q(I,JM1,k,IC) = q(I,JM,k,IC)
endif
C Initialize DQ
DO 420 j=1,JM
DO 420 i=1,IM
420 DQ(i,j,k) = q(i,j,k,IC)*delp1(i,j,k)
C E-W advective cross term
call xadv(IM,JM,j1,j2,q(1,1,k,IC),U(1,1,k),JS(k),JN(k),IML,
& wk(1,1,1))
do 430 j=1,JM
do 430 i=1,IM
430 wk(i,j,1) = q(i,j,k,IC) + 0.5*wk(i,j,1)
C N-S advective cross term
do 66 j=j1,j2
do 66 i=1,IM
jt = float(j) - V(i,j,k)
66 wk(i,j,2) = V(i,j,k) * (q(i,jt,k,IC) - q(i,jt+1,k,IC))
do 77 j=j1,j2
do 77 i=1,IM
77 wk(i,j,2) = q(i,j,k,IC) + 0.5*wk(i,j,2)
C****6***0*********0*********0*********0*********0*********0**********72
C compute flux in E-W direction
C Return flux contribution from TPCORE in FX1_TP array (bey, 9/28/00)
call xtp(IM,JM,IML,j1,j2,JN(k),JS(k),PU(1,1,k),DQ(1,1,k),
& wk(1,1,2),CRX(1,1,k),fx(1,1,k),xmass(1,1,k),IORD,
& fx1_tp(:,:,k))
C compute flux in N-S direction
C Return flux contribution from TPCORE in FY1_TP array (bey, 9/28/00)
call ytp(IM,JM,j1,j2,acosp,RCAP,DQ(1,1,k),wk(1,1,1),
& CRY(1,1,k),DG2,ymass(1,1,k),WK(1,1,3),wk(1,1,4),
& WK(1,1,5),WK(1,1,6),fy(1,1,k),JORD,
& fy1_tp(:,:,k))
C****6***0*********0*********0*********0*********0*********0**********72
if(ZCROSS) then
C qz is the horizontal advection modified value for input to the
C vertical transport operator FZPPM
C Note: DQ contains only first order upwind contribution.
do 88 j=1,JM
do 88 i=1,IM
88 qz(i,j,k) = DQ(i,j,k) / delp(i,j,k)
else
do 99 j=1,JM
do 99 i=1,IM
99 qz(i,j,k) = q(i,j,k,IC)
endif
2500 continue ! k-loop
C****6***0*********0*********0*********0*********0*********0**********72
C Compute fluxes in the vertical direction
C Return flux contribution from FZPPM in FZ1_TP for ND26 (bey, 9/28/00)
call FZPPM(qz,fz,IM,JM,NL,DQ,W,delp,KORD,fz1_tp)
C****6***0*********0*********0*********0*********0*********0**********72
if( MFCT ) then
C qlow is the low order "monotonic" solution
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all
CMIC$* shared(NL,im,jm,j1,jm1,qlow,DQ,delp2)
CMIC$* private(i,j,k)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K )
#endif
#endif
DO k=1,NL
DO 560 j=1,JM
DO 560 i=1,IM
560 qlow(i,j,k) = DQ(i,j,k) / delp2(i,j,k)
if(j1.ne.2) then
DO 561 i=1,IM
qlow(i, 2,k) = qlow(i, 1,k)
qlow(i,JM1,k) = qlow(i,JM,k)
561 CONTINUE
endif
enddo
C****6***0*********0*********0*********0*********0*********0**********72
call FCT3D(Q(1,1,1,IC),qlow,fx,fy,fz,IM,JM,NL,j1,j2,delp2,
& DPI,qz,wk,Qmax,Qmin,DG2,U,V,acosp,RCAP)
C Note: Q is destroyed!!!
C****6***0*********0*********0*********0*********0*********0**********72
ENDIF
C Final update
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* private(i,j,k,sum1,sum2)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K, SUM1, SUM2 )
#endif
#endif
do 101 k=1,NL
do 425 j=j1,j2
do 425 i=1,IM
DQ(i,j,k) = DQ(i,j,k) + fx(i,j,k) - fx(i+1,j,k)
& + (fy(i,j,k) - fy(i,j+1,k))*acosp(j)
& + fz(i,j,k) - fz(i,j,k+1)
425 continue
sum1 = fy(IM,j1 ,k)
sum2 = fy(IM,J2+1,k)
do i=1,IM-1
sum1 = sum1 + fy(i,j1 ,k)
sum2 = sum2 + fy(i,J2+1,k)
enddo
DQ(1, 1,k) = DQ(1, 1,k) - sum1*RCAP + fz(1, 1,k) - fz(1, 1,k+1)
DQ(1,JM,k) = DQ(1,JM,k) + sum2*RCAP + fz(1,JM,k) - fz(1,JM,k+1)
do i=2,IM
DQ(i, 1,k) = DQ(1, 1,k)
DQ(i,JM,k) = DQ(1,JM,k)
enddo
101 continue
!=================================================================
! bey, 6/20/00. for mass-flux diagnostic
! NOTE: DIAG_FLUX is not called within a parallel loop,
! so parallelization can be done within the subroutine
!=================================================================
CALL DIAG_FLUX( IC, FX, FX1_TP, FY, FY1_TP,
& FZ, FZ1_TP, NDT, ACOSP )
C****6***0*********0*********0*********0*********0*********0**********72
if(FILL) call qckxyz(DQ,DG2,IM,JM,NL,j1,j2,cosp,acosp,IC,NSTEP)
C****6***0*********0*********0*********0*********0*********0**********72
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all
CMIC$* shared(q,IC,NL,j1,im,jm,jm1,DQ,delp2)
CMIC$* private(i,j,k)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K )
#endif
#endif
DO k=1,NL
DO 447 j=1,JM
DO 447 i=1,IM
447 Q(i,j,k,IC) = DQ(i,j,k) / delp2(i,j,k)
if(j1.ne.2) then
DO 450 I=1,IM
Q(I, 2,k,IC) = Q(I, 1,k,IC)
Q(I,JM1,k,IC) = Q(I,JM,k,IC)
450 CONTINUE
endif
enddo
5000 continue
RETURN
END SUBROUTINE TPCORE
!------------------------------------------------------------------------------
subroutine cosa(cosp,cose,JM,PI,DP) 2
C****6***0*********0*********0*********0*********0*********0**********72
REAL cosp(*),cose(*),sine(JM)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
do 10 j=2,JM
ph5 = -0.5*PI + (float(j-1)-0.5)*DP
10 sine(j) = SIN(ph5)
do 80 J=2,JM-1
80 cosp(J) = (sine(j+1)-sine(j))/DP
cosp( 1) = 0.
cosp(JM) = 0.
C Define cosine at edges..
do 90 j=2,JM
90 cose(j) = 0.5 * (cosp(j-1)+cosp(j))
cose(1) = cose(2)
return
end subroutine cosa
!------------------------------------------------------------------------------
subroutine cosc(cosp,cose,JNP,PI,DP) 2
C****6***0*********0*********0*********0*********0*********0**********72
REAL cosp(*),cose(*)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
phi = -0.5*PI
do 55 j=2,JNP-1
phi = phi + DP
55 cosp(j) = cos(phi)
cosp( 1) = 0.
cosp(JNP) = 0.
C
do 66 j=2,JNP
cose(j) = 0.5*(cosp(j)+cosp(j-1))
66 CONTINUE
C
do 77 j=2,JNP-1
cosp(j) = 0.5*(cose(j)+cose(j+1))
77 CONTINUE
return
end subroutine cosc
!------------------------------------------------------------------------------
subroutine FCT3D(P,plow,fx,fy,fz,im,jm,km,j1,j2,delp,adx,ady, 1,2
& wk1,Qmax,Qmin,wkx,CRX,CRY,acosp,RCAP)
C****6***0*********0*********0*********0*********0*********0**********72
! Added to pass C-preprocessor switches (bmy, 3/9/01)
# include "define.h"
C MFCT Limiter
C plow: low order solution matrix
C P: current solution matrix
PARAMETER (esl = 1.E-30)
REAL P(IM,JM,km),CRX(IM,JM,km),CRY(IM,JM,km),plow(IM,JM,km),
& Qmax(IM,JM,km),Qmin(IM,JM,km),acosp(*),delp(im,jm,km),
& adx(IM,JM,km),ady(IM,JM,km),fx(IM+1,JM,km),
& fy(IM,JM,km),fz(im,jm,km+1),wk1(IM,JM,km),
& wkx(im,jm),wkn(im,jm)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
JM1 = JM-1
C Find local min/max of the low-order monotone solution
call hilo3D(P,im,jm,km,j1,j2,adx,ady,Qmax,Qmin,wkx,wkn)
call hilo3D(plow,im,jm,km,j1,j2,Qmax,Qmin,wk1,P,wkx,wkn)
C P is destroyed!
C GOTO 123
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* shared(im,j1,j2,km,CRX,CRY,adx,ady,Qmax,Qmin)
CMIC$* private(i,j,k,IT,JT,PS1,PS2,PN1,PN2)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K, IT, JT, PS1, PS2, PN1, PN2 )
#endif
#endif
DO 1000 k=1,km
do j=j1,j2
DO i=1,IM
IT = NINT( float(i) - CRX(i,j,k) )
C Wrap around in E-W
if(IT .lt. 1) then
IT = IM + IT
elseif(IT .GT. IM) then
IT = IT - IM
endif
JT = NINT( float(j) - CRY(i,j,k) )
Qmax(i,j,k) = max(Qmax(i,j,k), adx(IT,JT,k))
Qmin(i,j,k) = min(Qmin(i,j,k), ady(IT,JT,k))
enddo
enddo
C Poles:
PS1 = max(Qmax(1, 1,k), adx(1, 1,k))
PS2 = min(Qmin(1, 1,k), ady(1, 1,k))
PN1 = max(Qmax(1,JM,k), adx(1,JM,k))
PN2 = min(Qmin(1,JM,k), ady(1,JM,k))
DO i=1,IM
Qmax(i, 1,k) = PS1
Qmin(i, 1,k) = PS2
Qmax(i,JM,k) = PN1
Qmin(i,JM,k) = PN2
enddo
1000 continue
123 continue
C Flux Limiter
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* shared(adx,ady,fx,fy,fz,plow,Qmax,Qmin,delp)
CMIC$* private(wkx,wkn)
CMIC$* private(i,j,k,ain,aou,bin,bou,cin,cou,btop,bdon)
#else
!$OMP PARALLEL DO PRIVATE( WKX, WKN, I, J, K, AIN, AOU,
!$OMP+ BIN, BOU, CIN, COU, BTOP, BDON )
#endif
#endif
DO 2000 k=1,km
DO j=j1,j2
DO i=1,IM
if(fx(i,j,k) .gt. 0.) then
Ain = fx(i,j,k)
Aou = 0.
else
Ain = 0.
Aou = -fx(i,j,k)
endif
if(fx(i+1,j,k) .gt. 0.) then
Aou = Aou + fx(i+1,j,k)
else
Ain = Ain - fx(i+1,j,k)
endif
if(fy(i,j,k) .gt. 0.) then
Bin = fy(i,j,k)
Bou = 0.
else
Bin = 0.
Bou = -fy(i,j,k)
endif
if(fy(i,j+1,k) .gt. 0.) then
Bou = Bou + fy(i,j+1,k)
else
Bin = Bin - fy(i,j+1,k)
endif
if(fz(i,j,k) .gt. 0.) then
Cin = fz(i,j,k)
Cou = 0.
else
Cin = 0.
Cou = -fz(i,j,k)
endif
if(fz(i,j,k+1) .gt. 0.) then
Cou = Cou + fz(i,j,k+1)
else
Cin = Cin - fz(i,j,k+1)
endif
C****6***0*********0*********0*********0*********0*********0**********72
wkx(i,j) = Ain + Bin*acosp(j) + Cin
wkn(i,j) = Aou + Bou*acosp(j) + Cou
C****6***0*********0*********0*********0*********0*********0**********72
enddo
enddo
DO j=j1,j2
DO i=1,IM
adx(i,j,k) = delp(i,j,k)*(Qmax(i,j,k)-plow(i,j,k))/(wkx(i,j)+esl)
ady(i,j,k) = delp(i,j,k)*(plow(i,j,k)-Qmin(i,j,k))/(wkn(i,j)+esl)
enddo
enddo
C S Pole
Ain = 0.
Aou = 0.
DO i=1,IM
if(fy(i,j1,k).gt. 0.) then
Aou = Aou + fy(i,j1,k)
else
Ain = Ain + fy(i,j1,k)
endif
enddo
Ain = -Ain * RCAP
Aou = Aou * RCAP
C add vertical contribution...
i=1
j=1
if(fz(i,j,k) .gt. 0.) then
Cin = fz(i,j,k)
Cou = 0.
else
Cin = 0.
Cou = -fz(i,j,k)
endif
if(fz(i,j,k+1) .gt. 0.) then
Cou = Cou + fz(i,j,k+1)
else
Cin = Cin - fz(i,j,k+1)
endif
C****6***0*********0*********0*********0*********0*********0**********72
btop = delp(1,1,k)*(Qmax(1,1,k)-plow(1,1,k))/(Ain+Cin+esl)
bdon = delp(1,1,k)*(plow(1,1,k)-Qmin(1,1,k))/(Aou+Cou+esl)
C****6***0*********0*********0*********0*********0*********0**********72
DO i=1,IM
adx(i,j,k) = btop
ady(i,j,k) = bdon
enddo
C N Pole
J=JM
Ain = 0.
Aou = 0.
DO i=1,IM
if(fy(i,j2+1,k).gt. 0.) then
Ain = Ain + fy(i,j2+1,k)
else
Aou = Aou + fy(i,j2+1,k)
endif
enddo
Ain = Ain * RCAP
Aou = -Aou * RCAP
C add vertical contribution...
i=1
if(fz(i,j,k) .gt. 0.) then
Cin = fz(i,j,k)
Cou = 0.
else
Cin = 0.
Cou = -fz(i,j,k)
endif
if(fz(i,j,k+1) .gt. 0.) then
Cou = Cou + fz(i,j,k+1)
else
Cin = Cin - fz(i,j,k+1)
endif
C****6***0*********0*********0*********0*********0*********0**********72
btop = delp(1,j,k)*(Qmax(1,j,k)-plow(1,j,k))/(Ain+Cin+esl)
bdon = delp(1,j,k)*(plow(1,j,k)-Qmin(1,j,k))/(Aou+Cou+esl)
C****6***0*********0*********0*********0*********0*********0**********72
DO i=1,IM
adx(i,j,k) = btop
ady(i,j,k) = bdon
enddo
if(j1 .ne. 2) then
DO i=1,IM
C SP
adx(i,2,k) = adx(i,1,k)
ady(i,2,k) = ady(i,1,k)
C NP
adx(i,JM1,k) = adx(i,JM,k)
ady(i,JM1,k) = ady(i,JM,k)
enddo
endif
2000 continue
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* shared(fz,adx,ady,im,jm,km)
CMIC$* private(i,j,k)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K )
#endif
#endif
DO 3000 k=1,km
DO j=j1,j2
do i=2,IM
if(fx(i,j,k) .gt. 0.) then
fx(i,j,k) = min(1.,ady(i-1,j,k),adx(i,j,k))*fx(i,j,k)
else
fx(i,j,k) = min(1.,adx(i-1,j,k),ady(i,j,k))*fx(i,j,k)
endif
enddo
enddo
C For i=1
DO j=j1,j2
if(fx(1,j,k) .gt. 0.) then
fx(1,j,k) = min(1.,ady(IM,j,k),adx(1,j,k))*fx(1,j,k)
else
fx(1,j,k) = min(1.,adx(IM,j,k),ady(1,j,k))*fx(1,j,k)
endif
fx(IM+1,j,k) = fx(1,j,k)
enddo
do j=j1,j2+1
do i=1,IM
if(fy(i,j,k) .gt. 0.) then
fy(i,j,k) = min(1.,ady(i,j-1,k),adx(i,j,k))*fy(i,j,k)
else
fy(i,j,k) = min(1.,adx(i,j-1,k),ady(i,j,k))*fy(i,j,k)
endif
enddo
enddo
if(k .ne. 1) then
do j=1,jm
do i=1,im
if(fz(i,j,k) .gt. 0.) then
fz(i,j,k) = min(1.,ady(i,j,k-1),adx(i,j,k))*fz(i,j,k)
else
fz(i,j,k) = min(1.,adx(i,j,k-1),ady(i,j,k))*fz(i,j,k)
endif
enddo
enddo
endif
3000 continue
return
end subroutine fct3d
!------------------------------------------------------------------------------
subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny) 2
C****6***0*********0*********0*********0*********0*********0**********72
REAL q(IMR,*),qtmp(JNP,IMR)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
ipx = 0
C Copy & swap direction for vectorization.
do 25 i=1,imr
do 25 j=j1,j2
25 qtmp(j,i) = q(i,j)
C
do 55 i=2,imr-1
do 55 j=j1,j2
if(qtmp(j,i).lt.0.) then
ipx = 1
c west
d0 = max(0.,qtmp(j,i-1))
d1 = min(-qtmp(j,i),d0)
qtmp(j,i-1) = qtmp(j,i-1) - d1
qtmp(j,i) = qtmp(j,i) + d1
c east
d0 = max(0.,qtmp(j,i+1))
d2 = min(-qtmp(j,i),d0)
qtmp(j,i+1) = qtmp(j,i+1) - d2
qtmp(j,i) = qtmp(j,i) + d2 + tiny
endif
55 continue
c
i=1
do 65 j=j1,j2
if(qtmp(j,i).lt.0.) then
ipx = 1
c west
d0 = max(0.,qtmp(j,imr))
d1 = min(-qtmp(j,i),d0)
qtmp(j,imr) = qtmp(j,imr) - d1
qtmp(j,i) = qtmp(j,i) + d1
c east
d0 = max(0.,qtmp(j,i+1))
d2 = min(-qtmp(j,i),d0)
qtmp(j,i+1) = qtmp(j,i+1) - d2
c
qtmp(j,i) = qtmp(j,i) + d2 + tiny
endif
65 continue
i=IMR
do 75 j=j1,j2
if(qtmp(j,i).lt.0.) then
ipx = 1
c west
d0 = max(0.,qtmp(j,i-1))
d1 = min(-qtmp(j,i),d0)
qtmp(j,i-1) = qtmp(j,i-1) - d1
qtmp(j,i) = qtmp(j,i) + d1
c east
d0 = max(0.,qtmp(j,1))
d2 = min(-qtmp(j,i),d0)
qtmp(j,1) = qtmp(j,1) - d2
c
qtmp(j,i) = qtmp(j,i) + d2 + tiny
endif
75 continue
C
if(ipx.ne.0) then
do 85 j=j1,j2
do 85 i=1,imr
85 q(i,j) = qtmp(j,i)
else
C Pole
if(q(1,1).lt.0. or. q(1,JNP).lt.0.) ipx = 1
endif
return
end subroutine filew
!------------------------------------------------------------------------------
subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) 2
C****6***0*********0*********0*********0*********0*********0**********72
REAL q(IMR,*),cosp(*),acosp(*)
LOGICAL first
DATA first /.true./
SAVE cap1
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
if(first) then
DP = 4.*ATAN(1.)/float(JNP-1)
cap1 = IMR*(1.-COS((j1-1.5)*DP))/DP
first = .false.
endif
C
ipy = 0
do 55 j=j1+1,j2-1
DO 55 i=1,IMR
IF(q(i,j).LT.0.) THEN
ipy = 1
dq = - q(i,j)*cosp(j)
C North
dn = q(i,j+1)*cosp(j+1)
d0 = max(0.,dn)
d1 = min(dq,d0)
q(i,j+1) = (dn - d1)*acosp(j+1)
dq = dq - d1
C South
ds = q(i,j-1)*cosp(j-1)
d0 = max(0.,ds)
d2 = min(dq,d0)
q(i,j-1) = (ds - d2)*acosp(j-1)
q(i,j) = (d2 - dq)*acosp(j) + tiny
endif
55 continue
C
do i=1,imr
IF(q(i,j1).LT.0.) THEN
ipy = 1
dq = - q(i,j1)*cosp(j1)
C North
dn = q(i,j1+1)*cosp(j1+1)
d0 = max(0.,dn)
d1 = min(dq,d0)
q(i,j1+1) = (dn - d1)*acosp(j1+1)
q(i,j1) = (d1 - dq)*acosp(j1) + tiny
endif
enddo
C
j = j2
do i=1,imr
IF(q(i,j).LT.0.) THEN
ipy = 1
dq = - q(i,j)*cosp(j)
C South
ds = q(i,j-1)*cosp(j-1)
d0 = max(0.,ds)
d2 = min(dq,d0)
q(i,j-1) = (ds - d2)*acosp(j-1)
q(i,j) = (d2 - dq)*acosp(j) + tiny
endif
enddo
C
C Check Poles.
if(q(1,1).lt.0.) then
dq = q(1,1)*cap1/float(IMR)*acosp(j1)
do i=1,imr
q(i,1) = 0.
q(i,j1) = q(i,j1) + dq
if(q(i,j1).lt.0.) ipy = 1
enddo
endif
C
if(q(1,JNP).lt.0.) then
dq = q(1,JNP)*cap1/float(IMR)*acosp(j2)
do i=1,imr
q(i,JNP) = 0.
q(i,j2) = q(i,j2) + dq
if(q(i,j2).lt.0.) ipy = 1
enddo
endif
C
return
end subroutine filns
!------------------------------------------------------------------------------
subroutine fxppm(IMR,IML,UT,P,DC,fx1,fx2,IORD) 4,4
C****6***0*********0*********0*********0*********0*********0**********72
PARAMETER ( R3 = 1./3., R23 = 2./3. )
REAL UT(*),fx1(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1)
REAL AR(0:IMR),AL(0:IMR),A6(0:IMR),fx2(*)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
LMT = IORD - 3
C
DO 10 i=1,IMR
10 AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
C
do 20 i=1,IMR-1
20 AR(i) = AL(i+1)
AR(IMR) = AL(1)
C
do 30 i=1,IMR
30 A6(i) = 3.*(p(i)+p(i) - (AL(i)+AR(i)))
C
if(LMT.LE.2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
C
AL(0) = AL(IMR)
AR(0) = AR(IMR)
A6(0) = A6(IMR)
C
C Abs(UT(i)) < 1
DO i=1,IMR
IF(UT(i).GT.0.) then
fx1(i) = P(i-1)
fx2(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) +
& A6(i-1)*(1.-R23*UT(i)) )
else
fx1(i) = P(i)
fx2(i) = AL(i) - 0.5*UT(i)*(AR(i) - AL(i) +
& A6(i)*(1.+R23*UT(i)))
endif
enddo
C
DO i=1,IMR
fx2(i) = fx2(i) - fx1(i)
enddo
return
end subroutine fxppm
!------------------------------------------------------------------------------
subroutine fyppm(C,P,DC,fy1,fy2,IMR,JNP,j1,j2,A6,AR,AL,JORD) 3,3
C****6***0*********0*********0*********0*********0*********0**********72
PARAMETER ( R3 = 1./3., R23 = 2./3. )
REAL C(IMR,*),fy1(IMR,*),P(IMR,*),DC(IMR,*),fy2(IMR,JNP)
REAL AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
IMH = IMR / 2
JMR = JNP - 1
j11 = j1-1
IMJM1 = IMR*(J2-J1+2)
len = IMR*(J2-J1+3)
LMT = JORD - 3
C
DO 10 i=1,IMR*JMR
AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3
AR(i,1) = AL(i,2)
10 CONTINUE
C
C Poles:
C
DO i=1,IMH
AL(i,1) = AL(i+IMH,2)
AL(i+IMH,1) = AL(i,2)
C
AR(i,JNP) = AR(i+IMH,JMR)
AR(i+IMH,JNP) = AR(i,JMR)
enddo
C
do 30 i=1,len
30 A6(i,j11) = 3.*(p(i,j11)+p(i,j11) - (AL(i,j11)+AR(i,j11)))
C
if(LMT.le.2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11),
& AL(1,j11),P(1,j11),len,LMT)
C
DO 140 i=1,IMJM1
IF(C(i,j1).GT.0.) then
fy1(i,j1) = P(i,j11)
fy2(i,j1) = AR(i,j11) + 0.5*C(i,j1)*(AL(i,j11) - AR(i,j11) +
& A6(i,j11)*(1.-R23*C(i,j1)) )
else
fy1(i,j1) = P(i,j1)
fy2(i,j1) = AL(i,j1) - 0.5*C(i,j1)*(AR(i,j1) - AL(i,j1) +
& A6(i,j1)*(1.+R23*C(i,j1)))
endif
140 continue
c
DO i=1,IMJM1
fy2(i,j1) = fy2(i,j1) - fy1(i,j1)
ENDDO
return
end subroutine fyppm
!------------------------------------------------------------------------------
subroutine FZPPM(P,fz,IMR,JNP,NL,DQ,WZ,delp,KORD,fz1_tp) 2,6
C****6***0*********0*********0*********0*********0*********0**********72
! Added to pass C-preprocessor switches (bmy, 3/9/01)
# include "define.h"
PARAMETER ( R23 = 2./3., R3 = 1./3.)
REAL WZ(IMR,JNP,NL),P(IMR,JNP,NL),DQ(IMR,JNP,NL),
& fz(IMR,JNP,NL+1),delp(IMR,JNP,NL)
C local 2d arrays
REAL AR(IMR,NL),AL(IMR,NL),A6(IMR,NL),delq(IMR,NL),DC(IMR,NL)
! bey, 6/20/00. for mass-flux diagnostic
real fz1_tp(IMR,JNP,NL)
real lac
c real x, y, z
c real median
c median(x,y,z) = min(max(x,y), max(y,z), max(z,x))
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
km = NL
km1 = NL-1
LMT = max(KORD - 3, 0)
C find global min/max
! VMAX1D causes bus errors on SGI. Replace with F90 intrinsic
! functions "MAXVAL" and "MINVAL". These functions produced
! identical results as vmax1d in testing. "MAXVAL" and "MINVAL"
! should also execute more efficiently as well. (bmy, 4/24/00)
Tmax = MAXVAL( P(:,:,1) )
Tmin = MINVAL( P(:,:,1) )
Bmax = MAXVAL( P(:,:,NL) )
Bmin = MINVAL( P(:,:,NL) )
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* shared(LMT,Tmax,Tmin,Bmax,Bmin,JNP,IMR)
CMIC$* shared(fz,DQ,WZ,fz1_tp)
CMIC$* private(i,j,k,c1,c2,tmp,qmax,qmin,A1,A2,d1,d2,qm,dp,c3)
CMIC$* private(cmax,cmin,DC,delq,AR,AL,A6,CM,CP, qmp, lac)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K, C1, C2, TMP, QMAX, QMIN, A1, A2,
!$OMP+ D1, D2, QM, DP, C3, CMAX, CMIN, DC, DELQ,
!$OMP+ AR, AL, A6, CM, CP, QMP, LAC )
#endif
#endif
do 4000 j=1,JNP
do 500 k=2,km
do 500 i=1,IMR
500 A6(i,k) = delp(i,j,k-1) + delp(i,j,k)
do 1000 k=1,km1
do 1000 i=1,IMR
1000 delq(i,k) = P(i,j,k+1) - P(i,j,k)
DO 1220 k=2,km1
DO 1220 I=1,IMR
c1 = (delp(i,j,k-1)+0.5*delp(i,j,k))/A6(i,k+1)
c2 = (delp(i,j,k+1)+0.5*delp(i,j,k))/A6(i,k)
tmp = delp(i,j,k)*(c1*delq(i,k) + c2*delq(i,k-1))
& / (A6(i,k)+delp(i,j,k+1))
Qmax = max(P(i,j,k-1),P(i,j,k),P(i,j,k+1)) - P(i,j,k)
Qmin = P(i,j,k) - min(P(i,j,k-1),P(i,j,k),P(i,j,k+1))
DC(i,k) = sign(min(abs(tmp),Qmax,Qmin), tmp)
1220 CONTINUE
C****6***0*********0*********0*********0*********0*********0**********72
C Compute the first guess at cell interface
C First guesses are required to be continuous.
C****6***0*********0*********0*********0*********0*********0**********72
C Interior.
DO 12 k=3,km1
DO 12 i=1,IMR
c1 = delq(i,k-1)*delp(i,j,k-1) / A6(i,k)
A1 = A6(i,k-1) / (A6(i,k) + delp(i,j,k-1))
A2 = A6(i,k+1) / (A6(i,k) + delp(i,j,k))
AL(i,k) = P(i,j,k-1) + c1 + 2./(A6(i,k-1)+A6(i,k+1)) *
& ( delp(i,j,k )*(c1*(A1 - A2)+A2*DC(i,k-1)) -
& delp(i,j,k-1)*A1*DC(i,k ) )
12 CONTINUE
C Area preserving cubic with 2nd deriv. = 0 at the boundaries
C Top
DO 10 i=1,IMR
d1 = delp(i,j,1)
d2 = delp(i,j,2)
qm = (d2*P(i,j,1)+d1*P(i,j,2)) / (d1+d2)
dp = 2.*(P(i,j,2)-P(i,j,1)) / (d1+d2)
c1 = 4.*(AL(i,3)-qm-d2*dp) / ( d2*(2.*d2*d2+d1*(d2+3.*d1)) )
c3 = dp - 0.5*c1*(d2*(5.*d1+d2)-3.*d1**2)
AL(i,2) = qm - 0.25*c1*d1*d2*(d2+3.*d1)
AL(i,1) = d1*(2.*c1*d1**2-c3) + AL(i,2)
DC(i,1) = P(i,j,1) - AL(i,1)
C No over- and undershoot condition
AL(i,1) = max(Tmin,AL(i,1))
AL(i,1) = min(Tmax,AL(i,1))
Cmax = max(P(i,j,1), P(i,j,2))
Cmin = min(P(i,j,1), P(i,j,2))
AL(i,2) = max(Cmin,AL(i,2))
AL(i,2) = min(Cmax,AL(i,2))
10 continue
C Bottom
DO 15 i=1,IMR
d1 = delp(i,j,km )
d2 = delp(i,j,km1)
qm = (d2*P(i,j,km)+d1*P(i,j,km1)) / (d1+d2)
dp = 2.*(P(i,j,km1)-P(i,j,km)) / (d1+d2)
c1 = 4.*(AL(i,km1)-qm-d2*dp) / (d2*(2.*d2*d2+d1*(d2+3.*d1)))
c3 = dp - 0.5*c1*(d2*(5.*d1+d2)-3.*d1**2)
AL(i,km) = qm - 0.25*c1*d1*d2*(d2+3.*d1)
AR(i,km) = d1*(2.*c1*d1**2-c3) + AL(i,km)
DC(i,km) = AR(i,km) - P(i,j,km)
C No over- and undershoot condition
Cmax = max(P(i,j,km), P(i,j,km1))
Cmin = min(P(i,j,km), P(i,j,km1))
AL(i,km) = max(Cmin,AL(i,km))
AL(i,km) = min(Cmax,AL(i,km))
AR(i,km) = max(Bmin,AR(i,km))
AR(i,km) = min(Bmax,AR(i,km))
15 continue
do 20 k=1,km1
do 20 i=1,IMR
AR(i,k) = AL(i,k+1)
20 continue
C f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 )
C Top 2 layers
do k=1,2
do i=1,IMR
A6(i,k) = 3.*(P(i,j,k)+P(i,j,k) - (AL(i,k)+AR(i,k)))
enddo
call lmtppm(DC(1,k),A6(1,k),AR(1,k),AL(1,k),P(1,j,k),
& IMR,0)
enddo
C Interior.
if(LMT.LE.2) then
do k=3,NL-2
do i=1,IMR
A6(i,k) = 3.*(P(i,j,k)+P(i,j,k) - (AL(i,k)+AR(i,k)))
enddo
call lmtppm(DC(1,k),A6(1,k),AR(1,k),AL(1,k),P(1,j,k),
& IMR,LMT)
enddo
elseif(LMT .eq. 4) then
c****6***0*********0*********0*********0*********0*********0**********72
C Huynh's 2nd constraint
c****6***0*********0*********0*********0*********0*********0**********72
do k=2, NL-1
do i=1,imr
DC(i,k) = delq(i,k) - delq(i,k-1)
enddo
enddo
do k=3, NL-2
do i=1, imr
C Right edges
qmp = P(i,j,k) + 2.0*delq(i,k-1)
lac = P(i,j,k) + 1.5*DC(i,k-1) + 0.5*delq(i,k-1)
qmin = min(P(i,j,k), qmp, lac)
qmax = max(P(i,j,k), qmp, lac)
c AR(i,k) = median(AR(i,k), qmin, qmax)
AR(i,k) = min(max(AR(i,k), qmin), qmax)
C Left edges
qmp = P(i,j,k) - 2.0*delq(i,k)
lac = P(i,j,k) + 1.5*DC(i,k+1) - 0.5*delq(i,k)
qmin = min(P(i,j,k), qmp, lac)
qmax = max(P(i,j,k), qmp, lac)
c AL(i,k) = median(AL(i,k), qmin, qmax)
AL(i,k) = min(max(AL(i,k), qmin), qmax)
C Recompute A6
A6(i,k) = 3.*(2.*P(i,j,k) - (AR(i,k)+AL(i,k)))
enddo
enddo
endif
C Bottom 2 layers
do k=NL-1,NL
do i=1,IMR
A6(i,k) = 3.*(P(i,j,k)+P(i,j,k) - (AL(i,k)+AR(i,k)))
enddo
call lmtppm(DC(1,k),A6(1,k),AR(1,k),AL(1,k),P(1,j,k),
& IMR,0)
enddo
DO 140 k=2,NL
DO 140 i=1,IMR
IF(WZ(i,j,k-1).GT.0.) then
CM = WZ(i,j,k-1) / delp(i,j,k-1)
DC(i,k) = P(i,j,k-1)
fz(i,j,k) = AR(i,k-1)+0.5*CM*(AL(i,k-1)-AR(i,k-1)+
& A6(i,k-1)*(1.-R23*CM))
else
CP = WZ(i,j,k-1) / delp(i,j,k)
DC(i,k) = P(i,j,k)
fz(i,j,k) = AL(i,k)+0.5*CP*(AL(i,k)-AR(i,k)-
& A6(i,k)*(1.+R23*CP))
endif
140 continue
DO 250 k=2,NL
DO 250 i=1,IMR
fz(i,j,k) = WZ(i,j,k-1) * (fz(i,j,k) - DC(i,k))
DC(i,k) = WZ(i,j,k-1) * DC(i,k)
250 continue
do 350 i=1,IMR
fz(i,j, 1) = 0.
fz(i,j,NL+1) = 0.
DQ(i,j, 1) = DQ(i,j, 1) - DC(i, 2)
DQ(i,j,NL) = DQ(i,j,NL) + DC(i,NL)
fz1_tp(i,j,1) = 0. ! PHS
fz1_tp(i,j,NL) = DC(i,NL) ! PHS - flux b/w 1st and second layer
350 continue
!-----------------------------------------------------------------------------
! bey, 6/20/00. for mass-flux diagnostic, loop had to be extended
! do 360 k=2,km1
! do 360 i=1,IMR
!360 DQ(i,j,k) = DQ(i,j,k) + DC(i,k) - DC(i,k+1)
!-----------------------------------------------------------------------------
do k=2,km1
do i=1,IMR
DQ(i,j,k) = DQ(i,j,k) + DC(i,k) - DC(i,k+1)
! bey, 6/20/00. for mass-flux diagnostic
fz1_tp(i,j,k) = DC(i,k)
enddo
enddo
4000 continue
return
end subroutine fzppm
!------------------------------------------------------------------------------
subroutine hilo(q,im,jm,j1,j2,qmax,qmin,bt,bd) 2
C****6***0*********0*********0*********0*********0*********0**********72
REAL q(IM,JM),Qmax(IM,JM),Qmin(IM,JM),bt(IM,*),bd(IM,*)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C y-sweep
DO j=j1,j2
DO i=1,IM
bt(i,j) = max(q(i,j-1),q(i,j),q(i,j+1))
bd(i,j) = min(q(i,j-1),q(i,j),q(i,j+1))
enddo
enddo
C
C x-sweep
IM1 = IM-1
DO j=j1,j2
DO i=2,IM1
Qmax(i,j) = max(bt(i-1,j),bt(i,j),bt(i+1,j))
Qmin(i,j) = min(bd(i-1,j),bd(i,j),bd(i+1,j))
enddo
enddo
C
DO j=j1,j2
C i = 1
Qmax(1,j) = max(bt(IM,j),bt(1,j),bt(2,j))
Qmin(1,j) = min(bd(IM,j),bd(1,j),bd(2,j))
C i = IM
Qmax(IM,j) = max(bt(IM1,j),bt(IM,j),bt(1,j))
Qmin(IM,j) = min(bd(IM1,j),bd(IM,j),bd(1,j))
enddo
C
C N. Pole:
Pmax = q(1,JM)
Pmin = q(1,JM)
do i=1,IM
if(q(i,j2) .gt. Pmax) then
Pmax = q(i,j2)
elseif(q(i,j2) .lt. Pmin) then
Pmin = q(i,j2)
endif
enddo
C
do i=1,IM
Qmax(i,JM) = Pmax
Qmin(i,JM) = Pmin
enddo
C
C S. Pole:
Pmax = q(1,1)
Pmin = q(1,1)
do i=1,IM
if(q(i,j1) .gt. Pmax) then
Pmax = q(i,j1)
elseif(q(i,j1) .lt. Pmin) then
Pmin = q(i,j1)
endif
enddo
C
do i=1,IM
Qmax(i,1) = Pmax
Qmin(i,1) = Pmin
enddo
C
if(j1 .ne. 2) then
JM1 = JM-1
do i=1,IM
Qmax(i,2) = Qmax(i,1)
Qmin(i,2) = Qmin(i,1)
C
Qmax(i,JM1) = Qmax(i,JM)
Qmin(i,JM1) = Qmin(i,JM)
enddo
endif
return
end subroutine hilo
!------------------------------------------------------------------------------
subroutine hilo3D(P,im,jm,km,j1,j2,Pmax,Pmin,Qmax,Qmin,bt,bd) 2,2
C****6***0*********0*********0*********0*********0*********0**********72
! Added to pass C-preprocessor switches (bmy, 3/9/01)
# include "define.h"
REAL P(IM,JM,km),Pmax(IM,JM,km),Pmin(IM,JM,km),
& Qmax(IM,JM,km),Qmin(IM,JM,km),bt(im,jm),bd(im,jm)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* shared(P,Qmax,Qmin,im,jm,j1,j2)
CMIC$* private(k,bt,bd)
#else
!$OMP PARALLEL DO PRIVATE( K, BT, BD )
#endif
#endif
DO 1000 k=1,km
call hilo(P(1,1,k),im,jm,j1,j2,Qmax(1,1,k),Qmin(1,1,k),bt,bd)
1000 continue
km1 = km-1
km2 = km-2
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* shared(Pmax,Pmin,Qmax,Qmin,im,jm,km,km1,km2)
CMIC$* private(i,j)
#else
!$OMP PARALLEL DO PRIVATE( I, J )
#endif
#endif
DO 2000 j=1,jm
DO 2000 i=1,im
C k=1 and k=km
Pmax(i,j, 1) = max(Qmax(i,j, 2),Qmax(i,j, 1))
Pmin(i,j, 1) = min(Qmin(i,j, 2),Qmin(i,j, 1))
Pmax(i,j,km) = max(Qmax(i,j,km1),Qmax(i,j,km))
Pmin(i,j,km) = min(Qmin(i,j,km1),Qmin(i,j,km))
C k=2 and k=km1
Pmax(i,j, 2) = max(Qmax(i,j, 3),Pmax(i,j, 1))
Pmin(i,j, 2) = min(Qmin(i,j, 3),Pmin(i,j, 1))
Pmax(i,j,km1) = max(Qmax(i,j,km2),Pmax(i,j,km))
Pmin(i,j,km1) = min(Qmin(i,j,km2),Pmin(i,j,km))
2000 continue
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* shared(Pmax,Pmin,Qmax,Qmin,im,jm,km,km1,km2)
CMIC$* private(i,j,k)
#else
!$OMP PARALLEL DO PRIVATE( I, J, K )
#endif
#endif
DO 3000 k=3,km2
DO 3000 j=1,jm
DO 3000 i=1,im
Pmax(i,j,k) = max(Qmax(i,j,k-1),Qmax(i,j,k),Qmax(i,j,k+1))
Pmin(i,j,k) = min(Qmin(i,j,k-1),Qmin(i,j,k),Qmin(i,j,k+1))
3000 continue
return
end subroutine hilo3D
!------------------------------------------------------------------------------
subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT) 5
C****6***0*********0*********0*********0*********0*********0**********72
C
C A6 = CURVATURE OF THE TEST PARABOLA
C AR = RIGHT EDGE VALUE OF THE TEST PARABOLA
C AL = LEFT EDGE VALUE OF THE TEST PARABOLA
C DC = 0.5 * MISMATCH
C P = CELL-AVERAGED VALUE
C IM = VECTOR LENGTH
C
C OPTIONS:
C
C LMT = 0: FULL MONOTONICITY
C LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS)
C LMT = 2: POSITIVE-DEFINITE CONSTRAINT
C
PARAMETER ( R12 = 1./12. )
REAL A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
if(LMT.eq.0) then
C Full constraint
do 100 i=1,IM
if(DC(i).eq.0.) then
AR(i) = p(i)
AL(i) = p(i)
A6(i) = 0.
else
da1 = AR(i) - AL(i)
da2 = da1**2
A6DA = A6(i)*da1
if(A6DA .lt. -da2) then
A6(i) = 3.*(AL(i)-p(i))
AR(i) = AL(i) - A6(i)
elseif(A6DA .gt. da2) then
A6(i) = 3.*(AR(i)-p(i))
AL(i) = AR(i) - A6(i)
endif
endif
100 continue
elseif(LMT.eq.1) then
C Semi-monotonic constraint
do 150 i=1,IM
if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 150
if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
AR(i) = p(i)
AL(i) = p(i)
A6(i) = 0.
elseif(AR(i) .gt. AL(i)) then
A6(i) = 3.*(AL(i)-p(i))
AR(i) = AL(i) - A6(i)
else
A6(i) = 3.*(AR(i)-p(i))
AL(i) = AR(i) - A6(i)
endif
150 continue
elseif(LMT.eq.2) then
do 250 i=1,IM
if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 250
fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12
if(fmin.ge.0.) go to 250
if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
AR(i) = p(i)
AL(i) = p(i)
A6(i) = 0.
elseif(AR(i) .gt. AL(i)) then
A6(i) = 3.*(AL(i)-p(i))
AR(i) = AL(i) - A6(i)
else
A6(i) = 3.*(AR(i)-p(i))
AL(i) = AR(i) - A6(i)
endif
250 continue
endif
return
end subroutine lmtppm
!------------------------------------------------------------------------------
SUBROUTINE qckxyz(Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp,IC,NSTEP) 2,4
C****6***0*********0*********0*********0*********0*********0**********72
! Added to pass C-preprocessor switches (bmy, 3/9/01)
# include "define.h"
PARAMETER ( tiny = 1.E-30 )
PARAMETER ( kmax = 200 )
REAL Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)
integer IP(kmax)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
NLM1 = NLAY-1
C Do horizontal filling.
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* private(i,j,L,qtmp)
#else
!$OMP PARALLEL DO PRIVATE( I, J, L, QTMP )
#endif
#endif
do 1000 L=1,NLAY
call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ip(L),tiny)
if(ip(L).ne.0)
& call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ip(L),tiny)
1000 continue
ipz = 0
do L=1,NLAY
if(ip(L) .ne. 0) then
ipz = L
go to 111
endif
enddo
return
111 continue
if(ipz .eq. 0) return
if(ipz .eq. 1) then
lpz = 2
else
lpz = ipz
endif
C Do vertical filling.
#if defined( multitask )
#if defined( CRAY )
CMIC$ do all autoscope
CMIC$* private(i,j,L,qup,qly,dup)
#else
!$OMP PARALLEL DO PRIVATE( I, J, L, QUP, QLY, DUP )
#endif
#endif
do 2000 j=j1,j2
if(ipz .eq. 1) then
C Top layer
do i=1,IMR
if(Q(i,j,1).LT.0.) then
Q(i,j,2) = Q(i,j,2) + Q(i,j,1)
Q(i,j,1) = 0.
endif
enddo
endif
DO 225 L = lpz,NLM1
do i=1,IMR
IF( Q(i,j,L).LT.0.) THEN
C From above
qup = Q(i,j,L-1)
qly = -Q(i,j,L)
dup = min(qly,qup)
Q(i,j,L-1) = qup - dup
Q(i,j,L ) = dup-qly
C Below
Q(i,j,L+1) = Q(i,j,L+1) + Q(i,j,L)
Q(i,j,L) = 0.
ENDIF
enddo
225 CONTINUE
C BOTTOM LAYER
L = NLAY
do i=1,IMR
IF( Q(i,j,L).LT.0.) THEN
C From above
qup = Q(i,j,NLM1)
qly = -Q(i,j,L)
dup = min(qly,qup)
Q(i,j,NLM1) = qup - dup
C From "below" the surface.
Q(i,j,L) = 0.
ENDIF
enddo
2000 continue
RETURN
END SUBROUTINE qckxyz
!------------------------------------------------------------------------------
subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx) 2
C****6***0*********0*********0*********0*********0*********0**********72
REAL p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
JMR = JNP-1
do 1309 j=j1,j2
if(J.GT.JS .and. J.LT.JN) GO TO 1309
C
do i=1,IMR
qtmp(i) = p(i,j)
enddo
C
do i=-IML,0
qtmp(i) = p(IMR+i,j)
qtmp(IMR+1-i) = p(1-i,j)
enddo
C
DO i=1,IMR
iu = UA(i,j)
ru = UA(i,j) - iu
iiu = i-iu
if(UA(i,j).GE.0.) then
adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
else
adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu)-qtmp(iiu+1))
endif
enddo
do i=1,IMR
adx(i,j) = adx(i,j) - p(i,j)
enddo
1309 continue
C Eulerian upwind
do j=JS+1,JN-1
C
do i=1,IMR
qtmp(i) = p(i,j)
enddo
C
qtmp(0) = p(IMR,J)
qtmp(IMR+1) = p(1,J)
C
DO i=1,IMR
IP = i - UA(i,j)
adx(i,j) = UA(i,j)*(qtmp(ip)-qtmp(ip+1))
enddo
enddo
C
if(j1.ne.2) then
do i=1,IMR
adx(i, 2) = 0.
adx(i,JMR) = 0.
enddo
endif
C set cross term due to x-adv at the poles to zero.
do i=1,IMR
adx(i, 1) = 0.
adx(i,JNP) = 0.
enddo
return
end subroutine xadv
!------------------------------------------------------------------------------
subroutine xmist(IMR,IML,P,DC) 6
C****6***0*********0*********0*********0*********0*********0**********72
REAL P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
C 2nd order version.
C
do 10 i=1,IMR
tmp = 0.25*(p(i+1) - p(i-1))
Pmax = max(P(i-1), p(i), p(i+1)) - p(i)
Pmin = p(i) - min(P(i-1), p(i), p(i+1))
10 DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
return
end subroutine xmist
!------------------------------------------------------------------------------
subroutine xtp(im,jm,IML,j1,j2,JN,JS,PU,DQ,q,C,fx2,xmass,IORD, 4,10
& fx1_tp)
C****6***0*********0*********0*********0*********0*********0**********72
REAL C(im,*),DC(-IML:im+IML+1),xmass(im,jm),
& fx1(im+1),DQ(im,jm),qtmp(-IML:im+1+IML)
REAL PU(im,jm),q(im,jm)
real fx2(im+1,jm)
INTEGER isave(im)
! bey, 6/20/00. for mass-flux diagnostic
real fx1_tp(im,jm)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
IMP = im + 1
C
C van Leer at high latitudes
jvan = max(1,jm/20)
j1vl = j1+jvan
j2vl = j2-jvan
C
do 1310 j=j1,j2
C
do i=1,im
qtmp(i) = q(i,j)
enddo
C
if(j.ge.JN .or. j.le.JS) goto 2222
C****6***0*********0*********0*********0*********0*********0**********72
C *** Eulerian ***
C****6***0*********0*********0*********0*********0*********0**********72
C
qtmp(0) = q(im,J)
qtmp(-1) = q(im-1,J)
qtmp(IMP) = q(1,J)
qtmp(IMP+1) = q(2,J)
IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
do i=1,im
iu = float(i) - c(i,j)
fx1(i) = qtmp(iu)
enddo
C Zero high order contribution
DO i=1,im
fx2(i,j) = 0.
enddo
ELSE
call xmist(im,IML,Qtmp,DC)
DC(0) = DC(im)
C
if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then
DO i=1,im
iu = float(i) - c(i,j)
fx1(i ) = qtmp(iu)
fx2(i,j) = DC(iu)*(sign(1.,c(i,j))-c(i,j))
enddo
else
call fxppm(im,IML,C(1,j),Qtmp,DC,fx1,fx2(1,j),IORD)
endif
C
ENDIF
C
DO i=1,im
fx1(i ) = fx1(i )*xmass(i,j)
fx2(i,j) = fx2(i,j)*xmass(i,j)
enddo
C
goto 1309
C
C****6***0*********0*********0*********0*********0*********0**********72
C *** Conservative (flux-form) Semi-Lagrangian transport ***
C****6***0*********0*********0*********0*********0*********0**********72
2222 continue
C ghost zone for the western edge:
iuw = -c(1,j)
iuw = min(0, iuw)
do i=iuw, 0
qtmp(i) = q(im+i,j)
enddo
C ghost zone for the eastern edge:
iue = imp - c(im,j)
iue = max(imp, iue)
do i=imp, iue
qtmp(i) = q(i-im,j)
enddo
if(iord.eq.1 .or. j.eq.j1. or. j.eq.j2) then
do i=1,im
iu = c(i,j)
if(c(i,j) .le. 0.) then
itmp = i - iu
isave(i) = itmp - 1
else
itmp = i - iu - 1
isave(i) = itmp + 1
endif
fx1(i) = (c(i,j)-iu) * qtmp(itmp)
enddo
C Zero high order contribution
do i=1,im
fx2(i,j) = 0.
enddo
ELSE
call xmist(im,IML,qtmp,dc)
do i=iuw, 0
dc(i) = dc(im+i)
enddo
do i=imp, iue
dc(i) = dc(i-im)
enddo
do i=1,im
iu = c(i,j)
rut = c(i,j) - iu
if(c(i,j) .le. 0.) then
itmp = i - iu
isave(i) = itmp - 1
fx2(i,j) = -rut*dc(itmp)*(1.+rut)
else
itmp = i - iu - 1
isave(i) = itmp + 1
fx2(i,j) = rut*dc(itmp)*(1.-rut)
endif
fx1(i) = rut*qtmp(itmp)
enddo
ENDIF
do i=1,im
IF(c(i,j).GT.1.) then
CDIR$ NOVECTOR
do ist = isave(i),i-1
fx1(i) = fx1(i) + qtmp(ist)
enddo
elseIF(c(i,j).LT.-1.) then
CDIR$ NOVECTOR
do ist = i,isave(i)
fx1(i) = fx1(i) - qtmp(ist)
enddo
endif
enddo
CDIR$ VECTOR
do i=1,im
fx1(i) = PU(i,j)*fx1(i)
fx2(i,j) = PU(i,j)*fx2(i,j)
enddo
1309 fx1(IMP ) = fx1(1 )
fx2(IMP,j) = fx2(1,j)
C Update using low order fluxes.
DO i=1,im
DQ(i,j) = DQ(i,j) + fx1(i)-fx1(i+1)
! bey, 6/20/00. for mass-flux diagnostic
fx1_tp(i,j) = fx1(i)
enddo
1310 continue
return
end subroutine xtp
!------------------------------------------------------------------------------
subroutine ymist(IMR,JNP,j1,P,DC) 3
C****6***0*********0*********0*********0*********0*********0**********72
PARAMETER ( R24 = 1./24. )
REAL P(IMR,JNP),DC(IMR,JNP)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C
C 2nd order version for scalars
C
IMH = IMR / 2
JMR = JNP - 1
C
do 10 i=1,IMR*(JMR-1)
tmp = 0.25*(p(i,3) - p(i,1))
Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)
10 CONTINUE
C
C Poles:
C
if(j1.ne.2) then
do i=1,IMR
DC(i,1) = 0.
DC(i,JNP) = 0.
enddo
else
C Determine slopes in polar caps for scalars!
C
do 20 i=1,IMH
C South
tmp = 0.25*(p(i,2) - p(i+imh,2))
Pmax = max(p(i,2),p(i,1), p(i+imh,2)) - p(i,1)
Pmin = p(i,1) - min(p(i,2),p(i,1), p(i+imh,2))
DC(i,1)=sign(min(abs(tmp),Pmax,Pmin),tmp)
C North.
tmp = 0.25*(p(i+imh,JMR) - p(i,JMR))
Pmax = max(p(i+imh,JMR),p(i,jnp), p(i,JMR)) - p(i,JNP)
Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR))
DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp)
20 continue
C
C Scalars:
do 25 i=imh+1,IMR
DC(i, 1) = - DC(i-imh, 1)
DC(i,JNP) = - DC(i-imh,JNP)
25 continue
endif
return
end subroutine ymist
!------------------------------------------------------------------------------
subroutine ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,C,DC2 3,6
& ,ymass,fy1,A6,AR,AL,fy2,JORD,
& fy1_tp)
C****6***0*********0*********0*********0*********0*********0**********72
REAL P(IMR,JNP),C(IMR,JNP),ymass(IMR,JNP),fy2(IMR,JNP),
& DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP)
! bey, 6/20/00. for mass-flux diagnostic
REAL fy1_tp(IMR,JNP)
C============================================================================
C Cray NOBOUNDS directive will turn off all subscript bounds checking.
C This directive is also legal on SGI compilers (bmy, 4/24/00)
CDIR$ NOBOUNDS
C============================================================================
C Work array
REAL fy1(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
C
JMR = JNP - 1
len = IMR*(J2-J1+2)
if(JORD.eq.1) then
DO 1000 i=1,len
JT = float(J1) - C(i,J1)
1000 fy1(i,j1) = p(i,JT)
DO 1050 i=1,len
1050 fy2(i,j1) = 0.
else
call ymist(IMR,JNP,j1,P,DC2)
C
if(JORD.LE.0 .or. JORD.GE.3) then
call fyppm(C,P,DC2,fy1,fy2,IMR,JNP,j1,j2,A6,AR,AL,JORD)
else
DO 1200 i=1,len
JT = float(J1) - C(i,J1)
fy1(i,j1) = p(i,JT)
1200 fy2(i,j1) = (sign(1.,C(i,j1))-C(i,j1))*DC2(i,JT)
endif
endif
C
DO 1300 i=1,len
fy1(i,j1) = fy1(i,j1)*ymass(i,j1)
1300 fy2(i,j1) = fy2(i,j1)*ymass(i,j1)
C
!=============================================================================
! This loop had to be extended for the mass-flux diagnostics (bmy, 4/26/00)
! DO 1400 j=j1,j2
! DO 1400 i=1,IMR
!1400 DQ(i,j) = DQ(i,j) + (fy1(i,j) - fy1(i,j+1)) * acosp(j)
!=============================================================================
DO j=j1,j2
DO i=1,IMR
DQ(i,j) = DQ(i,j) + (fy1(i,j) - fy1(i,j+1)) * acosp(j)
! bey, 6/20/00. for mass-flux diagnostic
fy1_tp(i,j) = fy1(i,j)
ENDDO
ENDDO
C
C Poles
sum1 = fy1(IMR,j1 )
sum2 = fy1(IMR,J2+1)
do i=1,IMR-1
sum1 = sum1 + fy1(i,j1 )
sum2 = sum2 + fy1(i,J2+1)
enddo
C
sum1 = DQ(1, 1) - sum1 * RCAP
sum2 = DQ(1,JNP) + sum2 * RCAP
do i=1,IMR
DQ(i, 1) = sum1
DQ(i,JNP) = sum2
enddo
C
if(j1.ne.2) then
do i=1,IMR
DQ(i, 2) = sum1
DQ(i,JMR) = sum2
enddo
endif
return
end subroutine ytp
!------------------------------------------------------------------------------
SUBROUTINE PRESS_FIX( FX, FY, NDT, ACOSP, J1 ) 2,10
!
!******************************************************************************
! Subroutine PRESS_FIX is a wrapper for the Pressure fixer DYN0. PRESS_FIX
! takes the mass fluxes in pressure units and converts them to [kg air/s]
! using the correct geometry for TPCORE. (bdf, bmy, 10/11/01, 2/4/03)
!
! Arguments as Input:
! ============================================================================
! (1 ) FX (REAL*8 ) : E-W flux passed from TPCORE [mb/timestep]
! (2 ) FY (REAL*8 ) : N-S flux passed from TPCORE [mb/timestep]
! (3 ) NDT (INTEGER) : Dynamic timestep for TPCORE [s]
! (4 ) ACOSP (REAL*8 ) : Array of inverse cosines [unitless]
! (5 ) J1 (INTEGER) : TPCORE polar cap extent [# of boxes]
!
! NOTES:
! (1 ) Adapted from original code from LLNL. Added comments and F90 syntax
! for declarations. (bdf, bmy, 10/11/01)
! (2 ) For now, assumes that JGLOB=JJPAR, and DXYP(J) is equivalent to
! DXYP(J+J0). (bmy, 10/11/01)
! (3 ) Now declare DXYP as a local array, and initialize it with calls
! to routine GET_AREA_M2 of "grid_mod.f". Now use function GET_TS_DYN
! from "time_mod.f". Remove reference to CMN header file. (bmy, 2/4/03)
!******************************************************************************
!
! References to F90 modules
USE GRID_MOD, ONLY : GET_AREA_M2
USE TIME_MOD, ONLY : GET_TS_DYN
IMPLICIT NONE
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! Diagnostic switches
# include "CMN_GCTM" ! g0_100
! Arguments
INTEGER, INTENT(IN) :: NDT, J1
REAL*8, INTENT(IN) :: ACOSP(JJPAR)
REAL*8, INTENT(INOUT) :: FX(IIPAR,JJPAR,LLPAR)
REAL*8, INTENT(INOUT) :: FY(IIPAR,JJPAR,LLPAR)
! Local variables
INTEGER :: I, J, J2, K, K2, L
REAL*8 :: DTC, DTDYN, NSDYN, SUM1, SUM2
REAL*8 :: NP_FLUX(IIPAR,LLPAR)
REAL*8 :: SP_FLUX(IIPAR,LLPAR)
REAL*8 :: ALFA(IIPAR+1,JJPAR,LLPAR)
REAL*8 :: BETA(IIPAR,JJPAR+1,LLPAR)
REAL*8 :: GAMA(IIPAR,JJPAR,LLPAR+1)
REAL*8 :: UMFLX(IIPAR,JJPAR,LLPAR)
REAL*8 :: VMFLX(IIPAR,JJPAR,LLPAR)
! Local SAVEd variables
LOGICAL, SAVE :: FIRST = .TRUE.
REAL*8, SAVE :: DXYP(JJPAR)
!=================================================================
! PRESS_FIX begins here!
!
! K is the vertical index down from the atmosphere top downwards
! K2 is the vertical index up from the surface
!=================================================================
! Initialize arrays
ALFA = 0d0
BETA = 0d0
GAMA = 0d0
! NSDYN is the dynamic time step in seconds
NSDYN = GET_TS_DYN() * 60d0
! J2 is the south polar edge
J2 = JJPAR - J1 + 1
! DTDYN = double precision value for NDT, the dynamic timestep
DTDYN = DBLE( NDT )
! Save grid box surface areas [m2] into the local DXYP array
IF ( FIRST ) THEN
DO J = 1, JJPAR
DXYP(J) = GET_AREA_M2( J )
ENDDO
! Reset first-time flag
FIRST = .FALSE.
ENDIF
!=================================================================
! FX is the E-W mass flux from TPCORE in [mb/timestep].
! UMFLX is the mass flux in [kg air/s], which is what DYN0 needs.
!
! FY is the E-W mass flux from TPCORE in [mb/timestep].
! VMFLX is the mass flux in [kg air/s], which is what DYN0 needs.
!
! The unit conversion from [mb/timestep] to [kg air/s] is:
!
! mb | 100 Pa | 1 kg air | s^2 | step | DXYP m^2 kg air
! ------+--------+----------+-------+--------+---------- = -------
! step | mb | Pa m s^2 | 9.8 m | DTDYN s| s s
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K, K2, DTC )
DO K = 1, LLPAR
K2 = LLPAR - K + 1
! Compute UMFLX from FX
DO J = 1, JJPAR
DO I = 1, IIPAR
UMFLX(I,J,K2) = FX(I,J,K) * ( G0_100 * DXYP(J) ) / DTDYN
ENDDO
ENDDO
! Compute VMFLX from FY
DO I = 1, IIPAR
DO J = J1, J2+1
IF ( FY(I,J,K) .GE. 0 ) THEN
DTC = FY(I,J,K) * G0_100 * ACOSP(J) * DXYP(J) / DTDYN
ELSE
DTC = FY(I,J,K) * G0_100 * ACOSP(J-1)* DXYP(J-1) / DTDYN
ENDIF
VMFLX(I,J,K2) = DTC
ENDDO
ENDDO
!=================================================================
! TREATMENT OF THE POLES: 1
! copy ymass values strait into vmflx at poles for pressure fixer
!=================================================================
DO I = 1, IIPAR
VMFLX(I,1,K2) = FY(I,1,K)
VMFLX(I,J1-1,K2) = FY(I,J1-1,K)
VMFLX(I,JJPAR,K2) = FY(I,JJPAR,K)
ENDDO
ENDDO
!$OMP END PARALLEL DO
DO K = 1, LLPAR
!=================================================================
! TREATMENT OF THE POLES: 2
! North polar cap: J=1
!=================================================================
SUM1 = FY(IIPAR,J1,K)
DO I = 1, IIPAR-1
SUM1 = SUM1 + FY(I,J1,K)
ENDDO
! NORTH POLE FLUX IN KG.
DO I = 1, IIPAR
NP_FLUX(I,K) = SUM1 * G0_100 * ACOSP(1) * DXYP(1)
ENDDO
!==============================================================
! TREATMENT OF THE POLES: 3
! South polar cap: J=JJPAR
!==============================================================
SUM2 = FY(IIPAR,J2+1,K)
DO I = 1, IIPAR-1
SUM2 = SUM2 + FY(I,J2+1,K)
ENDDO
DO I = 1, IIPAR
SP_FLUX(I,K) = SUM2 * G0_100 * ACOSP(JJPAR) * DXYP(JJPAR)
ENDDO
ENDDO
!=================================================================
! Call DYN0 to fix the pressures
!=================================================================
CALL DYN0( NSDYN, J1, NP_FLUX, SP_FLUX,
& UMFLX, VMFLX, ALFA, BETA, GAMA )
!=================================================================
! ALFA is the E-W mass flux adjusted by DYN0 in [kg air/s]
! FX is the E-W mass flux for TPCORE in [mb/timestep].
!
! BETA is the N-S mass flux adjusted by DYN0 in [kg air/s]
! FY is the E-W mass flux for TPCORE in [mb/timestep].
!
! The unit conversion from to [kg air/s] to [mb/timestep] is:
!
! kg air | Pa m s^2 | 9.8 m | 1 | DTDYN s | mb mb
! --------+----------+-------+----------+---------+------- = ----
! s | 1 kg air | s^2 | DXYP m^2 | step | 100 Pa step
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K, K2 )
DO K = 1, LLPAR
K2 = LLPAR - K + 1
! Update FX from ALFA
DO J = 1, JJPAR
DO I = 1, IIPAR
FX(I,J,K) = ALFA(I,J,K2) * DTDYN / ( G0_100 * DXYP(J) )
ENDDO
ENDDO
! Update FY from BETA
DO I = 1, IIPAR
DO J = J1, J2+1
IF ( BETA(I,J,K) .GE. 0 ) THEN
FY(I,J,K) = BETA(I,J,K2) * DTDYN /
& ( G0_100 * ACOSP(J) * DXYP(J) )
ELSE
FY(I,J,K) = BETA(I,J,K2) * DTDYN /
& ( G0_100 * ACOSP(J-1) * DXYP(J-1) )
ENDIF
ENDDO
ENDDO
! Special treatment of BETA at the poles
DO I = 1, IIPAR
FY(I,1,K) = BETA(I,1,K2)
FY(I,J1-1,K) = BETA(I,J1-1,K2)
FY(I,JJPAR,K) = BETA(I,JJPAR,K2)
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE PRESS_FIX
!------------------------------------------------------------------------------
SUBROUTINE DYN0( DTWIND, J1, NP_FLUX, SP_FLUX, 2,16
& UMFLX, VMFLX, ALFA, BETA, GAMA )
!
!******************************************************************************
! Subroutine DYN0 is the pressure fixer for TPCORE. DYN0 readjusts the
! mass fluxes ALFA, BETA, GAMA, so that they are consistent with the
! met fields. (bdf, bmy, 10/11/01, 7/20/04)
!
! Arguments as Input:
! ============================================================================
! (1 ) DTWIND (REAL*8 ) : Time step between wind intervals [s]
! (2 ) J1 (INTEGER) : TPCORE polar cap width
! (4 ) NP_FLUX (REAL*8 ) : North polar flux (from PRESS_FIX) in [kg]
! (5 ) SP_FLUX (REAL*8 ) : South polar flux (from PRESS_FIX) in [kg]
! (6 ) UMFLX (REAL*8 ) : Wet air mass flux in E-W direction [kg air/s]
! (7 ) VMFLX (REAL*8 ) : Wet air mass flux in N-S direction [kg air/s]
! (8 ) ALFA (REAL*8 ) : Dry air mass flux in E-W direction [kg air/s]
! (9 ) BETA (REAL*8 ) : Dry air mass flux in N-S direction [kg air/s]
! (10) GAMA (REAL*8 ) : Dry air mass flux in up/down direction [kg air/s]
!
! Arguments as Output:
! ============================================================================
! (8 ) ALFA (REAL*8 ) : ALFA air mass, after pressure fix is applied
! (9 ) BETA (REAL*8 ) : BETA air mass, after pressure fix is applied
! (10) GAMA (REAL*8 ) : GAMA air mass, after pressure fix is applied
!
! NOTES:
! (1 ) Adapted from original code from LLNL. Added comments and F90 syntax
! for declarations. (bdf, bmy, 10/10/01)
! (2 ) For a global run (as we usually do in GEOS-CHEM) IM=ID=IIPAR and
! JM=JD=JJPAR. (bmy, 10/10/01)
! (3 ) For now, assumes that JGLOB=JJPAR, and DXYP(J) is equivalent to
! DXYP(J+J0). (bmy, 10/11/01)
! (4 ) Rename AD to AD_L so as not to conflict with the AD array in
! the header file "CMN" (bmy, 10/11/01)
! (5 ) Now reference PSC2 instead of PS from "dao_mod.f". Replace all
! instances of PS with PSC2. Updated comments. (bdf, bmy, 4/1/02)
! (6 ) Removed obsolete code from 4/1/02 (bdf, bmy, 8/22/02)
! (7 ) Now declare DXYP as a local array, and initialize it with calls
! to routine GET_AREA_M2 and GET_YOFFSET of "grid_mod.f". Now also
! references GET_BP from "pressure_mod.f" (bmy, 2/11/03)
! (8 ) Removed reference to CMN (bmy, 7/20/04)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : SPHU, PSC2, AIRDEN, AIRVOL
USE GRID_MOD, ONLY : GET_AREA_M2, GET_YOFFSET
USE PRESSURE_MOD, ONLY : GET_BP
IMPLICIT NONE
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: J1
REAL*8, INTENT(IN) :: DTWIND
REAL*8, INTENT(IN) :: NP_FLUX(IIPAR,LLPAR)
REAL*8, INTENT(IN) :: SP_FLUX(IIPAR,LLPAR)
REAL*8, INTENT(IN) :: UMFLX(IIPAR,JJPAR,LLPAR)
REAL*8, INTENT(IN) :: VMFLX(IIPAR,JJPAR,LLPAR)
REAL*8, INTENT(INOUT) :: ALFA(IIPAR+1,JJPAR,LLPAR)
REAL*8, INTENT(INOUT) :: BETA(IIPAR,JJPAR+1,LLPAR)
REAL*8, INTENT(INOUT) :: GAMA(IIPAR,JJPAR,LLPAR+1)
! Local variables
LOGICAL :: LSP, LNP, LEW
INTEGER :: IIX, JJX, KM, JB, JE, IEPZ, IMZ
INTEGER :: I, J, J2, K, L
REAL*8 :: ALFAX, UFILT, VFILT, PCTM8, G100
REAL*8 :: AIRQAV, AWE, SUMAD0, SUMAW0, AIRWET
REAL*8 :: AIRH2O, AIRQKG ,SUM1, SUMA, SUMP, SUMQ
REAL*8 :: SUMU, SUMV, SUMW, ZIMZ, ZDTW, G0
REAL*8 :: AD_L(IIPAR,JJPAR,LLPAR)
REAL*8 :: AIRD(IIPAR,JJPAR,LLPAR)
REAL*8 :: AIRNEW(IIPAR,JJPAR,LLPAR)
REAL*8 :: AIRX(IIPAR,JJPAR,LLPAR)
REAL*8 :: AX(IIPAR+1,JJPAR)
REAL*8 :: BX(IIPAR,JJPAR+1)
REAL*8 :: MERR(IIPAR,JJPAR)
REAL*8 :: PCTM(IIPAR,JJPAR)
REAL*8 :: PERR(IIPAR,JJPAR)
REAL*8 :: SPHU_KG(IIPAR,JJPAR,LLPAR)
REAL*8 :: SUMAQ(IIPAR,JJPAR)
REAL*8 :: XYB(IIPAR,JJPAR)
REAL*8 :: XYZB(IIPAR,JJPAR,LLPAR)
! Local saved variables
LOGICAL, SAVE :: FIRST = .TRUE.
REAL*8, SAVE :: DXYP(JJPAR)
REAL*8, SAVE :: DSIG(LLPAR)
!=================================================================
! DYN0 begins here!
!
! UNITS OF AIR MASS AND TRACER = (kg)
!
! Air mass (kg) is given by:
! area (m^2) * pressure thickness (Pa) / g0
!
! DXYP(J) = area of [I,J] [m^2] (is longitude-symmetric)
!
! PSC2(I,J) = surf pressure [Pa] averaged in extended zone.
! this is the surface pressure at the end of the
! current dynamic timestep, passed to TPCORE.
!
! SPHU_KG(I,J,K) = specific humidity of grid box
! [kg H2O/kg wet air] averaged in extended zone.
!
! AIRQKG(I,J) = Mass of H2O [kg] at each level
! = PS(I,J)) * SPHU_KG(I,J,K)
!
! AIRD(I,J,K) = dry-air mass [kg] in each box as calculated
! in CTM at the beginning of each time step,
! updated at end of DYN0.
!
! PCTM(I,J) = inferred wet-air (total) surf press [Pa] calc.
! in CTM (using SUMAQ & AIRD-X-NEW)
!
! AIRNEW(I,J,K) = new dry-air mass in each CTM box after
! horizontal divergence (ALFA+BETA) over time
! step DTWIND (sec)
!
! AIRX(I,J,K) = expected dry-air mass in each CTM box after
! calculating the vertical divergence (GAMA)
! (also used for GCM dry mass)
! = XYZA(I,J,K) + XYZB(I,J,K)*PCTM(I,J) - AIRQKG
!
! DTWIND = time step [s] that applies to the averaged
! wind fields (i.e., the time between successive
! pressures.
!
!-----------------------------------------------------------------
!
! Assume that we have "wet-air" mass fluxes across each boundary
!
! UMFLX(I,J,K) ==> [I,J,K] ==> UMFLX(I+1,J,K) [kg air/s]
! VMFLX(I,J,K) ==> [I,J,K] ==> VMFLX(I,J+1,K) [kg air/s]
!
! Convert to "dry-air" mass flux in/out of box using
! average Q at boundary
!
! ALFA(I,J,K) ==> [I,J,K] ==> ALFA(I+1,J,K) [kg air/s]
! BETA(I,J,K) ==> [I,J,K] ==> BETA(I,J+1,K) [kg air/s]
!
! Calculate convergence in each layer of dry air, compare with
! expected dry air mass (AIRX) and then calculate vertical
! dry-mass fluxes
!
! GAMA(I,J,K) ==> [I,J,K] ==> GAMA(I,J,K+1) [kg air/s]
!
! Horizontal pressure filter adjusts UMFLX & VMFLX to reduce
! error in [PCTM - PSC2]
!
! UMFLX + pressure filter ==> UMFLX#,
! VMFLX + filter ==> VMFLX# (temporary)
!
! The pressure filter does nearest neighbor flux
! (adjusting ALFA/BETA)
!
!-----------------------------------------------------------------
!
! Note that K->K+1 is downward (increasing pressure) and
! that boundaries:
! GAMA(I,J,1) = GAMA(I,J,KM+1) = 0 no flux across
! upper/lower boundaries
!
! BETA(I,1,K) = BETA(I,JJPAR+1,K) = 0 no flux at S & N poles
!
! ALFA(1,J,K) = ALFA(IIPAR+1,J,K) is NOT ZERO, but cyclic
!
! Dimensions for ALFA, BETA, GAMA are extended by +1 beyond grid
! to allow simple formulation of fluxes in/out of final grid box.
!
! GCM input UMFLX,VMFLX,PSG is ALWAYS of GLOBAL dimensions
! (IIPAR x JJPAR x LLPAR)
!
! Indices of ALFA, BETA, GAMA, SPHU_KG & PSC2 are always LOCAL
! (IIPAR x JJPAR x KM): FOR GEOS-CHEM, KM = LLPAR (bmy
!
! Indices of tracer (STT), and diagnostics are local
! (w.r.t. WINDOW. WINDOW calculations are defined by an
! offset and size
!
! I0 .ge.0 and IIPAR+I0 .le. IIPAR
! J0 .ge.0 and JJPAR+J0 .le. JJPAR
! K0 .ge.0 and KM+K0 .le. LLPAR
!
! The WINDOW calculation must allow for a boundary layer
! of grid boxes:
!
! IG(abs. coords) = IW(in window) + I0
! JG(abs. coords) = JW(in window) + J0
! KG(abs. coords) = KW(in window) + K0
!
! vertical window (NEW) allows for an upper boundary with flow
! across it and specified mixing ratio b.c.'s at KG = K0
!=================================================================
! First-time initialization
IF ( FIRST ) THEN
! Surface area [m2]
DO J = 1, JJPAR
DXYP(J) = GET_AREA_M2( J )
ENDDO
! Sigma-level thickness [unitless]
! Assumes we are using a pure-sigma grid
DO L = 1, LLPAR
DSIG(L) = GET_BP(L) - GET_BP(L+1)
ENDDO
! Reset first-time flag
FIRST = .FALSE.
ENDIF
! for tpcore poles.
J2 = JJPAR - J1 + 1
! geos code
G0 = 9.8d0
!=================================================================
! XYZB is the factor needed to get mass in kg of gridbox
! mass (kg) = XYZB (kg/mb) * P (mb)
!
! AD_L is the dry air mass in the grid box
!
! SPHU_KG is the water vapor [kg H2O/kg air]
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO L = 1, LLPAR
DO J = 1, JJPAR
DO I = 1, IIPAR
XYZB(I,J,L) = DSIG(L) * DXYP(J) * 1.d2 / G0
AD_L(I,J,L) = AIRDEN(L,I,J) * AIRVOL(I,J,L)
SPHU_KG(I,J,L) = SPHU(I,J,L) / 1000d0
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! XYB is the factor needed to get mass in kg of column
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J )
DO J = 1, JJPAR
DO I = 1, IIPAR
XYB(I,J) = SUM( XYZB(I,J,1:LLPAR) )
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! Define other variables
!=================================================================
G100 = 100.D0 / G0
ZDTW = 1.D0 / DTWIND
LSP = ( GET_YOFFSET() .EQ. 0 )
LNP = ( JJPAR + GET_YOFFSET() .EQ. JJPAR )
LEW = ( IIPAR .EQ. IIPAR )
!=================================================================
! Initialize ALFA with UMFLX and BETA with VMFLX
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO L = 1, LLPAR
DO J = 1, JJPAR
DO I = 1,IIPAR
ALFA(I,J,L) = UMFLX(I,J,L)
ENDDO
ALFA(IIPAR+1,J,L) = ALFA(1,J,L)
ENDDO
DO J = 2, JJPAR
DO I = 1, IIPAR
BETA(I,J,L) = VMFLX(I,J,L)
ENDDO
DO I = 1, IIPAR
BETA(I,1,L) = 0.D0
BETA(I,JJPAR+1,L) = 0.D0
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! SUMAQ(I,J): column integral of water (kg)
! Check on air mass
!=================================================================
SUMAD0 = 0.D0
SUMAW0 = 0.D0
DO J = 1, JJPAR
DO I = 1, IIPAR
SUMAQ(I,J) = 0.D0
DO K = 1, LLPAR
AIRWET = PSC2(I,J) * XYZB(I,J,K)
AIRH2O = SPHU_KG(I,J,K) * AIRWET
SUMAQ(I,J) = SUMAQ(I,J) + AIRH2O
SUMAD0 = SUMAD0 + AIRWET
SUMAW0 = SUMAW0 + AIRH2O
ENDDO
ENDDO
ENDDO
SUMAD0 = SUMAD0 - SUMAW0
!=================================================================
! Initialize AIRD, the dry-air mass [kg] in each box as calculated
! in CTM at the start of each time step, updated at end of DYN0.
!
! Compute AIRNEW, the new dry-air mass in each CTM box after
! horizontal divergence (ALFA+BETA) over time step DTWIND (sec)
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K )
DO K = 1, LLPAR
DO J = J1, J2
DO I = 1, IIPAR
AIRD(I,J,K) = AD_L(I,J,K)
AIRNEW(I,J,K) = AIRD(I,J,K) + DTWIND *
& ( ALFA(I,J,K) - ALFA(I+1,J,K) +
& BETA(I,J,K) - BETA(I,J+1,K) )
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! treatment of the poles for tpcore.
! j=2 and j=jjpar-1 don't have any airmass change.
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, K )
DO K = 1, LLPAR
! J=1
DO I = 1, IIPAR
AIRNEW(I,1,K) = AD_L(I,1,K) - NP_FLUX(I,K)
ENDDO
! J=JJPAR
DO I = 1, IIPAR
AIRNEW(I,JJPAR,K) = AD_L(I,JJPAR,K) + SP_FLUX(I,K)
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! Average AIRNEW at the South pole
!=================================================================
ZIMZ = 1.D0 / DBLE( IIPAR )
IF ( LSP ) THEN
JB = 2
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, K, SUMA )
DO K = 1, LLPAR
SUMA = SUM( AIRNEW(1:IIPAR,1,K) ) * ZIMZ
DO I = 1, IIPAR
AIRNEW(I,1,K) = SUMA
ENDDO
ENDDO
!$OMP END PARALLEL DO
ELSE
JB = 1
ENDIF
!=================================================================
! Average AIRNEW at the North pole
!=================================================================
IF ( LNP ) THEN
JE = JJPAR - 1
! poles, just average AIRNEW
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, K, SUMA )
DO K = 1, LLPAR
SUMA = SUM( AIRNEW(1:IIPAR,JJPAR,K) ) * ZIMZ
DO I = 1, IIPAR
AIRNEW(I,JJPAR,K) = SUMA
ENDDO
ENDDO
!$OMP END PARALLEL DO
ELSE
JE = JJPAR
ENDIF
!================================================================
! BEGIN FILTER of PRESSURE ERRORS
!
! Define the error in surface pressure PERR expected at end of
! time step filter by error in adjacent boxes, weight by areas,
! adjust ALFA & BETA
!
! PCTM(I,J) = new CTM wet-air column based on
! dry-air convergence (Pascals)
! PERR(I,J) = pressure-error between CTM-GCM at new time
! (before filter)
!================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J )
DO J = 1, JJPAR
DO I = 1, IIPAR
PCTM(I,J) = SUM( AIRNEW(I,J,:) ) / XYB(I,J)
! special case for j=2, jjpar-1 for tpcore pole configuration.
IF ( J .eq. 2 .OR. J .eq. JJPAR-1 ) THEN
PCTM(I,J) = PSC2(I,J)
ENDIF
PERR(I,J) = PCTM(I,J) - PSC2(I,J)
MERR(I,J) = PERR(I,J) * DXYP(J) * G100
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Call pressure filter
CALL PFILTR( MERR, AX, BX, DXYP, IIPAR, JJPAR,
& IIPAR, JJPAR, 1, LSP, LNP, LEW )
!=================================================================
! Calculate corrections to ALFA from the filtered AX
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, IIX, J, K, UFILT )
DO J = JB, JE
DO I = 1, IIPAR+1
IIX = MIN(I,IIPAR)
UFILT = AX(I,J) / ( XYB(IIX,J) * DTWIND )
DO K = 1, LLPAR
ALFA(I,J,K) = ALFA(I,J,K) + UFILT * XYZB(IIX,J,K)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! Calculate corrections to BETA from the filtered BX
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, JJX, K, VFILT )
DO J = 1, JJPAR+1
JJX = J
IF ( J+J .gt. JJPAR ) JJX = J - 1
DO I = 1, IIPAR
VFILT = BX(I,J) / ( XYB(I,JJX) * DTWIND )
DO K = 1, LLPAR
BETA(I,J,K) = BETA(I,J,K) + VFILT * XYZB(I,JJX,K)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! Calculate the corrected AIRNEW's & PCTM after P-filter:
! has changed ALFA+BETAs and ctm surface pressure (PCTM)
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K )
DO K = 1, LLPAR
DO J = 1, JJPAR
DO I = 1, IIPAR
AIRNEW(I,J,K) = AIRD(I,J,K) + DTWIND *
& ( ALFA(I,J,K) - ALFA(I+1,J,K) +
& BETA(I,J,K) - BETA(I,J+1,K) )
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! Average the adjusted AIRNEW at the South pole
!=================================================================
ZIMZ = 1.D0 / DBLE( IIPAR )
IF ( LSP ) THEN
JB = 2
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, K, SUMA )
DO K = 1, LLPAR
SUMA = SUM( AIRNEW(1:IIPAR,1,K ) ) * ZIMZ
DO I = 1, IIPAR
AIRNEW(I,1,K) = SUMA
ENDDO
ENDDO
!$OMP END PARALLEL DO
ELSE
JB = 1
ENDIF
!=================================================================
! Average the adjusted AIRNEW at the North pole
!=================================================================
IF ( LNP ) THEN
JE = JJPAR -1
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, K, SUMA )
DO K = 1, LLPAR
SUMA = SUM( AIRNEW(1:IIPAR,JJPAR,K) ) * ZIMZ
DO I = 1,IIPAR
AIRNEW(I,JJPAR,K) = SUMA
ENDDO
ENDDO
!$OMP END PARALLEL DO
ELSE
JE = JJPAR
ENDIF
!=================================================================
! END OF PRESSURE FILTER
!
! GAMA: redistribute the new dry-air mass consistent with the
! new CTM surface pressure, rigid upper b.c., no change in PCTM
!
! AIRX(I,J,K) = dry-air mass expected, based on PCTM
! PCTM(I,J) & PERR(I,J)
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K, PCTM8, AIRQKG )
DO J = 1, JJPAR
DO I = 1, IIPAR
PCTM8 = ( SUM( AIRNEW(I,J,:) ) + SUMAQ(I,J) ) / XYB(I,J)
PCTM(I,J) = PCTM8
PERR(I,J) = PCTM8 - PSC2(I,J)
DO K = 1, LLPAR
AIRQKG = SPHU_KG(I,J,K) * ( XYZB(I,J,K) * PSC2(I,J) )
AIRX(I,J,K) = PCTM8 * XYZB(I,J,K) - AIRQKG
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! GAMA from top down to be consistent with AIRX, AIRNEW not reset!
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K )
DO J = 1, JJPAR
DO I = 1, IIPAR
GAMA(I,J,LLPAR+1) = 0.D0
DO K = LLPAR, 2, -1
GAMA(I,J,K) = GAMA(I,J,K+1) - (AIRNEW(I,J,K) - AIRX(I,J,K))
ENDDO
! GAMA(I,J,1) will not be exactly ZERO, but it must be set so!
GAMA(I,J,1) = 0.D0
DO K = 2, LLPAR
GAMA(I,J,K) = GAMA(I,J,K) * ZDTW
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE DYN0
!------------------------------------------------------------------------------
SUBROUTINE PFILTR( MERR, ALFAX, BETAX, AXY, ID, JD, 2,4
& IM, JM, NITR, LSP, LNP, LEW )
!
!******************************************************************************
! Subroutine PFILTR applies the pressure-filter, the pressure
! between predicted Ps(CTM) and Ps(GCM). (bdf, bmy, 10/11/01)
!
! Arguments as Input:
! ============================================================================
! (1 ) MERR(ID,JD) (REAL*8 ) : mass error
! (2 ) ALFAX(ID+1,JD) (REAL*8 ) : perturbed ALFA by MERR
! (3 ) BETAX(ID,JD+1) (REAL*8 ) : perturbed BETA by MERR
! (4 ) AXY(ID,JD) (REAL*8 ) : area of grid box (I,J) in [m^2]
! (5-6) ID, JD (INTEGER) : "Global" array dimensions for lon, lat
! (7-8) IM, JM (INTEGER) : "Window" array dimensions for lon, lat
! (9 ) NITR (INTEGER) : number of iterations (NITR .LE. 4)
! (10 ) LSP (LOGICAL) : true if J=1 is S. POLE
! (11 ) LNP (LOGICAL) : true if J=JM is N. POLE
! (12 ) LEW (LOGICAL) : true if cyclic in W-E direction
! (i.e. if I=1 connects to I=IM)
!
! Arguments as Output:
! ============================================================================
! (1 ) MERR(ID,JD) (REAL*8 ) : adjusted mass error
! (2 ) ALFAX(ID+1,JD) (REAL*8 ) : adjusted ALFAX
! (3 ) BETAX(ID,JD+1) (REAL*8 ) : adjusted BETAX
!
! NOTES:
! (1 ) Adapted from original code from LLNL. Added comments and F90 syntax
! for declarations. (bdf, bmy, 10/1/01)
! (2 ) For a global run (as we usually do in GEOS-CHEM) IM=ID=IIPAR and
! JM=JD=JJPAR. (bmy, 10/11/01)
! (3 ) Removed IMEPZ -- we don't need this for GEOS-CHEM. (bmy, 10/18/01)
!******************************************************************************
!
IMPLICIT NONE
! Arguments
LOGICAL, INTENT(IN) :: LSP,LNP,LEW
INTEGER, INTENT(IN) :: ID, JD, IM, JM, NITR
REAL*8, INTENT(IN) :: AXY(JD)
REAL*8, INTENT(INOUT) :: MERR(ID,JD)
REAL*8, INTENT(INOUT) :: ALFAX(ID+1,JD)
REAL*8, INTENT(INOUT) :: BETAX(ID,JD+1)
! Local variables
LOGICAL :: LPOLE
INTEGER :: I, J, K
REAL*8 :: X0(ID,JD)
!=================================================================
! PFILTR begins here!
!=================================================================
! LPOLE is true if J=1 is the SOUTH POLE and J=JM is the NORTH POLE
! (this is the way GEOS-CHEM is set up, so LPOLE should be TRUE!)
LPOLE = ( LSP .AND. LNP )
! Zero ALFAX, BETAX, save MERR in X0
DO J = 1, JM
DO I = 1, IM
ALFAX(I,J) = 0.D0
BETAX(I,J) = 0.D0
X0(I,J) = MERR(I,J)
ENDDO
ALFAX(IM+1,J) = 0.D0
ENDDO
DO I = 1, IM
BETAX(I,JM+1) = 0.D0
ENDDO
!=================================================================
! Call LOCFLT to do the local filtering
!=================================================================
CALL LOCFLT( MERR, ALFAX, BETAX, AXY, ID, JD,
& IM, JM, 5, LSP, LNP, LEW )
!=================================================================
! Call POLFLT to do the pole filtering (if necessary)
!=================================================================
IF ( LPOLE ) THEN
CALL POLFLT( MERR, BETAX, AXY, 1.D0, ID, JD, IM, JM )
ENDIF
!=================================================================
! Compute mass error MERR and return
! MERR, ALFAX, and BETAX are now adjusted
!=================================================================
DO J = 1, JM
DO I = 1, IM
MERR(I,J) = X0(I,J) + ALFAX(I,J) - ALFAX(I+1,J)
& + BETAX(I,J) - BETAX(I,J+1)
ENDDO
ENDDO
! Return to calling program
END SUBROUTINE PFILTR
!------------------------------------------------------------------------------
SUBROUTINE LOCFLT( XERR, AX, BX, AXY, ID, JD, 2
& IM, JM, NITR, LSP, LNP, LEW )
!
!******************************************************************************
! Subroutine LOCFLT applies the pressure-filter to non-polar boxes.
! LOCFLT is called from subroutine PFILTR above (bdf, bmy, 10/11/01)
!
! Arguments as Input:
! ============================================================================
! (1 ) XERR(ID,JD) (REAL*8 ) : mass error
! (2 ) AX(ID+1,JD) (REAL*8 ) : perturbed ALFA by XERR
! (3 ) BX(ID,JD+1) (REAL*8 ) : perturbed BETA by XERR
! (4 ) AXY(ID,JD) (REAL*8 ) : area of grid box (I,J) in [m^2]
! (5-6) ID, JD (INTEGER) : "Global" array dimensions for lon, lat
! (7-8) IM, JM (INTEGER) : "Window" array dimensions for lon, lat
! (9 ) NITR (INTEGER) : number of iterations (NITR .LE. 4)
! (10 ) LSP (LOGICAL) : true if J=1 is S. POLE
! (11 ) LNP (LOGICAL) : true if J=JM is N. POLE
! (12 ) LEW (LOGICAL) : true if cyclic in W-E direction
! (i.e. if I=1 connects to I=IM)
!
! Arguments as Output:
! ============================================================================
! (1 ) XERR(ID,JD) (REAL*8 ) : adjusted mass error
! (2 ) AX(ID+1,JD) (REAL*8 ) : adjusted AX
! (3 ) BX(ID,JD+1) (REAL*8 ) : adjusted BX
!
! NOTES:
! (1 ) Adapted from original code from LLNL. Added comments and F90 syntax
! for declarations. (bdf, bmy, 10/11/01)
! (2 ) For a global run (as we usually do in GEOS-CHEM) IM=ID=IIPAR and
! JM=JD=JJPAR. (bmy, 10/11/01)
!******************************************************************************
!
IMPLICIT NONE
! Arguments
LOGICAL, INTENT(IN) :: LSP, LNP, LEW
INTEGER, INTENT(IN) :: ID, JD, IM, JM, NITR
REAL*8, INTENT(IN) :: AXY(JD)
REAL*8, INTENT(INOUT) :: XERR(ID,JD)
REAL*8, INTENT(INOUT) :: AX(ID+1,JD)
REAL*8, INTENT(INOUT) :: BX(ID,JD+1)
! Local variables
INTEGER :: I, IA, NAZ, J, J1, J2, NFLTR
REAL*8 :: SUMA, FNAZ8
REAL*8 :: X0(ID,JD)
!=================================================================
! LOCFLT begins here!
!
! Initialize corrective column mass flows (kg): AX->alfa, BX->beta
!=================================================================
DO J = 1, JM
DO I = 1, IM
X0(I,J) = XERR(I,J)
ENDDO
ENDDO
!=================================================================
! Iterate over mass-error filter
! accumulate corrections in AX & BX
!=================================================================
DO NFLTR = 1, NITR
!==============================================================
! calculate AX = E-W filter
!==============================================================
! Compute polar box limits
J1 = 1
J2 = JM
IF ( LSP ) J1 = 2
IF ( LNP ) J2 = JM - 1
! Loop over non-polar latitudes
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, FNAZ8 )
DO J = J1, J2
! Calculate pressure-filter E-W wind between boxes [I-1] & [I].
! Enhance filtered wind by size of EPZ, will redistribute
! later within
FNAZ8 = 0.125d0
DO I = 2, IM
AX(I,J) = AX(I,J) + FNAZ8 *(XERR(I-1,J) - XERR(I,J))
ENDDO
! calculate pressure-filter E-W wind at edges I=1 & I=IM+1
IF ( LEW ) THEN
AX(IM+1,J) = AX(IM+1,J) + FNAZ8 * (XERR(IM,J) -XERR(1,J))
AX(1,J) = AX(1,J) + FNAZ8 * (XERR(IM,J) -XERR(1,J))
ELSE
! WINDOW, assume zero error outside window
AX(1,J) = AX(1,J) - FNAZ8 * XERR(1,J)
AX(IM+1,J)= AX(IM+1,J) + FNAZ8 * XERR(IM,J)
ENDIF
ENDDO
!$OMP END PARALLEL DO
!==============================================================
! calculate BX = N-S filter, N-S wind between boxes [J-1] & [J]
!==============================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, FNAZ8 )
DO J = 3, JM-1
FNAZ8 = 0.25D0 * AXY(J) / ( AXY(J-1) + AXY(J) )
DO I = 1, IM
BX(I,J) = BX(I,J) + FNAZ8 * ( XERR(I,J-1) - XERR(I,J) )
ENDDO
ENDDO
!$OMP END PARALLEL DO
! enhance the filtering by factor of 2 ONLY into/out-of polar caps
FNAZ8 = 0.5D0 * AXY(2) / ( AXY(1) + AXY(2) )
! When LSP=TRUE then J=1 is SOUTH POLE
IF ( LSP ) THEN
DO I = 1, IM
BX(I,2) = BX(I,2) + FNAZ8 * (XERR(I,1) -XERR(I,2))
ENDDO
ELSE
DO I = 1, IM
BX(I,1)= BX(I,1) -0.5D0 *FNAZ8 * XERR(I,1)
BX(I,2)= BX(I,2) +0.5D0 *FNAZ8 * (XERR(I,1) - XERR(I,2))
ENDDO
ENDIF
FNAZ8 = 0.5D0 * AXY(JM) / ( AXY(JM-1) + AXY(JM) )
! When LNP=TRUE, then J=JM is NORTH POLE
IF ( LNP ) THEN
DO I = 1, IM
BX(I,JM) = BX(I,JM) +FNAZ8 *(XERR(I,JM-1) -XERR(I,JM))
ENDDO
ELSE
DO I = 1,IM
BX(I,JM+1)= BX(I,JM+1) + 0.5D0 *FNAZ8 * XERR(I,JM)
BX(I,JM) = BX(I,JM) + 0.5D0 *FNAZ8 *
& (XERR(I,JM-1) -XERR(I,JM))
ENDDO
ENDIF
!==============================================================
! need N-S flux across boundaries if window calculation
! (assume XERR=0 outside)
!
! JM for optimal matrix/looping, it would be best to
! define XERR=0 for an oversized array XERR(0:IM+1,0:JM+1)
! Update the mass error (XERR)
!==============================================================
DO J = 1, JM
DO I = 1, IM
XERR(I,J) = X0(I,J) + AX(I,J) - AX(I+1,J)
& + BX(I,J) - BX(I,J+1)
ENDDO
ENDDO
ENDDO ! NFLTR
! Return to calling program
END SUBROUTINE LOCFLT
!------------------------------------------------------------------------------
SUBROUTINE POLFLT( XERR, BX, AXY, COEF, ID, JD, IM, JM ) 2
!
!******************************************************************************
! Subroutine POLFLT applies the pressure-filter to polar boxes.
! POLFLT is called from subroutine PFILTR above (bdf, bmy, 10/10/01)
!
! Arguments as Input:
! ============================================================================
! (1 ) XERR(ID,JD) (REAL*8 ) : mass error
! (2 ) BX(ID,JD+1) (REAL*8 ) : perturbed BETA by XERR
! (3 ) AXY(ID,JD) (REAL*8 ) : area of grid box (I,J) in [m^2]
! (4 ) COEF (REAL*8 ) : Multiplicative coefficient ?????
! (5-6) ID, JD (INTEGER) : "Window" array dimensions for lon, lat
! (7-8) IM, JM (INTEGER) : "Global" array dimensions for lon, lat
!
! Arguments as Output:
! ============================================================================
! (1 ) XERR(ID,JD) (REAL*8 ) : adjusted mass error
! (2 ) BX(ID,JD+1) (REAL*8 ) : adjusted BX
!
! NOTES:
! (1 ) Adapted from original code from LLNL. Added comments and F90 syntax
! for declarations. (bdf, bmy, 10/10/01)
! (2 ) For a global run (as we usually do in GEOS-CHEM) IM=ID=IIPAR and
! JM=JD=JJPAR. (bmy, 10/20/01)
!******************************************************************************
!
IMPLICIT NONE
! Arguments
INTEGER, INTENT(IN) :: ID, JD, IM, JM
REAL*8, INTENT(IN) :: AXY(JD)
REAL*8, INTENT(IN) :: COEF
REAL*8, INTENT(INOUT) :: XERR(ID,JD)
REAL*8, INTENT(INOUT) :: BX(ID,JD+1)
! Local variables
INTEGER :: I, J
REAL*8 :: ERAV, BXJ(JD+1), TOTAL
!=================================================================
! POLFLT begins here!
!
! Initialize corrective column mass flows (kg): BXJ->beta
!=================================================================
DO I = 1, IM
! Initialize
ERAV = 0.D0
TOTAL = 0.D0
! Sum XERR in ERAV and sum AXY in TOTAL
DO J = 1, JM
ERAV = ERAV + XERR(I,J)
TOTAL = TOTAL + AXY(J)
ENDDO
! Compute area-weighted mass error total
ERAV = ERAV / TOTAL
! mass-error filter, make corrections in BX
BXJ(1) = 0.D0
DO J = 2, JM
BXJ(J) = BXJ(J-1) + XERR(I,J-1) - AXY(J-1) * ERAV
ENDDO
DO J = 2, JM
BX(I,J) = BX(I,J) + COEF * BXJ(J)
ENDDO
ENDDO ! I
! Update XERR
DO J = 1, JM
DO I = 1, IM
XERR(I,J) = XERR(I,J) + BX(I,J) - BX(I,J+1)
ENDDO
ENDDO
! Return to calling program
END SUBROUTINE POLFLT
!------------------------------------------------------------------------------
SUBROUTINE DIAG_FLUX( IC, FX, FX1_TP, FY, FY1_TP, 2,14
& FZ, FZ1_TP, NDT, ACOSP )
!
!******************************************************************************
! Subroutine DIAG_FLUX archives the mass fluxes in TPCORE version 7.1.
! (bey, bmy, 9/20/00, 7/21/05)
!
! Arguments as Input:
! ============================================================================
! (1 ) IC (INTEGER) : Current tracer #
! (2,3) FX,FX1_TP (REAL*8 ) : Flux into the west side of grid box (I,J,K)
! (4,5) FY,FY1_TP (REAL*8 ) : Flux into the south side of grid box (I,J,K)
! (6,7) FZ,FZ1_TP (REAL*8 ) : Flux into top of grid box (I,J,K)
! (8 ) NDT (INTEGER) : Dynamic timestep in seconds
! (9 ) ACOSP (INTEGER) : Inverse cosine at latitude (J)
!
! Included via header files:
! ============================================================================
! (1 ) DXYP (REAL*8) : Surface area of grid box [m2]
! (2 ) g0_100 (REAL*8) : The value 100 / 9.8
!
! Diagnostics archived:
! ============================================================================
! (1 ) ND24 : Eastward flux of tracer in kg/s
! (2 ) ND25 : Westward flux of tracer in kg/s
! (3 ) ND26 : Upward flux of tracer in kg/s
!
! NOTES:
! (1 ) Original code & algorithm is from Isabelle Bey, as installed in
! TPCORE v. 4.1 (1998, 1999)
! (2 ) DXYP is of dimension JGLOB, so reference it by DXYP(JREF),
! where JREF = J + J0. (bmy, 9/28/00)
! (3 ) Add parallel processor directives to do-loops (bmy, 9/29/00)
! (4 ) Archive CO budget array TCO for CO-OH run (bnd, bmy, 10/16/00)
! (5 ) Also archive X-trop flux for CH4 simulation in TCH4 (bmy, 1/17/01)
! (6 ) Added to "tpcore_mod.f" (bmy, 7/16/01)
! (7 ) Now replace DXYP(JREF) with routine GET_AREA_M2 of "grid_mod.f".
! Also remove all references to JREF. (bmy, 2/11/03)
! (8 ) Now references TCVV and ITS_A_CH4_SIM from "tracer_mod.f"
! (bmy, 7/20/04)
! (9 ) Remove references obsolete to CO-OH param code (bmy, 6/24/05)
! (10) Bug fix: FX should be dimensioned with IIPAR+1 and FZ should be
! dimensioned with LLPAR+1 (bmy, 7/21/05)
!******************************************************************************
!
! References to F90 modules
USE DIAG_MOD, ONLY : MASSFLEW, MASSFLNS, MASSFLUP
USE GLOBAL_CH4_MOD, ONLY : XNUMOL_CH4, TCH4
USE GRID_MOD, ONLY : GET_AREA_M2
USE TRACER_MOD, ONLY : ITS_A_CH4_SIM, TCVV
IMPLICIT NONE
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! Diagnostic switches
# include "CMN_GCTM" ! g0_100
! Arguments
INTEGER, INTENT(IN) :: IC, NDT
REAL*8, INTENT(IN) :: FX(IIPAR+1,JJPAR,LLPAR)
REAL*8, INTENT(IN) :: FX1_TP(IIPAR,JJPAR,LLPAR)
REAL*8, INTENT(IN) :: FY(IIPAR,JJPAR,LLPAR)
REAL*8, INTENT(IN) :: FY1_TP(IIPAR,JJPAR,LLPAR)
REAL*8, INTENT(IN) :: FZ(IIPAR,JJPAR,LLPAR+1)
REAL*8, INTENT(IN) :: FZ1_TP(IIPAR,JJPAR,LLPAR)
REAL*8, INTENT(IN) :: ACOSP(JJPAR)
! Local variables
INTEGER :: I, J, K, K2
REAL*8 :: DTC, DTDYN, AREA_M2
!=================================================================
! DIAG_FLUX begins here!
!
! FX, FX1_TP, FY, FY1_TP, FZ, FZ1_TP have units of [mb/timestep].
!
! To get tracer fluxes in kg/s :
! * (100./9.8) => kg/m2
! * DXYP(J)/(DTDYN * TCVV(IC)) => kg/s
!
! Direction of the fluxes :
! ----------------------------------------------------------------
! FX(I,J,K) => flux coming into the west edge of the box I
! (from I-1 to I).
! => a positive flux goes from west to east.
!
! FY(I,J,K) => flux coming into the south edge of the box J
! (from J to J-1).
! => a positive flux goes from south to north
! (from J-1 to J)
!
! FZ(I,J,K) => flux coming down into the box k.
! => a positive flux goes down.
!=================================================================
! DTDYN = double precision value for NDT, the dynamic timestep
DTDYN = DBLE( NDT )
!=================================================================
! ND24 Diagnostic: Eastward flux of tracer in [kg/s]
!=================================================================
IF ( ND24 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K, K2, AREA_M2, DTC )
DO K = 1, LLPAR
! K is the vertical index down from the atmosphere top downwards
! K2 is the vertical index up from the surface
K2 = LLPAR - K + 1
DO J = 1, JJPAR
! Grid box surface area [m2]
AREA_M2 = GET_AREA_M2( J )
DO I = 1, IIPAR
DTC = ( FX(I,J,K) + FX1_TP(I,J,K) ) *
& ( g0_100 * AREA_M2 ) /
& ( TCVV(IC) * DTDYN )
MASSFLEW(I,J,K2,IC) = MASSFLEW(I,J,K2,IC) + DTC
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
!=================================================================
! ND25 Diagnostic: Northward flux of tracer in [kg/s]
!=================================================================
IF ( ND25 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K, K2, AREA_M2, DTC )
DO K = 1, LLPAR
! K is the vertical index down from the atmosphere top downwards
! K2 is the vertical index up from the surface
K2 = LLPAR - K + 1
DO J = 1, JJPAR
! Grid box surface area [m2]
AREA_M2 = GET_AREA_M2( J )
DO I = 1, IIPAR
DTC = ( FY(I,J,K) + FY1_TP(I,J,K) ) *
& ( ACOSP(J) * g0_100 * AREA_M2 ) /
& ( TCVV(IC) * DTDYN )
! Contribution for CH4 run (bmy, 1/17/01)
IF ( ITS_A_CH4_SIM() ) THEN
TCH4(I,J,K,10) = TCH4(I,J,K,10) +
& ( DTC * DTDYN * XNUMOL_CH4 )
ENDIF
MASSFLNS(I,J,K2,IC) = MASSFLNS(I,J,K2,IC) + DTC
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
!=================================================================
! ND26 Diagnostic : Upward flux of tracer in [kg/s]
!=================================================================
IF ( ND26 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, K, K2, DTC )
DO K = 1, LLPAR
! K is the vertical index down from the atmosphere top downwards
! K2 is the vertical index up from the surface
K2 = LLPAR - K + 1
DO J = 1, JJPAR
! Grid box surface area [m2]
AREA_M2 = GET_AREA_M2( J )
DO I = 1, IIPAR
DTC = ( FZ(I,J,K) + FZ1_TP(I,J,K) ) *
& ( g0_100 * AREA_M2 ) /
& ( TCVV(IC) * DTDYN )
MASSFLUP(I,J,K2,IC) = MASSFLUP(I,J,K2,IC) + DTC
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE DIAG_FLUX
!------------------------------------------------------------------------------
END MODULE TPCORE_MOD