C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.35 2005/06/17 16:51:24 ce107 Exp $
C $Name: $
#include "FIZHI_OPTIONS.h"
subroutine MOISTIO (ndmoist,istrip,npcs,
. lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup,
. pz,plz,plze,dpres,pkht,pkl,uz,vz,tz,qz,bi,bj,ntracerin,ptracer,
. qqz,dumoist,dvmoist,dtmoist,dqmoist,cumfric,
. 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
c Input Variables
c ---------------
integer im,jm,lm
integer ndmoist,istrip,npcs
integer bi,bj,ntracerin,ptracer
integer lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup
_RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
_RL pkht(im,jm,lm+1),pkl(im,jm,lm)
_RL tz(im,jm,lm),qz(im,jm,lm,ntracerin)
_RL uz(im,jm,lm),vz(im,jm,lm)
_RL qqz(im,jm,lm)
_RL dumoist(im,jm,lm),dvmoist(im,jm,lm)
_RL dtmoist(im,jm,lm),dqmoist(im,jm,lm,ntracerin)
logical cumfric
_RL ptop
integer iras
_RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
integer nswcld,nswlz
_RL cldlsp_sw(im,jm,lm),cldras_sw(im,jm,lm)
_RL cldtot_sw(im,jm,lm),swlz(im,jm,lm)
integer nlwcld,nlwlz
_RL cldlsp_lw(im,jm,lm),cldras_lw(im,jm,lm)
_RL cldtot_lw(im,jm,lm),lwlz(im,jm,lm)
logical lpnt
integer myid
c Local Variables
c ---------------
integer ncrnd,nsecf
_RL fracqq, dum
integer snowcrit
parameter (fracqq = 0.1)
_RL one
parameter (one=1.)
_RL cldsr(im,jm,lm)
_RL srcld(istrip,lm)
_RL plev
_RL cldnow,cldlsp_mem,cldlsp,cldras_mem,cldras
_RL watnow,watmin,cldmin
_RL cldprs(im,jm),cldtmp(im,jm)
_RL cldhi (im,jm),cldlow(im,jm)
_RL cldmid(im,jm),totcld(im,jm)
_RL CLDLS(im,jm,lm) , CPEN(im,jm,lm)
_RL tmpimjm(im,jm)
_RL lsp_new(im,jm)
_RL conv_new(im,jm)
_RL snow_new(im,jm)
_RL qqcolmin(im,jm)
_RL qqcolmax(im,jm)
integer levpbl(im,jm)
c Gathered Arrays for Variable Cloud Base
c ---------------------------------------
_RL raincgath(im*jm)
_RL pigather(im*jm)
_RL thgather(im*jm,lm)
_RL shgather(im*jm,lm)
_RL pkzgather(im*jm,lm)
_RL pkegather(im*jm,lm+1)
_RL plzgather(im*jm,lm)
_RL plegather(im*jm,lm+1)
_RL dpgather(im*jm,lm)
_RL tmpgather(im*jm,lm)
_RL deltgather(im*jm,lm)
_RL delqgather(im*jm,lm)
_RL ugather(im*jm,lm,ntracerin+2-ptracer)
_RL delugather(im*jm,lm,ntracerin+2-ptracer)
_RL deltrnev(im*jm,lm)
_RL delqrnev(im*jm,lm)
integer nindeces(lm)
integer pblindex(im*jm)
integer levgather(im*jm)
c Stripped Arrays
c ---------------
_RL saveth (istrip,lm)
_RL saveq (istrip,lm)
_RL saveu (istrip,lm,ntracerin+2-ptracer)
_RL usubcl (istrip, ntracerin+2-ptracer)
_RL ple(istrip,lm+1)
_RL dp(istrip,lm)
_RL TL(ISTRIP,lm) , SHL(ISTRIP,lm)
_RL PL(ISTRIP,lm) , PLK(ISTRIP,lm)
_RL PLKE(ISTRIP,lm+1)
_RL TH(ISTRIP,lm) ,CVTH(ISTRIP,lm)
_RL CVQ(ISTRIP,lm)
_RL UL(ISTRIP,lm,ntracerin+2-ptracer)
_RL cvu(istrip,lm,ntracerin+2-ptracer)
_RL CLMAXO(ISTRIP,lm),CLBOTH(ISTRIP,lm)
_RL CLSBTH(ISTRIP,lm)
_RL TMP1(ISTRIP,lm), TMP2(ISTRIP,lm)
_RL TMP3(ISTRIP,lm), TMP4(ISTRIP,lm+1)
_RL TMP5(ISTRIP,lm+1)
integer ITMP1(ISTRIP,lm), ITMP2(ISTRIP,lm)
_RL PRECIP(ISTRIP), PCNET(ISTRIP)
_RL SP(ISTRIP), PREP(ISTRIP)
_RL PCPEN (ISTRIP,lm)
integer pbl(istrip),depths(lm)
_RL cldlz(istrip,lm), cldwater(im,jm,lm)
_RL rhfrac(istrip), rhmin, pup, ppbl, rhcrit(istrip,lm)
_RL offset, alpha, rasmax
logical first
logical lras
_RL clfrac (istrip,lm)
_RL cldmas (istrip,lm)
_RL detrain(istrip,lm)
_RL psubcld (istrip), psubcldg (im,jm)
_RL psubcld_cnt(istrip), psubcldgc(im,jm)
_RL rnd(lm/2)
DATA FIRST /.TRUE./
integer imstp,nsubcl,nlras
integer i,j,iloop,indx,indgath,l,nn,num,numdeps,nt
_RL tmstp,tminv,sday,grav,alhl,cp,elocp,gamfac
_RL rkappa,p0kappa,p0kinv,ptopkap,pcheck
_RL tice,getcon,pi
integer ntracer,ntracedim, ntracex
#ifdef ALLOW_DIAGNOSTICS
logical diagnostics_is_on
external
_RL tmpdiag(im,jm),tmpdiag2(im,jm)
#endif
C **********************************************************************
C **** INITIALIZATION ****
C **********************************************************************
C Add U and V components to tracer array for cumulus friction
if(cumfric) then
ntracer = ntracerin + 2
else
ntracer = ntracerin
endif
ntracedim= max(ntracer-ptracer,1)
ntracex= ntracer-ptracer
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.1 .and. bi.eq.1 ) 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 ----------------------------------------------------------------
indx = 0
do L = nsubmin,nltop,-1
do j = 1,jm
do i = 1,im
if(levpbl(i,j).eq.L) then
indx = indx + 1
pblindex(indx) = (j-1)*im + i
endif
enddo
enddo
enddo
do indx = 1,im*jm
levgather(indx) = levpbl(pblindex(indx),1)
pigather(indx) = pz(pblindex(indx),1)
pkegather(indx,lm+1) = pkht(pblindex(indx),1,lm+1)
plegather(indx,lm+1) = plze(pblindex(indx),1,lm+1)
enddo
do L = 1,lm
do indx = 1,im*jm
thgather(indx,L) = tz(pblindex(indx),1,L)
shgather(indx,L) = qz(pblindex(indx),1,L,1)
pkegather(indx,L) = pkht(pblindex(indx),1,L)
pkzgather(indx,L) = pkl(pblindex(indx),1,L)
plegather(indx,L) = plze(pblindex(indx),1,L)
plzgather(indx,L) = plz(pblindex(indx),1,L)
dpgather(indx,L) = dpres(pblindex(indx),1,L)
enddo
enddo
C General Tracers
C----------------
do nt = 1,ntracerin-ptracer
do L = 1,lm
do indx = 1,im*jm
ugather(indx,L,nt) = qz(pblindex(indx),1,L,nt+ptracer)
enddo
enddo
enddo
if(cumfric) then
C Cumulus Friction - load u and v wind components into tracer array
C------------------------------------------------------------------
do L = 1,lm
do indx = 1,im*jm
ugather(indx,L,ntracerin-ptracer+1) = uz(pblindex(indx),1,L)
ugather(indx,L,ntracerin-ptracer+2) = vz(pblindex(indx),1,L)
enddo
enddo
endif
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
dumoist(i,j,l) = 0.
dvmoist(i,j,l) = 0.
dtmoist(i,j,l) = 0.
do nt = 1,ntracerin
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),ntracedim,ntracex,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 **********************************************************************
call PASTE( cldmas,tmpgather,istrip,im*jm,lm,NN)
call PASTE(detrain,pkegather,istrip,im*jm,lm,NN)
call PASTE(psubcld ,psubcldg ,istrip,im*jm,1,NN)
call PASTE(psubcld_cnt,psubcldgc,istrip,im*jm,1,NN)
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,one,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 ------------------------------------------------------------
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)
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)
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)
CALL PASTE (pcnet,raincgath,ISTRIP,im*jm,1,NN)
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)
call BACK2GRD (raincgath,pblindex,raincgath,im*jm)
c Subcloud Layer Pressure
c -----------------------
call BACK2GRD (psubcldg ,pblindex,psubcldg ,im*jm)
call BACK2GRD (psubcldgc,pblindex,psubcldgc,im*jm)
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