/[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.12 - (hide annotations) (download)
Sun Jul 18 23:17:00 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.11: +9 -5 lines
Fix last issues: save statements, random number generator

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

  ViewVC Help
Powered by ViewVC 1.1.22