!
! $Id: pdyn2.F90,v 1.1 2005/10/18 20:10:17 mbhat Exp $
!
!---(pdyn2.f)-------generic CTM shell from UCIrvine (p-code 1.0, 1/31/95)
!---subroutines: DYN2U, DYN2V, DYN2W, QXZON, QZONX, QLIMIT, QVECT, FLIMIT
SUBROUTINE DYN2U(DTALFA) 4,3
!---WEST-EAST ADVECTION OF TRACE COMPOUNDS---V 10.0 (S.O.M.)
REAL*8 DTALFA
REAL*8 QM(99), & ! air mass in box at start, !**=REDEFINED after flow
& QU(100), & ! air mass to be moved (I-1) -> (I) in advective step
& QBC1, QBC2 ! tracer mixing ratio boundary conditions at (1)&(99)
LOGICAL LCYCLE ! .T. = cyclic b.c., .F. = mixing ratio b.c.
REAL*8 QTT(99) !** tracer mass in box (I)
REAL*8 QXT(99),QXX(99), & !** 1st & 2nd moments of tracer in X-direct
& QXY(99),QXZ(99), & !** coupled 1st moments including X component
& QYZ(99),QYT(99),QYY(99),QZT(99),QZZ(99), & !** other SOMs without X
& Q0F(99+1) !** computed tracer flux from (I-1) to (I)
# include "uci.h"
!xx$$
save
!xx$$
!-----------------------------------------------------------------------
DO 50 L=1,LM
DO 40 J=1,JM
IF(J+J0.EQ.1 .OR. J+J0.EQ.JMX) GOTO 40 !NO E-W advect @ Poles
IMB = IMRSLV(J+J0)
FIMB = IMB
N1N2 = (IMB+1)*(IMB+2)
NQ = IM/IMB
!---for each tracer calculate 'pipe-flow' in E-W direction
DO 30 N=1,NTRACE
!---
IF(IMB.EQ.1) THEN
!---transfer STT() & moments into local array for piped flow
DO I=1,NQ
QM(I) = MAX(AD(I,J,L),0.0)
QU(I) = ALFA(I,J,L)*DTALFA
QTT(I) = STT(I,J,L,N)
QXT(I) = SUT(I,J,L,N)
QYT(I) = SVT(I,J,L,N)
QZT(I) = SWT(I,J,L,N)
QXX(I) = SUU(I,J,L,N)
QYY(I) = SVV(I,J,L,N)
QZZ(I) = SWW(I,J,L,N)
QXY(I) = SUV(I,J,L,N)
QXZ(I) = SUW(I,J,L,N)
QYZ(I) = SVW(I,J,L,N)
ENDDO
ELSE
!---create NQ-EXTENDED ZONES near poles, assume equal masses (i.e. same Psurf)
DO IQ=1,NQ
I = 1 + IMB*(IQ-1)
QU(IQ) = ALFA(I,J,L)*DTALFA
!---combine IMB equal-mass boxes beginning at I:
CALL QXZON
( STT(I,J,L,N),SUT(I,J,L,N),SUU(I,J,L,N),SUV(I,J,L,N), &
& SUW(I,J,L,N),SVT(I,J,L,N),SVV(I,J,L,N),SWT(I,J,L,N), &
& SWW(I,J,L,N),SVW(I,J,L,N), AD(I,J,L), IMB, &
& QTT(IQ),QXT(IQ),QXX(IQ),QXY(IQ),QXZ(IQ),QYT(IQ),QYY(IQ), &
& QZT(IQ),QZZ(IQ),QYZ(IQ),QM(IQ))
ENDDO
ENDIF
!---WINDOWS: tracer in (IM+1,J,L) is eastern b.c,
!--- (IM+2,J,L) is western b.c, adjacent to (1,J,L)
IF(IM.LT.IMX) THEN
LCYCLE = .FALSE.
QBC1 = STT(IM+2,J,L,N)/AD(IM+2,J,L)
QBC2 = STT(IM+1,J,L,N)/AD(IM+1,J,L)
QU(NQ+1) = ALFA(IM+1,J,L)*DTALFA
ELSE
!---PERIODIC boundary conditions, (IM,J,L) connects to (1,J,L)
LCYCLE = .TRUE.
QBC1 = 0.0
QBC2 = 0.0
QU(NQ+1) = 0.0
ENDIF
!------------------------------------------------------
CALL QVECT(QM,QU,QBC1,QBC2,NQ,LCYCLE, LMTSOM, &
& QTT,QXT,QXX,QXY,QXZ,QYT,QYY,QZT,QZZ,QYZ, Q0F)
!------------------------------------------------------
IF(IMB.EQ.1) THEN
!---transfer Q__() back to STT() & moments
DO I=1,NQ
STT(I,J,L,N) = QTT(I)
SUT(I,J,L,N) = QXT(I)
SVT(I,J,L,N) = QYT(I)
SWT(I,J,L,N) = QZT(I)
SUU(I,J,L,N) = QXX(I)
SVV(I,J,L,N) = QYY(I)
SWW(I,J,L,N) = QZZ(I)
SUV(I,J,L,N) = QXY(I)
SUW(I,J,L,N) = QXZ(I)
SVW(I,J,L,N) = QYZ(I)
ENDDO
ELSE
!---remap EXTENDED ZONES near poles from Q__() back to STT():
DO IQ=1,NQ
I = 1 + IMB*(IQ-1)
CALL QZONX( QTT(IQ),QXT(IQ),QXX(IQ),QXY(IQ),QXZ(IQ), &
& QYT(IQ),QYY(IQ),QZT(IQ),QZZ(IQ),QYZ(IQ),QM(IQ), &
& STT(I,J,L,N),SUT(I,J,L,N),SUU(I,J,L,N),SUV(I,J,L,N), &
& SUW(I,J,L,N),SVT(I,J,L,N),SVV(I,J,L,N),SWT(I,J,L,N), &
& SWW(I,J,L,N),SVW(I,J,L,N), AD(I,J,L), IMB,LMTSOM)
!---simple fill for W->E flux across boundaries in extended zones
DO II=I+1,I+IMB-1
Q0F(II) = Q0F(I)
ENDDO
ENDDO
ENDIF
!---DIAGNOSTICS: flux across meridion, into (W->E) boxes ID18/ID19
IF(ND18.GT.0) AJL(J,L,ND18+N) = AJL(J,L,ND18+N) + Q0F(ID18)
IF(ND19.GT.0) AJL(J,L,ND19+N) = AJL(J,L,ND19+N) + Q0F(ID19)
30 CONTINUE
!---can NOW reset dry air mass in boxes along E-W pipe, as calculated in QVECT
IF(IMB.EQ.1) THEN
DO I=1,NQ
AD(I,J,L) = QM(I)
ENDDO
ELSE
DO IQ=1,NQ
QM(IQ) = QM(IQ)/REAL(IMB)
I = 1 + IMB*(IQ-1)
DO II=I,I+IMB-1
AD(II,J,L) = QM(IQ)
ENDDO
ENDDO
ENDIF
40 CONTINUE
50 CONTINUE
RETURN
END
SUBROUTINE DYN2V(DTBETA) 3,3
!---SOUTH-NORTH ADVECTION OF TRACE COMPOUNDS---V 10.0 (S.O.M.)
!xx DOUBLE PRECISION STTPOL, STTNP,STTSP
real*8 STTPOL, STTNP,STTSP
REAL*8 SWTPOL,SWWPOL,ADDPOL,SWTNP,SWTSP,SWWNP,SWWSP,ADDNP,ADDSP
REAL*8 DTBETA
REAL*8 QM(99), & ! air mass in box at start, !**=REDEFINED after flow
& QU(100), & ! air mass to be moved (I-1) -> (I) in advective step
& QBC1, QBC2 ! tracer mixing ratio boundary conditions at (1)&(99)
LOGICAL LCYCLE ! .T. = cyclic b.c., .F. = mixing ratio b.c.
REAL*8 QTT(99) !** tracer mass in box (I)
REAL*8 QXT(99),QXX(99), & !** 1st & 2nd moments of tracer in X-direct
& QXY(99),QXZ(99), & !** coupled 1st moments including X component
& QYZ(99),QYT(99),QYY(99),QZT(99),QZZ(99), & !** other SOMs without X
& Q0F(99+1) !** computed tracer flux from (I-1) to (I)
!xx DOUBLE PRECISION ETT ! tracer mass in extended zone
real*8 ETT ! tracer mass in extended zone
REAL*8 EXT,EXX,EXY,EXZ,EYZ,EYT,EYY,EZT,EZZ, & ! moments in extended zone
& EM ! dry air mass in extended zone
# include "uci.h"
!xx$$
save
!xx$$
!-----------------------------------------------------------------------
!---1st DO N-S ADVECTION IN MERIDIONAL STRIPS, 2nd SMOOTH OVER EXT.ZONES
!---assume BETA(I,J,L) has been smoothed in setups
!---WINDOWS: southern b.c is south pole if J+J0 begins at 1 (i.e., J0=0)
!--- northern b.c. is north pole if JM+J0 = JMX
!-----------------------------------------------------------------------
DO 60 L=1,LM
!---SETUP for poles
!---collect POLAR BOXES (pie-wedges), sum over I, and put sum in each box
!---do N-S transport with full POLAR CIRCLE for each meridion I, fix later
DO 20 J=1,JM,JM-1
IF(J+J0.EQ.1 .OR. J+J0.EQ.JMX) THEN
!---check if doing polar boxes, put dry air mass total in each pie-wedge
ADDPOL = 0.D0
DO I=1,IM
ADDPOL = ADDPOL + AD(I,J,L)
ENDDO
DO I=1,IM
AD(I,J,L) = ADDPOL
ENDDO
!---put tracer mass & moments total in each pie-wedge, zero all moments ex W
DO N=1,NTRACE
STTPOL = 0.0D0
SWTPOL = 0.0
SWWPOL = 0.0
DO I=1,IM
STTPOL = STTPOL + STT(I,J,L,N)
SWTPOL = SWTPOL + SWT(I,J,L,N)
SWWPOL = SWWPOL + SWW(I,J,L,N)
ENDDO
DO I=1,IM
STT(I,J,L,N) = STTPOL
SWT(I,J,L,N) = SWTPOL
SWW(I,J,L,N) = SWWPOL
SVT(I,J,L,N) = 0.0
SVV(I,J,L,N) = 0.0
SUT(I,J,L,N) = 0.0
SUU(I,J,L,N) = 0.0
SUV(I,J,L,N) = 0.0
SVW(I,J,L,N) = 0.0
SUW(I,J,L,N) = 0.0
ENDDO
ENDDO
ENDIF
20 CONTINUE
!---loop over each meridion:
DO 40 I=1,IM
!---store polar boxes for readjustment later (only used if J0=0, J0+JM=JMX)
ADDSP = AD(I,1,L)
ADDNP = AD(I,JM,L)
!---for each tracer calculate 'pipe-flow' in N-S direction
DO 30 N=1,NTRACE
STTSP = STT(I,1,L,N)
SWTSP = SWT(I,1,L,N)
SWWSP = SWW(I,1,L,N)
STTNP = STT(I,JM,L,N)
SWTNP = SWT(I,JM,L,N)
SWWNP = SWW(I,JM,L,N)
DO J=1,JM
QM(J) = MAX(AD(I,J,L),0.0)
QU(J) = BETA(I,J,L)*DTBETA
QTT(J) = STT(I,J,L,N)
QXT(J) = SUT(I,J,L,N)
QYT(J) = SVT(I,J,L,N)
QZT(J) = SWT(I,J,L,N)
QXX(J) = SUU(I,J,L,N)
QYY(J) = SVV(I,J,L,N)
QZZ(J) = SWW(I,J,L,N)
QXY(J) = SUV(I,J,L,N)
QXZ(J) = SUW(I,J,L,N)
QYZ(J) = SVW(I,J,L,N)
ENDDO
!----N-S flow is never periodic:
LCYCLE = .FALSE.
NQ = JM
!---WINDOWS: tracer in (I,JM+1,L) is northern b.c.
IF(J0+JM.EQ.JMX) THEN
QBC2 = 0.0
QU(NQ+1) = 0.0
ELSE
QBC2 = STT(I,JM+1,L,N)/AD(I,JM+1,L)
QU(NQ+1) = BETA(I,JM+1,L)*DTBETA
ENDIF
!---WNDOWS: tracer in (I,JM+2,L) is southern b.c, adjacent to (I,1,L)
IF(J0.EQ.0) THEN
QBC1 = 0.0
QU(1) = 0.0
ELSE
QBC1 = STT(I,JM+2,L,N)/AD(I,JM+2,L)
ENDIF
!---
CALL QVECT(QM,QU,QBC1,QBC2,NQ,LCYCLE, LMTSOM, &
& QTT,QYT,QYY,QYZ,QXY,QZT,QZZ,QXT,QXX,QXZ, Q0F)
!---
DO J=1,NQ
STT(I,J,L,N) = QTT(J)
SUT(I,J,L,N) = QXT(J)
SVT(I,J,L,N) = QYT(J)
SWT(I,J,L,N) = QZT(J)
SUU(I,J,L,N) = QXX(J)
SVV(I,J,L,N) = QYY(J)
SWW(I,J,L,N) = QZZ(J)
SUV(I,J,L,N) = QXY(J)
SUW(I,J,L,N) = QXZ(J)
SVW(I,J,L,N) = QYZ(J)
ENDDO
!---now fix POLAR BOXES by removing the (IM-1) pie wedges put in for transport
!---N.B. the remaining mass (STT & AD) may be negative, and must be left so
!--- they can be averaged later
IF(J0.EQ.0) THEN
STT(I,1,L,N) = STT(I,1,L,N) - STTSP*REAL(IM-1)/REAL(IM)
SWT(I,1,L,N) = SWT(I,1,L,N) - SWTSP*REAL(IM-1)/REAL(IM)
SWW(I,1,L,N) = SWW(I,1,L,N) - SWWSP*REAL(IM-1)/REAL(IM)
ENDIF
IF(J0+JM.EQ.JMX) THEN
STT(I,JM,L,N) = STT(I,JM,L,N) - STTNP*REAL(IM-1)/REAL(IM)
SWT(I,JM,L,N) = SWT(I,JM,L,N) - SWTNP*REAL(IM-1)/REAL(IM)
SWW(I,JM,L,N) = SWW(I,JM,L,N) - SWWNP*REAL(IM-1)/REAL(IM)
ENDIF
!---diagnostics of north_south flux
IF(ND10.GT.0) THEN
JMAX = MIN(JM+1,JMX)
DO J=1,JMAX
AJL(J,L,N+ND10) = AJL(J,L,N+ND10) + Q0F(J)
ENDDO
ENDIF
!---flux across latitude circle (north = +)
IF(ND51.GT.0) AIL(I,L,N+ND51) = AIL(I,L,N+ND51) + Q0F(JD51)
IF(ND52.GT.0) AIL(I,L,N+ND52) = AIL(I,L,N+ND52) + Q0F(JD52)
30 CONTINUE
!---reset dry air mass
DO J=1,NQ
AD(I,J,L) = QM(J)
ENDDO
IF(J0.EQ.0) THEN
AD(I,1,L) = AD(I,1,L) - ADDSP*REAL(IM-1)/REAL(IM)
ENDIF
IF(J0+JM.EQ.JMX) THEN
AD(I,JM,L) = AD(I,JM,L) - ADDNP*REAL(IM-1)/REAL(IM)
ENDIF
40 CONTINUE
60 CONTINUE
!-----------------------------------------------------------------------
!---now average both tracer & dry air mass over EXTENDED ZONES, incl POLES
DO 72 J=1,JM
IMB = IMRSLV(J+J0)
IF(IMB.EQ.1) GOTO 72
NQ = IM/IMB
DO 70 L=1,LM
DO N=1,NTRACE
DO IQ=1,NQ
I = 1 + IMB*(IQ-1)
!---put tracer into EXTENDED ZONES near poles (assumes equal-mass boxes)
!--- NQ-zones, each of size IMB
CALL QXZON( STT(I,J,L,N),SUT(I,J,L,N),SUU(I,J,L,N),SUV(I,J,L,N), &
& SUW(I,J,L,N),SVT(I,J,L,N),SVV(I,J,L,N),SWT(I,J,L,N), &
& SWW(I,J,L,N),SVW(I,J,L,N), AD(I,J,L), IMB, &
& ETT,EXT,EXX,EXY,EXZ,EYT,EYY,EZT,EZZ,EYZ,EM)
!---if complete zone, then zero out E-W moments (e.g., poles)
IF(IMB.EQ.IMX) THEN
EXT = 0.0
EXX = 0.0
EXY = 0.0
EXZ = 0.0
ENDIF
!---if polar zone, ALSO zero out N-S moments (later fix for N-S moments)
IF(J+J0.EQ.1 .OR. J+J0.EQ.JMX) THEN
EYT = 0.0
EYY = 0.0
EYZ = 0.0
ENDIF
!---remap EXTENDED ZONES to STT(): (includes redef of AD())
CALL QZONX( ETT,EXT,EXX,EXY,EXZ,EYT,EYY,EZT,EZZ,EYZ,EM, &
& STT(I,J,L,N),SUT(I,J,L,N),SUU(I,J,L,N),SUV(I,J,L,N), &
& SUW(I,J,L,N),SVT(I,J,L,N),SVV(I,J,L,N),SWT(I,J,L,N), &
& SWW(I,J,L,N),SVW(I,J,L,N), AD(I,J,L), IMB,LMTSOM)
ENDDO
ENDDO
70 CONTINUE
72 CONTINUE
RETURN
END
SUBROUTINE DYN2W(DTGAMA) 2,1
!---------VERTICAL ADVECTION OF TRACE COMPOUNDS---V 10.0 (S.O.M.)
REAL*8 DTGAMA
REAL*8 QM(99), & ! air mass in box at start, !**=REDEFINED after flow
& QU(100), & ! air mass to be moved (I-1) -> (I) in advective step
& QBC1, QBC2 ! tracer mixing ratio boundary conditions at (1)&(99)
LOGICAL LCYCLE ! .T. = cyclic b.c., .F. = mixing ratio b.c.
REAL*8 QTT(99) !** tracer mass in box (I)
REAL*8 QXT(99),QXX(99), & !** 1st & 2nd moments of tracer in X-direct
& QXY(99),QXZ(99), & !** coupled 1st moments including X component
& QYZ(99),QYT(99),QYY(99),QZT(99),QZZ(99), & !** other SOMs without X
& Q0F(99+1) !** computed tracer flux from (I-1) to (I)
# include "uci.h"
!xx$$
save
!xx$$
!-----------------------------------------------------------------------
!-----CALCULATE VERTICAL TRANSPORTS SEPARATELY FOR BOXES IN EXTENDED ZONE
!-----****ASSUME THAT GAMA() HAS BEEN AVERAGED OVER ZONE.****
!-----DO all pie-wedges in polar boxes (for simplicity)
!-----------------------------------------------------------------------
DO 50 J=1,JM
DO 40 I=1,IM
!---for each tracer calculate 'pipe-flow' in UP-DOWN direction
DO 30 N=1,NTRACE
!---
DO L=1,LM
QM(L) = MAX(AD(I,J,L),0.0)
QU(L) = GAMA(I,J,L)*DTGAMA
QTT(L) = STT(I,J,L,N)
QXT(L) = SUT(I,J,L,N)
QYT(L) = SVT(I,J,L,N)
QZT(L) = SWT(I,J,L,N)
QXX(L) = SUU(I,J,L,N)
QYY(L) = SVV(I,J,L,N)
QZZ(L) = SWW(I,J,L,N)
QXY(L) = SUV(I,J,L,N)
QXZ(L) = SUW(I,J,L,N)
QYZ(L) = SVW(I,J,L,N)
ENDDO
NQ = LM
!----vertical transport always with rigid boundaries:
LCYCLE = .FALSE.
QBC1 = 0.0
QBC2 = 0.0
QU(1) = 0.0
QU(NQ+1) = 0.0
!---
CALL QVECT(QM,QU,QBC1,QBC2,NQ,LCYCLE, LMTSOM, &
& QTT,QZT,QZZ,QXZ,QYZ,QXT,QXX,QYT,QYY,QXY, Q0F)
!---
DO L=1,LM
STT(I,J,L,N) = QTT(L)
SUT(I,J,L,N) = QXT(L)
SVT(I,J,L,N) = QYT(L)
SWT(I,J,L,N) = QZT(L)
SUU(I,J,L,N) = QXX(L)
SVV(I,J,L,N) = QYY(L)
SWW(I,J,L,N) = QZZ(L)
SUV(I,J,L,N) = QXY(L)
SUW(I,J,L,N) = QXZ(L)
SVW(I,J,L,N) = QYZ(L)
!---have dropped 'eddy' vertical flux. need better diagnostic
IF(ND13.GT.0) AJL(J,L,N+ND13) = AJL(J,L,N+ND13) + Q0F(L+1)
ENDDO
!---vertical transport thru TOP of layer (LD41/2) --> (LD41/2 + 1)
!--- need to fix vertical GAMA & DIAG's to mean: L=1 -> layer (1)
IF(ND41.GT.0) AIJ(I,J,N+ND41) = AIJ(I,J,N+ND41) + Q0F(LD41+1)
IF(ND42.GT.0) AIJ(I,J,N+ND42) = AIJ(I,J,N+ND42) + Q0F(LD42+1)
30 CONTINUE
!---have finished all tracers, can now reset mass in column (I,J)
DO L=1,LM
AD(I,J,L) = QM(L)
ENDDO
40 CONTINUE
50 CONTINUE
RETURN
END
!-----qxzon.f, qzonx.f
SUBROUTINE QXZON(QTT,QXT,QXX,QXY,QXZ,QYT,QYY,QZT,QZZ,QYZ,QM, NQ, & 2
& ETT,EXT,EXX,EXY,EXZ,EYT,EYY,EZT,EZZ,EYZ,EM)
!---Combines NQ equal mass boxes (QTT,QXT,QXX,...,QYX) aligned in X-direction
!--- into 1 Extended zone (ETT,EXT,EXX,...,EYZ)
IMPLICIT NONE
INTEGER NQ ! # of boxes in initial distribution to be combined
Real*8 QTT(NQ) ! tracer mass in box (I)
Real*8 QXT(NQ),QXX(NQ), & ! 1st & 2nd moments of tracer in X-direct
& QXY(NQ),QXZ(NQ), & ! coupled 1st moments including X component
& QYZ(NQ),QYT(NQ),QYY(NQ),QZT(NQ),QZZ(NQ), & ! other SOMs without X
& QM(NQ) ! dry air mass in box (I)
Real*8 ETT ! tracer mass in extended zone
Real*8 EXT,EXX,EXY,EXZ,EYZ,EYT,EYY,EZT,EZZ, & ! moments in extended zone
& EM ! dry air mass in extended zone
!---local variables
Real*8 SUMT,SUMF
Real*8 FLTNQ,FLTIQ
INTEGER N1N2,IQ
!xx$$
save
!xx$$
!---since assume equal-mass boxes, rescale STT() so that no spurious moments
!--- from variations in AIRD() with constant STT/AD
SUMT = 0.D0
SUMF = 0.D0
DO IQ=1,NQ
SUMT = SUMT + QTT(IQ)
IF (QM(IQ).GT.0.D0) QTT(IQ) = QTT(IQ)/QM(IQ)
SUMF = SUMF + QTT(IQ)
ENDDO
IF (SUMF.GT.0.D0) SUMT = SUMT/SUMF
DO IQ=1,NQ
QTT(IQ) = QTT(IQ)*SUMT
ENDDO
!---
FLTNQ = NQ
N1N2 = (NQ+1)*(NQ+2)
ETT = 0.0D0
EXT = 0.0
EXX = 0.0
EXY = 0.0
EXZ = 0.0
EYZ = 0.0
EYT = 0.0
EYY = 0.0
EZT = 0.0
EZZ = 0.0
EM = 0.0
!---
DO IQ=1,NQ
FLTIQ = REAL(IQ+IQ-NQ-1)
ETT = ETT + QTT(IQ)
EXT = EXT + QXT(IQ) + QTT(IQ)*FLTIQ*3.
EYT = EYT + QYT(IQ)
EZT = EZT + QZT(IQ)
EXX = EXX + QXX(IQ) + QXT(IQ)*FLTIQ*5. + &
& QTT(IQ)*REAL(5*(6*IQ*(IQ-NQ-1)+N1N2))
EYY = EYY + QYY(IQ)
EZZ = EZZ + QZZ(IQ)
EXY = EXY + QXY(IQ) + QYT(IQ)*FLTIQ*3.
EXZ = EXZ + QXZ(IQ) + QZT(IQ)*FLTIQ*3.
EYZ = EYZ + QYZ(IQ)
EM = EM + QM(IQ)
ENDDO
EXT = EXT/FLTNQ
EXX = EXX/(FLTNQ*FLTNQ)
EXY = EXY/FLTNQ
EXZ = EXZ/FLTNQ
RETURN
END
SUBROUTINE QZONX(ETT,EXT,EXX,EXY,EXZ,EYT,EYY,EZT,EZZ,EYZ,EM, & 2,1
& QTT,QXT,QXX,QXY,QXZ,QYT,QYY,QZT,QZZ,QYZ,QM, NQ,LIMTR)
!---Divide one Extended zone (ETT,EXT,EXX,...,EYZ) into equal-mass
!--- boxes (QTT,QXT,QXX,...,QYX) aligned in X-direction
IMPLICIT NONE
Integer NQ,LIMTR ! # of boxes in initial distribution to be combined & lim
Real*8 QTT(NQ) ! tracer mass in box (I)
Real*8 QXT(NQ),QXX(NQ), & ! 1st & 2nd moments of tracer in X-direct
& QXY(NQ),QXZ(NQ), & ! coupled 1st moments including X component
& QYZ(NQ),QYT(NQ),QYY(NQ),QZT(NQ),QZZ(NQ), & ! other SOMs without X
& QM(NQ) ! dry air mass in box (I)
Real*8 ETT ! tracer mass in extended zone
Real*8 EXT,EXX,EXY,EXZ,EYZ,EYT,EYY,EZT,EZZ, & ! moments in extended zone
& EM ! dry air mass in extended zone
!---local variables
Real*8 ETT0,EXT0,EXX0, FLTNQ8,FLTIQ,FLTIQ2
Real*8 FLTNQ, EXY0,EXZ0,EYZ0,EYT0,EYY0,EZT0,EZZ0,EM0
INTEGER N1N2,IQ
!xx$$
save
!xx$$
!---
FLTNQ = NQ
FLTNQ8= NQ
N1N2 = (NQ+1)*(NQ+2)
!---N.B. have not scaled Y,Z moments (EYT,...) to ETT, because not reversible
!---must call positive limiter to ensure So > 0 in sub-boxes
CALL QLIMIT(ETT,EXT,EXX,EXY,EXZ,EYT,EYY,EZT,EZZ,EYZ,EM,LIMTR)
ETT0 = ETT/FLTNQ8
EXT0 = EXT/(FLTNQ*FLTNQ)
EXX0 = EXX/(FLTNQ*FLTNQ*FLTNQ)
EXY0 = EXY/(FLTNQ*FLTNQ)
EXZ0 = EXZ/(FLTNQ*FLTNQ)
EYT0 = EYT/FLTNQ
EZT0 = EZT/FLTNQ
EYY0 = EYY/FLTNQ
EZZ0 = EZZ/FLTNQ
EYZ0 = EYZ/FLTNQ
EM0 = EM/FLTNQ
DO IQ=1,NQ
FLTIQ = (IQ+IQ-NQ-1)
FLTIQ2 = (6*IQ*(IQ-NQ-1)+N1N2)
QTT(IQ) = ETT0 + FLTIQ*EXT0 + FLTIQ2*EXX0
QXT(IQ) = EXT0 + 3.*FLTIQ*EXX0
QXX(IQ) = EXX0
QYT(IQ) = EYT0 + FLTIQ*EXY0
QZT(IQ) = EZT0 + FLTIQ*EXZ0
QXY(IQ) = EXY0
QXZ(IQ) = EXZ0
QYY(IQ) = EYY0
QZZ(IQ) = EZZ0
QYZ(IQ) = EYZ0
QM(IQ) = EM0
ENDDO
RETURN
END
SUBROUTINE QLIMIT(ETT,EXT,EXX,EXY,EXZ,EYT,EYY,EZT,EZZ,EYZ,EM,LIM) 2
!---Place limits on moments in box to attain some properties:
!--- LIM=0 place no constraints, allow negative tracer
!--- LIM=1 require positive definite along X-axis (averaged over Y,Z)
!--- this is the standard limiter on the original CTM
!--- LIM=2 force positive, monotonic along X-axis (averaged over Y,Z)
IMPLICIT NONE
INTEGER LIM
Real*8 ETT ! tracer mass in extended zone
Real*8 EXT,EXX, & ! 1st & 2nd moments of tracer in X-direct
& EXY,EXZ,EYZ,EYT,EYY,EZT,EZZ, & ! other moments
& EM ! dry air mass in extended zone
Real*8 Q0N,Q1N,Q2N,Q3N
!xx$$
save
!xx$$
!---
IF (LIM.EQ.0) GOTO 99
IF (LIM.EQ.1) THEN
!---positive definite, reset Sxx to allow maximum Sx (1.5*So)
Q0N = 0.0
IF(ETT.GT.1.0D-35) Q0N = ETT
Q3N = 1.499*Q0N
Q1N = MIN(+Q3N,MAX(-Q3N,EXT))
Q2N = MIN((Q0N+Q0N-0.334*ABS(Q1N)),MAX(ABS(Q1N)-Q0N,EXX))
EXY = MIN(+Q0N,MAX(-Q0N,EXY))
EXZ = MIN(+Q0N,MAX(-Q0N,EXZ))
EXT = Q1N
EXX = Q2N
!---NEW, restrict other moments to +-1.5*mass
! EYZ = MIN(+Q3N,MAX(-Q3N,EYZ))
! EYT = MIN(+Q3N,MAX(-Q3N,EYT))
! EYY = MIN(+Q3N,MAX(-Q3N,EYY))
! EZT = MIN(+Q3N,MAX(-Q3N,EZT))
! EZZ = MIN(+Q3N,MAX(-Q3N,EZZ))
ENDIF
!-----
IF (LIM.GE.2) THEN
Q0N = 0.0
IF(ETT.GT.1.0D-35) Q0N = ETT
Q3N = 1.50*Q0N
Q1N = EXT
Q2N = EXX
IF(Q2N.GT.0.0) THEN
!----do not allow reversal of curvature: Q2N > 0
Q2N = MIN(Q2N, ABS(Q1N)/3.0, 0.50*Q0N)
IF (Q1N.LT.0.0) THEN
Q1N = MAX(-(Q0N+Q2N), Q1N)
ELSE
Q1N = MIN(+(Q0N+Q2N), Q1N)
ENDIF
ELSE
!------Sxx < 0 = curved down at ends, allow if Sx < So
Q1N = MIN(+Q0N ,MAX(-Q0N, Q1N))
Q2N = MAX(Q2N, -Q0N+ABS(Q1N), -ABS(Q1N)/3.0)
ENDIF
EXY = MIN(+Q0N, MAX(-Q0N, EXY))
EXZ = MIN(+Q0N, MAX(-Q0N, EXZ))
EXT = Q1N
EXX = Q2N
!---NEW, restrict other moments to +-1.5*mass
EYZ = MIN(+Q3N,MAX(-Q3N,EYZ))
EYT = MIN(+Q3N,MAX(-Q3N,EYT))
EYY = MIN(+Q3N,MAX(-Q3N,EYY))
EZT = MIN(+Q3N,MAX(-Q3N,EZT))
EZZ = MIN(+Q3N,MAX(-Q3N,EZZ))
ENDIF
99 RETURN
END
!-----qvect.f: hybrid version for ZRcode tests with LIMTR=4 (mjp 9/94)
SUBROUTINE QVECT(QM,QU,QBC1,QBC2,NQ,LCYCLE, LIMTR, & 4,2
& QTT,QXT,QXX,QXY,QXZ,QYT,QYY,QZT,QZZ,QYZ, Q0F)
! Computes SOM advection for a single 'piped flow' of tracer in the X-direction
!-----LCYCLE=.T.= connect QM(NQ) to QM(1), assume QU cyclic: QU(NQ+1)=QU(1)
!----- =QU(1)=> BOX(1) =QU(2)=> ... =QU(NQ)=> BOX(NQ) =QU(1)=> BOX(1) ...
!-----LCYCLE=.F.= use mixing ratio b.c. at (1) and (NQ), need QU(NQ+1)
!----- (QBC1) =QU(1)=> BOX(1) =QU(2)=>... =QU(NQ)=> BOX(NQ) =QU(NQ+1)=> (QBC2)
!
!--- LIMTR=3 = monotonic internal, with min/max from connecting boxes,
!--- QXT, QXX adjusted only.
!--- =4 = bounds, min/max limits applied to ALL X-terms: QXY, QXZ
!--- ****not recommended, smooths out strat isolines too much***
!
IMPLICIT NONE
INTEGER NQ,LIMTR ! # of boxes, does not include boundary conditions&LIM
Real*8 QM(NQ), & ! air mass in box at start, !**=REDEFINED after flow
& QU(NQ+1), & ! air mass to be moved (I-1) -> (I) in advective step
& QBC1, QBC2 ! tracer mixing ratio boundary conditions at (1)&(NQ)
LOGICAL LCYCLE ! .T. = cyclic b.c., .F. = mixing ratio b.c.
Real*8 QTT(NQ) !** tracer mass in box (I)
Real*8 QXT(NQ),QXX(NQ), & !** 1st & 2nd moments of tracer in X-direct
& QXY(NQ),QXZ(NQ), & !** coupled 1st moments including X component
& QYZ(NQ),QYT(NQ),QYY(NQ),QZT(NQ),QZZ(NQ), & !** other SOMs without X
& Q0F(NQ+1) !** computed tracer flux from (I-1) to (I)
!-----local variables:
Real*8 FTT(99) !parallel to QTT
Real*8 FXT(99),FXX(99),FYT(99),FZT(99),FXY(99),FXZ(99),FYY(99), &
& FZZ(99),FYZ(99) !parallel to Q__s
Real*8 EPS,EPSQ,EPS1,EPS1Q,EPS0, QTT4
REAL*8 FTTMX(99),FTTMN(99),FTTMX0(99),FTTMN0(99)
REAL*8 FMAX,FMIN, F0,F1,F2,F0X,FSUM
real*8 arg1,arg2
INTEGER I ! counter
!xx$$
save
!xx$$
!
IF(LIMTR.GT.0) THEN
IF(LIMTR.LE.2) THEN
!----LIMTR = 1:2
!-----restrain moments (only positive tracer transport) BEFORE transport
DO I=1,NQ
CALL QLIMIT(QTT(I),QXT(I),QXX(I), &
& QXY(I),QXZ(I),QYT(I),QYY(I),QZT(I),QZZ(I),QYZ(I),QM(I),LIMTR)
ENDDO
ELSE
!----LIMTR = 3:4
!---calculate initial MIN/MAX of tracer mix ratio in each box, limit x-terms
DO I=1,NQ
QTT4 = QTT(I)
arg1=-1.5*QTT4
arg2=QXY(I)
arg2=MAX(arg1,arg2)
QXY(I) = MIN(+1.5*QTT4,arg2)
arg1=-1.5*QTT4
arg2=QXZ(I)
arg2=MAX(arg1,arg2)
QXZ(I) = MIN(+1.5*QTT4,arg2)
! option QYZ(I) = MIN(+1.5*QTT4,MAX(-1.5*QTT4, QYZ(I)))
! QYT(I) = MIN(+1.5*QTT4,MAX(-1.5*QTT4, QYT(I)))
! QYY(I) = MIN(+1.5*QTT4,MAX(-1.5*QTT4, QYY(I)))
! QZT(I) = MIN(+1.5*QTT4,MAX(-1.5*QTT4, QZT(I)))
! QZZ(I) = MIN(+1.5*QTT4,MAX(-1.5*QTT4, QZZ(I)))
F0 = QTT(I)/QM(I)
F1 = ABS(QXT(I)) /QM(I)
IF(LIMTR.EQ.4) &
& F1 = ( ABS(QXT(I)) + ABS(QXY(I)) + ABS(QXZ(I)) )/QM(I)
F2 = QXX(I)/QM(I)
FTTMN0(I) = F0 - F1 + F2
FTTMX0(I) = F0 + F1 + F2
!---NB this extra test, OR the write above, saves a zero divide that probably
!--- results from the optimization of the F0X calc before the IF test.
IF(ABS(F2).GT.0.0) THEN
IF(3.*ABS(F2).GT.F1) THEN ! local extrema in interval
F0X = F0 - F1*F1/(6.*F2) - 0.5*F2
FTTMN0(I) = MIN(FTTMN0(I),F0X)
FTTMX0(I) = MAX(FTTMX0(I),F0X)
ENDIF
ENDIF
FTTMN(I) = FTTMN0(I)
FTTMX(I) = FTTMX0(I)
ENDDO
ENDIF
ENDIF
!---clear moments/fluxes between boxes (I-1) & (I):
DO I=1,NQ+1
Q0F(I) = 0.0
FTT(I) = 0.0
FXT(I) = 0.0
FXX(I) = 0.0
FYT(I) = 0.0
FYY(I) = 0.0
FZT(I) = 0.0
FZZ(I) = 0.0
FXY(I) = 0.0
FXZ(I) = 0.0
FYZ(I) = 0.0
ENDDO
!---Calculate transport for all interior faces: I = 2, NQ
!---QU(I) = mass FLUX across boundary: BOX(I-1) <-> BOX(I)
DO 24 I=2,NQ
IF(QU(I).EQ.0.) GOTO 24
!---QU(I) = 0 ==> do nothing if no mass flow across boundary
IF(QU(I).LT.0.) GOTO 22
!---QU(I) > 0 ==> FLUX FROM (I-1) -> (I)
!---FLUX limited by CURRENT mass in box (I-1) (cant get blood from a stone)
IF(QM(I-1).LE.0.0) GOTO 24
EPS = QU(I)/QM(I-1)
arg1 = EPS
EPS = MIN(arg1, 1.0)
EPSQ = EPS*EPS
EPS1 = 1.0 - EPS
EPS1Q = EPS1*EPS1
!---COMPUTE moments of tracer in the parcel transferring from (I-1) to (I)
FTT(I) = EPS*(QTT(I-1) + EPS1*(QXT(I-1) + (EPS1-EPS)*QXX(I-1)))
IF(LIMTR.GT.0) FTT(I) = MIN(MAX(FTT(I), 0.D0), QTT(I-1))
FXT(I) = EPSQ*(QXT(I-1)+3.0*EPS1*QXX(I-1))
FXX(I) = EPS*EPSQ*QXX(I-1)
FYT(I) = EPS*(QYT(I-1) + EPS1*QXY(I-1))
FZT(I) = EPS*(QZT(I-1) + EPS1*QXZ(I-1))
FXY(I) = EPSQ*QXY(I-1)
FXZ(I) = EPSQ*QXZ(I-1)
FYY(I) = EPS*QYY(I-1)
FZZ(I) = EPS*QZZ(I-1)
FYZ(I) = EPS*QYZ(I-1)
!---READJUST air mass and tracer moments remaining in box (I-1)
QM(I-1) = QM(I-1) - QU(I)
QTT(I-1) = QTT(I-1) - FTT(I)
QXT(I-1) = EPS1Q*(QXT(I-1) - 3.0*EPS*QXX(I-1))
QXX(I-1) = EPS1*EPS1Q*QXX(I-1)
QYT(I-1) = QYT(I-1) - FYT(I)
QYY(I-1) = QYY(I-1) - FYY(I)
QZT(I-1) = QZT(I-1) - FZT(I)
QZZ(I-1) = QZZ(I-1) - FZZ(I)
QXY(I-1) = EPS1Q*QXY(I-1)
QXZ(I-1) = EPS1Q*QXZ(I-1)
QYZ(I-1) = QYZ(I-1) - FYZ(I)
!---diagnostic flux across boundary:
Q0F(I) = +FTT(I)
GOTO 24
22 CONTINUE
!---QU(I) < 0 ==> FLUX FROM (I) -> (I-1=I-1)
IF(QM(I).LE.0.0) GOTO 24
EPS = -QU(I)/QM(I)
arg1=EPS
EPS = MIN(arg1, 1.0)
EPSQ = EPS*EPS
EPS1 = 1.0 - EPS
EPS1Q = EPS1*EPS1
!---COMPUTE moments of tracer in the parcel transferring from (I) to (I-1)
FTT(I) = EPS*(QTT(I) - EPS1*(QXT(I) - (EPS1-EPS)*QXX(I)))
IF(LIMTR.GT.0) FTT(I) = MIN(MAX(FTT(I), 0.D0), QTT(I))
FXT(I) = EPSQ*(QXT(I)-3.0*EPS1*QXX(I))
FXX(I) = EPS*EPSQ*QXX(I)
FYT(I) = EPS*(QYT(I) - EPS1*QXY(I))
FZT(I) = EPS*(QZT(I) - EPS1*QXZ(I))
FXY(I) = EPSQ*QXY(I)
FXZ(I) = EPSQ*QXZ(I)
FYY(I) = EPS*QYY(I)
FZZ(I) = EPS*QZZ(I)
FYZ(I) = EPS*QYZ(I)
!---READJUST air mass and tracer moments remaining in box (I) if not b.c.
QM(I) = QM(I) + QU(I)
QTT(I) = QTT(I) - FTT(I)
QXT(I) = EPS1Q*(QXT(I) + 3.0*EPS*QXX(I))
QXX(I) = EPS1*EPS1Q*QXX(I)
QYT(I) = QYT(I) - FYT(I)
QYY(I) = QYY(I) - FYY(I)
QZT(I) = QZT(I) - FZT(I)
QZZ(I) = QZZ(I) - FZZ(I)
QXY(I) = EPS1Q*QXY(I)
QXZ(I) = EPS1Q*QXZ(I)
QYZ(I) = QYZ(I) - FYZ(I)
!---diagnostic flux across boundary:
Q0F(I) = -FTT(I)
24 CONTINUE
!---end loop over internal faces
!---B.C. calculate transport for edge faces (boundary conditions): 1 & NQ+1
IF (LCYCLE) THEN
QU(NQ+1) = QU(1)
!---Connect (NQ) with (1) using air mass flux = QU(1), structure as above
IF(QU(1).EQ.0.) GOTO 27
IF(QU(1).LT.0.) GOTO 25
IF(QM(NQ).LE.0.0) GOTO 27
!---QU(1)>0 flow from (NQ) to (1)
EPS = QU(1)/QM(NQ)
arg1=eps
EPS = MIN(arg1, 1.0)
EPSQ = EPS*EPS
EPS1 = 1.0 - EPS
EPS1Q = EPS1*EPS1
FTT(1) = EPS*(QTT(NQ) + EPS1*(QXT(NQ) + (EPS1-EPS)*QXX(NQ)))
IF(LIMTR.GT.0) FTT(1) = MIN(MAX(FTT(1), 0.D0), QTT(NQ))
FXT(1) = EPSQ*(QXT(NQ)+3.0*EPS1*QXX(NQ))
FXX(1) = EPS*EPSQ*QXX(NQ)
FYT(1) = EPS*(QYT(NQ) + EPS1*QXY(NQ))
FZT(1) = EPS*(QZT(NQ) + EPS1*QXZ(NQ))
FXY(1) = EPSQ*QXY(NQ)
FXZ(1) = EPSQ*QXZ(NQ)
FYY(1) = EPS*QYY(NQ)
FZZ(1) = EPS*QZZ(NQ)
FYZ(1) = EPS*QYZ(NQ)
!---READJUST air mass and tracer moments remaining in box (NQ)
QM(NQ) = QM(NQ) - QU(1)
QTT(NQ) = QTT(NQ) - FTT(1)
QXT(NQ) = EPS1Q*(QXT(NQ) - 3.0*EPS*QXX(NQ))
QXX(NQ) = EPS1*EPS1Q*QXX(NQ)
QYT(NQ) = QYT(NQ) - FYT(1)
QYY(NQ) = QYY(NQ) - FYY(1)
QZT(NQ) = QZT(NQ) - FZT(1)
QZZ(NQ) = QZZ(NQ) - FZZ(1)
QXY(NQ) = EPS1Q*QXY(NQ)
QXZ(NQ) = EPS1Q*QXZ(NQ)
QYZ(NQ) = QYZ(NQ) - FYZ(1)
!---flux diagnostic
Q0F(1) = +FTT(1)
GOTO 26
25 CONTINUE
IF(QM(1).LE.0.0) GOTO 26
!---flow from (1) to (NQ+1)
EPS = -QU(1)/QM(1)
arg1=eps
EPS = MIN(arg1, 1.0)
EPSQ = EPS*EPS
EPS1 = 1.0 - EPS
EPS1Q = EPS1*EPS1
FTT(1) = EPS*(QTT(1) - EPS1*(QXT(1) - (EPS1-EPS)*QXX(1)))
IF(LIMTR.GT.0) FTT(1) = MIN(MAX(FTT(1), 0.D0), QTT(1))
FXT(1) = EPSQ*(QXT(1)-3.0*EPS1*QXX(1))
FXX(1) = EPS*EPSQ*QXX(1)
FYT(1) = EPS*(QYT(1) - EPS1*QXY(1))
FZT(1) = EPS*(QZT(1) - EPS1*QXZ(1))
FXY(1) = EPSQ*QXY(1)
FXZ(1) = EPSQ*QXZ(1)
FYY(1) = EPS*QYY(1)
FZZ(1) = EPS*QZZ(1)
FYZ(1) = EPS*QYZ(1)
!---READJUST air mass and tracer moments remaining in box (1)
QM(1) = QM(1) + QU(1)
QTT(1) = QTT(1) - FTT(1)
QXT(1) = EPS1Q*(QXT(1) + 3.0*EPS*QXX(1))
QXX(1) = EPS1*EPS1Q*QXX(1)
QYT(1) = QYT(1) - FYT(1)
QYY(1) = QYY(1) - FYY(1)
QZT(1) = QZT(1) - FZT(1)
QZZ(1) = QZZ(1) - FZZ(1)
QXY(1) = EPS1Q*QXY(1)
QXZ(1) = EPS1Q*QXZ(1)
QYZ(1) = QYZ(1) - FYZ(1)
!---flux diagnostic
Q0F(1) = -FTT(1)
!----equate (0)-(1) fluxes with (NQ)-(NQ+1) fluxes
26 CONTINUE
FTT(NQ+1) = FTT(1)
FXT(NQ+1) = FXT(1)
FXX(NQ+1) = FXX(1)
FYT(NQ+1) = FYT(1)
FYY(NQ+1) = FYY(1)
FZT(NQ+1) = FZT(1)
FZZ(NQ+1) = FZZ(1)
FXY(NQ+1) = FXY(1)
FXZ(NQ+1) = FXZ(1)
FYZ(NQ+1) = FYZ(1)
Q0F(NQ+1) = Q0F(1)
27 CONTINUE
ELSE
!---LCYCLE = .F. connect QBC1 with (1) using air mass flux = QU(1)
!--- and QBC2 with (NQ) using air mass flux = QU(NQ+1)
IF(QU(1).EQ.0.) GOTO 28
IF(QU(1).GT.0.) THEN
!--- bring in air with QBC1 mixing ratio, no moments
FTT(1) = QBC1*QU(1)
FXT(1) = 0.0
FXX(1) = 0.0
FYT(1) = 0.0
FZT(1) = 0.0
FXY(1) = 0.0
FXZ(1) = 0.0
FYY(1) = 0.0
FZZ(1) = 0.0
FYZ(1) = 0.0
ELSE
!--- take air out of (1) and throw away
IF(QM(1).GT.0.0) THEN
EPS = -QU(1)/QM(1)
arg1=EPS
EPS = MIN(arg1, 1.0)
EPSQ = EPS*EPS
EPS1 = 1.0 - EPS
EPS1Q = EPS1*EPS1
FTT(1) = EPS*(QTT(1) - EPS1*(QXT(1) - (EPS1-EPS)*QXX(1)))
IF(LIMTR.GT.0) FTT(1) = MIN(MAX(FTT(1), 0.D0), QTT(1))
! FXT(1) = EPSQ*(QXT(1)-3.0*EPS1*QXX(1))
! FXX(1) = EPS*EPSQ*QXX(1)
FYT(1) = EPS*(QYT(1) - EPS1*QXY(1))
FZT(1) = EPS*(QZT(1) - EPS1*QXZ(1))
FYY(1) = EPS*QYY(1)
FZZ(1) = EPS*QZZ(1)
! FXY(1) = EPSQ*QXY(1)
! FXZ(1) = EPSQ*QXZ(1)
FYZ(1) = EPS*QYZ(1)
QM(1) = QM(1) + QU(1)
QTT(1) = QTT(1) - FTT(1)
QXT(1) = EPS1Q*(QXT(1) + 3.0*EPS*QXX(1))
QXX(1) = EPS1*EPS1Q*QXX(1)
QYT(1) = QYT(1) - FYT(1)
QZT(1) = QZT(1) - FZT(1)
QYY(1) = QYY(1) - FYY(1)
QZZ(1) = QZZ(1) - FZZ(1)
QXY(1) = EPS1Q*QXY(1)
QXZ(1) = EPS1Q*QXZ(1)
QYZ(1) = QYZ(1) - FYZ(1)
ENDIF
ENDIF
28 CONTINUE
IF(QU(NQ+1).EQ.0.) GOTO 29
IF(QU(NQ+1).GT.0.) THEN
!--- take air out of (NQ) and throw away
IF(QM(NQ).GT.0.0) THEN
EPS = QU(NQ+1)/QM(NQ)
arg1=EPS
EPS = MIN(arg1, 1.0)
EPSQ = EPS*EPS
EPS1 = 1.0 - EPS
EPS1Q = EPS1*EPS1
FTT(NQ+1) = EPS*(QTT(NQ) + EPS1*(QXT(NQ) + (EPS1-EPS)*QXX(NQ)))
IF(LIMTR.GT.0) FTT(NQ+1) = MIN(MAX(FTT(NQ+1), 0.D0), QTT(NQ))
! FXT(NQ+1) = EPSQ*(QXT(NQ)+3.0*EPS1*QXX(NQ))
! FXX(NQ+1) = EPS*EPSQ*QXX(NQ)
FYT(NQ+1) = EPS*(QYT(NQ) + EPS1*QXY(NQ))
FZT(NQ+1) = EPS*(QZT(NQ) + EPS1*QXZ(NQ))
FYY(NQ+1) = EPS*QYY(NQ)
FZZ(NQ+1) = EPS*QZZ(NQ)
! FXY(NQ+1) = EPSQ*QXY(NQ)
! FXZ(NQ+1) = EPSQ*QXZ(NQ)
FYZ(NQ+1) = EPS*QYZ(NQ)
QM(NQ) = QM(NQ) - QU(NQ+1)
QTT(NQ) = QTT(NQ) - FTT(NQ+1)
QXT(NQ) = EPS1Q*(QXT(NQ) - 3.0*EPS*QXX(NQ))
QXX(NQ) = EPS1*EPS1Q*QXX(NQ)
QYT(NQ) = QYT(NQ) - FYT(NQ+1)
QYY(NQ) = QYY(NQ) - FYY(NQ+1)
QZT(NQ) = QZT(NQ) - FZT(NQ+1)
QZZ(NQ) = QZZ(NQ) - FZZ(NQ+1)
QXY(NQ) = EPS1Q*QXY(NQ)
QXZ(NQ) = EPS1Q*QXZ(NQ)
QYZ(NQ) = QYZ(NQ) - FYZ(NQ+1)
ENDIF
ELSE
!--- bring in air with QBC2 mixing ratio, no moments
FTT(NQ+1) = -QBC2*QU(NQ+1)
FXT(NQ+1) = 0.0
FXX(NQ+1) = 0.0
FYT(NQ+1) = 0.0
FZT(NQ+1) = 0.0
FXY(NQ+1) = 0.0
FXZ(NQ+1) = 0.0
FYY(NQ+1) = 0.0
FZZ(NQ+1) = 0.0
FYZ(NQ+1) = 0.0
ENDIF
29 CONTINUE
ENDIF
!---end of separating pieces (LCYCLE if/then/else)
!---NOW can put back the intermediate 'zones' at the interfaces of the grid
DO 44 I=2,NQ
IF(QU(I).EQ.0.) GOTO 44
!---QU(I) = 0 ==> do nothing if no mass flow across boundary
!---QU(I) > 0 ==> intermediate transfer box (I) goes into (I)
IF(QU(I).GT.0.) THEN
QM(I) = QM(I) + QU(I)
EPS = QU(I)/QM(I)
EPS1 = 1.0 - EPS
EPS0 = EPS*QTT(I) - EPS1*FTT(I)
!------N.B. Be careful about order of following statements:
QXX(I) = EPS*EPS*FXX(I) + EPS1*EPS1*QXX(I) &
& + 5.*(EPS*EPS1*(QXT(I)-FXT(I)) - (EPS1-EPS)*EPS0)
QXT(I) = EPS*FXT(I) + EPS1*QXT(I) + 3.0*EPS0
QTT(I) = QTT(I) + FTT(I)
QXY(I) = EPS*FXY(I) + EPS1*QXY(I) + 3.*(-EPS1*FYT(I)+EPS*QYT(I))
QXZ(I) = EPS*FXZ(I) + EPS1*QXZ(I) + 3.*(-EPS1*FZT(I)+EPS*QZT(I))
QYT(I) = QYT(I) + FYT(I)
QZT(I) = QZT(I) + FZT(I)
QYY(I) = QYY(I) + FYY(I)
QZZ(I) = QZZ(I) + FZZ(I)
QYZ(I) = QYZ(I) + FYZ(I)
ELSE
!---QU(I) < 0 ==> intermediate transfer box (I) goes into (I-1)
QM(I-1) = QM(I-1) - QU(I)
EPS = -QU(I)/QM(I-1)
EPS1 = 1.0 - EPS
EPS0 = -EPS*QTT(I-1) + EPS1*FTT(I)
QXX(I-1) = EPS*EPS*FXX(I) + EPS1*EPS1*QXX(I-1) &
& + 5.*(EPS*EPS1*(-QXT(I-1)+FXT(I)) + (EPS1-EPS)*EPS0)
QXT(I-1) = EPS*FXT(I) + EPS1*QXT(I-1) + 3.0*EPS0
QTT(I-1) = QTT(I-1) + FTT(I)
QXY(I-1)=EPS*FXY(I)+EPS1*QXY(I-1) + 3.*(EPS1*FYT(I)-EPS*QYT(I-1))
QXZ(I-1)=EPS*FXZ(I)+EPS1*QXZ(I-1) + 3.*(EPS1*FZT(I)-EPS*QZT(I-1))
QYT(I-1) = QYT(I-1) + FYT(I)
QZT(I-1) = QZT(I-1) + FZT(I)
QYY(I-1) = QYY(I-1) + FYY(I)
QZZ(I-1) = QZZ(I-1) + FZZ(I)
QYZ(I-1) = QYZ(I-1) + FYZ(I)
ENDIF
44 CONTINUE
!---end loop over internal faces, now fix b.c.
!-----do boundary at (1) if transport out of box (1)
IF(QU(1).GT.0.0) THEN
QM(1) = QM(1) + QU(1)
EPS = QU(1)/QM(1)
EPS1 = 1.0 - EPS
EPS0 = EPS*QTT(1) - EPS1*FTT(1)
!---N.B. Be careful about order of following statements:
QXX(1) = EPS*EPS*FXX(1) + EPS1*EPS1*QXX(1) &
& + 5.*(EPS*EPS1*(QXT(1)-FXT(1)) - (EPS1-EPS)*EPS0)
QXT(1) = EPS*FXT(1) + EPS1*QXT(1) + 3.0*EPS0
QTT(1) = QTT(1) + FTT(1)
QXY(1) = EPS*FXY(1) + EPS1*QXY(1) + 3.*(-EPS1*FYT(1)+EPS*QYT(1))
QXZ(1) = EPS*FXZ(1) + EPS1*QXZ(1) + 3.*(-EPS1*FZT(1)+EPS*QZT(1))
QYT(1) = QYT(1) + FYT(1)
QZT(1) = QZT(1) + FZT(1)
QYY(1) = QYY(1) + FYY(1)
QZZ(1) = QZZ(1) + FZZ(1)
QYZ(1) = QYZ(1) + FYZ(1)
ENDIF
!-----do boundary at NQ+1 if transport out of box (NQ)
IF(QU(NQ+1).LT.0.0) THEN
QM(NQ) = QM(NQ) - QU(1)
EPS = -QU(1)/QM(NQ)
EPS1 = 1.0 - EPS
EPS0 = -EPS*QTT(NQ) + EPS1*FTT(NQ+1)
!---N.B. Be careful about order of following statements:
QXX(NQ) = EPS*EPS*FXX(NQ+1) + EPS1*EPS1*QXX(NQ) &
& + 5.*(EPS*EPS1*(-QXT(NQ)+FXT(NQ+1)) + (EPS1-EPS)*EPS0)
QXT(NQ) = EPS*FXT(NQ+1) + EPS1*QXT(NQ) + 3.0*EPS0
QTT(NQ) = QTT(NQ) + FTT(NQ+1)
QXY(NQ)=EPS*FXY(NQ+1)+EPS1*QXY(NQ)+3.*(EPS1*FYT(NQ+1)-EPS*QYT(NQ))
QXZ(NQ)=EPS*FXZ(NQ+1)+EPS1*QXZ(NQ)+3.*(EPS1*FZT(NQ+1)-EPS*QZT(NQ))
QYT(NQ) = QYT(NQ) + FYT(NQ+1)
QZT(NQ) = QZT(NQ) + FZT(NQ+1)
QYY(NQ) = QYY(NQ) + FYY(NQ+1)
QZZ(NQ) = QZZ(NQ) + FZZ(NQ+1)
QYZ(NQ) = QYZ(NQ) + FYZ(NQ+1)
ENDIF
IF(LIMTR.GE.3) THEN !special section for LIMTR=3:4
!---NOW APPLY MIN/MAX TRANSPORT BOXES:
DO 52 I=2,NQ
IF(QU(I).EQ.0.) GOTO 52
IF(QU(I).GT.0.) THEN
!---QU(I) > 0 ==> (I-1) into (I)
FTTMX(I) = MAX(FTTMX0(I-1),FTTMX(I))
FTTMN(I) = MIN(FTTMN0(I-1),FTTMN(I))
ELSE
!---QU(I) < 0 ==> (I) goes into (I-1)
FTTMX(I-1) = MAX(FTTMX(I-1),FTTMX0(I))
FTTMN(I-1) = MIN(FTTMN(I-1),FTTMN0(I))
ENDIF
52 CONTINUE
!---B.C.
IF (LCYCLE) THEN !CYCLIC, connect (1) to (NQ)
IF(QU(1).GT.0.) THEN
!---QU(1)>0 flow from (NQ) into (1)
FTTMX(1) = MAX(FTTMX0(NQ),FTTMX(1))
FTTMN(1) = MIN(FTTMN0(NQ),FTTMN(1))
ELSE
!---flow from (1) into (NQ)
FTTMX(NQ) = MAX(FTTMX0(1),FTTMX(NQ))
FTTMN(NQ) = MIN(FTTMN0(1),FTTMN(NQ))
ENDIF
ELSE !connect QBC1 with (1), QBC2 with (NQ)
IF(QU(1).GT.0.) THEN
FTTMX(1) = MAX(QBC1,FTTMX(1))
FTTMN(1) = MIN(QBC1,FTTMN(1))
ENDIF
IF(QU(NQ+1).LT.0.) THEN
FTTMX(NQ) = MAX(QBC2,FTTMX(NQ))
FTTMN(NQ) = MIN(QBC2,FTTMN(NQ))
ENDIF
ENDIF
!---now place limits on X-moments: QXY & QXZ only changed if LIMTR=4
DO I=1,NQ
FMIN = FTTMN(I)*QM(I)
FMAX = FTTMX(I)*QM(I)
FSUM = QTT(I)
CALL FLIMIT(FSUM,QXT(I),QXY(I),QXZ(I),QXX(I),FMIN,FMAX,LIMTR)
ENDDO
ENDIF !end of LIMTR=3:4 section
RETURN
END
SUBROUTINE FLIMIT(F0,F1,F1A,F1B,F2,FMIN,FMAX,LIM) 1
!---given moments/mass = mixing ratio moments: F0,F1,F2, and FMIN<FMAX
!---adjust the moments to be monotonic, within limits
!---the cross-moments (eg, F1A=x-y, F1B=x-z) are also limited!
IMPLICIT NONE
REAL*8 F0,F1,F1A,F1B,F2, FMIN,FMAX
INTEGER LIM
REAL*8 FFMAX,FFMIN,EXCESS,SLACK,FX,FX0
real*8 arg1,arg2
!xx$$
save
!xx$$
!---FX is positive slope (true in some direction), rescale F1,F1A,F1B later
FX = ABS(F1)
IF(LIM.EQ.4) FX = ABS(F1) + ABS(F1A) + ABS(F1B)
!---check on reasonableness of limits:
IF(FMAX.LT.F0 .OR. FMIN.GT.F0 .OR. FX.LE.0.0) THEN
F1 = 0.0
F1A= 0.0
F1B= 0.0
F2 = 0.0
GOTO 99
ENDIF
FX0 = FX
!---limit slope by max-min range
arg1=FX
arg2=0.5*(FMAX-FMIN)
FX = MIN(arg1,arg2)
!---force F2 to be monotonic within range (0,1), x(extrema) = 1/2 - 1/6*(b/c)
IF(F2.GT.0.0) THEN
arg1=F2
arg2=FX/3.
F2 = MIN(arg1,arg2)
ELSE
arg1=F2
arg2=-FX/3.
F2 = MAX(arg1,arg2)
ENDIF
!---check if internal Max exceeds specified (FMAX)
FFMAX = F0+FX+F2
IF(FFMAX.GT.FMAX) THEN
!---EXCESS = amount beyond max allowed, SLACK = amount that can be taken
!--- from FX without violating monotonicity
EXCESS = FFMAX - FMAX
arg2=FX-3.*ABS(F2)
arg1=EXCESS
SLACK = MIN(arg1, arg2)
IF(SLACK.GT.0.0) THEN
FX = FX - SLACK
!---force F2 to be monotonic again
IF(F2.GT.0.0) THEN
arg1=F2
arg2=FX/3.
F2 = MIN(arg1,arg2)
ELSE
arg1=F2
arg2=-FX/3.
F2 = MAX(arg1,arg2)
ENDIF
ENDIF
!---now correct FX & F2 simultaneously for remaining excess:
EXCESS = EXCESS - SLACK
IF(EXCESS.GT.0.0) THEN
IF(F2.GT.0.0) THEN
F2 = F2 - 0.25*EXCESS
FX = FX - 0.75*EXCESS
ELSE
F2 = F2 + 0.5*EXCESS
FX = FX - 1.5*EXCESS
ENDIF
ENDIF
ENDIF ! end of fix if FFMAX>FMAX
!---now repeat for Min value:
FFMIN = F0-FX+F2
IF(FFMIN.LT.FMIN) THEN
EXCESS = FMIN - FFMIN
arg1=EXCESS
arg2=FX-3.*ABS(F2)
SLACK = MIN(arg1, arg2)
IF(SLACK.GT.0.0) THEN
FX = FX - SLACK
!---force F2 to be monotonic again
IF(F2.GT.0.0) THEN
arg1=F2
arg2=FX/3.
F2 = MIN(arg1,arg2)
ELSE
arg1=F2
arg2=-FX/3.
F2 = MAX(arg1,arg2)
ENDIF
ENDIF
!---now correct FX & F2 simultaneously for remaining excess:
EXCESS = EXCESS - SLACK
IF(EXCESS.GT.0.0) THEN
IF(F2.GT.0.0) THEN
F2 = F2 - 0.5*EXCESS
FX = FX - 1.5*EXCESS
ELSE
F2 = F2 + 0.25*EXCESS
FX = FX - 0.75*EXCESS
ENDIF
ENDIF
ENDIF ! end of fix if FFMAX>FMAX
!-----rescale the F1's
F1 = F1 *(FX/FX0)
IF(LIM.EQ.4) THEN
F1A = F1A *(FX/FX0)
F1B = F1B *(FX/FX0)
ENDIF
99 RETURN
END