/[MITgcm]/MITgcm/pkg/fizhi/fizhi_moist.F
ViewVC logotype

Annotation of /MITgcm/pkg/fizhi/fizhi_moist.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.16 - (hide annotations) (download)
Tue Aug 10 15:13:47 2004 UTC (19 years, 10 months ago) by molod
Branch: MAIN
Changes since 1.15: +129 -103 lines 
Debugging
1 molod 1.16 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.15 2004/07/28 01:25:07 molod Exp $
2 molod 1.1 C $Name: $
3 molod 1.2
4 molod 1.13 #include "FIZHI_OPTIONS.h"
5 molod 1.5 subroutine moistio (ndmoist,istrip,npcs,
6 molod 1.4 . lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup,
7 molod 1.6 . pz,plz,plze,dpres,pkht,pkl,tz,qz,bi,bj,ntracer,ptracer,
8 molod 1.5 . qqz,dumoist,dvmoist,dtmoist,dqmoist,
9 molod 1.4 . im,jm,lm,ptop,
10 molod 1.1 . iras,rainlsp,rainconv,snowfall,
11     . nswcld,cldtot_sw,cldras_sw,cldlsp_sw,nswlz,swlz,
12     . nlwcld,cldtot_lw,cldras_lw,cldlsp_lw,nlwlz,lwlz,
13 molod 1.2 . lpnt,myid)
14 molod 1.1
15 molod 1.8 implicit none
16    
17 molod 1.2 #ifdef ALLOW_DIAGNOSTICS
18 molod 1.8 #include "SIZE.h"
19     #include "diagnostics_SIZE.h"
20 molod 1.1 #include "diagnostics.h"
21 molod 1.2 #endif
22 molod 1.1
23     c Input Variables
24     c ---------------
25 molod 1.8 integer im,jm,lm
26 molod 1.5 integer ndmoist,istrip,npcs
27 molod 1.8 integer bi,bj,ntracer,ptracer
28 molod 1.4 integer lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup
29 molod 1.13 _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
30     _RL pkht(im,jm,lm+1),pkl(im,jm,lm)
31     _RL tz(im,jm,lm),qz(im,jm,lm,ntracer)
32     _RL qqz(im,jm,lm)
33     _RL dumoist(im,jm,lm),dvmoist(im,jm,lm)
34     _RL dtmoist(im,jm,lm),dqmoist(im,jm,lm,ntracer)
35     _RL ptop
36 molod 1.5 integer iras
37 molod 1.13 _RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
38 molod 1.5 integer nswcld,nswlz
39 molod 1.13 _RL cldlsp_sw(im,jm,lm),cldras_sw(im,jm,lm)
40     _RL cldtot_sw(im,jm,lm),swlz(im,jm,lm)
41 molod 1.5 integer nlwcld,nlwlz
42 molod 1.13 _RL cldlsp_lw(im,jm,lm),cldras_lw(im,jm,lm)
43     _RL cldtot_lw(im,jm,lm),lwlz(im,jm,lm)
44 molod 1.5 logical lpnt
45     integer myid
46 molod 1.1
47     c Local Variables
48     c ---------------
49 molod 1.4 integer ncrnd,nsecf
50 molod 1.1
51 molod 1.13 _RL fracqq, dum
52 molod 1.4 integer snowcrit
53 molod 1.1 parameter (fracqq = 0.1)
54 molod 1.14 _RL one
55     parameter (one=1.)
56 molod 1.1
57 molod 1.13 _RL cldsr(im,jm,lm)
58     _RL srcld(istrip,lm)
59 molod 1.1
60 molod 1.13 _RL plev
61     _RL cldnow,cldlsp_mem,cldlsp,cldras_mem,cldras
62     _RL watnow,watmin,cldmin
63     _RL cldprs(im,jm),cldtmp(im,jm)
64     _RL cldhi (im,jm),cldlow(im,jm)
65     _RL cldmid(im,jm),totcld(im,jm)
66    
67     _RL CLDLS(im,jm,lm) , CPEN(im,jm,lm)
68     _RL tmpimjm(im,jm)
69     _RL lsp_new(im,jm)
70     _RL conv_new(im,jm)
71     _RL snow_new(im,jm)
72 molod 1.1
73 molod 1.13 _RL qqcolmin(im,jm)
74     _RL qqcolmax(im,jm)
75 molod 1.1 integer levpbl(im,jm)
76    
77     c Gathered Arrays for Variable Cloud Base
78     c ---------------------------------------
79 molod 1.13 _RL raincgath(im*jm)
80     _RL pigather(im*jm)
81     _RL thgather(im*jm,lm)
82     _RL shgather(im*jm,lm)
83     _RL pkzgather(im*jm,lm)
84     _RL pkegather(im*jm,lm+1)
85     _RL plzgather(im*jm,lm)
86     _RL plegather(im*jm,lm+1)
87     _RL dpgather(im*jm,lm)
88     _RL tmpgather(im*jm,lm)
89     _RL deltgather(im*jm,lm)
90     _RL delqgather(im*jm,lm)
91     _RL ugather(im*jm,lm,ntracer)
92     _RL delugather(im*jm,lm,ntracer)
93     _RL deltrnev(im*jm,lm)
94     _RL delqrnev(im*jm,lm)
95 molod 1.1
96     integer nindeces(lm)
97     integer pblindex(im*jm)
98     integer levgather(im*jm)
99    
100     c Stripped Arrays
101     c ---------------
102 molod 1.13 _RL saveth (istrip,lm)
103     _RL saveq (istrip,lm)
104     _RL saveu (istrip,lm,ntracer)
105     _RL usubcl (istrip, ntracer)
106    
107     _RL ple(istrip,lm+1)
108     _RL dp(istrip,lm)
109     _RL TL(ISTRIP,lm) , SHL(ISTRIP,lm)
110     _RL PL(ISTRIP,lm) , PLK(ISTRIP,lm)
111     _RL PLKE(ISTRIP,lm+1)
112     _RL TH(ISTRIP,lm) ,CVTH(ISTRIP,lm)
113     _RL CVQ(ISTRIP,lm)
114     _RL UL(ISTRIP,lm,ntracer)
115     _RL cvu(istrip,lm,ntracer)
116     _RL CLMAXO(ISTRIP,lm),CLBOTH(ISTRIP,lm)
117     _RL CLSBTH(ISTRIP,lm)
118     _RL TMP1(ISTRIP,lm), TMP2(ISTRIP,lm)
119     _RL TMP3(ISTRIP,lm), TMP4(ISTRIP,lm+1)
120     _RL TMP5(ISTRIP,lm+1)
121 molod 1.1 integer ITMP1(ISTRIP,lm), ITMP2(ISTRIP,lm)
122    
123 molod 1.13 _RL PRECIP(ISTRIP), PCNET(ISTRIP)
124     _RL SP(ISTRIP), PREP(ISTRIP)
125     _RL PCPEN (ISTRIP,lm)
126 molod 1.1 integer pbl(istrip),depths(lm)
127    
128 molod 1.13 _RL cldlz(istrip,lm), cldwater(im,jm,lm)
129     _RL rhfrac(istrip), rhmin, pup, ppbl, rhcrit(istrip,lm)
130     _RL offset, alpha, rasmax
131 molod 1.1
132     logical first
133     logical lras
134 molod 1.13 _RL clfrac (istrip,lm)
135     _RL cldmas (istrip,lm)
136     _RL detrain(istrip,lm)
137     _RL psubcld (istrip), psubcldg (im,jm)
138     _RL psubcld_cnt(istrip), psubcldgc(im,jm)
139     _RL rnd(lm/2)
140 molod 1.1 DATA FIRST /.TRUE./
141    
142 molod 1.4 integer imstp,nsubcl,nlras
143 molod 1.16 integer i,j,iloop,indx,indgath,l,nn,num,numdeps,nt
144 molod 1.13 _RL tmstp,tminv,sday,grav,alhl,cp,elocp,gamfac
145     _RL rkappa,p0kappa,p0kinv,ptopkap,pcheck
146     _RL tice,getcon,pi
147 molod 1.1
148     C **********************************************************************
149     C **** INITIALIZATION ****
150     C **********************************************************************
151    
152     IMSTP = nsecf(NDMOIST)
153     TMSTP = FLOAT(IMSTP)
154     TMINV = 1. / TMSTP
155    
156     C Minimum Large-Scale Cloud Fraction at rhcrit
157     alpha = 0.80
158 molod 1.2 C Difference in fraction between SR and LS Threshold
159 molod 1.1 offset = 0.10
160 molod 1.2 C Large-Scale Relative Humidity Threshold in PBL
161 molod 1.1 rhmin = 0.90
162     C Maximum Cloud Fraction associated with RAS
163     rasmax = 1.00
164    
165     nn = 3*3600.0/tmstp + 1
166     C Threshold for Cloud Fraction Memory
167     cldmin = rasmax*(1.0-tmstp/3600.)**nn
168     C Threshold for Cloud Liquid Water Memory
169     watmin = 1.0e-8
170    
171     SDAY = GETCON('SDAY')
172     GRAV = GETCON('GRAVITY')
173     ALHL = GETCON('LATENT HEAT COND')
174     CP = GETCON('CP')
175     ELOCP = GETCON('LATENT HEAT COND') / GETCON('CP')
176     GAMFAC = GETCON('LATENT HEAT COND') * GETCON('EPS') * ELOCP
177     . / GETCON('RGAS')
178     RKAPPA = GETCON('KAPPA')
179     P0KAPPA = 1000.0**RKAPPA
180     P0KINV = 1. / P0KAPPA
181     PTOPKAP = PTOP**RKAPPA
182     tice = getcon('FREEZING-POINT')
183     PI = 4.*atan(1.)
184    
185 molod 1.4 c Determine Total number of Random Clouds to Check
186 molod 1.1 c ---------------------------------------------
187     ncrnd = (lm-nltop+1)/2
188    
189 molod 1.15 if(first .and. myid.eq.1) then
190 molod 1.1 print *
191 molod 1.4 print *,'Top Level Allowed for Convection : ',nltop
192     print *,' Highest Sub-Cloud Level: ',nsubmax
193     print *,' Lowest Sub-Cloud Level: ',nsubmin
194 molod 1.1 print *,' Total Number of Random Clouds: ',ncrnd
195     print *
196     first = .false.
197     endif
198    
199     c And now find PBL depth - the level where qq = fracqq * qq at surface
200     c --------------------------------------------------------------------
201     do j = 1,jm
202     do i = 1,im
203     qqcolmin(i,j) = qqz(i,j,lm)*fracqq
204     qqcolmax(i,j) = qqz(i,j,lm)
205     levpbl(i,j) = lm+1
206     enddo
207     enddo
208    
209     DO L = lm-1,1,-1
210     DO j = 1,jm
211     DO i = 1,im
212     IF((qqz(i,j,l).gt.qqcolmax(i,j))
213     1 .and.(levpbl(i,j).eq.lm+1))then
214     qqcolmax(i,j) = qqz(i,j,l)
215     qqcolmin(i,j) = fracqq*qqcolmax(i,j)
216     endif
217     if((qqz(i,j,l).lt.qqcolmin(i,j))
218     1 .and.(levpbl(i,j).eq.lm+1))levpbl(i,j)=L+1
219     enddo
220     enddo
221     enddo
222    
223     do j = 1,jm
224     do i = 1,im
225     if(levpbl(i,j).gt.nsubmin) levpbl(i,j) = nsubmin
226     if(levpbl(i,j).lt.nsubmax) levpbl(i,j) = nsubmax
227     enddo
228     enddo
229    
230    
231     c Set up the array of indeces of subcloud levels for the gathering
232     c ----------------------------------------------------------------
233 molod 1.16 indx = 0
234 molod 1.1 do L = nsubmin,nltop,-1
235     do j = 1,jm
236     do i = 1,im
237     if(levpbl(i,j).eq.L) then
238 molod 1.16 indx = indx + 1
239     pblindex(indx) = (j-1)*im + i
240 molod 1.1 endif
241     enddo
242     enddo
243     enddo
244    
245 molod 1.16 do indx = 1,im*jm
246     levgather(indx) = levpbl(pblindex(indx),1)
247     pigather(indx) = pz(pblindex(indx),1)
248     pkegather(indx,lm+1) = pkht(pblindex(indx),1,lm+1)
249     plegather(indx,lm+1) = plze(pblindex(indx),1,lm+1)
250 molod 1.1 enddo
251    
252     do L = 1,lm
253 molod 1.16 do indx = 1,im*jm
254     thgather(indx,L) = tz(pblindex(indx),1,L)
255     shgather(indx,L) = qz(pblindex(indx),1,L,1)
256     pkegather(indx,L) = pkht(pblindex(indx),1,L)
257     pkzgather(indx,L) = pkl(pblindex(indx),1,L)
258     plegather(indx,L) = plze(pblindex(indx),1,L)
259     plzgather(indx,L) = plz(pblindex(indx),1,L)
260     dpgather(indx,L) = dpres(pblindex(indx),1,L)
261 molod 1.1 enddo
262     enddo
263     do nt = 1,ntracer-ptracer
264     do L = 1,lm
265 molod 1.16 do indx = 1,im*jm
266     ugather(indx,L,nt) = qz(pblindex(indx),1,L,nt+ptracer)
267 molod 1.1 enddo
268     enddo
269     enddo
270    
271     c bump the counter for number of calls to convection
272     c --------------------------------------------------
273     iras = iras + 1
274     if( iras.ge.1e9 ) iras = 1
275    
276     c select the 'random' cloud detrainment levels for RAS
277     c ----------------------------------------------------
278     call rndcloud(iras,ncrnd,rnd,myid)
279    
280     do l=1,lm
281     do j=1,jm
282     do i=1,im
283     dtmoist(i,j,l) = 0.
284     do nt = 1,ntracer
285     dqmoist(i,j,l,nt) = 0.
286     enddo
287     enddo
288     enddo
289     enddo
290    
291     C***********************************************************************
292     C **** LOOP OVER NPCS PEICES ****
293     C **********************************************************************
294    
295     DO 1000 NN = 1,NPCS
296    
297     C **********************************************************************
298     C **** VARIABLE INITIALIZATION ****
299     C **********************************************************************
300    
301     CALL STRIP ( pigather, SP ,im*jm,ISTRIP,1 ,NN )
302     CALL STRIP ( pkzgather, PLK ,im*jm,ISTRIP,lm,NN )
303 molod 1.6 CALL STRIP ( pkegather, PLKE ,im*jm,ISTRIP,lm+1,NN )
304     CALL STRIP ( plzgather, PL ,im*jm,ISTRIP,lm,NN )
305     CALL STRIP ( plegather, PLE ,im*jm,ISTRIP,lm+1,NN )
306     CALL STRIP ( dpgather, dp ,im*jm,ISTRIP,lm,NN )
307 molod 1.1 CALL STRIP ( thgather, TH ,im*jm,ISTRIP,lm,NN )
308     CALL STRIP ( shgather, SHL ,im*jm,ISTRIP,lm,NN )
309     CALL STRINT( levgather, pbl ,im*jm,ISTRIP,1 ,NN )
310    
311     do nt = 1,ntracer-ptracer
312     call strip ( ugather(1,1,nt), ul(1,1,nt),im*jm,istrip,lm,nn )
313     enddo
314    
315     C **********************************************************************
316     C **** SETUP FOR RAS CUMULUS PARAMETERIZATION ****
317     C **********************************************************************
318    
319     DO L = 1,lm
320     DO I = 1,ISTRIP
321     TH(I,L) = TH(I,L) * P0KAPPA
322     CLMAXO(I,L) = 0.
323     CLBOTH(I,L) = 0.
324     cldmas(I,L) = 0.
325     detrain(I,L) = 0.
326     ENDDO
327     ENDDO
328    
329     do L = 1,lm
330     depths(L) = 0
331     enddo
332    
333     numdeps = 0
334     do L = nsubmin,nltop,-1
335     nindeces(L) = 0
336     do i = 1,istrip
337     if(pbl(i).eq.L) nindeces(L) = nindeces(L) + 1
338     enddo
339     if(nindeces(L).gt.0) then
340     numdeps = numdeps + 1
341     depths(numdeps) = L
342     endif
343     enddo
344    
345    
346     C Initiate a do-loop around RAS for the number of different
347     C sub-cloud layer depths in this strip
348     C --If all subcloud depths are the same, execute loop once
349     C Otherwise loop over different subcloud layer depths
350    
351     num = 1
352     DO iloop = 1,numdeps
353    
354     nsubcl = depths(iloop)
355    
356     c Compute sub-cloud values for Temperature and Spec.Hum.
357     c ------------------------------------------------------
358     DO 600 I=num,num+nindeces(nsubcl)-1
359     TMP1(I,2) = 0.
360     TMP1(I,3) = 0.
361     600 CONTINUE
362    
363     NLRAS = NSUBCL - NLTOP + 1
364     DO 601 L=NSUBCL,lm
365     DO 602 I=num,num+nindeces(nsubcl)-1
366     TMP1(I,2) = TMP1(I,2) + (PLE(I,L+1)-PLE(I,L))*TH (I,L)/sp(i)
367     TMP1(I,3) = TMP1(I,3) + (PLE(I,L+1)-PLE(I,L))*SHL(I,L)/sp(i)
368     602 CONTINUE
369     601 CONTINUE
370     DO 603 I=num,num+nindeces(nsubcl)-1
371     TMP1(I,4) = 1. / ( (PLE(I,lm+1)-PLE(I,NSUBCL))/sp(I) )
372     TH(I,NSUBCL) = TMP1(I,2)*TMP1(I,4)
373     SHL(I,NSUBCL) = TMP1(I,3)*TMP1(I,4)
374     603 CONTINUE
375    
376     c Save initial value of tracers and compute sub-cloud value
377     c ---------------------------------------------------------
378     DO NT = 1,ntracer-ptracer
379     do L = 1,lm
380     do i = num,num+nindeces(nsubcl)-1
381     saveu(i,L,nt) = ul(i,L,nt)
382     enddo
383     enddo
384     DO I=num,num+nindeces(nsubcl)-1
385     TMP1(I,2) = 0.
386     ENDDO
387     DO L=NSUBCL,lm
388     DO I=num,num+nindeces(nsubcl)-1
389     TMP1(I,2) = TMP1(I,2)+(PLE(I,L+1)-PLE(I,L))*UL(I,L,NT)/sp(i)
390     ENDDO
391     ENDDO
392     DO I=num,num+nindeces(nsubcl)-1
393     UL(I,NSUBCL,NT) = TMP1(I,2)*TMP1(I,4)
394     usubcl(i,nt) = ul(i,nsubcl,nt)
395     ENDDO
396     ENDDO
397    
398     c Compute Pressure Arrays for RAS
399     c -------------------------------
400     DO 111 L=1,lm
401     DO 112 I=num,num+nindeces(nsubcl)-1
402     TMP4(I,L) = PLE(I,L)
403     112 CONTINUE
404     111 CONTINUE
405     DO I=num,num+nindeces(nsubcl)-1
406     TMP5(I,1) = PTOPKAP / P0KAPPA
407     ENDDO
408     DO L=2,lm
409     DO I=num,num+nindeces(nsubcl)-1
410 molod 1.6 TMP5(I,L) = PLKE(I,L)*P0KINV
411 molod 1.1 ENDDO
412     ENDDO
413     DO I=num,num+nindeces(nsubcl)-1
414     TMP4(I,lm+1) = PLE (I,lm+1)
415 molod 1.6 TMP5(I,lm+1) = PLKE(I,lm+1)*P0KINV
416 molod 1.1 ENDDO
417     DO 113 I=num,num+nindeces(nsubcl)-1
418     TMP4(I,NSUBCL+1) = PLE (I,lm+1)
419 molod 1.6 TMP5(I,NSUBCL+1) = PLKE(I,lm+1)*P0KINV
420 molod 1.1 113 CONTINUE
421    
422     do i=num,num+nindeces(nsubcl)-1
423     C Temperature at top of sub-cloud layer
424     tmp2(i,1) = TH(i,NSUBCL) * PLKE(i,NSUBCL)/P0KAPPA
425     C Pressure at top of sub-cloud layer
426     tmp2(i,2) = tmp4(i,nsubcl)
427     enddo
428    
429     C CHANGED THIS: no RH requirement for RAS
430     c call vqsat ( tmp2(num,1),tmp2(num,2),tmp2(num,3),
431     c . dum,.false.,nindeces(nsubcl) )
432     c do i=num,num+nindeces(nsubcl)-1
433     c rh = SHL(I,NSUBCL) / tmp2(i,3)
434     c if (rh .le. 0.85) then
435     c rhfrac(i) = 0.
436     c else if (rh .ge. 0.95) then
437     c rhfrac(i) = 1.
438     c else
439     c rhfrac(i) = (rh-0.85)*10.
440     c endif
441     c enddo
442     do i=num,num+nindeces(nsubcl)-1
443     rhfrac(i) = 1.
444     enddo
445    
446     C Compute RH threshold for Large-scale condensation
447     C Used in Slingo-Ritter clouds as well - define offset between SR and LS
448    
449     C Top level of atan func above this rh_threshold = rhmin
450     pup = 600.
451     do i=num,num+nindeces(nsubcl)-1
452 molod 1.2 do L = nsubcl, lm
453     rhcrit(i,L) = 1.
454     enddo
455     do L = 1, nsubcl-1
456 molod 1.6 pcheck = pl(i,L)
457 molod 1.2 if (pcheck .le. pup) then
458     rhcrit(i,L) = rhmin
459     else
460 molod 1.6 ppbl = pl(i,nsubcl)
461 molod 1.2 rhcrit(i,L) = rhmin + (1.-rhmin)/(19.) *
462     . ((atan( (2.*(pcheck-pup)/(ppbl-pup)-1.) *
463 molod 1.1 . tan(20.*pi/21.-0.5*pi) )
464     . + 0.5*pi) * 21./pi - 1.)
465 molod 1.2 endif
466     enddo
467 molod 1.1 enddo
468    
469     c Save Initial Values of Temperature and Specific Humidity
470     c --------------------------------------------------------
471     do L = 1,lm
472     do i = num,num+nindeces(nsubcl)-1
473     saveth(i,L) = th (i,L)
474     saveq (i,L) = shl(i,L)
475     PCPEN (i,L) = 0.
476     CLFRAC(i,L) = 0.
477     enddo
478     enddo
479    
480     CALL RAS ( NN,istrip,nindeces(nsubcl),NLRAS,NLTOP,lm,TMSTP
481     1, UL(num,1,1),ntracer-ptracer,TH(num,NLTOP),SHL(num,NLTOP)
482     2, TMP4(num,NLTOP), TMP5(num,NLTOP),rnd, ncrnd, PCPEN(num,NLTOP)
483     3, CLBOTH(num,NLTOP), CLFRAC(num,NLTOP)
484     4, cldmas(num,nltop), detrain(num,nltop)
485     8, cp,grav,rkappa,alhl,rhfrac(num),rasmax )
486    
487     c Compute Diagnostic CLDMAS in RAS Subcloud Layers
488     c ------------------------------------------------
489     do L=nsubcl,lm
490     do I=num,num+nindeces(nsubcl)-1
491 molod 1.6 dum = dp(i,L)/(ple(i,lm+1)-ple(i,nsubcl))
492 molod 1.1 cldmas(i,L) = cldmas(i,L-1) - dum*cldmas(i,nsubcl-1)
493     enddo
494     enddo
495    
496     c Update Theta and Moisture due to RAS
497     c ------------------------------------
498     DO L=1,nsubcl
499     DO I=num,num+nindeces(nsubcl)-1
500     CVTH(I,L) = (TH (I,L) - saveth(i,l))
501     CVQ (I,L) = (SHL(I,L) - saveq (i,l))
502     ENDDO
503     ENDDO
504     DO L=nsubcl+1,lm
505     DO I=num,num+nindeces(nsubcl)-1
506     CVTH(I,L) = cvth(i,nsubcl)
507     CVQ (I,L) = cvq (i,nsubcl)
508     ENDDO
509     ENDDO
510    
511     DO L=nsubcl+1,lm
512     DO I=num,num+nindeces(nsubcl)-1
513     TH (I,L) = saveth(i,l) + cvth(i,l)
514     SHL(I,L) = saveq (i,l) + cvq (i,l)
515     ENDDO
516     ENDDO
517     DO L=1,lm
518     DO I=num,num+nindeces(nsubcl)-1
519     CVTH(I,L) = CVTH(I,L) *P0KINV*SP(I)*tminv
520     CVQ (I,L) = CVQ (I,L) *SP(I)*tminv
521     ENDDO
522     ENDDO
523    
524     c Compute Tracer Tendency due to RAS
525     c ----------------------------------
526     do nt = 1,ntracer-ptracer
527     DO L=1,nsubcl-1
528     DO I=num,num+nindeces(nsubcl)-1
529     CVU(I,L,nt) = ( UL(I,L,nt)-saveu(i,l,nt) )*sp(i)*tminv
530     ENDDO
531     ENDDO
532     DO L=nsubcl,lm
533     DO I=num,num+nindeces(nsubcl)-1
534     if( usubcl(i,nt).ne.0.0 ) then
535     cvu(i,L,nt) = ( ul(i,nsubcl,nt)-usubcl(i,nt) ) *
536     . ( saveu(i,L,nt)/usubcl(i,nt) )*sp(i)*tminv
537     else
538     cvu(i,L,nt) = 0.0
539     endif
540     ENDDO
541     ENDDO
542     enddo
543    
544     c Compute Diagnostic PSUBCLD (Subcloud Layer Pressure)
545     c ----------------------------------------------------
546     do i=num,num+nindeces(nsubcl)-1
547     lras = .false.
548     do L=nltop,nsubcl
549     if( cvq(i,L).ne.0.0 ) lras = .true.
550     enddo
551     psubcld (i) = 0.0
552     psubcld_cnt(i) = 0.0
553     if( lras ) then
554     psubcld (i) = sp(i)+ptop-ple(i,nsubcl)
555     psubcld_cnt(i) = 1.0
556     endif
557     enddo
558    
559    
560     C End of subcloud layer depth loop (iloop)
561    
562     num = num+nindeces(nsubcl)
563    
564     ENDDO
565    
566     C **********************************************************************
567     C **** TENDENCY UPDATES ****
568     C **** (Keep 'Gathered' tendencies in 'gather' arrays now) ****
569     C **********************************************************************
570    
571     call paste( CVTH,deltgather,istrip,im*jm,lm,NN )
572     call paste( CVQ,delqgather,istrip,im*jm,lm,NN )
573     do nt = 1,ntracer-ptracer
574     call paste( CVU(1,1,nt),delugather(1,1,nt),istrip,im*jm,lm,NN )
575     enddo
576    
577     C **********************************************************************
578     C And now paste some arrays for filling diagnostics
579     C (use pkegather to hold detrainment and tmpgather for cloud mass flux)
580     C **********************************************************************
581    
582     if(icldmas .gt.0) call paste( cldmas,tmpgather,istrip,im*jm,lm,NN)
583     if(idtrain .gt.0) call paste(detrain,pkegather,istrip,im*jm,lm,NN)
584     if(ipsubcld.gt.0) then
585     call paste(psubcld ,psubcldg ,istrip,im*jm,1,NN)
586     call paste(psubcld_cnt,psubcldgc,istrip,im*jm,1,NN)
587     endif
588    
589     C *********************************************************************
590     C **** RE-EVAPORATION OF PENETRATING CONVECTIVE RAIN ****
591     C *********************************************************************
592    
593     CALL STRIP ( thgather,TH ,im*jm,ISTRIP,lm,NN)
594     CALL STRIP ( shgather,SHL,im*jm,ISTRIP,lm,NN)
595     DO L=1,lm
596     DO I=1,ISTRIP
597     TH(I,L) = TH(I,L) + CVTH(I,L)*tmstp/SP(I)
598     SHL(I,L) = SHL(I,L) + CVQ(I,L)*tmstp/SP(I)
599     TL(I,L) = TH(I,L)*PLK(I,L)
600     saveth(I,L) = th(I,L)
601     saveq (I,L) = SHL(I,L)
602     ENDDO
603     ENDDO
604    
605 molod 1.6 CALL RNEVP (NN,ISTRIP,lm,TL,SHL,PCPEN,PL,CLFRAC,SP,DP,PLKE,
606 molod 1.1 . PLK,TH,TMP1,TMP2,TMP3,ITMP1,ITMP2,PCNET,PRECIP,
607 molod 1.14 . CLSBTH,TMSTP,one,cp,grav,alhl,gamfac,cldlz,rhcrit,offset,alpha)
608 molod 1.1
609     C **********************************************************************
610     C **** TENDENCY UPDATES ****
611     C **********************************************************************
612    
613     DO L=1,lm
614    
615     DO I =1,ISTRIP
616     TMP1(I,L) = sp(i) * (SHL(I,L)-saveq(I,L)) * tminv
617     ENDDO
618     CALL PSTBMP(TMP1(1,L),delqgather(1,L),ISTRIP,im*jm,1,NN)
619    
620     DO I =1,ISTRIP
621     TMP1(I,L) = sp(i) * ((TL(I,L)/PLK(I,L))-saveth(i,l)) * tminv
622     ENDDO
623     CALL PSTBMP(TMP1(1,L),deltgather(1,L),ISTRIP,im*jm,1,NN)
624    
625     C Paste rain evap tendencies into arrays for diagnostic output
626     c ------------------------------------------------------------
627     if(idtls.gt.0)then
628     DO I =1,ISTRIP
629     TMP1(I,L) = ((TL(I,L)/PLK(I,L))-saveth(i,l))*plk(i,l)*sday*tminv
630     ENDDO
631     call paste(tmp1(1,L),deltrnev(1,L),istrip,im*jm,1,NN)
632     endif
633    
634     if(idqls.gt.0)then
635     DO I =1,ISTRIP
636     TMP1(I,L) = (SHL(I,L)-saveq(I,L)) * 1000. * sday * tminv
637     ENDDO
638     call paste(tmp1(1,L),delqrnev(1,L),istrip,im*jm,1,NN)
639     endif
640    
641     ENDDO
642    
643     C *********************************************************************
644     C Add Non-Precipitating Clouds where the relative
645     C humidity is less than 100%
646     C Apply Cloud Top Entrainment Instability
647     C *********************************************************************
648    
649     do L=1,lm
650     do i=1,istrip
651     srcld(i,L) = -clsbth(i,L)
652     enddo
653     enddo
654    
655     call srclouds (saveth,saveq,plk,pl,plke,clsbth,cldlz,istrip,lm,
656     . rhcrit,offset,alpha)
657    
658     do L=1,lm
659     do i=1,istrip
660     srcld(i,L) = srcld(i,L)+clsbth(i,L)
661     enddo
662     enddo
663    
664     C *********************************************************************
665     C **** PASTE CLOUD AMOUNTS ****
666     C *********************************************************************
667    
668     call paste ( srcld, cldsr,istrip,im*jm,lm,nn )
669     call paste ( cldlz,cldwater,istrip,im*jm,lm,nn )
670     call paste ( clsbth, cldls,istrip,im*jm,lm,nn )
671     call paste ( clboth, cpen ,istrip,im*jm,lm,nn )
672    
673     c compute Total Accumulated Precip for Landsurface Model
674     c ------------------------------------------------------
675     do i = 1,istrip
676     C Initialize Rainlsp, Rainconv and Snowfall
677     tmp1(i,1) = 0.0
678     tmp1(i,2) = 0.0
679     tmp1(i,3) = 0.0
680     enddo
681    
682     do i = 1,istrip
683     prep(i) = PRECIP(I) + PCNET(I)
684     tmp1(i,1) = PRECIP(I)
685     tmp1(i,2) = pcnet(i)
686     enddo
687     c
688     c check whether there is snow
689     c-------------------------------------------------------
690     c snow algorthm:
691     c if temperature profile from the surface level to 700 mb
692     c uniformaly c below zero, then precipitation (total) is
693     c snowfall. Else there is no snow.
694     c-------------------------------------------------------
695    
696     do i = 1,istrip
697     snowcrit=0
698     do l=lup,lm
699     if (saveth(i,l)*plk(i,l).le. tice ) then
700     snowcrit=snowcrit+1
701     endif
702     enddo
703     if (snowcrit .eq. (lm-lup+1)) then
704     tmp1(i,3) = prep(i)
705     tmp1(i,1)=0.0
706     tmp1(i,2)=0.0
707     endif
708     enddo
709    
710     CALL paste (tmp1(1,1), lsp_new,ISTRIP,im*jm,1,NN)
711     CALL paste (tmp1(1,2),conv_new,ISTRIP,im*jm,1,NN)
712     CALL paste (tmp1(1,3),snow_new,ISTRIP,im*jm,1,NN)
713    
714     if(iprecon.gt.0) then
715     CALL paste (pcnet,raincgath,ISTRIP,im*jm,1,NN)
716     endif
717    
718     C *********************************************************************
719     C **** End Major Stripped Region ****
720     C *********************************************************************
721    
722     1000 CONTINUE
723    
724     C Large Scale Rainfall, Conv rain, and snowfall
725     c ---------------------------------------------
726     call back2grd ( lsp_new,pblindex, lsp_new,im*jm)
727     call back2grd (conv_new,pblindex,conv_new,im*jm)
728     call back2grd (snow_new,pblindex,snow_new,im*jm)
729    
730     if(iprecon.gt.0) then
731     call back2grd (raincgath,pblindex,raincgath,im*jm)
732     endif
733    
734     c Subcloud Layer Pressure
735     c -----------------------
736     if(ipsubcld.gt.0) then
737     call back2grd (psubcldg ,pblindex,psubcldg ,im*jm)
738     call back2grd (psubcldgc,pblindex,psubcldgc,im*jm)
739     endif
740    
741     do L = 1,lm
742     C Delta theta,q, convective, max and ls clouds
743     c --------------------------------------------
744     call back2grd (deltgather(1,L),pblindex, dtmoist(1,1,L) ,im*jm)
745     call back2grd (delqgather(1,L),pblindex, dqmoist(1,1,L,1),im*jm)
746     call back2grd ( cpen(1,1,L),pblindex, cpen(1,1,L) ,im*jm)
747     call back2grd ( cldls(1,1,L),pblindex, cldls(1,1,L) ,im*jm)
748     call back2grd (cldwater(1,1,L),pblindex,cldwater(1,1,L) ,im*jm)
749     call back2grd ( pkzgather(1,L),pblindex, pkzgather(1,L) ,im*jm)
750    
751     C Diagnostics:
752     c ------------
753     if(icldmas.gt.0)call back2grd(tmpgather(1,L),pblindex,
754     . tmpgather(1,L),im*jm)
755     if(idtrain.gt.0)call back2grd(pkegather(1,L),pblindex,
756     . pkegather(1,L),im*jm)
757     if(idtls.gt.0)call back2grd(deltrnev(1,L),pblindex,
758     . deltrnev(1,L),im*jm)
759     if(idqls.gt.0)call back2grd(delqrnev(1,L),pblindex,
760     . delqrnev(1,L),im*jm)
761     if(icldnp.gt.0)call back2grd(cldsr(1,1,L),pblindex,
762     . cldsr(1,1,L),im*jm)
763     enddo
764    
765     c Tracers
766     c -------
767     do nt = 1,ntracer-ptracer
768     do L = 1,lm
769     call back2grd (delugather(1,L,nt),pblindex,
770     . dqmoist(1,1,L,ptracer+nt),im*jm)
771     enddo
772     enddo
773    
774    
775     C **********************************************************************
776     C BUMP DIAGNOSTICS
777     C **********************************************************************
778    
779 molod 1.16
780 molod 1.1 c Sub-Cloud Layer
781     c -------------------------
782     if( ipsubcld.ne.0 ) then
783     do j = 1,jm
784     do i = 1,im
785 molod 1.3 qdiag(i,j,ipsubcld,bi,bj) = qdiag(i,j,ipsubcld,bi,bj) +
786     . psubcldg (i,j)
787     qdiag(i,j,ipsubcldc,bi,bj) = qdiag(i,j,ipsubcldc,bi,bj) +
788     . psubcldgc(i,j)
789 molod 1.1 enddo
790     enddo
791     endif
792    
793     c Non-Precipitating Cloud Fraction
794     c --------------------------------
795     if( icldnp.ne.0 ) then
796     do L = 1,lm
797     do j = 1,jm
798     do i = 1,im
799 molod 1.3 qdiag(i,j,icldnp+L-1,bi,bj) = qdiag(i,j,icldnp+L-1,bi,bj) +
800     . cldsr(i,j,L)
801 molod 1.1 enddo
802     enddo
803     enddo
804     ncldnp = ncldnp + 1
805     endif
806    
807     c Moist Processes Heating Rate
808     c ----------------------------
809     if(imoistt.gt.0) then
810     do L = 1,lm
811 molod 1.16 do j = 1,jm
812     do i = 1,im
813     indgath = (j-1)*im + i
814     qdiag(i,j,imoistt+L-1,bi,bj) = qdiag(i,j,imoistt+L-1,bi,bj) +
815     . (dtmoist(i,j,L)*sday*pkzgather(indgath,L)/pz(i,j))
816     enddo
817 molod 1.1 enddo
818     enddo
819     endif
820    
821     c Moist Processes Moistening Rate
822     c -------------------------------
823     if(imoistq.gt.0) then
824     do L = 1,lm
825     do j = 1,jm
826     do i = 1,im
827 molod 1.3 qdiag(i,j,imoistq+L-1,bi,bj) = qdiag(i,j,imoistq+L-1,bi,bj) +
828 molod 1.1 . (dqmoist(i,j,L,1)*sday*1000.0/pz(i,j))
829     enddo
830     enddo
831     enddo
832     endif
833    
834     c Cloud Mass Flux
835     c ---------------
836     if(icldmas.gt.0) then
837     do L = 1,lm
838 molod 1.16 do j = 1,jm
839     do i = 1,im
840     indgath = (j-1)*im + i
841     qdiag(i,j,icldmas+L-1,bi,bj) = qdiag(i,j,icldmas+L-1,bi,bj) +
842     . tmpgather(indgath,L)
843     enddo
844 molod 1.1 enddo
845     enddo
846     endif
847    
848     c Detrained Cloud Mass Flux
849     c -------------------------
850     if(idtrain.gt.0) then
851     do L = 1,lm
852 molod 1.16 do j = 1,jm
853     do i = 1,im
854     indgath = (j-1)*im + i
855     qdiag(i,j,idtrain+L-1,bi,bj) = qdiag(i,j,idtrain+L-1,bi,bj) +
856     . pkegather(indgath,L)
857     enddo
858 molod 1.1 enddo
859     enddo
860     endif
861    
862     c Grid-Scale Condensational Heating Rate
863     c --------------------------------------
864     if(idtls.gt.0) then
865     do L = 1,lm
866 molod 1.16 do j = 1,jm
867     do i = 1,im
868     indgath = (j-1)*im + i
869     qdiag(i,j,idtls+L-1,bi,bj) = qdiag(i,j,idtls+L-1,bi,bj) +
870     . deltrnev(indgath,L)
871     enddo
872 molod 1.1 enddo
873     enddo
874     endif
875    
876     c Grid-Scale Condensational Moistening Rate
877     c -----------------------------------------
878     if(idqls.gt.0) then
879     do L = 1,lm
880 molod 1.16 do j = 1,jm
881     do i = 1,im
882     indgath = (j-1)*im + i
883     qdiag(i,j,idqls+L-1,bi,bj) = qdiag(i,j,idqls+L-1,bi,bj) +
884     . delqrnev(indgath,L)
885     enddo
886 molod 1.1 enddo
887     enddo
888     endif
889    
890     c Total Precipitation
891     c -------------------
892     if(ipreacc.gt.0) then
893     do j = 1,jm
894     do i = 1,im
895 molod 1.3 qdiag(i,j,ipreacc,bi,bj) = qdiag(i,j,ipreacc,bi,bj)
896 molod 1.1 . + ( lsp_new(I,j)
897     . + snow_new(I,j)
898     . + conv_new(i,j) ) *sday*tminv
899     enddo
900     enddo
901     endif
902    
903     c Convective Precipitation
904     c ------------------------
905     if(iprecon.gt.0) then
906 molod 1.16 do j = 1,jm
907     do i = 1,im
908     indgath = (j-1)*im + i
909     qdiag(i,j,iprecon,bi,bj) = qdiag(i,j,iprecon,bi,bj) +
910     . raincgath(indgath)*sday*tminv
911     enddo
912 molod 1.1 enddo
913     endif
914    
915     C **********************************************************************
916     C **** Fill Rainfall and Snowfall Arrays for Land Surface Model ****
917     C **** Note: Precip Rates work when DT(turb)<DT(moist) ****
918     C **********************************************************************
919    
920     do j = 1,jm
921     do i = 1,im
922     rainlsp (i,j) = rainlsp (i,j) + lsp_new(i,j)*tminv
923     rainconv(i,j) = rainconv(i,j) + conv_new(i,j)*tminv
924     snowfall(i,j) = snowfall(i,j) + snow_new(i,j)*tminv
925     enddo
926     enddo
927    
928     C **********************************************************************
929     C *** Compute Time-averaged Quantities for Radiation ***
930     C *** CPEN => Cloud Fraction from RAS ***
931     C *** CLDLS => Cloud Fraction from RNEVP ***
932     C **********************************************************************
933    
934     do j = 1,jm
935     do i = 1,im
936     cldhi (i,j) = 0.
937     cldmid(i,j) = 0.
938     cldlow(i,j) = 0.
939     cldtmp(i,j) = 0.
940     cldprs(i,j) = 0.
941     tmpimjm(i,j) = 0.
942     enddo
943     enddo
944    
945     c Set Moist-Process Memory Coefficient
946     c ------------------------------------
947     cldras_mem = 1.0-tmstp/ 3600.0
948     cldlsp_mem = 1.0-tmstp/(3600.0*3)
949    
950     do L = 1,lm
951     do i = 1,im*jm
952 molod 1.6 plev = pl(i,L)
953 molod 1.1
954     c Compute Time-averaged Cloud and Water Amounts for Longwave Radiation
955     c --------------------------------------------------------------------
956     watnow = cldwater(i,1,L)
957     if( plev.le.500.0 ) then
958     cldras = min( max( cldras_lw(i,1,L)*cldras_mem,cpen(i,1,L)),1.0)
959     else
960     cldras = 0.0
961     endif
962     cldlsp = min( max( cldlsp_lw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0)
963    
964     if( cldras.lt.cldmin ) cldras = 0.0
965     if( cldlsp.lt.cldmin ) cldlsp = 0.0
966    
967     cldnow = max( cldlsp,cldras )
968    
969     lwlz(i,1,L) = ( nlwlz*lwlz(i,1,L) + watnow)/(nlwlz +1)
970     cldtot_lw(i,1,L) = (nlwcld*cldtot_lw(i,1,L) + cldnow)/(nlwcld+1)
971     cldlsp_lw(i,1,L) = (nlwcld*cldlsp_lw(i,1,L) + cldlsp)/(nlwcld+1)
972     cldras_lw(i,1,L) = (nlwcld*cldras_lw(i,1,L) + cldras)/(nlwcld+1)
973    
974    
975     c Compute Time-averaged Cloud and Water Amounts for Shortwave Radiation
976     c ---------------------------------------------------------------------
977     watnow = cldwater(i,1,L)
978     if( plev.le.500.0 ) then
979     cldras = min( max(cldras_sw(i,1,L)*cldras_mem, cpen(i,1,L)),1.0)
980     else
981     cldras = 0.0
982     endif
983     cldlsp = min( max(cldlsp_sw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0)
984    
985     if( cldras.lt.cldmin ) cldras = 0.0
986     if( cldlsp.lt.cldmin ) cldlsp = 0.0
987    
988     cldnow = max( cldlsp,cldras )
989    
990     swlz(i,1,L) = ( nswlz*swlz(i,1,L) + watnow)/(nswlz +1)
991     cldtot_sw(i,1,L) = (nswcld*cldtot_sw(i,1,L) + cldnow)/(nswcld+1)
992     cldlsp_sw(i,1,L) = (nswcld*cldlsp_sw(i,1,L) + cldlsp)/(nswcld+1)
993     cldras_sw(i,1,L) = (nswcld*cldras_sw(i,1,L) + cldras)/(nswcld+1)
994    
995    
996     c Compute Instantaneous Low-Mid-High Maximum Overlap Cloud Fractions
997     c ----------------------------------------------------------------------
998    
999     if( L.lt.midlevel ) cldhi (i,1) = max( cldnow,cldhi (i,1) )
1000     if( L.ge.midlevel .and.
1001     . L.lt.lowlevel ) cldmid(i,1) = max( cldnow,cldmid(i,1) )
1002     if( L.ge.lowlevel ) cldlow(i,1) = max( cldnow,cldlow(i,1) )
1003    
1004     c Compute Cloud-Top Temperature and Pressure
1005     c ------------------------------------------
1006     cldtmp(i,1) = cldtmp(i,1) + cldnow*pkzgather(i,L)
1007     . * ( tz(i,1,L) + dtmoist(i,1,L)*tmstp/pz(i,1) )
1008     cldprs(i,1) = cldprs(i,1) + cldnow*plev
1009     tmpimjm(i,1) = tmpimjm(i,1) + cldnow
1010    
1011     enddo
1012     enddo
1013    
1014     c Compute Instantanious Total 2-D Cloud Fraction
1015     c ----------------------------------------------
1016     do j = 1,jm
1017     do i = 1,im
1018     totcld(i,j) = 1.0 - (1.-cldhi (i,j))
1019     . * (1.-cldmid(i,j))
1020     . * (1.-cldlow(i,j))
1021     enddo
1022     enddo
1023    
1024    
1025     C **********************************************************************
1026     C *** Fill Cloud Top Pressure and Temperature Diagnostic ***
1027     C **********************************************************************
1028    
1029     if(icldtmp.gt.0) then
1030     do j = 1,jm
1031     do i = 1,im
1032     if( cldtmp(i,j).gt.0.0 ) then
1033 molod 1.3 qdiag(i,j,icldtmp,bi,bj) = qdiag(i,j,icldtmp,bi,bj) +
1034 molod 1.1 . cldtmp(i,j)*totcld(i,j)/tmpimjm(i,j)
1035 molod 1.3 qdiag(i,j,icttcnt,bi,bj) = qdiag(i,j,icttcnt,bi,bj) +
1036     . totcld(i,j)
1037 molod 1.1 endif
1038     enddo
1039     enddo
1040     endif
1041    
1042     if(icldprs.gt.0) then
1043     do j = 1,jm
1044     do i = 1,im
1045     if( cldprs(i,j).gt.0.0 ) then
1046 molod 1.3 qdiag(i,j,icldprs,bi,bj) = qdiag(i,j,icldprs,bi,bj) +
1047 molod 1.1 . cldprs(i,j)*totcld(i,j)/tmpimjm(i,j)
1048 molod 1.3 qdiag(i,j,ictpcnt,bi,bj) = qdiag(i,j,ictpcnt,bi,bj) +
1049     . totcld(i,j)
1050 molod 1.1 endif
1051     enddo
1052     enddo
1053     endif
1054    
1055     C **********************************************************************
1056     C **** INCREMENT COUNTERS ****
1057     C **********************************************************************
1058    
1059     nlwlz = nlwlz + 1
1060     nswlz = nswlz + 1
1061    
1062     nlwcld = nlwcld + 1
1063     nswcld = nswcld + 1
1064    
1065 molod 1.16 #ifdef ALLOW_DIAGNOSTICS
1066     if( (bi.eq.1) .and. (bj.eq.1) ) then
1067 molod 1.1 nmoistt = nmoistt + 1
1068     nmoistq = nmoistq + 1
1069     npreacc = npreacc + 1
1070     nprecon = nprecon + 1
1071    
1072     ncldmas = ncldmas + 1
1073     ndtrain = ndtrain + 1
1074    
1075     ndtls = ndtls + 1
1076     ndqls = ndqls + 1
1077 molod 1.16 endif
1078     #endif
1079 molod 1.1
1080     RETURN
1081     END
1082 molod 1.16 SUBROUTINE RAS( NN, LNG, LENC, K, NLTOP, nlayr, DT
1083 molod 1.1 *, UOI, ntracer, POI, QOI, PRS, PRJ, rnd, ncrnd
1084     *, RAINS, CLN, CLF, cldmas, detrain
1085     *, cp,grav,rkappa,alhl,rhfrac,rasmax )
1086     C
1087     C*********************************************************************
1088     C********************* SUBROUTINE RAS *****************************
1089     C********************** 16 MARCH 1988 ******************************
1090     C*********************************************************************
1091     C
1092 molod 1.8 implicit none
1093    
1094 molod 1.9 C Argument List
1095 molod 1.16 integer nn,lng,lenc,k,nltop,nlayr
1096 molod 1.1 integer ntracer
1097 molod 1.9 integer ncrnd
1098 molod 1.13 _RL dt
1099 molod 1.16 _RL UOI(lng,nlayr,ntracer), POI(lng,K)
1100     _RL QOI(lng,K), PRS(lng,K+1), PRJ(lng,K+1)
1101 molod 1.13 _RL rnd(ncrnd)
1102 molod 1.16 _RL RAINS(lng,K), CLN(lng,K), CLF(lng,K)
1103     _RL cldmas(lng,K), detrain(lng,K)
1104     _RL cp,grav,rkappa,alhl,rhfrac(lng),rasmax
1105 molod 1.9
1106     C Local Variables
1107 molod 1.16 _RL TCU(lng,K), QCU(lng,K)
1108     _RL ucu(lng,K,ntracer)
1109     _RL ALF(lng,K), BET(lng,K), GAM(lng,K)
1110     *, ETA(lng,K), HOI(lng,K)
1111     *, PRH(lng,K), PRI(lng,K)
1112     _RL HST(lng,K), QOL(lng,K), GMH(lng,K)
1113    
1114     _RL TX1(lng), TX2(lng), TX3(lng), TX4(lng), TX5(lng)
1115     *, TX6(lng), TX7(lng), TX8(lng), TX9(lng)
1116     *, TX11(lng), TX12(lng), TX13(lng), TX14(lng,ntracer)
1117     *, TX15(lng)
1118     *, WFN(lng)
1119     integer IA1(lng), IA2(lng), IA3(lng)
1120     _RL cloudn(lng), pcu(lng)
1121 molod 1.1
1122 molod 1.8 integer krmin,icm
1123 molod 1.13 _RL rknob, cmb2pa
1124 molod 1.8 PARAMETER (KRMIN=01)
1125     PARAMETER (ICM=1000)
1126     PARAMETER (CMB2PA=100.0)
1127     PARAMETER (rknob = 10.)
1128 molod 1.9
1129     integer IC(ICM), IRND(icm)
1130 molod 1.16 _RL cmass(lng,K)
1131 molod 1.9 LOGICAL SETRAS
1132    
1133     integer i,L,nc,ib,nt
1134 molod 1.8 integer km1,kp1,kprv,kcr,kfx,ncmx
1135 molod 1.13 _RL p00, crtmsf, frac, rasblf
1136 molod 1.8
1137     do L = 1,k
1138     do I = 1,LENC
1139     rains(i,l) = 0.
1140     enddo
1141     enddo
1142 molod 1.1
1143     p00 = 1000.
1144     crtmsf = 0.
1145    
1146     C The numerator here is the fraction of the subcloud layer mass flux
1147     C allowed to entrain into the cloud
1148    
1149     CCC FRAC = 1./dt
1150     FRAC = 0.5/dt
1151    
1152     KM1 = K - 1
1153     KP1 = K + 1
1154     C we want the ras adjustment time scale to be one hour (indep of dt)
1155     RASBLF = 1./3600.
1156     C
1157     KPRV = KM1
1158     C Removed KRMAX parameter
1159     KCR = MIN(KM1,nlayr-2)
1160     KFX = KM1 - KCR
1161     NCMX = KFX + NCRND
1162     C
1163     IF (KFX .GT. 0) THEN
1164     DO NC=1,KFX
1165     IC(NC) = K - NC
1166     ENDDO
1167     ENDIF
1168     C
1169     IF (NCRND .GT. 0) THEN
1170     DO I=1,ncrnd
1171     IRND(I) = (RND(I)-0.0005)*(KCR-KRMIN+1)
1172     IRND(I) = IRND(I) + KRMIN
1173     ENDDO
1174     C
1175     DO NC=1,NCRND
1176     IC(KFX+NC) = IRND(NC)
1177     ENDDO
1178     ENDIF
1179     C
1180     DO 100 NC=1,NCMX
1181     C
1182     IF (NC .EQ. 1 ) THEN
1183     SETRAS = .TRUE.
1184     ELSE
1185     SETRAS = .FALSE.
1186     ENDIF
1187     IB = IC(NC)
1188    
1189 molod 1.16 print *,' Calling cloud for cloud ',nc,' det at ',ic(nc)
1190    
1191 molod 1.1 c Initialize Cloud Fraction Array
1192     c -------------------------------
1193     do i = 1,lenc
1194     cloudn(i) = 0.0
1195     enddo
1196    
1197 molod 1.16 CALL CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IB, RASBLF,SETRAS,FRAC
1198 molod 1.1 *, CP, ALHL, RKAPPA, GRAV, P00, CRTMSF
1199     *, POI, QOI, UOI, Ntracer, PRS, PRJ
1200     *, PCU, CLOUDN, TCU, QCU, UCU, CMASS
1201     *, ALF, BET, GAM, PRH, PRI, HOI, ETA
1202     *, HST, QOL, GMH
1203     *, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9
1204     *, WFN, TX11, TX12, TX13, TX14, TX15
1205     *, IA1,IA2,IA3,rhfrac)
1206    
1207     C Compute fraction of grid box into which rain re-evap occurs (clf)
1208     c -----------------------------------------------------------------
1209     do i = 1,lenc
1210    
1211     c mass in detrainment layer
1212     c -------------------------
1213     tx1(i) = cmb2pa * (prs(i,ib+1) - prs(i,ib))/(grav*dt)
1214    
1215     c ratio of detraining cloud mass to mass in detrainment layer
1216     c -----------------------------------------------------------
1217     tx1(i) = rhfrac(i)*rknob * cmass(i,ib) / tx1(i)
1218     if(cmass(i,K).gt.0.) clf(i,ib) = clf(i,ib) + tx1(i)
1219     if( clf(i,ib).gt.1.) clf(i,ib) = 1.
1220     enddo
1221    
1222     c Compute Total Cloud Mass Flux
1223     c *****************************
1224     do L=ib,k
1225     do i=1,lenc
1226     cmass(i,L) = rhfrac(i)*cmass(i,L) * dt
1227     enddo
1228     enddo
1229    
1230     do L=ib,k
1231     do i=1,lenc
1232     cldmas(i,L) = cldmas(i,L) + cmass(i,L)
1233     enddo
1234     enddo
1235    
1236     do i=1,lenc
1237     detrain(i,ib) = detrain(i,ib) + cmass(i,ib)
1238     enddo
1239    
1240     DO L=IB,K
1241     DO I=1,LENC
1242     POI(I,L) = POI(I,L) + TCU(I,L) * DT * rhfrac(i)
1243     QOI(I,L) = QOI(I,L) + QCU(I,L) * DT * rhfrac(i)
1244     ENDDO
1245     ENDDO
1246     DO NT=1,Ntracer
1247     DO L=IB,K
1248     DO I=1,LENC
1249     UOI(I,L+nltop-1,NT)=UOI(I,L+nltop-1,NT)+UCU(I,L,NT)*DT*rhfrac(i)
1250     ENDDO
1251     ENDDO
1252     ENDDO
1253     DO I=1,LENC
1254     rains(I,ib) = rains(I,ib) + PCU(I)*dt * rhfrac(i)
1255     ENDDO
1256    
1257     100 CONTINUE
1258    
1259     c Fill Convective Cloud Fractions based on 3-D Rain Amounts
1260     c ---------------------------------------------------------
1261     do L=k-1,1,-1
1262     do i=1,lenc
1263     tx1(i) = 100*(prs(i,L+1)-prs(i,L))/grav
1264     cln(i,L) = min(1600*rains(i,L)/tx1(i),rasmax )
1265     enddo
1266     enddo
1267    
1268     RETURN
1269     END
1270     subroutine rndcloud (iras,nrnd,rnd,myid)
1271     implicit none
1272     integer n,iras,nrnd,myid
1273 molod 1.13 _RL random_numbx
1274     _RL rnd(nrnd)
1275 molod 1.1 integer irm
1276     parameter (irm = 1000)
1277 molod 1.13 _RL random(irm)
1278 molod 1.16 integer i,mcheck,numrand,iseed,indx
1279 molod 1.1 logical first
1280     data first /.true./
1281     integer iras0
1282     data iras0 /0/
1283     save random, iras0
1284    
1285     if(nrnd.eq.0.)then
1286     do i = 1,nrnd
1287     rnd(i) = 0
1288     enddo
1289 molod 1.15 if(first .and. myid.eq.1) print *,' NO RANDOM CLOUDS IN RAS '
1290 molod 1.1 go to 100
1291     endif
1292    
1293     mcheck = mod(iras-1,irm/nrnd)
1294    
1295     c First Time In From a Continuing RESTART (IRAS.GT.1) or Reading a New RESTART
1296     c ----------------------------------------------------------------------------
1297     if( first.and.(iras.gt.1) .or. iras.ne.iras0+1 )then
1298 molod 1.15 if( myid.eq.1 ) print *, 'Recreating Rand Numb Array in RNDCLOUD'
1299     if( myid.eq.1 ) print *, 'IRAS: ',iras,' IRAS0: ',iras0
1300 molod 1.1 numrand = mod(iras,irm/nrnd) * nrnd
1301     iseed = iras * nrnd - numrand
1302     call random_seedx(iseed)
1303     do i = 1,irm
1304 molod 1.12 random(i) = random_numbx(iseed)
1305 molod 1.1 enddo
1306 molod 1.16 indx = (iras-1)*nrnd
1307 molod 1.1
1308     c Multiple Time In But have Used Up all 1000 numbers (MCHECK.EQ.0)
1309     c ----------------------------------------------------------------
1310     else if (mcheck.eq.0) then
1311     iseed = (iras-1)*nrnd
1312     call random_seedx(iseed)
1313     do i = 1,irm
1314 molod 1.12 random(i) = random_numbx(iseed)
1315 molod 1.1 enddo
1316 molod 1.16 indx = iseed
1317 molod 1.1
1318     c Multiple Time In But have NOT Used Up all 1000 numbers (MCHECK.NE.0)
1319     c --------------------------------------------------------------------
1320     else
1321 molod 1.16 indx = (iras-1)*nrnd
1322 molod 1.1 endif
1323    
1324 molod 1.16 indx = mod(indx,irm)
1325     if( indx+nrnd.gt.1000 ) indx=1000-nrnd
1326 molod 1.1
1327     do n = 1,nrnd
1328 molod 1.16 rnd(n) = random(indx+n)
1329 molod 1.1 enddo
1330    
1331     100 continue
1332     first = .false.
1333     iras0 = iras
1334     return
1335     end
1336 molod 1.12 function random_numbx(iseed)
1337 molod 1.1 implicit none
1338 molod 1.12 integer iseed
1339     real *8 seed,port_rand
1340 molod 1.13 _RL random_numbx
1341 molod 1.11 random_numbx = 0
1342 molod 1.10 #ifdef CRAY
1343 molod 1.13 _RL ranf
1344 molod 1.1 random_numbx = ranf()
1345 molod 1.12 #else
1346 molod 1.10 #ifdef SGI
1347 molod 1.13 _RL rand
1348 molod 1.1 random_numbx = rand()
1349     #endif
1350 molod 1.12 random_numbx = port_rand(seed)
1351     #endif
1352 molod 1.1 return
1353     end
1354     subroutine random_seedx (iseed)
1355     implicit none
1356     integer iseed
1357 molod 1.10 #ifdef CRAY
1358 molod 1.1 call ranset (iseed)
1359     #endif
1360 molod 1.10 #ifdef SGI
1361 molod 1.1 integer*4 seed
1362     seed = iseed
1363     call srand (seed)
1364     #endif
1365     return
1366     end
1367 molod 1.16 SUBROUTINE CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IC, RASALF
1368 molod 1.1 *, SETRAS, FRAC
1369     *, CP, ALHL, RKAP, GRAV, P00, CRTMSF
1370     *, POI, QOI, UOI, Ntracer, PRS, PRJ
1371     *, PCU, CLN, TCU, QCU, UCU, CMASS
1372     *, ALF, BET, GAM, PRH, PRI, HOL, ETA
1373     *, HST, QOL, GMH
1374     *, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, ALM
1375     *, WFN, AKM, QS1, CLF, UHT, WLQ
1376     *, IA, I1, I2,rhfrac)
1377     C
1378     C*********************************************************************
1379     C******************** Relaxed Arakawa-Schubert ***********************
1380     C********************* Plug Compatible Version **********************
1381     C************************ SUBROUTINE CLOUD ***************************
1382     C************************* 23 JULY 1992 ***************************
1383     C*********************************************************************
1384     C*********************************************************************
1385     C*********************************************************************
1386     C************************** Developed By *****************************
1387     C************************** *****************************
1388     C************************ Shrinivas Moorthi **************************
1389     C************************ and **************************
1390     C************************ Max J. Suarez *****************************
1391     C************************ *****************************
1392     C******************** Laboratory for Atmospheres *********************
1393     C****************** NASA/GSFC, Greenbelt, MD 20771 *******************
1394     C*********************************************************************
1395     C*********************************************************************
1396     C
1397     C The calculations of Moorthi and Suarez (1992, MWR) are
1398     C contained in the CLOUD routine.
1399     C It is probably advisable, at least initially, to treat CLOUD
1400     C as a black box that computes the single cloud adjustments. RAS,
1401     C on the other hand, can be tailored to each GCMs configuration
1402     C (ie, number and placement of levels, nature of boundary layer,
1403     C time step and frequency with which RAS is called).
1404     C
1405     C
1406     C Input:
1407     C ------
1408     C
1409 molod 1.16 C lng : The inner dimension of update and input arrays.
1410 molod 1.1 C
1411     C LENC : The run: the number of soundings processes in a single call.
1412 molod 1.16 C RAS works on the first LENC of the lng soundings
1413 molod 1.1 C passed. This allows working on pieces of the world
1414     C say for multitasking, without declaring temporary arrays
1415     C and copying the data to and from them. This is an f77
1416     C version. An F90 version would have to allow more
1417     C flexibility in the argument declarations. Obviously
1418 molod 1.16 C (LENC<=lng).
1419 molod 1.1 C
1420     C K : Number of vertical layers (increasing downwards).
1421     C Need not be the same as the number of layers in the
1422     C GCM, since it is the outer dimension. The bottom layer
1423     C (K) is the subcloud layer.
1424     C
1425     C IC : Detrainment level to check for presence of convection
1426     C
1427     C RASALF : Relaxation parameter (< 1.) for present cloud-type
1428     C
1429     C SETRAS : Logical parameter to control re-calculation of
1430     C saturation specific humidity and mid level P**kappa
1431     C
1432     C FRAC : Fraction of the PBL (layer K) mass allowed to be used
1433     C by a cloud-type in time DT
1434     C
1435     C CP : Specific heat at constant pressure
1436     C
1437     C ALHL : Latent Heat of condensation
1438     C
1439     C RKAP : R/Cp, where R is the gas constant
1440     C
1441     C GRAV : Acceleration due to gravity
1442     C
1443     C P00 : A reference pressure in hPa, useually 1000 hPa
1444     C
1445     C CRTMSF : Critical value of mass flux above which cloudiness at
1446     C the detrainment layer of that cloud-type is assumed.
1447     C Affects only cloudiness calculation.
1448     C
1449 molod 1.16 C POI : 2D array of dimension (lng,K) containing potential
1450 molod 1.1 C temperature. Updated but not initialized by RAS.
1451     C
1452 molod 1.16 C QOI : 2D array of dimension (lng,K) containing specific
1453 molod 1.1 C humidity. Updated but not initialized by RAS.
1454     C
1455 molod 1.16 C UOI : 3D array of dimension (lng,K,NTRACER) containing tracers
1456 molod 1.1 C Updated but not initialized by RAS.
1457     C
1458 molod 1.16 C PRS : 2D array of dimension (lng,K+1) containing pressure
1459 molod 1.1 C in hPa at the interfaces of K-layers from top of the
1460     C atmosphere to the bottom. Not modified.
1461     C
1462 molod 1.16 C PRJ : 2D array of dimension (lng,K+1) containing (PRS/P00) **
1463 molod 1.1 C RKAP. i.e. Exner function at layer edges. Not modified.
1464     C
1465 molod 1.16 C rhfrac : 1D array of dimension (lng) containing a rel.hum. scaling
1466 molod 1.1 C fraction. Not modified.
1467     C
1468     C Output:
1469     C -------
1470     C
1471 molod 1.16 C PCU : 1D array of length lng containing accumulated
1472 molod 1.1 C precipitation in mm/sec.
1473     C
1474 molod 1.16 C CLN : 2D array of dimension (lng,K) containing cloudiness
1475 molod 1.1 C Note: CLN is bumped but NOT initialized
1476     C
1477 molod 1.16 C TCU : 2D array of dimension (lng,K) containing accumulated
1478 molod 1.1 C convective heating (K/sec).
1479     C
1480 molod 1.16 C QCU : 2D array of dimension (lng,K) containing accumulated
1481 molod 1.1 C convective drying (kg/kg/sec).
1482     C
1483 molod 1.16 C CMASS : 2D array of dimension (lng,K) containing the
1484 molod 1.1 C cloud mass flux (kg/sec). Filled from cloud top
1485     C to base.
1486     C
1487     C Temporaries:
1488     C
1489     C ALF, BET, GAM, ETA, PRH, PRI, HOI, HST, QOL, GMH are temporary
1490     C 2D real arrays of dimension of at least (LENC,K) where LENC is
1491     C the horizontal dimension over which convection is invoked.
1492     C
1493     C
1494     C TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9, AKM, QS1, CLF, UHT
1495     C VHT, WLQ WFN are temporary real arrays of length at least LENC
1496     C
1497     C IA, I1, and I2 are temporary integer arrays of length LENC
1498     C
1499     C
1500     C************************************************************************
1501 molod 1.9 implicit none
1502     C Argument List declarations
1503 molod 1.16 integer nn,lng,LENC,K,NLTOP,nlayr,ic,ntracer
1504 molod 1.13 _RL rasalf
1505 molod 1.9 LOGICAL SETRAS
1506 molod 1.13 _RL frac, cp, alhl, rkap, grav, p00, crtmsf
1507 molod 1.16 _RL POI(lng,K),QOI(lng,K),PRS(lng,K+1),PRJ(lng,K+1)
1508     _RL uoi(lng,nlayr,ntracer)
1509     _RL PCU(LENC), CLN(lng)
1510     _RL TCU(lng,K), QCU(lng,K), ucu(lng,k,ntracer), CMASS(lng,K)
1511     _RL ALF(lng,K), BET(lng,K), GAM(lng,K), PRH(lng,K), PRI(lng,K)
1512 molod 1.13 _RL HOL(LENC,K), ETA(LENC,K), HST(LENC,K), QOL(LENC,K)
1513     _RL GMH(LENC,K)
1514     _RL TX1(LENC), TX2(LENC), TX3(LENC), TX4(LENC)
1515     _RL TX5(LENC), TX6(LENC), TX7(LENC), TX8(LENC)
1516     _RL ALM(LENC), WFN(LENC), AKM(LENC), QS1(LENC)
1517     _RL WLQ(LENC), CLF(LENC)
1518 molod 1.16 _RL uht(lng,ntracer)
1519 molod 1.9 integer IA(LENC), I1(LENC),I2(LENC)
1520 molod 1.16 _RL rhfrac(lng)
1521 molod 1.1
1522 molod 1.9 C Local Variables
1523 molod 1.13 _RL daylen,half,one,zero,cmb2pa,rhmax
1524 molod 1.1 PARAMETER (DAYLEN=86400.0, HALF=0.5, ONE=1.0, ZERO=0.0)
1525     PARAMETER (CMB2PA=100.0)
1526     PARAMETER (RHMAX=0.9999)
1527 molod 1.13 _RL rkapp1,onebcp,albcp,onebg,cpbg,twobal
1528 molod 1.1 C
1529 molod 1.9 integer nt,km1,ic1,i,L,len1,len2,isav,len11,ii
1530 molod 1.16 integer lena,lena1,lenb
1531     _RL tem,tem1
1532 molod 1.1
1533     c Explicit Inline Directives
1534     c --------------------------
1535 molod 1.10 #ifdef CRAY
1536     #ifdef f77
1537 molod 1.1 cfpp$ expand (qsat)
1538     #endif
1539     #endif
1540    
1541     RKAPP1 = 1.0 + RKAP
1542     ONEBCP = 1.0 / CP
1543     ALBCP = ALHL * ONEBCP
1544     ONEBG = 1.0 / GRAV
1545     CPBG = CP * ONEBG
1546     TWOBAL = 2.0 / ALHL
1547     C
1548     KM1 = K - 1
1549     IC1 = IC + 1
1550     C
1551 molod 1.9 C SETTING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T.
1552 molod 1.1 C
1553    
1554     IF (SETRAS) THEN
1555    
1556     DO 2050 L=1,K
1557     DO 2030 I=1,LENC
1558     PRH(I,L) = (PRJ(I,L+1)*PRS(I,L+1) - PRJ(I,L)*PRS(I,L))
1559     * / ((PRS(I,L+1)-PRS(I,L)) * RKAPP1)
1560     2030 CONTINUE
1561     2050 CONTINUE
1562    
1563     DO 2070 L=1,K
1564     DO 2060 I=1,LENC
1565     TX5(I) = POI(I,L) * PRH(I,L)
1566     TX1(I) = (PRS(I,L) + PRS(I,L+1)) * 0.5
1567     TX3(I) = TX5(I)
1568     CALL QSAT(TX3(I), TX1(I), TX2(I), TX4(I), .TRUE.)
1569     ALF(I,L) = TX2(I) - TX4(I) * TX5(I)
1570     BET(I,L) = TX4(I) * PRH(I,L)
1571     GAM(I,L) = 1.0 / ((1.0 + TX4(I)*ALBCP) * PRH(I,L))
1572     PRI(I,L) = (CP/CMB2PA) / (PRS(I,L+1) - PRS(I,L))
1573     2060 CONTINUE
1574     2070 CONTINUE
1575    
1576     ENDIF
1577     C
1578     C
1579     DO 10 L=1,K
1580 molod 1.16 DO 10 I=1,lng
1581 molod 1.1 TCU(I,L) = 0.0
1582     QCU(I,L) = 0.0
1583     CMASS(I,L) = 0.0
1584     10 CONTINUE
1585    
1586     do nt = 1,ntracer
1587     do L=1,K
1588     do I=1,LENC
1589     ucu(I,L,nt) = 0.0
1590     enddo
1591     enddo
1592     enddo
1593     C
1594     DO 30 I=1,LENC
1595     TX1(I) = PRJ(I,K+1) * POI(I,K)
1596     QS1(I) = ALF(I,K) + BET(I,K)*POI(I,K)
1597     QOL(I,K) = MIN(QS1(I)*RHMAX,QOI(I,K))
1598    
1599     HOL(I,K) = TX1(I)*CP + QOL(I,K)*ALHL
1600     ETA(I,K) = ZERO
1601     TX2(I) = (PRJ(I,K+1) - PRJ(I,K)) * POI(I,K) * CP
1602     30 CONTINUE
1603     C
1604     IF (IC .LT. KM1) THEN
1605     DO 3703 L=KM1,IC1,-1
1606     DO 50 I=1,LENC
1607     QS1(I) = ALF(I,L) + BET(I,L)*POI(I,L)
1608     QOL(I,L) = MIN(QS1(I)*RHMAX,QOI(I,L))
1609     C
1610     TEM1 = TX2(I) + PRJ(I,L+1) * POI(I,L) * CP
1611     HOL(I,L) = TEM1 + QOL(I,L )* ALHL
1612     HST(I,L) = TEM1 + QS1(I) * ALHL
1613    
1614     TX1(I) = (PRJ(I,L+1) - PRJ(I,L)) * POI(I,L)
1615     ETA(I,L) = ETA(I,L+1) + TX1(I)*CPBG
1616     TX2(I) = TX2(I) + TX1(I)*CP
1617     50 CONTINUE
1618     C
1619     3703 CONTINUE
1620     ENDIF
1621    
1622    
1623     DO 70 I=1,LENC
1624     HOL(I,IC) = TX2(I)
1625     QS1(I) = ALF(I,IC) + BET(I,IC)*POI(I,IC)
1626     QOL(I,IC) = MIN(QS1(I)*RHMAX,QOI(I,IC))
1627     c
1628     TEM1 = TX2(I) + PRJ(I,IC1) * POI(I,IC) * CP
1629     HOL(I,IC) = TEM1 + QOL(I,IC) * ALHL
1630     HST(I,IC) = TEM1 + QS1(I) * ALHL
1631     C
1632     TX3(I ) = (PRJ(I,IC1) - PRH(I,IC)) * POI(I,IC)
1633     ETA(I,IC) = ETA(I,IC1) + CPBG * TX3(I)
1634     70 CONTINUE
1635     C
1636     DO 130 I=1,LENC
1637     TX2(I) = HOL(I,K) - HST(I,IC)
1638     TX1(I) = ZERO
1639    
1640     130 CONTINUE
1641     C
1642     C ENTRAINMENT PARAMETER
1643     C
1644     DO 160 L=IC,KM1
1645     DO 160 I=1,LENC
1646     TX1(I) = TX1(I) + (HST(I,IC) - HOL(I,L)) * (ETA(I,L) - ETA(I,L+1))
1647     160 CONTINUE
1648     C
1649     LEN1 = 0
1650     LEN2 = 0
1651     ISAV = 0
1652     DO 195 I=1,LENC
1653     IF (TX1(I) .GT. ZERO .AND. TX2(I) .GT. ZERO
1654     . .AND. rhfrac(i).ne.0.0 ) THEN
1655     LEN1 = LEN1 + 1
1656     IA(LEN1) = I
1657     ALM(LEN1) = TX2(I) / TX1(I)
1658     ENDIF
1659     195 CONTINUE
1660     C
1661     LEN2 = LEN1
1662     if (IC1 .lt. K) then
1663     DO 196 I=1,LENC
1664     IF (TX2(I) .LE. 0.0 .AND. (HOL(I,K) .GT. HST(I,IC1))
1665     . .AND. rhfrac(i).ne.0.0 ) THEN
1666     LEN2 = LEN2 + 1
1667     IA(LEN2) = I
1668     ALM(LEN2) = 0.0
1669     ENDIF
1670     196 CONTINUE
1671     endif
1672     C
1673     IF (LEN2 .EQ. 0) THEN
1674     DO 5010 I=1,LENC*K
1675     HST(I,1) = 0.0
1676     QOL(I,1) = 0.0
1677     5010 CONTINUE
1678     DO 5020 I=1,LENC
1679     PCU(I) = 0.0
1680     5020 CONTINUE
1681     RETURN
1682     ENDIF
1683     LEN11 = LEN1 + 1
1684     C
1685     C NORMALIZED MASSFLUX
1686     C
1687     DO 250 I=1,LEN2
1688     ETA(I,K) = 1.0
1689     II = IA(I)
1690     TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
1691     TX4(I) = PRS(II,K)
1692     250 CONTINUE
1693     C
1694     DO 252 I=LEN11,LEN2
1695     WFN(I) = 0.0
1696     II = IA(I)
1697     IF (HST(II,IC1) .LT. HST(II,IC)) THEN
1698     TX6(I) = (HST(II,IC1)-HOL(II,K))/(HST(II,IC1)-HST(II,IC))
1699     ELSE
1700     TX6(I) = 0.0
1701     ENDIF
1702     TX2(I) = 0.5 * (PRS(II,IC1)+PRS(II,IC1+1)) * (1.0-TX6(I))
1703     * + TX2(I) * TX6(I)
1704     252 CONTINUE
1705     C
1706     CALL ACRITN(LEN2, TX2, TX4, TX3)
1707     C
1708     DO 260 L=KM1,IC,-1
1709     DO 255 I=1,LEN2
1710     TX1(I) = ETA(IA(I),L)
1711     255 CONTINUE
1712     DO 260 I=1,LEN2
1713     ETA(I,L) = 1.0 + ALM(I) * TX1(I)
1714     260 CONTINUE
1715     C
1716     C CLOUD WORKFUNCTION
1717     C
1718     IF (LEN1 .GT. 0) THEN
1719     DO 270 I=1,LEN1
1720     II = IA(I)
1721     WFN(I) = - GAM(II,IC) * (PRJ(II,IC1) - PRH(II,IC))
1722     * * HST(II,IC) * ETA(I,IC1)
1723     270 CONTINUE
1724     ENDIF
1725     C
1726     DO 290 I=1,LEN2
1727     II = IA(I)
1728     TX1(I) = HOL(II,K)
1729     290 CONTINUE
1730     C
1731     IF (IC1 .LE. KM1) THEN
1732    
1733     DO 380 L=KM1,IC1,-1
1734     DO 380 I=1,LEN2
1735     II = IA(I)
1736     TEM = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * HOL(II,L)
1737     C
1738     PCU(I) = PRJ(II,L+1) - PRH(II,L)
1739     TEM1 = ETA(I,L+1) * PCU(I)
1740     TX1(I) = TX1(I)*PCU(I)
1741     C
1742     PCU(I) = PRH(II,L) - PRJ(II,L)
1743     TEM1 = (TEM1 + ETA(I,L) * PCU(I)) * HST(II,L)
1744     TX1(I) = TX1(I) + TEM*PCU(I)
1745     C
1746     WFN(I) = WFN(I) + (TX1(I) - TEM1) * GAM(II,L)
1747     TX1(I) = TEM
1748     380 CONTINUE
1749     ENDIF
1750     C
1751     LENA = 0
1752     IF (LEN1 .GT. 0) THEN
1753     DO 512 I=1,LEN1
1754     II = IA(I)
1755     WFN(I) = WFN(I) + TX1(I) * GAM(II,IC)*(PRJ(II,IC1)-PRH(II,IC))
1756     * - TX3(I)
1757     IF (WFN(I) .GT. 0.0) THEN
1758     LENA = LENA + 1
1759     I1(LENA) = IA(I)
1760     I2(LENA) = I
1761     TX1(LENA) = WFN(I)
1762     TX2(LENA) = QS1(IA(I))
1763     TX6(LENA) = 1.0
1764     ENDIF
1765     512 CONTINUE
1766     ENDIF
1767     LENB = LENA
1768     DO 515 I=LEN11,LEN2
1769     WFN(I) = WFN(I) - TX3(I)
1770     IF (WFN(I) .GT. 0.0 .AND. TX6(I) .GT. 0.0) THEN
1771     LENB = LENB + 1
1772     I1(LENB) = IA(I)
1773     I2(LENB) = I
1774     TX1(LENB) = WFN(I)
1775     TX2(LENB) = QS1(IA(I))
1776     TX4(LENB) = TX6(I)
1777     ENDIF
1778     515 CONTINUE
1779     C
1780     IF (LENB .LE. 0) THEN
1781     DO 5030 I=1,LENC*K
1782     HST(I,1) = 0.0
1783     QOL(I,1) = 0.0
1784     5030 CONTINUE
1785     DO 5040 I=1,LENC
1786     PCU(I) = 0.0
1787     5040 CONTINUE
1788     RETURN
1789     ENDIF
1790    
1791     C
1792     DO 516 I=1,LENB
1793     WFN(I) = TX1(I)
1794     QS1(I) = TX2(I)
1795     516 CONTINUE
1796     C
1797     DO 520 L=IC,K
1798     DO 517 I=1,LENB
1799     TX1(I) = ETA(I2(I),L)
1800     517 CONTINUE
1801     DO 520 I=1,LENB
1802     ETA(I,L) = TX1(I)
1803     520 CONTINUE
1804     C
1805     LENA1 = LENA + 1
1806     C
1807     DO 510 I=1,LENA
1808     II = I1(I)
1809     TX8(I) = HST(II,IC) - HOL(II,IC)
1810     510 CONTINUE
1811     DO 530 I=LENA1,LENB
1812     II = I1(I)
1813     TX6(I) = TX4(I)
1814     TEM = TX6(I) * (HOL(II,IC)-HOL(II,IC1)) + HOL(II,IC1)
1815     TX8(I) = HOL(II,K) - TEM
1816    
1817     TEM1 = TX6(I) * (QOL(II,IC)-QOL(II,IC1)) + QOL(II,IC1)
1818     TX5(I) = TEM - TEM1 * ALHL
1819     QS1(I) = TEM1 + TX8(I)*(ONE/ALHL)
1820     TX3(I) = HOL(II,IC)
1821     530 CONTINUE
1822     C
1823     C
1824     DO 620 I=1,LENB
1825     II = I1(I)
1826     WLQ(I) = QOL(II,K) - QS1(I) * ETA(I,IC)
1827     TX7(I) = HOL(II,K)
1828     620 CONTINUE
1829     DO NT=1,Ntracer
1830     DO 621 I=1,LENB
1831     II = I1(I)
1832     UHT(I,NT) = UOI(II,K+nltop-1,NT)-UOI(II,IC+nltop-1,NT) * ETA(I,IC)
1833     621 CONTINUE
1834     ENDDO
1835     C
1836     DO 635 L=KM1,IC,-1
1837     DO 630 I=1,LENB
1838     II = I1(I)
1839     TEM = ETA(I,L) - ETA(I,L+1)
1840     WLQ(I) = WLQ(I) + TEM * QOL(II,L)
1841     630 CONTINUE
1842     635 CONTINUE
1843     DO NT=1,Ntracer
1844     DO L=KM1,IC,-1
1845     DO I=1,LENB
1846     II = I1(I)
1847     TEM = ETA(I,L) - ETA(I,L+1)
1848     UHT(I,NT) = UHT(I,NT) + TEM * UOI(II,L+nltop-1,NT)
1849     ENDDO
1850     ENDDO
1851     ENDDO
1852     C
1853     C CALCULATE GS AND PART OF AKM (THAT REQUIRES ETA)
1854     C
1855     DO 690 I=1,LENB
1856     II = I1(I)
1857     c TX7(I) = HOL(II,K)
1858     TEM = (POI(II,KM1) - POI(II,K)) / (PRH(II,K) - PRH(II,KM1))
1859     HOL(I,K) = TEM * (PRJ(II,K)-PRH(II,KM1))*PRH(II,K)*PRI(II,K)
1860     HOL(I,KM1) = TEM * (PRH(II,K)-PRJ(II,K))*PRH(II,KM1)*PRI(II,KM1)
1861     AKM(I) = ZERO
1862     TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
1863     690 CONTINUE
1864    
1865     IF (IC1 .LE. KM1) THEN
1866     DO 750 L=KM1,IC1,-1
1867     DO 750 I=1,LENB
1868     II = I1(I)
1869     TEM = (POI(II,L-1) - POI(II,L)) * ETA(I,L)
1870     * / (PRH(II,L) - PRH(II,L-1))
1871     C
1872     HOL(I,L) = TEM * (PRJ(II,L)-PRH(II,L-1)) * PRH(II,L)
1873     * * PRI(II,L) + HOL(I,L)
1874     HOL(I,L-1) = TEM * (PRH(II,L)-PRJ(II,L)) * PRH(II,L-1)
1875     * * PRI(II,L-1)
1876     C
1877     AKM(I) = AKM(I) - HOL(I,L)
1878     * * (ETA(I,L) * (PRH(II,L)-PRJ(II,L)) +
1879     * ETA(I,L+1) * (PRJ(II,L+1)-PRH(II,L))) / PRH(II,L)
1880     750 CONTINUE
1881     ENDIF
1882     C
1883     C
1884     CALL RNCL(LENB, TX2, TX1, CLF)
1885    
1886     DO 770 I=1,LENB
1887     TX2(I) = (ONE - TX1(I)) * WLQ(I)
1888     WLQ(I) = TX1(I) * WLQ(I)
1889     C
1890     TX1(I) = HOL(I,IC)
1891     770 CONTINUE
1892     DO 790 I=LENA1, LENB
1893     II = I1(I)
1894     TX1(I) = TX1(I) + (TX5(I)-TX3(I)+QOL(II,IC)*ALHL)*(PRI(II,IC)/CP)
1895     790 CONTINUE
1896    
1897     DO 800 I=1,LENB
1898     HOL(I,IC) = TX1(I) - TX2(I) * ALBCP * PRI(I1(I),IC)
1899     800 CONTINUE
1900    
1901     IF (LENA .GT. 0) THEN
1902     DO 810 I=1,LENA
1903     II = I1(I)
1904     AKM(I) = AKM(I) - ETA(I,IC1) * (PRJ(II,IC1) - PRH(II,IC))
1905     * * TX1(I) / PRH(II,IC)
1906     810 CONTINUE
1907     ENDIF
1908     c
1909     C CALCULATE GH
1910     C
1911     DO 830 I=1,LENB
1912     II = I1(I)
1913     TX3(I) = QOL(II,KM1) - QOL(II,K)
1914     GMH(I,K) = HOL(I,K) + TX3(I) * PRI(II,K) * (ALBCP)
1915    
1916     AKM(I) = AKM(I) + GAM(II,KM1)*(PRJ(II,K)-PRH(II,KM1))
1917     * * GMH(I,K)
1918     TX3(I) = zero
1919     830 CONTINUE
1920     C
1921     IF (IC1 .LE. KM1) THEN
1922     DO 840 L=KM1,IC1,-1
1923     DO 840 I=1,LENB
1924     II = I1(I)
1925     TX2(I) = TX3(I)
1926     TX3(I) = (QOL(II,L-1) - QOL(II,L)) * ETA(I,L)
1927     TX2(I) = TX2(I) + TX3(I)
1928     C
1929     GMH(I,L) = HOL(I,L) + TX2(I) * PRI(II,L) * (ALBCP*HALF)
1930     840 CONTINUE
1931     C
1932     C
1933     ENDIF
1934     DO 850 I=LENA1,LENB
1935     TX3(I) = TX3(I) + TWOBAL
1936     * * (TX7(I) - TX8(I) - TX5(I) - QOL(I1(I),IC)*ALHL)
1937     850 CONTINUE
1938     DO 860 I=1,LENB
1939     GMH(I,IC) = TX1(I) + PRI(I1(I),IC) * ONEBCP
1940     * * (TX3(I)*(ALHL*HALF) + ETA(I,IC) * TX8(I))
1941     860 CONTINUE
1942     C
1943     C CALCULATE HC PART OF AKM
1944     C
1945     IF (IC1 .LE. KM1) THEN
1946     DO 870 I=1,LENB
1947     TX1(I) = GMH(I,K)
1948     870 CONTINUE
1949     DO 3725 L=KM1,IC1,-1
1950     DO 880 I=1,LENB
1951     II = I1(I)
1952     TX1(I) = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * GMH(I,L)
1953     TX2(I) = GAM(II,L-1) * (PRJ(II,L) - PRH(II,L-1))
1954     880 CONTINUE
1955     C
1956     IF (L .EQ. IC1) THEN
1957     DO 890 I=LENA1,LENB
1958     TX2(I) = ZERO
1959     890 CONTINUE
1960     ENDIF
1961     DO 900 I=1,LENB
1962     II = I1(I)
1963     AKM(I) = AKM(I) + TX1(I) *
1964     * (TX2(I) + GAM(II,L)*(PRH(II,L)-PRJ(II,L)))
1965     900 CONTINUE
1966     3725 CONTINUE
1967     ENDIF
1968     C
1969     DO 920 I=LENA1,LENB
1970     II = I1(I)
1971     TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
1972     * + 0.5*(PRS(II,IC+2) - PRS(II,IC)) * (ONE-TX6(I))
1973     c
1974     TX1(I) = PRS(II,IC1)
1975     TX5(I) = 0.5 * (PRS(II,IC1) + PRS(II,IC+2))
1976     C
1977     IF ((TX2(I) .GE. TX1(I)) .AND. (TX2(I) .LT. TX5(I))) THEN
1978     TX6(I) = ONE - (TX2(I) - TX1(I)) / (TX5(I) - TX1(I))
1979     C
1980     TEM = PRI(II,IC1) / PRI(II,IC)
1981     HOL(I,IC1) = HOL(I,IC1) + HOL(I,IC) * TEM
1982     HOL(I,IC) = ZERO
1983     C
1984     GMH(I,IC1) = GMH(I,IC1) + GMH(I,IC) * TEM
1985     GMH(I,IC) = ZERO
1986     ELSEIF (TX2(I) .LT. TX1(I)) THEN
1987     TX6(I) = 1.0
1988     ELSE
1989     TX6(I) = 0.0
1990     ENDIF
1991     920 CONTINUE
1992     C
1993     C
1994     DO I=1,LENC
1995     PCU(I) = 0.0
1996     ENDDO
1997    
1998     DO 970 I=1,LENB
1999     II = I1(I)
2000     IF (AKM(I) .LT. ZERO .AND. WLQ(I) .GE. 0.0) THEN
2001     WFN(I) = - TX6(I) * WFN(I) * RASALF / AKM(I)
2002     ELSE
2003     WFN(I) = ZERO
2004     ENDIF
2005     TEM = (PRS(II,K+1)-PRS(II,K))*(CMB2PA*FRAC)
2006     WFN(I) = MIN(WFN(I), TEM)
2007     C
2008     C compute cloud amount
2009     C
2010     CC TX1(I) = CLN(II)
2011     CC IF (WFN(I) .GT. CRTMSF) TX1(I) = TX1(I) + CLF(I)
2012     CC IF (TX1(I) .GT. ONE) TX1(I) = ONE
2013     C
2014     C PRECIPITATION
2015     C
2016     PCU(II) = WLQ(I) * WFN(I) * ONEBG
2017     C
2018     C CUMULUS FRICTION AT THE BOTTOM LAYER
2019     C
2020     TX4(I) = WFN(I) * (1.0/ALHL)
2021     TX5(I) = WFN(I) * ONEBCP
2022     970 CONTINUE
2023     C
2024     C compute cloud mass flux for diagnostic output
2025     C
2026     DO L = IC,K
2027     DO I=1,LENB
2028     II = I1(I)
2029     if(L.lt.K)then
2030     CMASS(II,L) = ETA(I,L+1) * WFN(I) * ONEBG
2031     else
2032     CMASS(II,L) = WFN(I) * ONEBG
2033     endif
2034     ENDDO
2035     ENDDO
2036     C
2037     CC DO 975 I=1,LENB
2038     CC II = I1(I)
2039     CC CLN(II) = TX1(I)
2040     CC975 CONTINUE
2041     C
2042     C THETA AND Q CHANGE DUE TO CLOUD TYPE IC
2043     C
2044    
2045     c TEMA = 0.0
2046     c TEMB = 0.0
2047     DO 990 L=IC,K
2048     DO 980 I=1,LENB
2049     II = I1(I)
2050     TEM = (GMH(I,L) - HOL(I,L)) * TX4(I)
2051     TEM1 = HOL(I,L) * TX5(I)
2052     C
2053     TCU(II,L) = TEM1 / PRH(II,L)
2054     QCU(II,L) = TEM
2055     980 CONTINUE
2056    
2057     c I = I1(IP1)
2058     c
2059     c TEM = (PRS(I,L+1)-PRS(I,L)) * (ONEBG*100.0)
2060     c TEMA = TEMA + TCU(I,L) * PRH(I,L) * TEM * (CP/ALHL)
2061     c TEMB = TEMB + QCU(I,L) * TEM
2062     C
2063     990 CONTINUE
2064     C
2065     c Compute Tracer Tendencies
2066     c -------------------------
2067     do nt = 1,ntracer
2068    
2069     c Tracer Tendency at the Bottom Layer
2070     c -----------------------------------
2071     DO 995 I=1,LENB
2072     II = I1(I)
2073     TEM = half*TX5(I) * PRI(II,K)
2074     TX1(I) = (UOI(II,KM1+nltop-1,nt) - UOI(II,K+nltop-1,nt))
2075     ucu(II,K,nt) = TEM * TX1(I)
2076     995 CONTINUE
2077    
2078     c Tracer Tendency at all other Levels
2079     c -----------------------------------
2080     DO 1020 L=KM1,IC1,-1
2081     DO 1010 I=1,LENB
2082     II = I1(I)
2083     TEM = half*TX5(I) * PRI(II,L)
2084     TEM1 = TX1(I)
2085     TX1(I) = (UOI(II,L-1+nltop-1,nt)-UOI(II,L+nltop-1,nt)) * ETA(I,L)
2086     TX3(I) = (TX1(I) + TEM1) * TEM
2087     1010 CONTINUE
2088     DO 1020 I=1,LENB
2089     II = I1(I)
2090     ucu(II,L,nt) = TX3(I)
2091     1020 CONTINUE
2092    
2093     DO 1030 I=1,LENB
2094     II = I1(I)
2095     IF (TX6(I) .GE. 1.0) THEN
2096     TEM = half*TX5(I) * PRI(II,IC)
2097     ELSE
2098     TEM = 0.0
2099     ENDIF
2100     TX1(I) = (TX1(I) + UHT(I,nt) + UHT(I,nt)) * TEM
2101     1030 CONTINUE
2102     DO 1040 I=1,LENB
2103     II = I1(I)
2104     ucu(II,IC,nt) = TX1(I)
2105     1040 CONTINUE
2106    
2107     enddo
2108     C
2109     C PENETRATIVE CONVECTION CALCULATION OVER
2110     C
2111    
2112     RETURN
2113     END
2114 molod 1.16 SUBROUTINE RNCL(lng, PL, RNO, CLF)
2115 molod 1.1 C
2116     C*********************************************************************
2117     C********************** Relaxed Arakawa-Schubert *********************
2118     C************************ SUBROUTINE RNCL ************************
2119     C**************************** 23 July 1992 ***************************
2120     C*********************************************************************
2121 molod 1.9 implicit none
2122     C Argument List declarations
2123 molod 1.16 integer lng
2124     _RL PL(lng), RNO(lng), CLF(lng)
2125 molod 1.1
2126 molod 1.9 C Local Variables
2127 molod 1.13 _RL p5,p8,pt8,pt2,pfac,p4,p6,p7,p9,cucld,cfac
2128 molod 1.1 PARAMETER (P5=500.0, P8=800.0, PT8=0.8, PT2=0.2)
2129     PARAMETER (PFAC=PT2/(P8-P5))
2130     PARAMETER (P4=400.0, P6=401.0)
2131     PARAMETER (P7=700.0, P9=900.0)
2132     PARAMETER (CUCLD=0.5,CFAC=CUCLD/(P6-P4))
2133 molod 1.9
2134     integer i
2135 molod 1.1 C
2136 molod 1.16 DO 10 I=1,lng
2137 molod 1.1 rno(i) = 1.0
2138     ccc if( pl(i).le.400.0 ) rno(i) = max( 0.75, 1.0-0.0025*(400.0-pl(i)) )
2139    
2140     ccc IF ( PL(I).GE.P7 .AND. PL(I).LE.P9 ) THEN
2141     ccc RNO(I) = ((P9-PL(I))/(P9-P7)) **2
2142     ccc ELSE IF (PL(I).GT.P9) THEN
2143     ccc RNO(I) = 0.
2144     ccc ENDIF
2145    
2146     CLF(I) = CUCLD
2147     C
2148     CARIESIF (PL(I) .GE. P5 .AND. PL(I) .LE. P8) THEN
2149     CARIES RNO(I) = (P8-PL(I))*PFAC + PT8
2150     CARIESELSEIF (PL(I) .GT. P8 ) THEN
2151     CARIES RNO(I) = PT8
2152     CARIESENDIF
2153     CARIES
2154     IF (PL(I) .GE. P4 .AND. PL(I) .LE. P6) THEN
2155     CLF(I) = (P6-PL(I))*CFAC
2156     ELSEIF (PL(I) .GT. P6 ) THEN
2157     CLF(I) = 0.0
2158     ENDIF
2159     10 CONTINUE
2160     C
2161     RETURN
2162     END
2163 molod 1.16 SUBROUTINE ACRITN ( lng,PL,PLB,ACR )
2164 molod 1.1
2165     C*********************************************************************
2166     C********************** Relaxed Arakawa-Schubert *********************
2167     C************************** SUBROUTINE ACRIT *********************
2168     C****************** modified August 28, 1996 L.Takacs ************
2169     C**** *****
2170     C**** Note: Data obtained from January Mean After-Analysis *****
2171     C**** from 4x5 46-layer GEOS Assimilation *****
2172     C**** *****
2173     C*********************************************************************
2174 molod 1.9 implicit none
2175     C Argument List declarations
2176 molod 1.16 integer lng
2177     _RL PL(lng), PLB(lng), ACR(lng)
2178 molod 1.1
2179 molod 1.9 C Local variables
2180     integer lma
2181 molod 1.1 parameter (lma=18)
2182 molod 1.13 _RL p(lma)
2183     _RL a(lma)
2184 molod 1.9 integer i,L
2185 molod 1.13 _RL temp
2186 molod 1.1
2187     data p / 93.81, 111.65, 133.46, 157.80, 186.51,
2188     . 219.88, 257.40, 301.21, 352.49, 409.76,
2189     . 471.59, 535.04, 603.33, 672.79, 741.12,
2190     . 812.52, 875.31, 930.20/
2191    
2192     data a / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364,
2193     . 1.01971, .65846, .45867, .38687, .31002,
2194     . .25574, .20347, .17254, .15260, .16756,
2195     . .09916, .10360, .05880/
2196    
2197    
2198     do L=1,lma-1
2199 molod 1.16 do i=1,lng
2200 molod 1.1 if( pl(i).ge.p(L) .and.
2201     . pl(i).le.p(L+1)) then
2202     temp = ( pl(i)-p(L) )/( p(L+1)-p(L) )
2203     acr(i) = a(L+1)*temp + a(L)*(1-temp)
2204     endif
2205     enddo
2206     enddo
2207    
2208 molod 1.16 do i=1,lng
2209 molod 1.1 if( pl(i).lt.p(1) ) acr(i) = a(1)
2210     if( pl(i).gt.p(lma) ) acr(i) = a(lma)
2211     enddo
2212    
2213 molod 1.16 do i=1,lng
2214 molod 1.1 acr(i) = acr(i) * (plb(i)-pl(i))
2215     enddo
2216    
2217     RETURN
2218     END
2219 molod 1.6 SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DP,PLKE,
2220 molod 1.1 1 PLK,TH,TEMP1,TEMP2,TEMP3,ITMP1,ITMP2,RCON,RLAR,CLSBTH,tmscl,
2221     2 tmfrc,cp,gravity,alhl,gamfac,cldlz,RHCRIT,offset,alpha)
2222    
2223 molod 1.9 implicit none
2224     C Argument List declarations
2225     integer nn,irun,nlay
2226 molod 1.13 _RL TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY),
2227 molod 1.9 . PL(IRUN,NLAY),CLFRAC(IRUN,NLAY),SP(IRUN),TEMP1(IRUN,NLAY),
2228     . TEMP2(IRUN,NLAY),PLKE(IRUN,NLAY+1),
2229     . RCON(IRUN),RLAR(IRUN),DP(IRUN,NLAY),PLK(IRUN,NLAY),TH(IRUN,NLAY),
2230     . TEMP3(IRUN,NLAY)
2231     integer ITMP1(IRUN,NLAY),ITMP2(IRUN,NLAY)
2232 molod 1.13 _RL CLSBTH(IRUN,NLAY)
2233     _RL tmscl,tmfrc,cp,gravity,alhl,gamfac,offset,alpha
2234     _RL cldlz(irun,nlay)
2235     _RL rhcrit(irun,nlay)
2236 molod 1.9 C
2237     C Local Variables
2238 molod 1.13 _RL zm1p04,zero,two89,zp44,zp01,half,zp578,one,thousand,z3600
2239     _RL zp1,zp001
2240 molod 1.1 PARAMETER (ZM1P04 = -1.04E-4 )
2241     PARAMETER (ZERO = 0.)
2242     PARAMETER (TWO89= 2.89E-5)
2243     PARAMETER ( ZP44= 0.44)
2244     PARAMETER ( ZP01= 0.01)
2245     PARAMETER ( ZP1 = 0.1 )
2246     PARAMETER ( ZP001= 0.001)
2247     PARAMETER ( HALF= 0.5)
2248     PARAMETER ( ZP578 = 0.578 )
2249     PARAMETER ( ONE = 1.)
2250     PARAMETER ( THOUSAND = 1000.)
2251     PARAMETER ( Z3600 = 3600.)
2252     C
2253 molod 1.13 _RL EVP9(IRUN,NLAY)
2254     _RL water(irun),crystal(irun)
2255     _RL watevap(irun),iceevap(irun)
2256     _RL fracwat,fracice, tice,rh,fact,dum
2257     _RL rainmax(irun)
2258     _RL getcon,rphf,elocp,cpog,relax
2259     _RL exparg,arearat,rpow
2260 molod 1.9
2261     integer i,L,n,nlaym1,irnlay,irnlm1
2262 molod 1.1
2263     c Explicit Inline Directives
2264     c --------------------------
2265 molod 1.10 #ifdef CRAY
2266     #ifdef f77
2267 molod 1.1 cfpp$ expand (qsat)
2268     #endif
2269     #endif
2270    
2271     tice = getcon('FREEZING-POINT')
2272    
2273     fracwat = 0.70
2274     fracice = 0.01
2275    
2276     NLAYM1 = NLAY - 1
2277     IRNLAY = IRUN*NLAY
2278     IRNLM1 = IRUN*(NLAY-1)
2279    
2280     RPHF = Z3600/tmscl
2281    
2282     ELOCP = alhl/cp
2283     CPOG = cp/gravity
2284    
2285     DO I = 1,IRUN
2286     RLAR(I) = 0.
2287     water(i) = 0.
2288     crystal(i) = 0.
2289     ENDDO
2290    
2291     do L = 1,nlay
2292     do i = 1,irun
2293     EVP9(i,L) = 0.
2294     TEMP1(i,L) = 0.
2295     TEMP2(i,L) = 0.
2296     TEMP3(i,L) = 0.
2297     CLSBTH(i,L) = 0.
2298     cldlz(i,L) = 0.
2299     enddo
2300     enddo
2301    
2302     C RHO(ZERO) / RHO FOR TERMINAL VELOCITY APPROX.
2303     c ---------------------------------------------
2304     DO L = 1,NLAY
2305     DO I = 1,IRUN
2306     TEMP2(I,L) = PL(I,L)*ZP001
2307     TEMP2(I,L) = SQRT(TEMP2(I,L))
2308     ENDDO
2309     ENDDO
2310    
2311     C INVERSE OF MASS IN EACH LAYER
2312     c -----------------------------
2313     DO L = 1,NLAY
2314     DO I = 1,IRUN
2315 molod 1.6 TEMP3(I,L) = GRAVITY*ZP01 / DP(I,L)
2316 molod 1.1 ENDDO
2317     ENDDO
2318    
2319     C DO LOOP FOR MOISTURE EVAPORATION ABILITY AND CONVEC EVAPORATION.
2320     c ----------------------------------------------------------------
2321     DO 100 L=1,NLAY
2322    
2323     DO I = 1,IRUN
2324     TEMP1(I,3) = TL(I,L)
2325     TEMP1(I,4) = QL(I,L)
2326     ENDDO
2327    
2328     DO 50 N=1,2
2329     IF(N.EQ.1)RELAX=HALF
2330     IF(N.GT.1)RELAX=ONE
2331    
2332     DO I = 1,IRUN
2333     call qsat ( temp1(i,3),pl(i,L),temp1(i,2),temp1(i,6),.true. )
2334     TEMP1(I,5)=TEMP1(I,2)-TEMP1(I,4)
2335     TEMP1(I,6)=TEMP1(I,6)*ELOCP
2336     TEMP1(I,5)=TEMP1(I,5)/(ONE+TEMP1(I,6))
2337     TEMP1(I,4)=TEMP1(I,4)+TEMP1(I,5)*RELAX
2338     TEMP1(I,3)=TEMP1(I,3)-TEMP1(I,5)*ELOCP*RELAX
2339     ENDDO
2340     50 CONTINUE
2341    
2342     DO I = 1,IRUN
2343     EVP9(I,L) = (TEMP1(I,4) - QL(I,L))/TEMP3(I,L)
2344     C convective detrained water
2345     cldlz(i,L) = rain(i,L)*temp3(i,L)
2346     if( tl(i,L).gt.tice-20.) then
2347     water(i) = water(i) + rain(i,L)
2348     else
2349     crystal(i) = crystal(i) + rain(i,L)
2350     endif
2351     ENDDO
2352    
2353     C**********************************************************************
2354     C FOR CONVECTIVE PRECIP, FIND THE "EVAPORATION EFFICIENCY" USING *
2355     C KESSLERS PARAMETERIZATION *
2356     C**********************************************************************
2357    
2358     DO 20 I=1,IRUN
2359    
2360     iceevap(i) = 0.
2361     watevap(i) = 0.
2362    
2363     if( (evp9(i,L).gt.0.) .and. (crystal(i).gt.0.) ) then
2364     iceevap(I) = EVP9(I,L)*fracice
2365     IF(iceevap(i).GE.crystal(i)) iceevap(i) = crystal(i)
2366     EVP9(I,L)=EVP9(I,L)-iceevap(I)
2367     crystal(i) = crystal(i) - iceevap(i)
2368     endif
2369    
2370     C and now warm precipitate
2371     if( (evp9(i,L).gt.0.) .and. (water(i).gt.0.) ) then
2372     exparg = ZM1P04*tmscl*((water(i)*RPHF*TEMP2(I,L))**ZP578)
2373     AREARAT = ONE-(EXP(EXPARG))
2374     watevap(I) = EVP9(I,L)*AREARAT*fracwat
2375     IF(watevap(I).GE.water(i)) watevap(I) = water(i)
2376     EVP9(I,L)=EVP9(I,L)-watevap(I)
2377     water(i) = water(i) - watevap(i)
2378     endif
2379    
2380     QL(I,L) = QL(I,L)+(iceevap(i)+watevap(i))*TEMP3(I,L)
2381     TL(I,L) = TL(I,L)-(iceevap(i)+watevap(i))*TEMP3(I,L)*ELOCP
2382    
2383     20 CONTINUE
2384    
2385     100 CONTINUE
2386    
2387     do i = 1,irun
2388     rcon(i) = water(i) + crystal(i)
2389     enddo
2390    
2391     C**********************************************************************
2392     C Large Scale Precip
2393     C**********************************************************************
2394    
2395     DO 200 L=1,NLAY
2396     DO I = 1,IRUN
2397     rainmax(i) = rhcrit(i,L)*evp9(i,L) +
2398     . ql(i,L)*(rhcrit(i,L)-1.)/temp3(i,L)
2399    
2400     if (rainmax(i).LE.0.0) then
2401     call qsat( tl(i,L),pl(i,L),rh,dum,.false.)
2402     rh = ql(i,L)/rh
2403    
2404     if( rhcrit(i,L).eq.1.0 ) then
2405     fact = 1.0
2406     else
2407     fact = min( 1.0, alpha + (1.0-alpha)*( rh-rhcrit(i,L)) /
2408     1 (1.0-rhcrit(i,L)) )
2409     endif
2410    
2411     C Do not allow clouds above 10 mb
2412     if( pl(i,L).ge.10.0 ) CLSBTH(I,L) = fact
2413     RLAR(I) = RLAR(I)-rainmax(I)
2414     QL(I,L) = QL(I,L)+rainmax(I)*TEMP3(I,L)
2415     TL(I,L) = TL(I,L)-rainmax(I)*TEMP3(I,L)*ELOCP
2416     C Large-scale water
2417     cldlz(i,L) = cldlz(i,L) - rainmax(i)*temp3(i,L)
2418     ENDIF
2419     ENDDO
2420    
2421     DO I=1,IRUN
2422     IF((RLAR(I).GT.0.0).AND.(rainmax(I).GT.0.0))THEN
2423     RPOW=(RLAR(I)*RPHF*TEMP2(I,L))**ZP578
2424     EXPARG = ZM1P04*tmscl*RPOW
2425     AREARAT = ONE-(EXP(EXPARG))
2426     TEMP1(I,7) = rainmax(I)*AREARAT
2427     IF(TEMP1(I,7).GE.RLAR(I)) TEMP1(I,7) = RLAR(I)
2428     RLAR(I) = RLAR(I)-TEMP1(I,7)
2429     QL(I,L) = QL(I,L)+TEMP1(I,7)*TEMP3(I,L)
2430     TL(I,L) = TL(I,L)-TEMP1(I,7)*TEMP3(I,L)*ELOCP
2431     ENDIF
2432     ENDDO
2433    
2434     200 CONTINUE
2435    
2436     RETURN
2437     END
2438     subroutine srclouds (th,q,plk,pl,plke,cloud,cldwat,irun,irise,
2439     1 rhc,offset,alpha)
2440     C***********************************************************************
2441     C
2442     C PURPOSE:
2443     C ========
2444     C Compute non-precipitating cloud fractions
2445     C based on Slingo and Ritter (1985).
2446     C Remove cloudiness where conditionally unstable.
2447     C
2448     C INPUT:
2449     C ======
2450     C th ......... Potential Temperature (irun,irise)
2451     C q .......... Specific Humidity (irun,irise)
2452     C plk ........ P**Kappa at mid-layer (irun,irise)
2453     C pl ......... Pressure at mid-layer (irun,irise)
2454     C plke ....... P**Kappa at edge (irun,irise+1)
2455     C irun ....... Horizontal dimension
2456     C irise ...... Vertical dimension
2457     C
2458     C OUTPUT:
2459     C =======
2460     C cloud ...... Cloud Fraction (irun,irise)
2461     C
2462     C***********************************************************************
2463    
2464     implicit none
2465     integer irun,irise
2466    
2467 molod 1.13 _RL th(irun,irise)
2468     _RL q(irun,irise)
2469     _RL plk(irun,irise)
2470     _RL pl(irun,irise)
2471     _RL plke(irun,irise+1)
2472    
2473     _RL cloud(irun,irise)
2474     _RL cldwat(irun,irise)
2475     _RL qs(irun,irise)
2476    
2477     _RL cp, alhl, getcon, akap
2478     _RL ratio, temp, elocp
2479     _RL rhcrit,rh,dum
2480 molod 1.9 integer i,L
2481 molod 1.1
2482 molod 1.13 _RL rhc(irun,irise)
2483     _RL offset,alpha
2484 molod 1.1
2485     c Explicit Inline Directives
2486     c --------------------------
2487 molod 1.10 #ifdef CRAY
2488     #ifdef f77
2489 molod 1.1 cfpp$ expand (qsat)
2490     #endif
2491     #endif
2492    
2493     cp = getcon('CP')
2494     alhl = getcon('LATENT HEAT COND')
2495     elocp = alhl/cp
2496     akap = getcon('KAPPA')
2497    
2498     do L = 1,irise
2499     do i = 1,irun
2500     temp = th(i,L)*plk(i,L)
2501     call qsat ( temp,pl(i,L),qs(i,L),dum,.false. )
2502     enddo
2503     enddo
2504    
2505     do L = 2,irise
2506     do i = 1,irun
2507     rh = q(i,L)/qs(i,L)
2508    
2509     rhcrit = rhc(i,L) - offset
2510     ratio = alpha*(rh-rhcrit)/offset
2511    
2512     if(cloud(i,L).eq. 0.0 .and. ratio.gt.0.0 ) then
2513     cloud(i,L) = min( ratio,1.0 )
2514     endif
2515    
2516     enddo
2517     enddo
2518    
2519     c Reduce clouds from conditionally unstable layer
2520     c -----------------------------------------------
2521     call ctei ( th,q,cloud,cldwat,pl,plk,plke,irun,irise )
2522    
2523     return
2524     end
2525    
2526     subroutine ctei ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm )
2527     implicit none
2528     integer im,lm
2529 molod 1.13 _RL th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm)
2530     _RL plk(im,lm),pl(im,lm),cldfrc(im,lm)
2531 molod 1.1 integer i,L
2532 molod 1.13 _RL getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq
2533     _RL k,krd,kmm,f
2534 molod 1.1
2535     cp = getcon('CP')
2536     alhl = getcon('LATENT HEAT COND')
2537     cpoel = cp/alhl
2538     elocp = alhl/cp
2539    
2540     do L=lm,2,-1
2541     do i=1,im
2542     dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1)
2543     if( dq.eq.0.0 ) dq = 1.0e-20
2544     k = 1.0 + cpoel*plke(i,L)*( th(i,L)-th(i,L-1) ) / dq
2545    
2546     t = th(i,L)*plk(i,L)
2547     p = pl(i,L)
2548     call qsat ( t,p,qs,dqsdt,.true. )
2549    
2550     krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t )
2551    
2552     kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t )
2553     . / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt )
2554    
2555     s = ( (k-krd)/(kmm-krd) )
2556     f = 1.0 - min( 1.0, max(0.0,1.0-exp(-s)) )
2557    
2558     cldfrc(i,L) = cldfrc(i,L)*f
2559     cldwat(i,L) = cldwat(i,L)*f
2560    
2561     enddo
2562     enddo
2563    
2564     return
2565     end
2566    
2567     subroutine back2grd(gathered,indeces,scattered,irun)
2568     implicit none
2569     integer i,irun,indeces(irun)
2570 molod 1.13 _RL gathered(irun),scattered(irun)
2571     _RL temp(irun)
2572 molod 1.1 do i = 1,irun
2573     temp(indeces(i)) = gathered(i)
2574     enddo
2575     do i = 1,irun
2576     scattered(i) = temp(i)
2577     enddo
2578     return
2579     end
  ViewVC Help
Powered by ViewVC 1.1.22