C**** QUSDEF.f
MODULE QUSDEF 15
!@sum QUSDEF contains definitions for manipulating moments
!@auth Maxwell Kelley
integer, parameter ::
& nmom=9,
& mx =1, my =2, mz =3,
& mxx=4, myy=5, mzz=6,
& mxy=7, mzx=8, myz=9
C**** some useful vector subscripts:
C**** moments with a vertical component
integer, dimension(4), parameter ::
& zmoms=(/MZ,MZZ,MYZ,MZX/)
C**** moments with no vertical component
integer, dimension(5), parameter ::
& xymoms=(/MX,MY,MXX,MXY,MYY/)
C**** moments with a horizontal component
integer, dimension(7), parameter ::
& ihmoms=(/MX,MY,MXX,MYY,MXY,MYZ,MZX/)
C**** moments with no horizontal component
integer, dimension(2), parameter ::
& zomoms=(/MZ,MZZ/)
C**** x-x, x-y, x-z switches
integer, dimension(nmom), parameter ::
& xdir=(/mx,my,mz,mxx,myy,mzz,mxy,myz,mzx/)
integer, dimension(nmom), parameter ::
& ydir=(/my,mx,mz,myy,mxx,mzz,mxy,mzx,myz/)
integer, dimension(nmom), parameter ::
& zdir=(/mz,my,mx,mzz,myy,mxx,myz,mxy,mzx/)
!@dbparam prather_limits forces +ve sub-grid scale profiles (default=0)
integer :: prather_limits = 0.
END MODULE QUSDEF
subroutine adv1d(s,smom, f,fmom, mass,dm, nx,qlimit,stride,dir 8,2
* ,ierr,nerr)
!@sum adv1d implements the quadratic upstream scheme in one dimension
!@auth G. Russell, modified by Maxwell Kelley
c--------------------------------------------------------------
c adv1d advects tracers in x-direction using the qus
c the order of the moments in dir is: x,y,z,xx,yy,zz,xy,yz,zx
c--------------------------------------------------------------
use QUSDEF
, only : nmom,prather_limits
implicit none
!
!@var s mean tracer amount (kg or J)
!@var smom qus tracer moments (kg or J)
!@var f tracer flux (diagnostic output) (kg or J)
!@var fmom tracer moment flux (diagnostic output) (kg or J)
!@var mass mass field (kg)
!@var dm mass flux (kg)
!@var nx length of 1D vector
!@var qlimit true if negative tracer is to be avoided
!@var stride spacing in s array between elements of relevant 1D array
!@var dir direction switch (equals one of xdir ydir or zdir)
!@var ierr, nerr error codes
integer, intent(in) :: nx,stride
logical, intent(in) :: qlimit
REAL*8, dimension(nx) :: dm, f
REAL*8, dimension(nx*stride) :: s,mass
REAL*8, dimension(nmom,nx*stride) :: smom
REAL*8, dimension(nmom,nx) :: fmom
integer, dimension(nmom) :: dir
integer :: mx,my,mz,mxx,myy,mzz,mxy,myz,mzx
integer :: n,np1,nm1,nn,ns
integer,intent(out) :: ierr,nerr
REAL*8 :: fracm,frac1,bymnew,mnew,dm2,tmp
! qlimit variables
REAL*8 :: an, anm1, fn, fnm1, sn, sxn, sxxn
!
ierr=0 ; nerr=0
mx = dir(1)
my = dir(2)
mz = dir(3)
mxx = dir(4)
myy = dir(5)
mzz = dir(6)
mxy = dir(7)
myz = dir(8)
mzx = dir(9)
c-----------------------------------------------------------
! calculate tracer mass flux f
c-----------------------------------------------------------
n = nx
do np1=1,nx
if(dm(n).lt.0.) then ! air mass flux is negative
nn=np1
frac1=+1.
else ! air mass flux is positive
nn=n
frac1=-1.
endif
ns=1+(nn-1)*stride
fracm=dm(n)/mass(ns)
if(mass(ns).le.0.d0) fracm=0.d0
frac1=fracm+frac1
f(n)=fracm*( & (frac1+fracm)*smom(mxx,ns)))
! temporary storage of fracm in fx, to be used below
fmom(mx,n)=fracm
! temporary storage of frac1 in fxx, to be used below
fmom(mxx,n)=frac1
!
n = np1
enddo
if(qlimit) then
nm1 = nx
do n=1,nx
ns=1+(n-1)*stride
an = fmom(mx,n) ! reading fracm which was stored in fx
anm1 = fmom(mx,nm1)
fn = f(n)
fnm1 = f(nm1)
sn = s(ns)
sxn = smom(mx,ns)
sxxn = smom(mxx,ns)
call limitq(anm1,an,fnm1,fn,sn,sxn,sxxn,ierr)
if (ierr.gt.0) then
nerr=n
if (ierr.eq.2) return
end if
f(n) = fn
f(nm1) = fnm1
smom(mx,ns) = sxn
smom(mxx,ns) = sxxn
nm1 = n
enddo
endif
c--------------------------------------------------------------------
! calculate tracer fluxes of slopes and curvatures
c--------------------------------------------------------------------
n = nx
do np1=1,nx
if(dm(n).lt.0.) then ! air mass flux is negative
nn=np1
else ! air mass flux is positive
nn=n
endif
ns=1+(nn-1)*stride
! retrieving fracm, which was stored in fx
fracm=fmom(mx,n)
! retrieving frac1, which was stored in fxx
frac1=fmom(mxx,n)
!
fmom(mx,n)=dm(n)*(fracm*fracm*(smom(mx,ns)
& -3.*frac1*smom(mxx,ns))-3.*f(n))
fmom(mxx,n)=dm(n)*(dm(n)*fracm**3 *smom(mxx,ns)
& -5.*(dm(n)*f(n)+fmom(mx,n)))
! cross moments
fmom(my,n) = fracm*(smom(my,ns)-frac1*smom(mxy,ns))
fmom(mxy,n) = dm(n)*(fracm*fracm*smom(mxy,ns)-3.*fmom(my,n))
fmom(mz,n) = fracm*(smom(mz,ns)-frac1*smom(mzx,ns))
fmom(mzx,n) = dm(n)*(fracm*fracm*smom(mzx,ns)-3.*fmom(mz,n))
fmom(myy,n) = fracm*smom(myy,ns)
fmom(mzz,n) = fracm*smom(mzz,ns)
fmom(myz,n) = fracm*smom(myz,ns)
n = np1
enddo
c-------------------------------------------------------------------
c update tracer mass, moments of tracer mass, air mass distribution
c-------------------------------------------------------------------
nm1 = nx
do n=1,nx
ns=1+(n-1)*stride
tmp=mass(ns)+dm(nm1)
mnew=tmp-dm(n)
! mnew=mass(ns)+dm(nm1)-dm(n)
bymnew = 1./mnew
dm2=dm(nm1)+dm(n)
tmp=s(ns)+f(nm1)
s(ns)=tmp-f(n)
! s(ns)=s(ns)+f(nm1)-f(n)
!
smom(mx,ns)=(smom(mx,ns)*mass(ns)-3.*(-dm2*s(ns)
& +mass(ns)*( smom(mxx,ns) = (smom(mxx,ns)*mass(ns)*mass(ns)
& +2.5*s(ns)*(mass(ns)*mass(ns)-mnew*mnew-3.*dm2*dm2)
& +5.*(mass(ns)*(mass(ns)*( & -fmom(mx,n))+dm2*smom(mx,ns)*mnew)
& +(fmom(mxx,nm1)-fmom(mxx,n))) * (bymnew*bymnew)
! cross moments
smom(my,ns)=smom(my,ns)+fmom(my,nm1)-fmom(my,n)
smom(mxy,ns)=(smom(mxy,ns)*mass(ns)-3.*(-dm2*smom(my,ns) +
& mass(ns)*(fmom(my,nm1)+fmom(my,n))) +
& (fmom(mxy,nm1)-fmom(mxy,n)))*bymnew
smom(mz,ns)=smom(mz,ns)+fmom(mz,nm1)-fmom(mz,n)
smom(mzx,ns)=(smom(mzx,ns)*mass(ns)-3.*(-dm2*smom(mz,ns) +
& mass(ns)*(fmom(mz,nm1)+fmom(mz,n))) +
& (fmom(mzx,nm1)-fmom(mzx,n)))*bymnew
!
smom(myy,ns)=smom(myy,ns)+fmom(myy,nm1)-fmom(myy,n)
smom(mzz,ns)=smom(mzz,ns)+fmom(mzz,nm1)-fmom(mzz,n)
smom(myz,ns)=smom(myz,ns)+fmom(myz,nm1)-fmom(myz,n)
!
c------------------------------------------------------------------
mass(ns) = mnew
if(mass(ns).le.0.) then
s(ns)=0.
smom(:,ns)=0.
endif
if (qlimit .and. prather_limits.eq.1) then ! force Prather limits
smom(mx,ns)=min(1.5*s(ns),max(-1.5*s(ns),smom(mx,ns)))
smom(mxx,ns)=min(2.*s(ns)-abs(smom(mx,ns))/3.,max(abs(smom(mx
* ,ns))-s(ns),smom(mxx,ns)))
smom(mxy,ns)=min( smom(mzx,ns)=min( end if
c-----------------------------------------------------------------
nm1 = n
enddo
return
end subroutine adv1d
subroutine limitq(anm1,an,fnm1,fn,sn,sx,sxx,ierr) 1
!@sum limitq adjusts moments to maintain non-neg. tracer means/fluxes
!@auth G. Russell, modified by Maxwell Kelley
implicit none
REAL*8 :: anm1,an,fnm1,fn,sn,sx,sxx
c local variables
REAL*8 :: sl,sc,sr, frl,frl1, frr,frr1, gamma,g13ab,
& fr,fr1, fsign,su,sd
integer, intent(out) ::ierr
ierr=0
c****
c**** modify the tracer moments so that the tracer mass in each
c**** division is non-negative
c****
c**** no air leaving the box
if(anm1.ge.0. .and. an.le.0.) return
c**** air is leaving through both the left and right edges
if(anm1.lt.0. .and. an.gt.0.) then
sl = -fnm1
sr = +fn
sc = sn - sl
sc = sc - sr
c**** all three divisions are non-negative
if(sl.ge.0. .and. sr.ge.0. .and. sc.ge.0.) return
c**** first check for the cases when only one of the three is negative
frl = anm1
frl1 = frl+1.
frr = an
frr1 = frr-1.
if(sl.ge.0. .and. sr.ge.0.) then ! center division
gamma = 1.+(frl-frr)
g13ab = gamma*gamma - 1. + 3.*(frl+frr)**2
sxx = sxx - sc*10.*g13ab /
& (gamma*(12.*(frl+frr)**2 + 5.*g13ab*g13ab))
sx = sx + sc*12.*(frl+frr) /
& (gamma*(12.*(frl+frr)**2 + 5.*g13ab*g13ab))
sl = -frl*(sn-frl1*(sx-(frl+frl1)*sxx))
sr = sn-sl
else if(sr.ge.0.) then ! leftmost division
sxx = sxx + sl*(frl+frl1)/(frl*frl1*(.6d0+(frl+frl1)**2))
sx = sn/frl1 + (frl+frl1)*sxx
sr = frr*(sn-frr1*(sx-(frr+frr1)*sxx))
sl = 0.
else if(sl.ge.0.) then ! rightmost division
sxx = sxx - sr*(frr+frr1)/(frr*frr1*(.6d0+(frr+frr1)**2))
sx = sn/frr1 + (frr+frr1)*sxx
sl = -frl*(sn-frl1*(sx-(frl+frl1)*sxx))
sr = 0.
endif
sc = sn - sl
sc = sc - sr
c check for the cases where two of the three divisions are nonpositive
c these cases arise due to adjustments made when only one of the three
c divisions was negative
if(sl.le.0. .and. sr.le.0.) then ! two outer divisions
gamma = 1.+(frl-frr)
sxx = sn*(1.+gamma)/(2.*gamma*frl1*frr1)
sx = sn*(frl+frr)*(1.+2.*gamma)/(2.*gamma*frl1*frr1)
sl = 0.
sr = 0.
else if(sl.le.0. .and. sc.le.0.) then ! center/left divs
sxx = sn/(2.*frr*frl1)
sx = sn*(frl+frr+.5)/(frr*frl1)
sl = 0.
sr = sn
else if(sr.le.0. .and. sc.le.0.) then ! center/right divs
sxx = sn/(2.*frl*frr1)
sx = sn*(frl+frr-.5)/(frl*frr1)
sl = sn
sr = 0.
endif
fnm1 = -sl
fn = +sr
else
c**** air is leaving only through one edge
if(an.gt.0.) then ! right edge
fr=an
sd=fn
fsign=-1.
else ! left edge
fr=anm1
sd=-fnm1
fsign=1.
endif
if(abs(fr).gt.1.) then
c**** give warnings if fractional mass loss > 1
ierr=1
write(6,*) "limitq warning: abs(a)>1",fr,sd
write(6,*) "limitq input: anm1,an,fnm1,fn,sn,sx,sxx",anm1
* ,an,fnm1,fn,sn,sx,sxx
c**** only stop if net tracer mass is negative
if (sn+fnm1-fn.lt.0) then
ierr=2
write(6,*) "limitq error: new sn < 0",sn,sn+fnm1-fn
return
end if
end if
su = sn-sd
if(sd.ge.0. .and. su.ge.0.) return
fr1=fr+fsign
if(sd.lt.0.) then
c**** downstream division is negative
sxx = sxx +fsign*sd*(fr+fr1)/(fr*fr1*(.6d0+(fr+fr1)**2))
sx = sn/fr1 + (fr+fr1)*sxx
su = sn
else
c**** upstream division is negative
sxx = sxx -fsign*su*(fr+fr1)/(fr*fr1*(.6d0+(fr+fr1)**2))
sx = sn/fr + (fr+fr1)*sxx
su = 0.
endif
sd = sn - su
if(an.gt.0.) then
fn=sd
else
fnm1=-sd
endif
endif
return
end subroutine limitq
SUBROUTINE CTMIX (RM,RMOM,FMAIR,FMIX,FRAT) 2,1
!@sum CTMIX Cloud top mixing of tracer moments (incl. q,t) from CONDSE
!@auth Gary Russell, Jean Lerner, Gavin Schmidt
!@ver 1.0
!@var RM,RX,RY,RZ,RXX,RYY,RZZ,RXY,RYZ,RZX mean and moments of tracer
USE QUSDEF
IMPLICIT NONE
REAL*8, DIMENSION(2) :: RM
REAL*8, DIMENSION(NMOM,2) :: RMOM
REAL*8 FMIX !@var FMIX fraction of lower box mixed up
REAL*8 FRAT !@var FRAT fraction of upper box mixed down
REAL*8 FMAIR !@var FMAIR mass of air mixed
REAL*8 RTEMP !@var RTEMP dummy variable
RTEMP = RM(1)*(1.-FMIX)+FRAT*RM(2)
RM(2) = RM(2)*(1.-FRAT)+FMIX*RM(1)
RM(1) = RTEMP
C X
RTEMP = RMOM(MX ,1)*(1.-FMIX)+FRAT*RMOM(MX ,2)
RMOM(MX ,2) = RMOM(MX ,2)*(1.-FRAT)+FMIX*RMOM(MX ,1)
RMOM(MX ,1) = RTEMP
C Y
RTEMP = RMOM(MY ,1)*(1.-FMIX)+FRAT*RMOM(MY ,2)
RMOM(MY ,2) = RMOM(MY ,2)*(1.-FRAT)+FMIX*RMOM(MY ,1)
RMOM(MY ,1) = RTEMP
C XX
RTEMP = RMOM(MXX,1)*(1.-FMIX)+FRAT*RMOM(MXX,2)
RMOM(MXX,2) = RMOM(MXX,2)*(1.-FRAT)+FMIX*RMOM(MXX,1)
RMOM(MXX,1) = RTEMP
C YY
RTEMP = RMOM(MYY,1)*(1.-FMIX)+FRAT*RMOM(MYY,2)
RMOM(MYY,2) = RMOM(MYY,2)*(1.-FRAT)+FMIX*RMOM(MYY,1)
RMOM(MYY,1) = RTEMP
C XY
RTEMP = RMOM(MXY,1)*(1.-FMIX)+FRAT*RMOM(MXY,2)
RMOM(MXY,2) = RMOM(MXY,2)*(1.-FRAT)+FMIX*RMOM(MXY,1)
RMOM(MXY,1) = RTEMP
C Z MOMENTS
RMOM(ZMOMS,:) = RMOM(ZMOMS,:)*(1.-FMAIR)
END SUBROUTINE CTMIX