C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.10 2004/07/16 20:08:08 molod Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
subroutine moistio (ndmoist,istrip,npcs,
. lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup,
. pz,plz,plze,dpres,pkht,pkl,tz,qz,bi,bj,ntracer,ptracer,
. qqz,dumoist,dvmoist,dtmoist,dqmoist,
. im,jm,lm,ptop,
. iras,rainlsp,rainconv,snowfall,
. nswcld,cldtot_sw,cldras_sw,cldlsp_sw,nswlz,swlz,
. nlwcld,cldtot_lw,cldras_lw,cldlsp_lw,nlwlz,lwlz,
. lpnt,myid)
implicit none
#ifdef ALLOW_DIAGNOSTICS
#include "SIZE.h"
#include "diagnostics_SIZE.h"
#include "diagnostics.h"
#endif
c Input Variables
c ---------------
integer im,jm,lm
integer ndmoist,istrip,npcs
integer bi,bj,ntracer,ptracer
integer lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup
real pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
real pkht(im,jm,lm+1),pkl(im,jm,lm)
real tz(im,jm,lm),qz(im,jm,lm,ntracer)
real qqz(im,jm,lm)
real dumoist(im,jm,lm),dvmoist(im,jm,lm)
real dtmoist(im,jm,lm),dqmoist(im,jm,lm,ntracer)
real ptop
integer iras
real rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
integer nswcld,nswlz
real cldlsp_sw(im,jm,lm),cldras_sw(im,jm,lm)
real cldtot_sw(im,jm,lm),swlz(im,jm,lm)
integer nlwcld,nlwlz
real cldlsp_lw(im,jm,lm),cldras_lw(im,jm,lm)
real cldtot_lw(im,jm,lm),lwlz(im,jm,lm)
logical lpnt
integer myid
c Local Variables
c ---------------
integer ncrnd,nsecf
real fracqq, dum
integer snowcrit
parameter (fracqq = 0.1)
real cldsr(im,jm,lm)
real srcld(istrip,lm)
real plev
real cldnow,cldlsp_mem,cldlsp,cldras_mem,cldras
real watnow,watmin,cldmin
real cldprs(im,jm),cldtmp(im,jm)
real cldhi (im,jm),cldlow(im,jm)
real cldmid(im,jm),totcld(im,jm)
real CLDLS(im,jm,lm) , CPEN(im,jm,lm)
real tmpimjm(im,jm)
real lsp_new(im,jm)
real conv_new(im,jm)
real snow_new(im,jm)
real qqcolmin(im,jm)
real qqcolmax(im,jm)
integer levpbl(im,jm)
c Gathered Arrays for Variable Cloud Base
c ---------------------------------------
real raincgath(im*jm)
real pigather(im*jm)
real thgather(im*jm,lm)
real shgather(im*jm,lm)
real pkzgather(im*jm,lm)
real pkegather(im*jm,lm+1)
real plzgather(im*jm,lm)
real plegather(im*jm,lm+1)
real dpgather(im*jm,lm)
real tmpgather(im*jm,lm)
real deltgather(im*jm,lm)
real delqgather(im*jm,lm)
real ugather(im*jm,lm,ntracer)
real delugather(im*jm,lm,ntracer)
real deltrnev(im*jm,lm)
real delqrnev(im*jm,lm)
integer nindeces(lm)
integer pblindex(im*jm)
integer levgather(im*jm)
c Stripped Arrays
c ---------------
real saveth (istrip,lm)
real saveq (istrip,lm)
real saveu (istrip,lm,ntracer)
real usubcl (istrip, ntracer)
real ple(istrip,lm+1)
real dp(istrip,lm)
real TL(ISTRIP,lm) , SHL(ISTRIP,lm)
real PL(ISTRIP,lm) , PLK(ISTRIP,lm)
real PLKE(ISTRIP,lm+1)
real TH(ISTRIP,lm) ,CVTH(ISTRIP,lm)
real CVQ(ISTRIP,lm)
real UL(ISTRIP,lm,ntracer)
real cvu(istrip,lm,ntracer)
real CLMAXO(ISTRIP,lm),CLBOTH(ISTRIP,lm)
real CLSBTH(ISTRIP,lm)
real TMP1(ISTRIP,lm), TMP2(ISTRIP,lm)
real TMP3(ISTRIP,lm), TMP4(ISTRIP,lm+1)
real TMP5(ISTRIP,lm+1)
integer ITMP1(ISTRIP,lm), ITMP2(ISTRIP,lm)
real PRECIP(ISTRIP), PCNET(ISTRIP)
real SP(ISTRIP), PREP(ISTRIP)
real PCPEN (ISTRIP,lm)
integer pbl(istrip),depths(lm)
real cldlz(istrip,lm), cldwater(im,jm,lm)
real rhfrac(istrip), rhmin, pup, ppbl, rhcrit(istrip,lm)
real offset, alpha, rasmax
logical first
logical lras
real clfrac (istrip,lm)
real cldmas (istrip,lm)
real detrain(istrip,lm)
real psubcld (istrip), psubcldg (im,jm)
real psubcld_cnt(istrip), psubcldgc(im,jm)
real rnd(lm/2)
DATA FIRST /.TRUE./
integer imstp,nsubcl,nlras
integer i,j,iloop,index,l,nn,num,numdeps,nt
real tmstp,tminv,sday,grav,alhl,cp,elocp,gamfac
real rkappa,p0kappa,p0kinv,ptopkap,pcheck
real tice,getcon,pi
C **********************************************************************
C **** INITIALIZATION ****
C **********************************************************************
IMSTP = nsecf(NDMOIST)
TMSTP = FLOAT(IMSTP)
TMINV = 1. / TMSTP
C Minimum Large-Scale Cloud Fraction at rhcrit
alpha = 0.80
C Difference in fraction between SR and LS Threshold
offset = 0.10
C Large-Scale Relative Humidity Threshold in PBL
rhmin = 0.90
C Maximum Cloud Fraction associated with RAS
rasmax = 1.00
nn = 3*3600.0/tmstp + 1
C Threshold for Cloud Fraction Memory
cldmin = rasmax*(1.0-tmstp/3600.)**nn
C Threshold for Cloud Liquid Water Memory
watmin = 1.0e-8
SDAY = GETCON('SDAY')
GRAV = GETCON('GRAVITY')
ALHL = GETCON('LATENT HEAT COND')
CP = GETCON('CP')
ELOCP = GETCON('LATENT HEAT COND') / GETCON('CP')
GAMFAC = GETCON('LATENT HEAT COND') * GETCON('EPS') * ELOCP
. / GETCON('RGAS')
RKAPPA = GETCON('KAPPA')
P0KAPPA = 1000.0**RKAPPA
P0KINV = 1. / P0KAPPA
PTOPKAP = PTOP**RKAPPA
tice = getcon('FREEZING-POINT')
PI = 4.*atan(1.)
c Determine Total number of Random Clouds to Check
c ---------------------------------------------
ncrnd = (lm-nltop+1)/2
if(first .and. myid.eq.0) then
print *
print *,'Top Level Allowed for Convection : ',nltop
print *,' Highest Sub-Cloud Level: ',nsubmax
print *,' Lowest Sub-Cloud Level: ',nsubmin
print *,' Total Number of Random Clouds: ',ncrnd
print *
first = .false.
endif
c And now find PBL depth - the level where qq = fracqq * qq at surface
c --------------------------------------------------------------------
do j = 1,jm
do i = 1,im
qqcolmin(i,j) = qqz(i,j,lm)*fracqq
qqcolmax(i,j) = qqz(i,j,lm)
levpbl(i,j) = lm+1
enddo
enddo
DO L = lm-1,1,-1
DO j = 1,jm
DO i = 1,im
IF((qqz(i,j,l).gt.qqcolmax(i,j))
1 .and.(levpbl(i,j).eq.lm+1))then
qqcolmax(i,j) = qqz(i,j,l)
qqcolmin(i,j) = fracqq*qqcolmax(i,j)
endif
if((qqz(i,j,l).lt.qqcolmin(i,j))
1 .and.(levpbl(i,j).eq.lm+1))levpbl(i,j)=L+1
enddo
enddo
enddo
do j = 1,jm
do i = 1,im
if(levpbl(i,j).gt.nsubmin) levpbl(i,j) = nsubmin
if(levpbl(i,j).lt.nsubmax) levpbl(i,j) = nsubmax
enddo
enddo
c Set up the array of indeces of subcloud levels for the gathering
c ----------------------------------------------------------------
index = 0
do L = nsubmin,nltop,-1
do j = 1,jm
do i = 1,im
if(levpbl(i,j).eq.L) then
index = index + 1
pblindex(index) = (j-1)*im + i
endif
enddo
enddo
enddo
do index = 1,im*jm
levgather(index) = levpbl(pblindex(index),1)
pigather(index) = pz(pblindex(index),1)
pkegather(index,lm+1) = pkht(pblindex(index),1,lm+1)
plegather(index,lm+1) = plze(pblindex(index),1,lm+1)
enddo
do L = 1,lm
do index = 1,im*jm
thgather(index,L) = tz(pblindex(index),1,L)
shgather(index,L) = qz(pblindex(index),1,L,1)
pkegather(index,L) = pkht(pblindex(index),1,L)
pkzgather(index,L) = pkl(pblindex(index),1,L)
plegather(index,L) = plze(pblindex(index),1,L)
plzgather(index,L) = plz(pblindex(index),1,L)
dpgather(index,L) = dpres(pblindex(index),1,L)
enddo
enddo
do nt = 1,ntracer-ptracer
do L = 1,lm
do index = 1,im*jm
ugather(index,L,nt) = qz(pblindex(index),1,L,nt+ptracer)
enddo
enddo
enddo
c bump the counter for number of calls to convection
c --------------------------------------------------
iras = iras + 1
if( iras.ge.1e9 ) iras = 1
c select the 'random' cloud detrainment levels for RAS
c ----------------------------------------------------
call rndcloud(iras,ncrnd,rnd,myid)
do l=1,lm
do j=1,jm
do i=1,im
dtmoist(i,j,l) = 0.
do nt = 1,ntracer
dqmoist(i,j,l,nt) = 0.
enddo
enddo
enddo
enddo
C***********************************************************************
C **** LOOP OVER NPCS PEICES ****
C **********************************************************************
DO 1000 NN = 1,NPCS
C **********************************************************************
C **** VARIABLE INITIALIZATION ****
C **********************************************************************
CALL STRIP ( pigather, SP ,im*jm,ISTRIP,1 ,NN )
CALL STRIP ( pkzgather, PLK ,im*jm,ISTRIP,lm,NN )
CALL STRIP ( pkegather, PLKE ,im*jm,ISTRIP,lm+1,NN )
CALL STRIP ( plzgather, PL ,im*jm,ISTRIP,lm,NN )
CALL STRIP ( plegather, PLE ,im*jm,ISTRIP,lm+1,NN )
CALL STRIP ( dpgather, dp ,im*jm,ISTRIP,lm,NN )
CALL STRIP ( thgather, TH ,im*jm,ISTRIP,lm,NN )
CALL STRIP ( shgather, SHL ,im*jm,ISTRIP,lm,NN )
CALL STRINT( levgather, pbl ,im*jm,ISTRIP,1 ,NN )
do nt = 1,ntracer-ptracer
call strip ( ugather(1,1,nt), ul(1,1,nt),im*jm,istrip,lm,nn )
enddo
C **********************************************************************
C **** SETUP FOR RAS CUMULUS PARAMETERIZATION ****
C **********************************************************************
DO L = 1,lm
DO I = 1,ISTRIP
TH(I,L) = TH(I,L) * P0KAPPA
CLMAXO(I,L) = 0.
CLBOTH(I,L) = 0.
cldmas(I,L) = 0.
detrain(I,L) = 0.
ENDDO
ENDDO
do L = 1,lm
depths(L) = 0
enddo
numdeps = 0
do L = nsubmin,nltop,-1
nindeces(L) = 0
do i = 1,istrip
if(pbl(i).eq.L) nindeces(L) = nindeces(L) + 1
enddo
if(nindeces(L).gt.0) then
numdeps = numdeps + 1
depths(numdeps) = L
endif
enddo
C Initiate a do-loop around RAS for the number of different
C sub-cloud layer depths in this strip
C --If all subcloud depths are the same, execute loop once
C Otherwise loop over different subcloud layer depths
num = 1
DO iloop = 1,numdeps
nsubcl = depths(iloop)
c Compute sub-cloud values for Temperature and Spec.Hum.
c ------------------------------------------------------
DO 600 I=num,num+nindeces(nsubcl)-1
TMP1(I,2) = 0.
TMP1(I,3) = 0.
600 CONTINUE
NLRAS = NSUBCL - NLTOP + 1
DO 601 L=NSUBCL,lm
DO 602 I=num,num+nindeces(nsubcl)-1
TMP1(I,2) = TMP1(I,2) + (PLE(I,L+1)-PLE(I,L))*TH (I,L)/sp(i)
TMP1(I,3) = TMP1(I,3) + (PLE(I,L+1)-PLE(I,L))*SHL(I,L)/sp(i)
602 CONTINUE
601 CONTINUE
DO 603 I=num,num+nindeces(nsubcl)-1
TMP1(I,4) = 1. / ( (PLE(I,lm+1)-PLE(I,NSUBCL))/sp(I) )
TH(I,NSUBCL) = TMP1(I,2)*TMP1(I,4)
SHL(I,NSUBCL) = TMP1(I,3)*TMP1(I,4)
603 CONTINUE
c Save initial value of tracers and compute sub-cloud value
c ---------------------------------------------------------
DO NT = 1,ntracer-ptracer
do L = 1,lm
do i = num,num+nindeces(nsubcl)-1
saveu(i,L,nt) = ul(i,L,nt)
enddo
enddo
DO I=num,num+nindeces(nsubcl)-1
TMP1(I,2) = 0.
ENDDO
DO L=NSUBCL,lm
DO I=num,num+nindeces(nsubcl)-1
TMP1(I,2) = TMP1(I,2)+(PLE(I,L+1)-PLE(I,L))*UL(I,L,NT)/sp(i)
ENDDO
ENDDO
DO I=num,num+nindeces(nsubcl)-1
UL(I,NSUBCL,NT) = TMP1(I,2)*TMP1(I,4)
usubcl(i,nt) = ul(i,nsubcl,nt)
ENDDO
ENDDO
c Compute Pressure Arrays for RAS
c -------------------------------
DO 111 L=1,lm
DO 112 I=num,num+nindeces(nsubcl)-1
TMP4(I,L) = PLE(I,L)
112 CONTINUE
111 CONTINUE
DO I=num,num+nindeces(nsubcl)-1
TMP5(I,1) = PTOPKAP / P0KAPPA
ENDDO
DO L=2,lm
DO I=num,num+nindeces(nsubcl)-1
TMP5(I,L) = PLKE(I,L)*P0KINV
ENDDO
ENDDO
DO I=num,num+nindeces(nsubcl)-1
TMP4(I,lm+1) = PLE (I,lm+1)
TMP5(I,lm+1) = PLKE(I,lm+1)*P0KINV
ENDDO
DO 113 I=num,num+nindeces(nsubcl)-1
TMP4(I,NSUBCL+1) = PLE (I,lm+1)
TMP5(I,NSUBCL+1) = PLKE(I,lm+1)*P0KINV
113 CONTINUE
do i=num,num+nindeces(nsubcl)-1
C Temperature at top of sub-cloud layer
tmp2(i,1) = TH(i,NSUBCL) * PLKE(i,NSUBCL)/P0KAPPA
C Pressure at top of sub-cloud layer
tmp2(i,2) = tmp4(i,nsubcl)
enddo
C CHANGED THIS: no RH requirement for RAS
c call vqsat ( tmp2(num,1),tmp2(num,2),tmp2(num,3),
c . dum,.false.,nindeces(nsubcl) )
c do i=num,num+nindeces(nsubcl)-1
c rh = SHL(I,NSUBCL) / tmp2(i,3)
c if (rh .le. 0.85) then
c rhfrac(i) = 0.
c else if (rh .ge. 0.95) then
c rhfrac(i) = 1.
c else
c rhfrac(i) = (rh-0.85)*10.
c endif
c enddo
do i=num,num+nindeces(nsubcl)-1
rhfrac(i) = 1.
enddo
C Compute RH threshold for Large-scale condensation
C Used in Slingo-Ritter clouds as well - define offset between SR and LS
C Top level of atan func above this rh_threshold = rhmin
pup = 600.
do i=num,num+nindeces(nsubcl)-1
do L = nsubcl, lm
rhcrit(i,L) = 1.
enddo
do L = 1, nsubcl-1
pcheck = pl(i,L)
if (pcheck .le. pup) then
rhcrit(i,L) = rhmin
else
ppbl = pl(i,nsubcl)
rhcrit(i,L) = rhmin + (1.-rhmin)/(19.) *
. ((atan( (2.*(pcheck-pup)/(ppbl-pup)-1.) *
. tan(20.*pi/21.-0.5*pi) )
. + 0.5*pi) * 21./pi - 1.)
endif
enddo
enddo
c Save Initial Values of Temperature and Specific Humidity
c --------------------------------------------------------
do L = 1,lm
do i = num,num+nindeces(nsubcl)-1
saveth(i,L) = th (i,L)
saveq (i,L) = shl(i,L)
PCPEN (i,L) = 0.
CLFRAC(i,L) = 0.
enddo
enddo
CALL RAS ( NN,istrip,nindeces(nsubcl),NLRAS,NLTOP,lm,TMSTP
1, UL(num,1,1),ntracer-ptracer,TH(num,NLTOP),SHL(num,NLTOP)
2, TMP4(num,NLTOP), TMP5(num,NLTOP),rnd, ncrnd, PCPEN(num,NLTOP)
3, CLBOTH(num,NLTOP), CLFRAC(num,NLTOP)
4, cldmas(num,nltop), detrain(num,nltop)
8, cp,grav,rkappa,alhl,rhfrac(num),rasmax )
c Compute Diagnostic CLDMAS in RAS Subcloud Layers
c ------------------------------------------------
do L=nsubcl,lm
do I=num,num+nindeces(nsubcl)-1
dum = dp(i,L)/(ple(i,lm+1)-ple(i,nsubcl))
cldmas(i,L) = cldmas(i,L-1) - dum*cldmas(i,nsubcl-1)
enddo
enddo
c Update Theta and Moisture due to RAS
c ------------------------------------
DO L=1,nsubcl
DO I=num,num+nindeces(nsubcl)-1
CVTH(I,L) = (TH (I,L) - saveth(i,l))
CVQ (I,L) = (SHL(I,L) - saveq (i,l))
ENDDO
ENDDO
DO L=nsubcl+1,lm
DO I=num,num+nindeces(nsubcl)-1
CVTH(I,L) = cvth(i,nsubcl)
CVQ (I,L) = cvq (i,nsubcl)
ENDDO
ENDDO
DO L=nsubcl+1,lm
DO I=num,num+nindeces(nsubcl)-1
TH (I,L) = saveth(i,l) + cvth(i,l)
SHL(I,L) = saveq (i,l) + cvq (i,l)
ENDDO
ENDDO
DO L=1,lm
DO I=num,num+nindeces(nsubcl)-1
CVTH(I,L) = CVTH(I,L) *P0KINV*SP(I)*tminv
CVQ (I,L) = CVQ (I,L) *SP(I)*tminv
ENDDO
ENDDO
c Compute Tracer Tendency due to RAS
c ----------------------------------
do nt = 1,ntracer-ptracer
DO L=1,nsubcl-1
DO I=num,num+nindeces(nsubcl)-1
CVU(I,L,nt) = ( UL(I,L,nt)-saveu(i,l,nt) )*sp(i)*tminv
ENDDO
ENDDO
DO L=nsubcl,lm
DO I=num,num+nindeces(nsubcl)-1
if( usubcl(i,nt).ne.0.0 ) then
cvu(i,L,nt) = ( ul(i,nsubcl,nt)-usubcl(i,nt) ) *
. ( saveu(i,L,nt)/usubcl(i,nt) )*sp(i)*tminv
else
cvu(i,L,nt) = 0.0
endif
ENDDO
ENDDO
enddo
c Compute Diagnostic PSUBCLD (Subcloud Layer Pressure)
c ----------------------------------------------------
do i=num,num+nindeces(nsubcl)-1
lras = .false.
do L=nltop,nsubcl
if( cvq(i,L).ne.0.0 ) lras = .true.
enddo
psubcld (i) = 0.0
psubcld_cnt(i) = 0.0
if( lras ) then
psubcld (i) = sp(i)+ptop-ple(i,nsubcl)
psubcld_cnt(i) = 1.0
endif
enddo
C End of subcloud layer depth loop (iloop)
num = num+nindeces(nsubcl)
ENDDO
C **********************************************************************
C **** TENDENCY UPDATES ****
C **** (Keep 'Gathered' tendencies in 'gather' arrays now) ****
C **********************************************************************
call paste( CVTH,deltgather,istrip,im*jm,lm,NN )
call paste( CVQ,delqgather,istrip,im*jm,lm,NN )
do nt = 1,ntracer-ptracer
call paste( CVU(1,1,nt),delugather(1,1,nt),istrip,im*jm,lm,NN )
enddo
C **********************************************************************
C And now paste some arrays for filling diagnostics
C (use pkegather to hold detrainment and tmpgather for cloud mass flux)
C **********************************************************************
if(icldmas .gt.0) call paste( cldmas,tmpgather,istrip,im*jm,lm,NN)
if(idtrain .gt.0) call paste(detrain,pkegather,istrip,im*jm,lm,NN)
if(ipsubcld.gt.0) then
call paste(psubcld ,psubcldg ,istrip,im*jm,1,NN)
call paste(psubcld_cnt,psubcldgc,istrip,im*jm,1,NN)
endif
C *********************************************************************
C **** RE-EVAPORATION OF PENETRATING CONVECTIVE RAIN ****
C *********************************************************************
CALL STRIP ( thgather,TH ,im*jm,ISTRIP,lm,NN)
CALL STRIP ( shgather,SHL,im*jm,ISTRIP,lm,NN)
DO L=1,lm
DO I=1,ISTRIP
TH(I,L) = TH(I,L) + CVTH(I,L)*tmstp/SP(I)
SHL(I,L) = SHL(I,L) + CVQ(I,L)*tmstp/SP(I)
TL(I,L) = TH(I,L)*PLK(I,L)
saveth(I,L) = th(I,L)
saveq (I,L) = SHL(I,L)
ENDDO
ENDDO
CALL RNEVP (NN,ISTRIP,lm,TL,SHL,PCPEN,PL,CLFRAC,SP,DP,PLKE,
. PLK,TH,TMP1,TMP2,TMP3,ITMP1,ITMP2,PCNET,PRECIP,
. CLSBTH,TMSTP,1.,cp,grav,alhl,gamfac,cldlz,rhcrit,offset,alpha)
C **********************************************************************
C **** TENDENCY UPDATES ****
C **********************************************************************
DO L=1,lm
DO I =1,ISTRIP
TMP1(I,L) = sp(i) * (SHL(I,L)-saveq(I,L)) * tminv
ENDDO
CALL PSTBMP(TMP1(1,L),delqgather(1,L),ISTRIP,im*jm,1,NN)
DO I =1,ISTRIP
TMP1(I,L) = sp(i) * ((TL(I,L)/PLK(I,L))-saveth(i,l)) * tminv
ENDDO
CALL PSTBMP(TMP1(1,L),deltgather(1,L),ISTRIP,im*jm,1,NN)
C Paste rain evap tendencies into arrays for diagnostic output
c ------------------------------------------------------------
if(idtls.gt.0)then
DO I =1,ISTRIP
TMP1(I,L) = ((TL(I,L)/PLK(I,L))-saveth(i,l))*plk(i,l)*sday*tminv
ENDDO
call paste(tmp1(1,L),deltrnev(1,L),istrip,im*jm,1,NN)
endif
if(idqls.gt.0)then
DO I =1,ISTRIP
TMP1(I,L) = (SHL(I,L)-saveq(I,L)) * 1000. * sday * tminv
ENDDO
call paste(tmp1(1,L),delqrnev(1,L),istrip,im*jm,1,NN)
endif
ENDDO
C *********************************************************************
C Add Non-Precipitating Clouds where the relative
C humidity is less than 100%
C Apply Cloud Top Entrainment Instability
C *********************************************************************
do L=1,lm
do i=1,istrip
srcld(i,L) = -clsbth(i,L)
enddo
enddo
call srclouds (saveth,saveq,plk,pl,plke,clsbth,cldlz,istrip,lm,
. rhcrit,offset,alpha)
do L=1,lm
do i=1,istrip
srcld(i,L) = srcld(i,L)+clsbth(i,L)
enddo
enddo
C *********************************************************************
C **** PASTE CLOUD AMOUNTS ****
C *********************************************************************
call paste ( srcld, cldsr,istrip,im*jm,lm,nn )
call paste ( cldlz,cldwater,istrip,im*jm,lm,nn )
call paste ( clsbth, cldls,istrip,im*jm,lm,nn )
call paste ( clboth, cpen ,istrip,im*jm,lm,nn )
c compute Total Accumulated Precip for Landsurface Model
c ------------------------------------------------------
do i = 1,istrip
C Initialize Rainlsp, Rainconv and Snowfall
tmp1(i,1) = 0.0
tmp1(i,2) = 0.0
tmp1(i,3) = 0.0
enddo
do i = 1,istrip
prep(i) = PRECIP(I) + PCNET(I)
tmp1(i,1) = PRECIP(I)
tmp1(i,2) = pcnet(i)
enddo
c
c check whether there is snow
c-------------------------------------------------------
c snow algorthm:
c if temperature profile from the surface level to 700 mb
c uniformaly c below zero, then precipitation (total) is
c snowfall. Else there is no snow.
c-------------------------------------------------------
do i = 1,istrip
snowcrit=0
do l=lup,lm
if (saveth(i,l)*plk(i,l).le. tice ) then
snowcrit=snowcrit+1
endif
enddo
if (snowcrit .eq. (lm-lup+1)) then
tmp1(i,3) = prep(i)
tmp1(i,1)=0.0
tmp1(i,2)=0.0
endif
enddo
CALL paste (tmp1(1,1), lsp_new,ISTRIP,im*jm,1,NN)
CALL paste (tmp1(1,2),conv_new,ISTRIP,im*jm,1,NN)
CALL paste (tmp1(1,3),snow_new,ISTRIP,im*jm,1,NN)
if(iprecon.gt.0) then
CALL paste (pcnet,raincgath,ISTRIP,im*jm,1,NN)
endif
C *********************************************************************
C **** End Major Stripped Region ****
C *********************************************************************
1000 CONTINUE
C Large Scale Rainfall, Conv rain, and snowfall
c ---------------------------------------------
call back2grd ( lsp_new,pblindex, lsp_new,im*jm)
call back2grd (conv_new,pblindex,conv_new,im*jm)
call back2grd (snow_new,pblindex,snow_new,im*jm)
if(iprecon.gt.0) then
call back2grd (raincgath,pblindex,raincgath,im*jm)
endif
c Subcloud Layer Pressure
c -----------------------
if(ipsubcld.gt.0) then
call back2grd (psubcldg ,pblindex,psubcldg ,im*jm)
call back2grd (psubcldgc,pblindex,psubcldgc,im*jm)
endif
do L = 1,lm
C Delta theta,q, convective, max and ls clouds
c --------------------------------------------
call back2grd (deltgather(1,L),pblindex, dtmoist(1,1,L) ,im*jm)
call back2grd (delqgather(1,L),pblindex, dqmoist(1,1,L,1),im*jm)
call back2grd ( cpen(1,1,L),pblindex, cpen(1,1,L) ,im*jm)
call back2grd ( cldls(1,1,L),pblindex, cldls(1,1,L) ,im*jm)
call back2grd (cldwater(1,1,L),pblindex,cldwater(1,1,L) ,im*jm)
call back2grd ( pkzgather(1,L),pblindex, pkzgather(1,L) ,im*jm)
C Diagnostics:
c ------------
if(icldmas.gt.0)call back2grd(tmpgather(1,L),pblindex,
. tmpgather(1,L),im*jm)
if(idtrain.gt.0)call back2grd(pkegather(1,L),pblindex,
. pkegather(1,L),im*jm)
if(idtls.gt.0)call back2grd(deltrnev(1,L),pblindex,
. deltrnev(1,L),im*jm)
if(idqls.gt.0)call back2grd(delqrnev(1,L),pblindex,
. delqrnev(1,L),im*jm)
if(icldnp.gt.0)call back2grd(cldsr(1,1,L),pblindex,
. cldsr(1,1,L),im*jm)
enddo
c Tracers
c -------
do nt = 1,ntracer-ptracer
do L = 1,lm
call back2grd (delugather(1,L,nt),pblindex,
. dqmoist(1,1,L,ptracer+nt),im*jm)
enddo
enddo
C **********************************************************************
C BUMP DIAGNOSTICS
C **********************************************************************
c Sub-Cloud Layer
c -------------------------
if( ipsubcld.ne.0 ) then
do j = 1,jm
do i = 1,im
qdiag(i,j,ipsubcld,bi,bj) = qdiag(i,j,ipsubcld,bi,bj) +
. psubcldg (i,j)
qdiag(i,j,ipsubcldc,bi,bj) = qdiag(i,j,ipsubcldc,bi,bj) +
. psubcldgc(i,j)
enddo
enddo
endif
c Non-Precipitating Cloud Fraction
c --------------------------------
if( icldnp.ne.0 ) then
do L = 1,lm
do j = 1,jm
do i = 1,im
qdiag(i,j,icldnp+L-1,bi,bj) = qdiag(i,j,icldnp+L-1,bi,bj) +
. cldsr(i,j,L)
enddo
enddo
enddo
ncldnp = ncldnp + 1
endif
c Moist Processes Heating Rate
c ----------------------------
if(imoistt.gt.0) then
do L = 1,lm
do i = 1,im*jm
qdiag(i,1,imoistt+L-1,bi,bj) = qdiag(i,1,imoistt+L-1,bi,bj) +
. (dtmoist(i,1,L)*sday*pkzgather(i,L)/pz(i,1))
enddo
enddo
endif
c Moist Processes Moistening Rate
c -------------------------------
if(imoistq.gt.0) then
do L = 1,lm
do j = 1,jm
do i = 1,im
qdiag(i,j,imoistq+L-1,bi,bj) = qdiag(i,j,imoistq+L-1,bi,bj) +
. (dqmoist(i,j,L,1)*sday*1000.0/pz(i,j))
enddo
enddo
enddo
endif
c Cloud Mass Flux
c ---------------
if(icldmas.gt.0) then
do L = 1,lm
do i = 1,im*jm
qdiag(i,1,icldmas+L-1,bi,bj) = qdiag(i,1,icldmas+L-1,bi,bj) +
. tmpgather(i,L)
enddo
enddo
endif
c Detrained Cloud Mass Flux
c -------------------------
if(idtrain.gt.0) then
do L = 1,lm
do i = 1,im*jm
qdiag(i,1,idtrain+L-1,bi,bj) = qdiag(i,1,idtrain+L-1,bi,bj) +
. pkegather(i,L)
enddo
enddo
endif
c Grid-Scale Condensational Heating Rate
c --------------------------------------
if(idtls.gt.0) then
do L = 1,lm
do i = 1,im*jm
qdiag(i,1,idtls+L-1,bi,bj) = qdiag(i,1,idtls+L-1,bi,bj) +
. deltrnev(i,L)
enddo
enddo
endif
c Grid-Scale Condensational Moistening Rate
c -----------------------------------------
if(idqls.gt.0) then
do L = 1,lm
do i = 1,im*jm
qdiag(i,1,idqls+L-1,bi,bj) = qdiag(i,1,idqls+L-1,bi,bj) +
. delqrnev(i,L)
enddo
enddo
endif
c Total Precipitation
c -------------------
if(ipreacc.gt.0) then
do j = 1,jm
do i = 1,im
qdiag(i,j,ipreacc,bi,bj) = qdiag(i,j,ipreacc,bi,bj)
. + ( lsp_new(I,j)
. + snow_new(I,j)
. + conv_new(i,j) ) *sday*tminv
enddo
enddo
endif
c Convective Precipitation
c ------------------------
if(iprecon.gt.0) then
do i = 1,im*jm
qdiag(i,1,iprecon,bi,bj) = qdiag(i,1,iprecon,bi,bj) +
. raincgath(i)*sday*tminv
enddo
endif
C **********************************************************************
C **** Fill Rainfall and Snowfall Arrays for Land Surface Model ****
C **** Note: Precip Rates work when DT(turb)
Cloud Fraction from RAS ***
C *** CLDLS => Cloud Fraction from RNEVP ***
C **********************************************************************
do j = 1,jm
do i = 1,im
cldhi (i,j) = 0.
cldmid(i,j) = 0.
cldlow(i,j) = 0.
cldtmp(i,j) = 0.
cldprs(i,j) = 0.
tmpimjm(i,j) = 0.
enddo
enddo
c Set Moist-Process Memory Coefficient
c ------------------------------------
cldras_mem = 1.0-tmstp/ 3600.0
cldlsp_mem = 1.0-tmstp/(3600.0*3)
do L = 1,lm
do i = 1,im*jm
plev = pl(i,L)
c Compute Time-averaged Cloud and Water Amounts for Longwave Radiation
c --------------------------------------------------------------------
watnow = cldwater(i,1,L)
if( plev.le.500.0 ) then
cldras = min( max( cldras_lw(i,1,L)*cldras_mem,cpen(i,1,L)),1.0)
else
cldras = 0.0
endif
cldlsp = min( max( cldlsp_lw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0)
if( cldras.lt.cldmin ) cldras = 0.0
if( cldlsp.lt.cldmin ) cldlsp = 0.0
cldnow = max( cldlsp,cldras )
lwlz(i,1,L) = ( nlwlz*lwlz(i,1,L) + watnow)/(nlwlz +1)
cldtot_lw(i,1,L) = (nlwcld*cldtot_lw(i,1,L) + cldnow)/(nlwcld+1)
cldlsp_lw(i,1,L) = (nlwcld*cldlsp_lw(i,1,L) + cldlsp)/(nlwcld+1)
cldras_lw(i,1,L) = (nlwcld*cldras_lw(i,1,L) + cldras)/(nlwcld+1)
c Compute Time-averaged Cloud and Water Amounts for Shortwave Radiation
c ---------------------------------------------------------------------
watnow = cldwater(i,1,L)
if( plev.le.500.0 ) then
cldras = min( max(cldras_sw(i,1,L)*cldras_mem, cpen(i,1,L)),1.0)
else
cldras = 0.0
endif
cldlsp = min( max(cldlsp_sw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0)
if( cldras.lt.cldmin ) cldras = 0.0
if( cldlsp.lt.cldmin ) cldlsp = 0.0
cldnow = max( cldlsp,cldras )
swlz(i,1,L) = ( nswlz*swlz(i,1,L) + watnow)/(nswlz +1)
cldtot_sw(i,1,L) = (nswcld*cldtot_sw(i,1,L) + cldnow)/(nswcld+1)
cldlsp_sw(i,1,L) = (nswcld*cldlsp_sw(i,1,L) + cldlsp)/(nswcld+1)
cldras_sw(i,1,L) = (nswcld*cldras_sw(i,1,L) + cldras)/(nswcld+1)
c Compute Instantaneous Low-Mid-High Maximum Overlap Cloud Fractions
c ----------------------------------------------------------------------
if( L.lt.midlevel ) cldhi (i,1) = max( cldnow,cldhi (i,1) )
if( L.ge.midlevel .and.
. L.lt.lowlevel ) cldmid(i,1) = max( cldnow,cldmid(i,1) )
if( L.ge.lowlevel ) cldlow(i,1) = max( cldnow,cldlow(i,1) )
c Compute Cloud-Top Temperature and Pressure
c ------------------------------------------
cldtmp(i,1) = cldtmp(i,1) + cldnow*pkzgather(i,L)
. * ( tz(i,1,L) + dtmoist(i,1,L)*tmstp/pz(i,1) )
cldprs(i,1) = cldprs(i,1) + cldnow*plev
tmpimjm(i,1) = tmpimjm(i,1) + cldnow
enddo
enddo
c Compute Instantanious Total 2-D Cloud Fraction
c ----------------------------------------------
do j = 1,jm
do i = 1,im
totcld(i,j) = 1.0 - (1.-cldhi (i,j))
. * (1.-cldmid(i,j))
. * (1.-cldlow(i,j))
enddo
enddo
C **********************************************************************
C *** Fill Cloud Top Pressure and Temperature Diagnostic ***
C **********************************************************************
if(icldtmp.gt.0) then
do j = 1,jm
do i = 1,im
if( cldtmp(i,j).gt.0.0 ) then
qdiag(i,j,icldtmp,bi,bj) = qdiag(i,j,icldtmp,bi,bj) +
. cldtmp(i,j)*totcld(i,j)/tmpimjm(i,j)
qdiag(i,j,icttcnt,bi,bj) = qdiag(i,j,icttcnt,bi,bj) +
. totcld(i,j)
endif
enddo
enddo
endif
if(icldprs.gt.0) then
do j = 1,jm
do i = 1,im
if( cldprs(i,j).gt.0.0 ) then
qdiag(i,j,icldprs,bi,bj) = qdiag(i,j,icldprs,bi,bj) +
. cldprs(i,j)*totcld(i,j)/tmpimjm(i,j)
qdiag(i,j,ictpcnt,bi,bj) = qdiag(i,j,ictpcnt,bi,bj) +
. totcld(i,j)
endif
enddo
enddo
endif
C **********************************************************************
C **** INCREMENT COUNTERS ****
C **********************************************************************
nlwlz = nlwlz + 1
nswlz = nswlz + 1
nlwcld = nlwcld + 1
nswcld = nswcld + 1
nmoistt = nmoistt + 1
nmoistq = nmoistq + 1
npreacc = npreacc + 1
nprecon = nprecon + 1
ncldmas = ncldmas + 1
ndtrain = ndtrain + 1
ndtls = ndtls + 1
ndqls = ndqls + 1
RETURN
END
SUBROUTINE RAS( NN, LEN, LENC, K, NLTOP, nlayr, DT
*, UOI, ntracer, POI, QOI, PRS, PRJ, rnd, ncrnd
*, RAINS, CLN, CLF, cldmas, detrain
*, cp,grav,rkappa,alhl,rhfrac,rasmax )
C
C*********************************************************************
C********************* SUBROUTINE RAS *****************************
C********************** 16 MARCH 1988 ******************************
C*********************************************************************
C
implicit none
C Argument List
integer nn,len,lenc,k,nltop,nlayr
integer ntracer
integer ncrnd
real dt
real UOI(len,nlayr,ntracer), POI(len,K)
real QOI(len,K), PRS(len,K+1), PRJ(len,K+1)
real rnd(ncrnd)
real RAINS(len,K), CLN(len,K), CLF(len,K)
real cldmas(len,K), detrain(len,K)
real cp,grav,rkappa,alhl,rhfrac(len),rasmax
C Local Variables
real TCU(len,K), QCU(len,K)
real ucu(len,K,ntracer)
real ALF(len,K), BET(len,K), GAM(len,K)
*, ETA(len,K), HOI(len,K)
*, PRH(len,K), PRI(len,K)
real HST(len,K), QOL(len,K), GMH(len,K)
real TX1(len), TX2(len), TX3(len), TX4(len), TX5(len)
*, TX6(len), TX7(len), TX8(len), TX9(len)
*, TX11(len), TX12(len), TX13(len), TX14(len,ntracer)
*, TX15(len)
*, WFN(len)
integer IA1(len), IA2(len), IA3(len)
real cloudn(len), pcu(len)
integer krmin,icm
real rknob, cmb2pa
PARAMETER (KRMIN=01)
PARAMETER (ICM=1000)
PARAMETER (CMB2PA=100.0)
PARAMETER (rknob = 10.)
integer IC(ICM), IRND(icm)
real cmass(len,K)
LOGICAL SETRAS
integer i,L,nc,ib,nt
integer km1,kp1,kprv,kcr,kfx,ncmx
real p00, crtmsf, frac, rasblf
do L = 1,k
do I = 1,LENC
rains(i,l) = 0.
enddo
enddo
p00 = 1000.
crtmsf = 0.
C The numerator here is the fraction of the subcloud layer mass flux
C allowed to entrain into the cloud
CCC FRAC = 1./dt
FRAC = 0.5/dt
KM1 = K - 1
KP1 = K + 1
C we want the ras adjustment time scale to be one hour (indep of dt)
RASBLF = 1./3600.
C
KPRV = KM1
C Removed KRMAX parameter
KCR = MIN(KM1,nlayr-2)
KFX = KM1 - KCR
NCMX = KFX + NCRND
C
IF (KFX .GT. 0) THEN
DO NC=1,KFX
IC(NC) = K - NC
ENDDO
ENDIF
C
IF (NCRND .GT. 0) THEN
DO I=1,ncrnd
IRND(I) = (RND(I)-0.0005)*(KCR-KRMIN+1)
IRND(I) = IRND(I) + KRMIN
ENDDO
C
DO NC=1,NCRND
IC(KFX+NC) = IRND(NC)
ENDDO
ENDIF
C
DO 100 NC=1,NCMX
C
IF (NC .EQ. 1 ) THEN
SETRAS = .TRUE.
ELSE
SETRAS = .FALSE.
ENDIF
IB = IC(NC)
c Initialize Cloud Fraction Array
c -------------------------------
do i = 1,lenc
cloudn(i) = 0.0
enddo
CALL CLOUD(nn,LEN, LENC, K, NLTOP, nlayr, IB, RASBLF,SETRAS,FRAC
*, CP, ALHL, RKAPPA, GRAV, P00, CRTMSF
*, POI, QOI, UOI, Ntracer, PRS, PRJ
*, PCU, CLOUDN, TCU, QCU, UCU, CMASS
*, ALF, BET, GAM, PRH, PRI, HOI, ETA
*, HST, QOL, GMH
*, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9
*, WFN, TX11, TX12, TX13, TX14, TX15
*, IA1,IA2,IA3,rhfrac)
C Compute fraction of grid box into which rain re-evap occurs (clf)
c -----------------------------------------------------------------
do i = 1,lenc
c mass in detrainment layer
c -------------------------
tx1(i) = cmb2pa * (prs(i,ib+1) - prs(i,ib))/(grav*dt)
c ratio of detraining cloud mass to mass in detrainment layer
c -----------------------------------------------------------
tx1(i) = rhfrac(i)*rknob * cmass(i,ib) / tx1(i)
if(cmass(i,K).gt.0.) clf(i,ib) = clf(i,ib) + tx1(i)
if( clf(i,ib).gt.1.) clf(i,ib) = 1.
enddo
c Compute Total Cloud Mass Flux
c *****************************
do L=ib,k
do i=1,lenc
cmass(i,L) = rhfrac(i)*cmass(i,L) * dt
enddo
enddo
do L=ib,k
do i=1,lenc
cldmas(i,L) = cldmas(i,L) + cmass(i,L)
enddo
enddo
do i=1,lenc
detrain(i,ib) = detrain(i,ib) + cmass(i,ib)
enddo
DO L=IB,K
DO I=1,LENC
POI(I,L) = POI(I,L) + TCU(I,L) * DT * rhfrac(i)
QOI(I,L) = QOI(I,L) + QCU(I,L) * DT * rhfrac(i)
ENDDO
ENDDO
DO NT=1,Ntracer
DO L=IB,K
DO I=1,LENC
UOI(I,L+nltop-1,NT)=UOI(I,L+nltop-1,NT)+UCU(I,L,NT)*DT*rhfrac(i)
ENDDO
ENDDO
ENDDO
DO I=1,LENC
rains(I,ib) = rains(I,ib) + PCU(I)*dt * rhfrac(i)
ENDDO
100 CONTINUE
c Fill Convective Cloud Fractions based on 3-D Rain Amounts
c ---------------------------------------------------------
do L=k-1,1,-1
do i=1,lenc
tx1(i) = 100*(prs(i,L+1)-prs(i,L))/grav
cln(i,L) = min(1600*rains(i,L)/tx1(i),rasmax )
enddo
enddo
RETURN
END
subroutine rndcloud (iras,nrnd,rnd,myid)
implicit none
integer n,iras,nrnd,myid
real random_numbx
real rnd(nrnd)
integer irm
parameter (irm = 1000)
real random(irm)
integer i,mcheck,numrand,iseed,index
logical first
data first /.true./
integer iras0
data iras0 /0/
save random, iras0
if(nrnd.eq.0.)then
do i = 1,nrnd
rnd(i) = 0
enddo
if(first .and. myid.eq.0) print *,' NO RANDOM CLOUDS IN RAS '
go to 100
endif
mcheck = mod(iras-1,irm/nrnd)
c First Time In From a Continuing RESTART (IRAS.GT.1) or Reading a New RESTART
c ----------------------------------------------------------------------------
if( first.and.(iras.gt.1) .or. iras.ne.iras0+1 )then
if( myid.eq.0 ) print *, 'Recreating Rand Numb Array in RNDCLOUD'
if( myid.eq.0 ) print *, 'IRAS: ',iras,' IRAS0: ',iras0
numrand = mod(iras,irm/nrnd) * nrnd
iseed = iras * nrnd - numrand
call random_seedx(iseed)
do i = 1,irm
random(i) = random_numbx()
enddo
index = (iras-1)*nrnd
c Multiple Time In But have Used Up all 1000 numbers (MCHECK.EQ.0)
c ----------------------------------------------------------------
else if (mcheck.eq.0) then
iseed = (iras-1)*nrnd
call random_seedx(iseed)
do i = 1,irm
random(i) = random_numbx()
enddo
index = iseed
c Multiple Time In But have NOT Used Up all 1000 numbers (MCHECK.NE.0)
c --------------------------------------------------------------------
else
index = (iras-1)*nrnd
endif
index = mod(index,irm)
if( index+nrnd.gt.1000 ) index=1000-nrnd
do n = 1,nrnd
rnd(n) = random(index+n)
enddo
100 continue
first = .false.
iras0 = iras
return
end
function random_numbx()
implicit none
real random_numbx
#ifdef CRAY
real ranf
random_numbx = ranf()
#endif
#ifdef SGI
real rand
random_numbx = rand()
#endif
return
end
subroutine random_seedx (iseed)
implicit none
integer iseed
#ifdef CRAY
call ranset (iseed)
#endif
#ifdef SGI
integer*4 seed
seed = iseed
call srand (seed)
#endif
return
end
SUBROUTINE CLOUD(nn,LEN, LENC, K, NLTOP, nlayr, IC, RASALF
*, SETRAS, FRAC
*, CP, ALHL, RKAP, GRAV, P00, CRTMSF
*, POI, QOI, UOI, Ntracer, PRS, PRJ
*, PCU, CLN, TCU, QCU, UCU, CMASS
*, ALF, BET, GAM, PRH, PRI, HOL, ETA
*, HST, QOL, GMH
*, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, ALM
*, WFN, AKM, QS1, CLF, UHT, WLQ
*, IA, I1, I2,rhfrac)
C
C*********************************************************************
C******************** Relaxed Arakawa-Schubert ***********************
C********************* Plug Compatible Version **********************
C************************ SUBROUTINE CLOUD ***************************
C************************* 23 JULY 1992 ***************************
C*********************************************************************
C*********************************************************************
C*********************************************************************
C************************** Developed By *****************************
C************************** *****************************
C************************ Shrinivas Moorthi **************************
C************************ and **************************
C************************ Max J. Suarez *****************************
C************************ *****************************
C******************** Laboratory for Atmospheres *********************
C****************** NASA/GSFC, Greenbelt, MD 20771 *******************
C*********************************************************************
C*********************************************************************
C
C The calculations of Moorthi and Suarez (1992, MWR) are
C contained in the CLOUD routine.
C It is probably advisable, at least initially, to treat CLOUD
C as a black box that computes the single cloud adjustments. RAS,
C on the other hand, can be tailored to each GCMs configuration
C (ie, number and placement of levels, nature of boundary layer,
C time step and frequency with which RAS is called).
C
C
C Input:
C ------
C
C LEN : The inner dimension of update and input arrays.
C
C LENC : The run: the number of soundings processes in a single call.
C RAS works on the first LENC of the LEN soundings
C passed. This allows working on pieces of the world
C say for multitasking, without declaring temporary arrays
C and copying the data to and from them. This is an f77
C version. An F90 version would have to allow more
C flexibility in the argument declarations. Obviously
C (LENC<=LEN).
C
C K : Number of vertical layers (increasing downwards).
C Need not be the same as the number of layers in the
C GCM, since it is the outer dimension. The bottom layer
C (K) is the subcloud layer.
C
C IC : Detrainment level to check for presence of convection
C
C RASALF : Relaxation parameter (< 1.) for present cloud-type
C
C SETRAS : Logical parameter to control re-calculation of
C saturation specific humidity and mid level P**kappa
C
C FRAC : Fraction of the PBL (layer K) mass allowed to be used
C by a cloud-type in time DT
C
C CP : Specific heat at constant pressure
C
C ALHL : Latent Heat of condensation
C
C RKAP : R/Cp, where R is the gas constant
C
C GRAV : Acceleration due to gravity
C
C P00 : A reference pressure in hPa, useually 1000 hPa
C
C CRTMSF : Critical value of mass flux above which cloudiness at
C the detrainment layer of that cloud-type is assumed.
C Affects only cloudiness calculation.
C
C POI : 2D array of dimension (LEN,K) containing potential
C temperature. Updated but not initialized by RAS.
C
C QOI : 2D array of dimension (LEN,K) containing specific
C humidity. Updated but not initialized by RAS.
C
C UOI : 3D array of dimension (LEN,K,NTRACER) containing tracers
C Updated but not initialized by RAS.
C
C PRS : 2D array of dimension (LEN,K+1) containing pressure
C in hPa at the interfaces of K-layers from top of the
C atmosphere to the bottom. Not modified.
C
C PRJ : 2D array of dimension (LEN,K+1) containing (PRS/P00) **
C RKAP. i.e. Exner function at layer edges. Not modified.
C
C rhfrac : 1D array of dimension (LEN) containing a rel.hum. scaling
C fraction. Not modified.
C
C Output:
C -------
C
C PCU : 1D array of length LEN containing accumulated
C precipitation in mm/sec.
C
C CLN : 2D array of dimension (LEN,K) containing cloudiness
C Note: CLN is bumped but NOT initialized
C
C TCU : 2D array of dimension (LEN,K) containing accumulated
C convective heating (K/sec).
C
C QCU : 2D array of dimension (LEN,K) containing accumulated
C convective drying (kg/kg/sec).
C
C CMASS : 2D array of dimension (LEN,K) containing the
C cloud mass flux (kg/sec). Filled from cloud top
C to base.
C
C Temporaries:
C
C ALF, BET, GAM, ETA, PRH, PRI, HOI, HST, QOL, GMH are temporary
C 2D real arrays of dimension of at least (LENC,K) where LENC is
C the horizontal dimension over which convection is invoked.
C
C
C TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9, AKM, QS1, CLF, UHT
C VHT, WLQ WFN are temporary real arrays of length at least LENC
C
C IA, I1, and I2 are temporary integer arrays of length LENC
C
C
C************************************************************************
implicit none
C Argument List declarations
integer nn,LEN,LENC,K,NLTOP,nlayr,ic,ntracer
real rasalf
LOGICAL SETRAS
real frac, cp, alhl, rkap, grav, p00, crtmsf
real POI(LEN,K),QOI(LEN,K),PRS(LEN,K+1),PRJ(LEN,K+1)
real uoi(len,nlayr,ntracer)
real PCU(LENC), CLN(LEN)
real TCU(LEN,K), QCU(LEN,K), ucu(len,k,ntracer), CMASS(LEN,K)
real ALF(LEN,K), BET(LEN,K), GAM(LEN,K), PRH(LEN,K), PRI(LEN,K)
real HOL(LENC,K), ETA(LENC,K), HST(LENC,K), QOL(LENC,K)
real GMH(LENC,K)
real TX1(LENC), TX2(LENC), TX3(LENC), TX4(LENC)
real TX5(LENC), TX6(LENC), TX7(LENC), TX8(LENC)
real ALM(LENC), WFN(LENC), AKM(LENC), QS1(LENC)
real WLQ(LENC), CLF(LENC)
real uht(len,ntracer)
integer IA(LENC), I1(LENC),I2(LENC)
real rhfrac(len)
C Local Variables
real daylen,half,one,zero,cmb2pa,rhmax
PARAMETER (DAYLEN=86400.0, HALF=0.5, ONE=1.0, ZERO=0.0)
PARAMETER (CMB2PA=100.0)
PARAMETER (RHMAX=0.9999)
real rkapp1,onebcp,albcp,onebg,cpbg,twobal
C
integer nt,km1,ic1,i,L,len1,len2,isav,len11,ii
integer lena,lena1,lenb,tem,tem1
c Explicit Inline Directives
c --------------------------
#ifdef CRAY
#ifdef f77
cfpp$ expand (qsat)
#endif
#endif
RKAPP1 = 1.0 + RKAP
ONEBCP = 1.0 / CP
ALBCP = ALHL * ONEBCP
ONEBG = 1.0 / GRAV
CPBG = CP * ONEBG
TWOBAL = 2.0 / ALHL
C
KM1 = K - 1
IC1 = IC + 1
C
C SETTING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T.
C
IF (SETRAS) THEN
DO 2050 L=1,K
DO 2030 I=1,LENC
PRH(I,L) = (PRJ(I,L+1)*PRS(I,L+1) - PRJ(I,L)*PRS(I,L))
* / ((PRS(I,L+1)-PRS(I,L)) * RKAPP1)
2030 CONTINUE
2050 CONTINUE
DO 2070 L=1,K
DO 2060 I=1,LENC
TX5(I) = POI(I,L) * PRH(I,L)
TX1(I) = (PRS(I,L) + PRS(I,L+1)) * 0.5
TX3(I) = TX5(I)
CALL QSAT(TX3(I), TX1(I), TX2(I), TX4(I), .TRUE.)
ALF(I,L) = TX2(I) - TX4(I) * TX5(I)
BET(I,L) = TX4(I) * PRH(I,L)
GAM(I,L) = 1.0 / ((1.0 + TX4(I)*ALBCP) * PRH(I,L))
PRI(I,L) = (CP/CMB2PA) / (PRS(I,L+1) - PRS(I,L))
2060 CONTINUE
2070 CONTINUE
ENDIF
C
C
DO 10 L=1,K
DO 10 I=1,LEN
TCU(I,L) = 0.0
QCU(I,L) = 0.0
CMASS(I,L) = 0.0
10 CONTINUE
do nt = 1,ntracer
do L=1,K
do I=1,LENC
ucu(I,L,nt) = 0.0
enddo
enddo
enddo
C
DO 30 I=1,LENC
TX1(I) = PRJ(I,K+1) * POI(I,K)
QS1(I) = ALF(I,K) + BET(I,K)*POI(I,K)
QOL(I,K) = MIN(QS1(I)*RHMAX,QOI(I,K))
HOL(I,K) = TX1(I)*CP + QOL(I,K)*ALHL
ETA(I,K) = ZERO
TX2(I) = (PRJ(I,K+1) - PRJ(I,K)) * POI(I,K) * CP
30 CONTINUE
C
IF (IC .LT. KM1) THEN
DO 3703 L=KM1,IC1,-1
DO 50 I=1,LENC
QS1(I) = ALF(I,L) + BET(I,L)*POI(I,L)
QOL(I,L) = MIN(QS1(I)*RHMAX,QOI(I,L))
C
TEM1 = TX2(I) + PRJ(I,L+1) * POI(I,L) * CP
HOL(I,L) = TEM1 + QOL(I,L )* ALHL
HST(I,L) = TEM1 + QS1(I) * ALHL
TX1(I) = (PRJ(I,L+1) - PRJ(I,L)) * POI(I,L)
ETA(I,L) = ETA(I,L+1) + TX1(I)*CPBG
TX2(I) = TX2(I) + TX1(I)*CP
50 CONTINUE
C
3703 CONTINUE
ENDIF
DO 70 I=1,LENC
HOL(I,IC) = TX2(I)
QS1(I) = ALF(I,IC) + BET(I,IC)*POI(I,IC)
QOL(I,IC) = MIN(QS1(I)*RHMAX,QOI(I,IC))
c
TEM1 = TX2(I) + PRJ(I,IC1) * POI(I,IC) * CP
HOL(I,IC) = TEM1 + QOL(I,IC) * ALHL
HST(I,IC) = TEM1 + QS1(I) * ALHL
C
TX3(I ) = (PRJ(I,IC1) - PRH(I,IC)) * POI(I,IC)
ETA(I,IC) = ETA(I,IC1) + CPBG * TX3(I)
70 CONTINUE
C
DO 130 I=1,LENC
TX2(I) = HOL(I,K) - HST(I,IC)
TX1(I) = ZERO
130 CONTINUE
C
C ENTRAINMENT PARAMETER
C
DO 160 L=IC,KM1
DO 160 I=1,LENC
TX1(I) = TX1(I) + (HST(I,IC) - HOL(I,L)) * (ETA(I,L) - ETA(I,L+1))
160 CONTINUE
C
LEN1 = 0
LEN2 = 0
ISAV = 0
DO 195 I=1,LENC
IF (TX1(I) .GT. ZERO .AND. TX2(I) .GT. ZERO
. .AND. rhfrac(i).ne.0.0 ) THEN
LEN1 = LEN1 + 1
IA(LEN1) = I
ALM(LEN1) = TX2(I) / TX1(I)
ENDIF
195 CONTINUE
C
LEN2 = LEN1
if (IC1 .lt. K) then
DO 196 I=1,LENC
IF (TX2(I) .LE. 0.0 .AND. (HOL(I,K) .GT. HST(I,IC1))
. .AND. rhfrac(i).ne.0.0 ) THEN
LEN2 = LEN2 + 1
IA(LEN2) = I
ALM(LEN2) = 0.0
ENDIF
196 CONTINUE
endif
C
IF (LEN2 .EQ. 0) THEN
DO 5010 I=1,LENC*K
HST(I,1) = 0.0
QOL(I,1) = 0.0
5010 CONTINUE
DO 5020 I=1,LENC
PCU(I) = 0.0
5020 CONTINUE
RETURN
ENDIF
LEN11 = LEN1 + 1
C
C NORMALIZED MASSFLUX
C
DO 250 I=1,LEN2
ETA(I,K) = 1.0
II = IA(I)
TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
TX4(I) = PRS(II,K)
250 CONTINUE
C
DO 252 I=LEN11,LEN2
WFN(I) = 0.0
II = IA(I)
IF (HST(II,IC1) .LT. HST(II,IC)) THEN
TX6(I) = (HST(II,IC1)-HOL(II,K))/(HST(II,IC1)-HST(II,IC))
ELSE
TX6(I) = 0.0
ENDIF
TX2(I) = 0.5 * (PRS(II,IC1)+PRS(II,IC1+1)) * (1.0-TX6(I))
* + TX2(I) * TX6(I)
252 CONTINUE
C
CALL ACRITN(LEN2, TX2, TX4, TX3)
C
DO 260 L=KM1,IC,-1
DO 255 I=1,LEN2
TX1(I) = ETA(IA(I),L)
255 CONTINUE
DO 260 I=1,LEN2
ETA(I,L) = 1.0 + ALM(I) * TX1(I)
260 CONTINUE
C
C CLOUD WORKFUNCTION
C
IF (LEN1 .GT. 0) THEN
DO 270 I=1,LEN1
II = IA(I)
WFN(I) = - GAM(II,IC) * (PRJ(II,IC1) - PRH(II,IC))
* * HST(II,IC) * ETA(I,IC1)
270 CONTINUE
ENDIF
C
DO 290 I=1,LEN2
II = IA(I)
TX1(I) = HOL(II,K)
290 CONTINUE
C
IF (IC1 .LE. KM1) THEN
DO 380 L=KM1,IC1,-1
DO 380 I=1,LEN2
II = IA(I)
TEM = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * HOL(II,L)
C
PCU(I) = PRJ(II,L+1) - PRH(II,L)
TEM1 = ETA(I,L+1) * PCU(I)
TX1(I) = TX1(I)*PCU(I)
C
PCU(I) = PRH(II,L) - PRJ(II,L)
TEM1 = (TEM1 + ETA(I,L) * PCU(I)) * HST(II,L)
TX1(I) = TX1(I) + TEM*PCU(I)
C
WFN(I) = WFN(I) + (TX1(I) - TEM1) * GAM(II,L)
TX1(I) = TEM
380 CONTINUE
ENDIF
C
LENA = 0
IF (LEN1 .GT. 0) THEN
DO 512 I=1,LEN1
II = IA(I)
WFN(I) = WFN(I) + TX1(I) * GAM(II,IC)*(PRJ(II,IC1)-PRH(II,IC))
* - TX3(I)
IF (WFN(I) .GT. 0.0) THEN
LENA = LENA + 1
I1(LENA) = IA(I)
I2(LENA) = I
TX1(LENA) = WFN(I)
TX2(LENA) = QS1(IA(I))
TX6(LENA) = 1.0
ENDIF
512 CONTINUE
ENDIF
LENB = LENA
DO 515 I=LEN11,LEN2
WFN(I) = WFN(I) - TX3(I)
IF (WFN(I) .GT. 0.0 .AND. TX6(I) .GT. 0.0) THEN
LENB = LENB + 1
I1(LENB) = IA(I)
I2(LENB) = I
TX1(LENB) = WFN(I)
TX2(LENB) = QS1(IA(I))
TX4(LENB) = TX6(I)
ENDIF
515 CONTINUE
C
IF (LENB .LE. 0) THEN
DO 5030 I=1,LENC*K
HST(I,1) = 0.0
QOL(I,1) = 0.0
5030 CONTINUE
DO 5040 I=1,LENC
PCU(I) = 0.0
5040 CONTINUE
RETURN
ENDIF
C
DO 516 I=1,LENB
WFN(I) = TX1(I)
QS1(I) = TX2(I)
516 CONTINUE
C
DO 520 L=IC,K
DO 517 I=1,LENB
TX1(I) = ETA(I2(I),L)
517 CONTINUE
DO 520 I=1,LENB
ETA(I,L) = TX1(I)
520 CONTINUE
C
LENA1 = LENA + 1
C
DO 510 I=1,LENA
II = I1(I)
TX8(I) = HST(II,IC) - HOL(II,IC)
510 CONTINUE
DO 530 I=LENA1,LENB
II = I1(I)
TX6(I) = TX4(I)
TEM = TX6(I) * (HOL(II,IC)-HOL(II,IC1)) + HOL(II,IC1)
TX8(I) = HOL(II,K) - TEM
TEM1 = TX6(I) * (QOL(II,IC)-QOL(II,IC1)) + QOL(II,IC1)
TX5(I) = TEM - TEM1 * ALHL
QS1(I) = TEM1 + TX8(I)*(ONE/ALHL)
TX3(I) = HOL(II,IC)
530 CONTINUE
C
C
DO 620 I=1,LENB
II = I1(I)
WLQ(I) = QOL(II,K) - QS1(I) * ETA(I,IC)
TX7(I) = HOL(II,K)
620 CONTINUE
DO NT=1,Ntracer
DO 621 I=1,LENB
II = I1(I)
UHT(I,NT) = UOI(II,K+nltop-1,NT)-UOI(II,IC+nltop-1,NT) * ETA(I,IC)
621 CONTINUE
ENDDO
C
DO 635 L=KM1,IC,-1
DO 630 I=1,LENB
II = I1(I)
TEM = ETA(I,L) - ETA(I,L+1)
WLQ(I) = WLQ(I) + TEM * QOL(II,L)
630 CONTINUE
635 CONTINUE
DO NT=1,Ntracer
DO L=KM1,IC,-1
DO I=1,LENB
II = I1(I)
TEM = ETA(I,L) - ETA(I,L+1)
UHT(I,NT) = UHT(I,NT) + TEM * UOI(II,L+nltop-1,NT)
ENDDO
ENDDO
ENDDO
C
C CALCULATE GS AND PART OF AKM (THAT REQUIRES ETA)
C
DO 690 I=1,LENB
II = I1(I)
c TX7(I) = HOL(II,K)
TEM = (POI(II,KM1) - POI(II,K)) / (PRH(II,K) - PRH(II,KM1))
HOL(I,K) = TEM * (PRJ(II,K)-PRH(II,KM1))*PRH(II,K)*PRI(II,K)
HOL(I,KM1) = TEM * (PRH(II,K)-PRJ(II,K))*PRH(II,KM1)*PRI(II,KM1)
AKM(I) = ZERO
TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
690 CONTINUE
IF (IC1 .LE. KM1) THEN
DO 750 L=KM1,IC1,-1
DO 750 I=1,LENB
II = I1(I)
TEM = (POI(II,L-1) - POI(II,L)) * ETA(I,L)
* / (PRH(II,L) - PRH(II,L-1))
C
HOL(I,L) = TEM * (PRJ(II,L)-PRH(II,L-1)) * PRH(II,L)
* * PRI(II,L) + HOL(I,L)
HOL(I,L-1) = TEM * (PRH(II,L)-PRJ(II,L)) * PRH(II,L-1)
* * PRI(II,L-1)
C
AKM(I) = AKM(I) - HOL(I,L)
* * (ETA(I,L) * (PRH(II,L)-PRJ(II,L)) +
* ETA(I,L+1) * (PRJ(II,L+1)-PRH(II,L))) / PRH(II,L)
750 CONTINUE
ENDIF
C
C
CALL RNCL(LENB, TX2, TX1, CLF)
DO 770 I=1,LENB
TX2(I) = (ONE - TX1(I)) * WLQ(I)
WLQ(I) = TX1(I) * WLQ(I)
C
TX1(I) = HOL(I,IC)
770 CONTINUE
DO 790 I=LENA1, LENB
II = I1(I)
TX1(I) = TX1(I) + (TX5(I)-TX3(I)+QOL(II,IC)*ALHL)*(PRI(II,IC)/CP)
790 CONTINUE
DO 800 I=1,LENB
HOL(I,IC) = TX1(I) - TX2(I) * ALBCP * PRI(I1(I),IC)
800 CONTINUE
IF (LENA .GT. 0) THEN
DO 810 I=1,LENA
II = I1(I)
AKM(I) = AKM(I) - ETA(I,IC1) * (PRJ(II,IC1) - PRH(II,IC))
* * TX1(I) / PRH(II,IC)
810 CONTINUE
ENDIF
c
C CALCULATE GH
C
DO 830 I=1,LENB
II = I1(I)
TX3(I) = QOL(II,KM1) - QOL(II,K)
GMH(I,K) = HOL(I,K) + TX3(I) * PRI(II,K) * (ALBCP)
AKM(I) = AKM(I) + GAM(II,KM1)*(PRJ(II,K)-PRH(II,KM1))
* * GMH(I,K)
TX3(I) = zero
830 CONTINUE
C
IF (IC1 .LE. KM1) THEN
DO 840 L=KM1,IC1,-1
DO 840 I=1,LENB
II = I1(I)
TX2(I) = TX3(I)
TX3(I) = (QOL(II,L-1) - QOL(II,L)) * ETA(I,L)
TX2(I) = TX2(I) + TX3(I)
C
GMH(I,L) = HOL(I,L) + TX2(I) * PRI(II,L) * (ALBCP*HALF)
840 CONTINUE
C
C
ENDIF
DO 850 I=LENA1,LENB
TX3(I) = TX3(I) + TWOBAL
* * (TX7(I) - TX8(I) - TX5(I) - QOL(I1(I),IC)*ALHL)
850 CONTINUE
DO 860 I=1,LENB
GMH(I,IC) = TX1(I) + PRI(I1(I),IC) * ONEBCP
* * (TX3(I)*(ALHL*HALF) + ETA(I,IC) * TX8(I))
860 CONTINUE
C
C CALCULATE HC PART OF AKM
C
IF (IC1 .LE. KM1) THEN
DO 870 I=1,LENB
TX1(I) = GMH(I,K)
870 CONTINUE
DO 3725 L=KM1,IC1,-1
DO 880 I=1,LENB
II = I1(I)
TX1(I) = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * GMH(I,L)
TX2(I) = GAM(II,L-1) * (PRJ(II,L) - PRH(II,L-1))
880 CONTINUE
C
IF (L .EQ. IC1) THEN
DO 890 I=LENA1,LENB
TX2(I) = ZERO
890 CONTINUE
ENDIF
DO 900 I=1,LENB
II = I1(I)
AKM(I) = AKM(I) + TX1(I) *
* (TX2(I) + GAM(II,L)*(PRH(II,L)-PRJ(II,L)))
900 CONTINUE
3725 CONTINUE
ENDIF
C
DO 920 I=LENA1,LENB
II = I1(I)
TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
* + 0.5*(PRS(II,IC+2) - PRS(II,IC)) * (ONE-TX6(I))
c
TX1(I) = PRS(II,IC1)
TX5(I) = 0.5 * (PRS(II,IC1) + PRS(II,IC+2))
C
IF ((TX2(I) .GE. TX1(I)) .AND. (TX2(I) .LT. TX5(I))) THEN
TX6(I) = ONE - (TX2(I) - TX1(I)) / (TX5(I) - TX1(I))
C
TEM = PRI(II,IC1) / PRI(II,IC)
HOL(I,IC1) = HOL(I,IC1) + HOL(I,IC) * TEM
HOL(I,IC) = ZERO
C
GMH(I,IC1) = GMH(I,IC1) + GMH(I,IC) * TEM
GMH(I,IC) = ZERO
ELSEIF (TX2(I) .LT. TX1(I)) THEN
TX6(I) = 1.0
ELSE
TX6(I) = 0.0
ENDIF
920 CONTINUE
C
C
DO I=1,LENC
PCU(I) = 0.0
ENDDO
DO 970 I=1,LENB
II = I1(I)
IF (AKM(I) .LT. ZERO .AND. WLQ(I) .GE. 0.0) THEN
WFN(I) = - TX6(I) * WFN(I) * RASALF / AKM(I)
ELSE
WFN(I) = ZERO
ENDIF
TEM = (PRS(II,K+1)-PRS(II,K))*(CMB2PA*FRAC)
WFN(I) = MIN(WFN(I), TEM)
C
C compute cloud amount
C
CC TX1(I) = CLN(II)
CC IF (WFN(I) .GT. CRTMSF) TX1(I) = TX1(I) + CLF(I)
CC IF (TX1(I) .GT. ONE) TX1(I) = ONE
C
C PRECIPITATION
C
PCU(II) = WLQ(I) * WFN(I) * ONEBG
C
C CUMULUS FRICTION AT THE BOTTOM LAYER
C
TX4(I) = WFN(I) * (1.0/ALHL)
TX5(I) = WFN(I) * ONEBCP
970 CONTINUE
C
C compute cloud mass flux for diagnostic output
C
DO L = IC,K
DO I=1,LENB
II = I1(I)
if(L.lt.K)then
CMASS(II,L) = ETA(I,L+1) * WFN(I) * ONEBG
else
CMASS(II,L) = WFN(I) * ONEBG
endif
ENDDO
ENDDO
C
CC DO 975 I=1,LENB
CC II = I1(I)
CC CLN(II) = TX1(I)
CC975 CONTINUE
C
C THETA AND Q CHANGE DUE TO CLOUD TYPE IC
C
c TEMA = 0.0
c TEMB = 0.0
DO 990 L=IC,K
DO 980 I=1,LENB
II = I1(I)
TEM = (GMH(I,L) - HOL(I,L)) * TX4(I)
TEM1 = HOL(I,L) * TX5(I)
C
TCU(II,L) = TEM1 / PRH(II,L)
QCU(II,L) = TEM
980 CONTINUE
c I = I1(IP1)
c
c TEM = (PRS(I,L+1)-PRS(I,L)) * (ONEBG*100.0)
c TEMA = TEMA + TCU(I,L) * PRH(I,L) * TEM * (CP/ALHL)
c TEMB = TEMB + QCU(I,L) * TEM
C
990 CONTINUE
C
c Compute Tracer Tendencies
c -------------------------
do nt = 1,ntracer
c Tracer Tendency at the Bottom Layer
c -----------------------------------
DO 995 I=1,LENB
II = I1(I)
TEM = half*TX5(I) * PRI(II,K)
TX1(I) = (UOI(II,KM1+nltop-1,nt) - UOI(II,K+nltop-1,nt))
ucu(II,K,nt) = TEM * TX1(I)
995 CONTINUE
c Tracer Tendency at all other Levels
c -----------------------------------
DO 1020 L=KM1,IC1,-1
DO 1010 I=1,LENB
II = I1(I)
TEM = half*TX5(I) * PRI(II,L)
TEM1 = TX1(I)
TX1(I) = (UOI(II,L-1+nltop-1,nt)-UOI(II,L+nltop-1,nt)) * ETA(I,L)
TX3(I) = (TX1(I) + TEM1) * TEM
1010 CONTINUE
DO 1020 I=1,LENB
II = I1(I)
ucu(II,L,nt) = TX3(I)
1020 CONTINUE
DO 1030 I=1,LENB
II = I1(I)
IF (TX6(I) .GE. 1.0) THEN
TEM = half*TX5(I) * PRI(II,IC)
ELSE
TEM = 0.0
ENDIF
TX1(I) = (TX1(I) + UHT(I,nt) + UHT(I,nt)) * TEM
1030 CONTINUE
DO 1040 I=1,LENB
II = I1(I)
ucu(II,IC,nt) = TX1(I)
1040 CONTINUE
enddo
C
C PENETRATIVE CONVECTION CALCULATION OVER
C
RETURN
END
SUBROUTINE RNCL(LEN, PL, RNO, CLF)
C
C*********************************************************************
C********************** Relaxed Arakawa-Schubert *********************
C************************ SUBROUTINE RNCL ************************
C**************************** 23 July 1992 ***************************
C*********************************************************************
implicit none
C Argument List declarations
integer len
real PL(LEN), RNO(LEN), CLF(LEN)
C Local Variables
real p5,p8,pt8,pt2,pfac,p4,p6,p7,p9,cucld,cfac
PARAMETER (P5=500.0, P8=800.0, PT8=0.8, PT2=0.2)
PARAMETER (PFAC=PT2/(P8-P5))
PARAMETER (P4=400.0, P6=401.0)
PARAMETER (P7=700.0, P9=900.0)
PARAMETER (CUCLD=0.5,CFAC=CUCLD/(P6-P4))
integer i
C
DO 10 I=1,LEN
rno(i) = 1.0
ccc if( pl(i).le.400.0 ) rno(i) = max( 0.75, 1.0-0.0025*(400.0-pl(i)) )
ccc IF ( PL(I).GE.P7 .AND. PL(I).LE.P9 ) THEN
ccc RNO(I) = ((P9-PL(I))/(P9-P7)) **2
ccc ELSE IF (PL(I).GT.P9) THEN
ccc RNO(I) = 0.
ccc ENDIF
CLF(I) = CUCLD
C
CARIESIF (PL(I) .GE. P5 .AND. PL(I) .LE. P8) THEN
CARIES RNO(I) = (P8-PL(I))*PFAC + PT8
CARIESELSEIF (PL(I) .GT. P8 ) THEN
CARIES RNO(I) = PT8
CARIESENDIF
CARIES
IF (PL(I) .GE. P4 .AND. PL(I) .LE. P6) THEN
CLF(I) = (P6-PL(I))*CFAC
ELSEIF (PL(I) .GT. P6 ) THEN
CLF(I) = 0.0
ENDIF
10 CONTINUE
C
RETURN
END
SUBROUTINE ACRITN ( LEN,PL,PLB,ACR )
C*********************************************************************
C********************** Relaxed Arakawa-Schubert *********************
C************************** SUBROUTINE ACRIT *********************
C****************** modified August 28, 1996 L.Takacs ************
C**** *****
C**** Note: Data obtained from January Mean After-Analysis *****
C**** from 4x5 46-layer GEOS Assimilation *****
C**** *****
C*********************************************************************
implicit none
C Argument List declarations
integer len
real PL(LEN), PLB(LEN), ACR(LEN)
C Local variables
integer lma
parameter (lma=18)
real p(lma)
real a(lma)
integer i,L
real temp
data p / 93.81, 111.65, 133.46, 157.80, 186.51,
. 219.88, 257.40, 301.21, 352.49, 409.76,
. 471.59, 535.04, 603.33, 672.79, 741.12,
. 812.52, 875.31, 930.20/
data a / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364,
. 1.01971, .65846, .45867, .38687, .31002,
. .25574, .20347, .17254, .15260, .16756,
. .09916, .10360, .05880/
do L=1,lma-1
do i=1,len
if( pl(i).ge.p(L) .and.
. pl(i).le.p(L+1)) then
temp = ( pl(i)-p(L) )/( p(L+1)-p(L) )
acr(i) = a(L+1)*temp + a(L)*(1-temp)
endif
enddo
enddo
do i=1,len
if( pl(i).lt.p(1) ) acr(i) = a(1)
if( pl(i).gt.p(lma) ) acr(i) = a(lma)
enddo
do i=1,len
acr(i) = acr(i) * (plb(i)-pl(i))
enddo
RETURN
END
SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DP,PLKE,
1 PLK,TH,TEMP1,TEMP2,TEMP3,ITMP1,ITMP2,RCON,RLAR,CLSBTH,tmscl,
2 tmfrc,cp,gravity,alhl,gamfac,cldlz,RHCRIT,offset,alpha)
implicit none
C Argument List declarations
integer nn,irun,nlay
real TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY),
. PL(IRUN,NLAY),CLFRAC(IRUN,NLAY),SP(IRUN),TEMP1(IRUN,NLAY),
. TEMP2(IRUN,NLAY),PLKE(IRUN,NLAY+1),
. RCON(IRUN),RLAR(IRUN),DP(IRUN,NLAY),PLK(IRUN,NLAY),TH(IRUN,NLAY),
. TEMP3(IRUN,NLAY)
integer ITMP1(IRUN,NLAY),ITMP2(IRUN,NLAY)
real CLSBTH(IRUN,NLAY)
real tmscl,tmfrc,cp,gravity,alhl,gamfac,offset,alpha
real cldlz(irun,nlay)
real rhcrit(irun,nlay)
C
C Local Variables
real zm1p04,zero,two89,zp44,zp01,half,zp578,one,thousand,z3600
real zp1,zp001
PARAMETER (ZM1P04 = -1.04E-4 )
PARAMETER (ZERO = 0.)
PARAMETER (TWO89= 2.89E-5)
PARAMETER ( ZP44= 0.44)
PARAMETER ( ZP01= 0.01)
PARAMETER ( ZP1 = 0.1 )
PARAMETER ( ZP001= 0.001)
PARAMETER ( HALF= 0.5)
PARAMETER ( ZP578 = 0.578 )
PARAMETER ( ONE = 1.)
PARAMETER ( THOUSAND = 1000.)
PARAMETER ( Z3600 = 3600.)
C
real EVP9(IRUN,NLAY)
real water(irun),crystal(irun)
real watevap(irun),iceevap(irun)
real fracwat,fracice, tice,rh,fact,dum
real rainmax(irun)
real getcon,rphf,elocp,cpog,relax
real exparg,arearat,rpow
integer i,L,n,nlaym1,irnlay,irnlm1
c Explicit Inline Directives
c --------------------------
#ifdef CRAY
#ifdef f77
cfpp$ expand (qsat)
#endif
#endif
tice = getcon('FREEZING-POINT')
fracwat = 0.70
fracice = 0.01
NLAYM1 = NLAY - 1
IRNLAY = IRUN*NLAY
IRNLM1 = IRUN*(NLAY-1)
RPHF = Z3600/tmscl
ELOCP = alhl/cp
CPOG = cp/gravity
DO I = 1,IRUN
RLAR(I) = 0.
water(i) = 0.
crystal(i) = 0.
ENDDO
do L = 1,nlay
do i = 1,irun
EVP9(i,L) = 0.
TEMP1(i,L) = 0.
TEMP2(i,L) = 0.
TEMP3(i,L) = 0.
CLSBTH(i,L) = 0.
cldlz(i,L) = 0.
enddo
enddo
C RHO(ZERO) / RHO FOR TERMINAL VELOCITY APPROX.
c ---------------------------------------------
DO L = 1,NLAY
DO I = 1,IRUN
TEMP2(I,L) = PL(I,L)*ZP001
TEMP2(I,L) = SQRT(TEMP2(I,L))
ENDDO
ENDDO
C INVERSE OF MASS IN EACH LAYER
c -----------------------------
DO L = 1,NLAY
DO I = 1,IRUN
TEMP3(I,L) = GRAVITY*ZP01 / DP(I,L)
ENDDO
ENDDO
C DO LOOP FOR MOISTURE EVAPORATION ABILITY AND CONVEC EVAPORATION.
c ----------------------------------------------------------------
DO 100 L=1,NLAY
DO I = 1,IRUN
TEMP1(I,3) = TL(I,L)
TEMP1(I,4) = QL(I,L)
ENDDO
DO 50 N=1,2
IF(N.EQ.1)RELAX=HALF
IF(N.GT.1)RELAX=ONE
DO I = 1,IRUN
call qsat ( temp1(i,3),pl(i,L),temp1(i,2),temp1(i,6),.true. )
TEMP1(I,5)=TEMP1(I,2)-TEMP1(I,4)
TEMP1(I,6)=TEMP1(I,6)*ELOCP
TEMP1(I,5)=TEMP1(I,5)/(ONE+TEMP1(I,6))
TEMP1(I,4)=TEMP1(I,4)+TEMP1(I,5)*RELAX
TEMP1(I,3)=TEMP1(I,3)-TEMP1(I,5)*ELOCP*RELAX
ENDDO
50 CONTINUE
DO I = 1,IRUN
EVP9(I,L) = (TEMP1(I,4) - QL(I,L))/TEMP3(I,L)
C convective detrained water
cldlz(i,L) = rain(i,L)*temp3(i,L)
if( tl(i,L).gt.tice-20.) then
water(i) = water(i) + rain(i,L)
else
crystal(i) = crystal(i) + rain(i,L)
endif
ENDDO
C**********************************************************************
C FOR CONVECTIVE PRECIP, FIND THE "EVAPORATION EFFICIENCY" USING *
C KESSLERS PARAMETERIZATION *
C**********************************************************************
DO 20 I=1,IRUN
iceevap(i) = 0.
watevap(i) = 0.
if( (evp9(i,L).gt.0.) .and. (crystal(i).gt.0.) ) then
iceevap(I) = EVP9(I,L)*fracice
IF(iceevap(i).GE.crystal(i)) iceevap(i) = crystal(i)
EVP9(I,L)=EVP9(I,L)-iceevap(I)
crystal(i) = crystal(i) - iceevap(i)
endif
C and now warm precipitate
if( (evp9(i,L).gt.0.) .and. (water(i).gt.0.) ) then
exparg = ZM1P04*tmscl*((water(i)*RPHF*TEMP2(I,L))**ZP578)
AREARAT = ONE-(EXP(EXPARG))
watevap(I) = EVP9(I,L)*AREARAT*fracwat
IF(watevap(I).GE.water(i)) watevap(I) = water(i)
EVP9(I,L)=EVP9(I,L)-watevap(I)
water(i) = water(i) - watevap(i)
endif
QL(I,L) = QL(I,L)+(iceevap(i)+watevap(i))*TEMP3(I,L)
TL(I,L) = TL(I,L)-(iceevap(i)+watevap(i))*TEMP3(I,L)*ELOCP
20 CONTINUE
100 CONTINUE
do i = 1,irun
rcon(i) = water(i) + crystal(i)
enddo
C**********************************************************************
C Large Scale Precip
C**********************************************************************
DO 200 L=1,NLAY
DO I = 1,IRUN
rainmax(i) = rhcrit(i,L)*evp9(i,L) +
. ql(i,L)*(rhcrit(i,L)-1.)/temp3(i,L)
if (rainmax(i).LE.0.0) then
call qsat( tl(i,L),pl(i,L),rh,dum,.false.)
rh = ql(i,L)/rh
if( rhcrit(i,L).eq.1.0 ) then
fact = 1.0
else
fact = min( 1.0, alpha + (1.0-alpha)*( rh-rhcrit(i,L)) /
1 (1.0-rhcrit(i,L)) )
endif
C Do not allow clouds above 10 mb
if( pl(i,L).ge.10.0 ) CLSBTH(I,L) = fact
RLAR(I) = RLAR(I)-rainmax(I)
QL(I,L) = QL(I,L)+rainmax(I)*TEMP3(I,L)
TL(I,L) = TL(I,L)-rainmax(I)*TEMP3(I,L)*ELOCP
C Large-scale water
cldlz(i,L) = cldlz(i,L) - rainmax(i)*temp3(i,L)
ENDIF
ENDDO
DO I=1,IRUN
IF((RLAR(I).GT.0.0).AND.(rainmax(I).GT.0.0))THEN
RPOW=(RLAR(I)*RPHF*TEMP2(I,L))**ZP578
EXPARG = ZM1P04*tmscl*RPOW
AREARAT = ONE-(EXP(EXPARG))
TEMP1(I,7) = rainmax(I)*AREARAT
IF(TEMP1(I,7).GE.RLAR(I)) TEMP1(I,7) = RLAR(I)
RLAR(I) = RLAR(I)-TEMP1(I,7)
QL(I,L) = QL(I,L)+TEMP1(I,7)*TEMP3(I,L)
TL(I,L) = TL(I,L)-TEMP1(I,7)*TEMP3(I,L)*ELOCP
ENDIF
ENDDO
200 CONTINUE
RETURN
END
subroutine srclouds (th,q,plk,pl,plke,cloud,cldwat,irun,irise,
1 rhc,offset,alpha)
C***********************************************************************
C
C PURPOSE:
C ========
C Compute non-precipitating cloud fractions
C based on Slingo and Ritter (1985).
C Remove cloudiness where conditionally unstable.
C
C INPUT:
C ======
C th ......... Potential Temperature (irun,irise)
C q .......... Specific Humidity (irun,irise)
C plk ........ P**Kappa at mid-layer (irun,irise)
C pl ......... Pressure at mid-layer (irun,irise)
C plke ....... P**Kappa at edge (irun,irise+1)
C irun ....... Horizontal dimension
C irise ...... Vertical dimension
C
C OUTPUT:
C =======
C cloud ...... Cloud Fraction (irun,irise)
C
C***********************************************************************
implicit none
integer irun,irise
real th(irun,irise)
real q(irun,irise)
real plk(irun,irise)
real pl(irun,irise)
real plke(irun,irise+1)
real cloud(irun,irise)
real cldwat(irun,irise)
real qs(irun,irise)
real cp, alhl, getcon, akap
real ratio, temp, elocp
real rhcrit,rh,dum
integer i,L
real rhc(irun,irise)
real offset,alpha
c Explicit Inline Directives
c --------------------------
#ifdef CRAY
#ifdef f77
cfpp$ expand (qsat)
#endif
#endif
cp = getcon('CP')
alhl = getcon('LATENT HEAT COND')
elocp = alhl/cp
akap = getcon('KAPPA')
do L = 1,irise
do i = 1,irun
temp = th(i,L)*plk(i,L)
call qsat ( temp,pl(i,L),qs(i,L),dum,.false. )
enddo
enddo
do L = 2,irise
do i = 1,irun
rh = q(i,L)/qs(i,L)
rhcrit = rhc(i,L) - offset
ratio = alpha*(rh-rhcrit)/offset
if(cloud(i,L).eq. 0.0 .and. ratio.gt.0.0 ) then
cloud(i,L) = min( ratio,1.0 )
endif
enddo
enddo
c Reduce clouds from conditionally unstable layer
c -----------------------------------------------
call ctei ( th,q,cloud,cldwat,pl,plk,plke,irun,irise )
return
end
subroutine ctei ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm )
implicit none
integer im,lm
real th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm)
real plk(im,lm),pl(im,lm),cldfrc(im,lm)
integer i,L
real getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq
real k,krd,kmm,f
cp = getcon('CP')
alhl = getcon('LATENT HEAT COND')
cpoel = cp/alhl
elocp = alhl/cp
do L=lm,2,-1
do i=1,im
dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1)
if( dq.eq.0.0 ) dq = 1.0e-20
k = 1.0 + cpoel*plke(i,L)*( th(i,L)-th(i,L-1) ) / dq
t = th(i,L)*plk(i,L)
p = pl(i,L)
call qsat ( t,p,qs,dqsdt,.true. )
krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t )
kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t )
. / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt )
s = ( (k-krd)/(kmm-krd) )
f = 1.0 - min( 1.0, max(0.0,1.0-exp(-s)) )
cldfrc(i,L) = cldfrc(i,L)*f
cldwat(i,L) = cldwat(i,L)*f
enddo
enddo
return
end
subroutine back2grd(gathered,indeces,scattered,irun)
implicit none
integer i,irun,indeces(irun)
real gathered(irun),scattered(irun)
real temp(irun)
do i = 1,irun
temp(indeces(i)) = gathered(i)
enddo
do i = 1,irun
scattered(i) = temp(i)
enddo
return
end