/[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.15 - (hide annotations) (download)
Wed Jul 28 01:25:07 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.14: +5 -5 lines 
debugging
1 molod 1.15 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.14 2004/07/26 19:51:08 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.1 integer i,j,iloop,index,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     index = 0
234     do L = nsubmin,nltop,-1
235     do j = 1,jm
236     do i = 1,im
237     if(levpbl(i,j).eq.L) then
238     index = index + 1
239     pblindex(index) = (j-1)*im + i
240     endif
241     enddo
242     enddo
243     enddo
244    
245     do index = 1,im*jm
246     levgather(index) = levpbl(pblindex(index),1)
247     pigather(index) = pz(pblindex(index),1)
248 molod 1.6 pkegather(index,lm+1) = pkht(pblindex(index),1,lm+1)
249 molod 1.8 plegather(index,lm+1) = plze(pblindex(index),1,lm+1)
250 molod 1.1 enddo
251    
252     do L = 1,lm
253     do index = 1,im*jm
254 molod 1.6 thgather(index,L) = tz(pblindex(index),1,L)
255     shgather(index,L) = qz(pblindex(index),1,L,1)
256 molod 1.1 pkegather(index,L) = pkht(pblindex(index),1,L)
257 molod 1.6 pkzgather(index,L) = pkl(pblindex(index),1,L)
258     plegather(index,L) = plze(pblindex(index),1,L)
259     plzgather(index,L) = plz(pblindex(index),1,L)
260     dpgather(index,L) = dpres(pblindex(index),1,L)
261 molod 1.1 enddo
262     enddo
263     do nt = 1,ntracer-ptracer
264     do L = 1,lm
265     do index = 1,im*jm
266     ugather(index,L,nt) = qz(pblindex(index),1,L,nt+ptracer)
267     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     c Sub-Cloud Layer
780     c -------------------------
781     if( ipsubcld.ne.0 ) then
782     do j = 1,jm
783     do i = 1,im
784 molod 1.3 qdiag(i,j,ipsubcld,bi,bj) = qdiag(i,j,ipsubcld,bi,bj) +
785     . psubcldg (i,j)
786     qdiag(i,j,ipsubcldc,bi,bj) = qdiag(i,j,ipsubcldc,bi,bj) +
787     . psubcldgc(i,j)
788 molod 1.1 enddo
789     enddo
790     endif
791    
792     c Non-Precipitating Cloud Fraction
793     c --------------------------------
794     if( icldnp.ne.0 ) then
795     do L = 1,lm
796     do j = 1,jm
797     do i = 1,im
798 molod 1.3 qdiag(i,j,icldnp+L-1,bi,bj) = qdiag(i,j,icldnp+L-1,bi,bj) +
799     . cldsr(i,j,L)
800 molod 1.1 enddo
801     enddo
802     enddo
803     ncldnp = ncldnp + 1
804     endif
805    
806     c Moist Processes Heating Rate
807     c ----------------------------
808     if(imoistt.gt.0) then
809     do L = 1,lm
810     do i = 1,im*jm
811 molod 1.3 qdiag(i,1,imoistt+L-1,bi,bj) = qdiag(i,1,imoistt+L-1,bi,bj) +
812 molod 1.1 . (dtmoist(i,1,L)*sday*pkzgather(i,L)/pz(i,1))
813     enddo
814     enddo
815     endif
816    
817     c Moist Processes Moistening Rate
818     c -------------------------------
819     if(imoistq.gt.0) then
820     do L = 1,lm
821     do j = 1,jm
822     do i = 1,im
823 molod 1.3 qdiag(i,j,imoistq+L-1,bi,bj) = qdiag(i,j,imoistq+L-1,bi,bj) +
824 molod 1.1 . (dqmoist(i,j,L,1)*sday*1000.0/pz(i,j))
825     enddo
826     enddo
827     enddo
828     endif
829    
830     c Cloud Mass Flux
831     c ---------------
832     if(icldmas.gt.0) then
833     do L = 1,lm
834     do i = 1,im*jm
835 molod 1.3 qdiag(i,1,icldmas+L-1,bi,bj) = qdiag(i,1,icldmas+L-1,bi,bj) +
836     . tmpgather(i,L)
837 molod 1.1 enddo
838     enddo
839     endif
840    
841     c Detrained Cloud Mass Flux
842     c -------------------------
843     if(idtrain.gt.0) then
844     do L = 1,lm
845     do i = 1,im*jm
846 molod 1.3 qdiag(i,1,idtrain+L-1,bi,bj) = qdiag(i,1,idtrain+L-1,bi,bj) +
847     . pkegather(i,L)
848 molod 1.1 enddo
849     enddo
850     endif
851    
852     c Grid-Scale Condensational Heating Rate
853     c --------------------------------------
854     if(idtls.gt.0) then
855     do L = 1,lm
856     do i = 1,im*jm
857 molod 1.3 qdiag(i,1,idtls+L-1,bi,bj) = qdiag(i,1,idtls+L-1,bi,bj) +
858     . deltrnev(i,L)
859 molod 1.1 enddo
860     enddo
861     endif
862    
863     c Grid-Scale Condensational Moistening Rate
864     c -----------------------------------------
865     if(idqls.gt.0) then
866     do L = 1,lm
867     do i = 1,im*jm
868 molod 1.3 qdiag(i,1,idqls+L-1,bi,bj) = qdiag(i,1,idqls+L-1,bi,bj) +
869     . delqrnev(i,L)
870 molod 1.1 enddo
871     enddo
872     endif
873    
874     c Total Precipitation
875     c -------------------
876     if(ipreacc.gt.0) then
877     do j = 1,jm
878     do i = 1,im
879 molod 1.3 qdiag(i,j,ipreacc,bi,bj) = qdiag(i,j,ipreacc,bi,bj)
880 molod 1.1 . + ( lsp_new(I,j)
881     . + snow_new(I,j)
882     . + conv_new(i,j) ) *sday*tminv
883     enddo
884     enddo
885     endif
886    
887     c Convective Precipitation
888     c ------------------------
889     if(iprecon.gt.0) then
890     do i = 1,im*jm
891 molod 1.3 qdiag(i,1,iprecon,bi,bj) = qdiag(i,1,iprecon,bi,bj) +
892     . raincgath(i)*sday*tminv
893 molod 1.1 enddo
894     endif
895    
896     C **********************************************************************
897     C **** Fill Rainfall and Snowfall Arrays for Land Surface Model ****
898     C **** Note: Precip Rates work when DT(turb)<DT(moist) ****
899     C **********************************************************************
900    
901     do j = 1,jm
902     do i = 1,im
903     rainlsp (i,j) = rainlsp (i,j) + lsp_new(i,j)*tminv
904     rainconv(i,j) = rainconv(i,j) + conv_new(i,j)*tminv
905     snowfall(i,j) = snowfall(i,j) + snow_new(i,j)*tminv
906     enddo
907     enddo
908    
909     C **********************************************************************
910     C *** Compute Time-averaged Quantities for Radiation ***
911     C *** CPEN => Cloud Fraction from RAS ***
912     C *** CLDLS => Cloud Fraction from RNEVP ***
913     C **********************************************************************
914    
915     do j = 1,jm
916     do i = 1,im
917     cldhi (i,j) = 0.
918     cldmid(i,j) = 0.
919     cldlow(i,j) = 0.
920     cldtmp(i,j) = 0.
921     cldprs(i,j) = 0.
922     tmpimjm(i,j) = 0.
923     enddo
924     enddo
925    
926     c Set Moist-Process Memory Coefficient
927     c ------------------------------------
928     cldras_mem = 1.0-tmstp/ 3600.0
929     cldlsp_mem = 1.0-tmstp/(3600.0*3)
930    
931     do L = 1,lm
932     do i = 1,im*jm
933 molod 1.6 plev = pl(i,L)
934 molod 1.1
935     c Compute Time-averaged Cloud and Water Amounts for Longwave Radiation
936     c --------------------------------------------------------------------
937     watnow = cldwater(i,1,L)
938     if( plev.le.500.0 ) then
939     cldras = min( max( cldras_lw(i,1,L)*cldras_mem,cpen(i,1,L)),1.0)
940     else
941     cldras = 0.0
942     endif
943     cldlsp = min( max( cldlsp_lw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0)
944    
945     if( cldras.lt.cldmin ) cldras = 0.0
946     if( cldlsp.lt.cldmin ) cldlsp = 0.0
947    
948     cldnow = max( cldlsp,cldras )
949    
950     lwlz(i,1,L) = ( nlwlz*lwlz(i,1,L) + watnow)/(nlwlz +1)
951     cldtot_lw(i,1,L) = (nlwcld*cldtot_lw(i,1,L) + cldnow)/(nlwcld+1)
952     cldlsp_lw(i,1,L) = (nlwcld*cldlsp_lw(i,1,L) + cldlsp)/(nlwcld+1)
953     cldras_lw(i,1,L) = (nlwcld*cldras_lw(i,1,L) + cldras)/(nlwcld+1)
954    
955    
956     c Compute Time-averaged Cloud and Water Amounts for Shortwave Radiation
957     c ---------------------------------------------------------------------
958     watnow = cldwater(i,1,L)
959     if( plev.le.500.0 ) then
960     cldras = min( max(cldras_sw(i,1,L)*cldras_mem, cpen(i,1,L)),1.0)
961     else
962     cldras = 0.0
963     endif
964     cldlsp = min( max(cldlsp_sw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0)
965    
966     if( cldras.lt.cldmin ) cldras = 0.0
967     if( cldlsp.lt.cldmin ) cldlsp = 0.0
968    
969     cldnow = max( cldlsp,cldras )
970    
971     swlz(i,1,L) = ( nswlz*swlz(i,1,L) + watnow)/(nswlz +1)
972     cldtot_sw(i,1,L) = (nswcld*cldtot_sw(i,1,L) + cldnow)/(nswcld+1)
973     cldlsp_sw(i,1,L) = (nswcld*cldlsp_sw(i,1,L) + cldlsp)/(nswcld+1)
974     cldras_sw(i,1,L) = (nswcld*cldras_sw(i,1,L) + cldras)/(nswcld+1)
975    
976    
977     c Compute Instantaneous Low-Mid-High Maximum Overlap Cloud Fractions
978     c ----------------------------------------------------------------------
979    
980     if( L.lt.midlevel ) cldhi (i,1) = max( cldnow,cldhi (i,1) )
981     if( L.ge.midlevel .and.
982     . L.lt.lowlevel ) cldmid(i,1) = max( cldnow,cldmid(i,1) )
983     if( L.ge.lowlevel ) cldlow(i,1) = max( cldnow,cldlow(i,1) )
984    
985     c Compute Cloud-Top Temperature and Pressure
986     c ------------------------------------------
987     cldtmp(i,1) = cldtmp(i,1) + cldnow*pkzgather(i,L)
988     . * ( tz(i,1,L) + dtmoist(i,1,L)*tmstp/pz(i,1) )
989     cldprs(i,1) = cldprs(i,1) + cldnow*plev
990     tmpimjm(i,1) = tmpimjm(i,1) + cldnow
991    
992     enddo
993     enddo
994    
995     c Compute Instantanious Total 2-D Cloud Fraction
996     c ----------------------------------------------
997     do j = 1,jm
998     do i = 1,im
999     totcld(i,j) = 1.0 - (1.-cldhi (i,j))
1000     . * (1.-cldmid(i,j))
1001     . * (1.-cldlow(i,j))
1002     enddo
1003     enddo
1004    
1005    
1006     C **********************************************************************
1007     C *** Fill Cloud Top Pressure and Temperature Diagnostic ***
1008     C **********************************************************************
1009    
1010     if(icldtmp.gt.0) then
1011     do j = 1,jm
1012     do i = 1,im
1013     if( cldtmp(i,j).gt.0.0 ) then
1014 molod 1.3 qdiag(i,j,icldtmp,bi,bj) = qdiag(i,j,icldtmp,bi,bj) +
1015 molod 1.1 . cldtmp(i,j)*totcld(i,j)/tmpimjm(i,j)
1016 molod 1.3 qdiag(i,j,icttcnt,bi,bj) = qdiag(i,j,icttcnt,bi,bj) +
1017     . totcld(i,j)
1018 molod 1.1 endif
1019     enddo
1020     enddo
1021     endif
1022    
1023     if(icldprs.gt.0) then
1024     do j = 1,jm
1025     do i = 1,im
1026     if( cldprs(i,j).gt.0.0 ) then
1027 molod 1.3 qdiag(i,j,icldprs,bi,bj) = qdiag(i,j,icldprs,bi,bj) +
1028 molod 1.1 . cldprs(i,j)*totcld(i,j)/tmpimjm(i,j)
1029 molod 1.3 qdiag(i,j,ictpcnt,bi,bj) = qdiag(i,j,ictpcnt,bi,bj) +
1030     . totcld(i,j)
1031 molod 1.1 endif
1032     enddo
1033     enddo
1034     endif
1035    
1036     C **********************************************************************
1037     C **** INCREMENT COUNTERS ****
1038     C **********************************************************************
1039    
1040     nlwlz = nlwlz + 1
1041     nswlz = nswlz + 1
1042    
1043     nlwcld = nlwcld + 1
1044     nswcld = nswcld + 1
1045    
1046     nmoistt = nmoistt + 1
1047     nmoistq = nmoistq + 1
1048     npreacc = npreacc + 1
1049     nprecon = nprecon + 1
1050    
1051     ncldmas = ncldmas + 1
1052     ndtrain = ndtrain + 1
1053    
1054     ndtls = ndtls + 1
1055     ndqls = ndqls + 1
1056    
1057     RETURN
1058     END
1059     SUBROUTINE RAS( NN, LEN, LENC, K, NLTOP, nlayr, DT
1060     *, UOI, ntracer, POI, QOI, PRS, PRJ, rnd, ncrnd
1061     *, RAINS, CLN, CLF, cldmas, detrain
1062     *, cp,grav,rkappa,alhl,rhfrac,rasmax )
1063     C
1064     C*********************************************************************
1065     C********************* SUBROUTINE RAS *****************************
1066     C********************** 16 MARCH 1988 ******************************
1067     C*********************************************************************
1068     C
1069 molod 1.8 implicit none
1070    
1071 molod 1.9 C Argument List
1072     integer nn,len,lenc,k,nltop,nlayr
1073 molod 1.1 integer ntracer
1074 molod 1.9 integer ncrnd
1075 molod 1.13 _RL dt
1076     _RL UOI(len,nlayr,ntracer), POI(len,K)
1077     _RL QOI(len,K), PRS(len,K+1), PRJ(len,K+1)
1078     _RL rnd(ncrnd)
1079     _RL RAINS(len,K), CLN(len,K), CLF(len,K)
1080     _RL cldmas(len,K), detrain(len,K)
1081     _RL cp,grav,rkappa,alhl,rhfrac(len),rasmax
1082 molod 1.9
1083     C Local Variables
1084 molod 1.13 _RL TCU(len,K), QCU(len,K)
1085     _RL ucu(len,K,ntracer)
1086     _RL ALF(len,K), BET(len,K), GAM(len,K)
1087 molod 1.1 *, ETA(len,K), HOI(len,K)
1088     *, PRH(len,K), PRI(len,K)
1089 molod 1.13 _RL HST(len,K), QOL(len,K), GMH(len,K)
1090 molod 1.1
1091 molod 1.13 _RL TX1(len), TX2(len), TX3(len), TX4(len), TX5(len)
1092 molod 1.1 *, TX6(len), TX7(len), TX8(len), TX9(len)
1093     *, TX11(len), TX12(len), TX13(len), TX14(len,ntracer)
1094 molod 1.9 *, TX15(len)
1095     *, WFN(len)
1096     integer IA1(len), IA2(len), IA3(len)
1097 molod 1.13 _RL cloudn(len), pcu(len)
1098 molod 1.1
1099 molod 1.8 integer krmin,icm
1100 molod 1.13 _RL rknob, cmb2pa
1101 molod 1.8 PARAMETER (KRMIN=01)
1102     PARAMETER (ICM=1000)
1103     PARAMETER (CMB2PA=100.0)
1104     PARAMETER (rknob = 10.)
1105 molod 1.9
1106     integer IC(ICM), IRND(icm)
1107 molod 1.13 _RL cmass(len,K)
1108 molod 1.9 LOGICAL SETRAS
1109    
1110     integer i,L,nc,ib,nt
1111 molod 1.8 integer km1,kp1,kprv,kcr,kfx,ncmx
1112 molod 1.13 _RL p00, crtmsf, frac, rasblf
1113 molod 1.8
1114     do L = 1,k
1115     do I = 1,LENC
1116     rains(i,l) = 0.
1117     enddo
1118     enddo
1119 molod 1.1
1120     p00 = 1000.
1121     crtmsf = 0.
1122    
1123     C The numerator here is the fraction of the subcloud layer mass flux
1124     C allowed to entrain into the cloud
1125    
1126     CCC FRAC = 1./dt
1127     FRAC = 0.5/dt
1128    
1129     KM1 = K - 1
1130     KP1 = K + 1
1131     C we want the ras adjustment time scale to be one hour (indep of dt)
1132     RASBLF = 1./3600.
1133     C
1134     KPRV = KM1
1135     C Removed KRMAX parameter
1136     KCR = MIN(KM1,nlayr-2)
1137     KFX = KM1 - KCR
1138     NCMX = KFX + NCRND
1139     C
1140     IF (KFX .GT. 0) THEN
1141     DO NC=1,KFX
1142     IC(NC) = K - NC
1143     ENDDO
1144     ENDIF
1145     C
1146     IF (NCRND .GT. 0) THEN
1147     DO I=1,ncrnd
1148     IRND(I) = (RND(I)-0.0005)*(KCR-KRMIN+1)
1149     IRND(I) = IRND(I) + KRMIN
1150     ENDDO
1151     C
1152     DO NC=1,NCRND
1153     IC(KFX+NC) = IRND(NC)
1154     ENDDO
1155     ENDIF
1156     C
1157     DO 100 NC=1,NCMX
1158     C
1159     IF (NC .EQ. 1 ) THEN
1160     SETRAS = .TRUE.
1161     ELSE
1162     SETRAS = .FALSE.
1163     ENDIF
1164     IB = IC(NC)
1165    
1166     c Initialize Cloud Fraction Array
1167     c -------------------------------
1168     do i = 1,lenc
1169     cloudn(i) = 0.0
1170     enddo
1171    
1172     CALL CLOUD(nn,LEN, LENC, K, NLTOP, nlayr, IB, RASBLF,SETRAS,FRAC
1173     *, CP, ALHL, RKAPPA, GRAV, P00, CRTMSF
1174     *, POI, QOI, UOI, Ntracer, PRS, PRJ
1175     *, PCU, CLOUDN, TCU, QCU, UCU, CMASS
1176     *, ALF, BET, GAM, PRH, PRI, HOI, ETA
1177     *, HST, QOL, GMH
1178     *, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9
1179     *, WFN, TX11, TX12, TX13, TX14, TX15
1180     *, IA1,IA2,IA3,rhfrac)
1181    
1182     C Compute fraction of grid box into which rain re-evap occurs (clf)
1183     c -----------------------------------------------------------------
1184     do i = 1,lenc
1185    
1186     c mass in detrainment layer
1187     c -------------------------
1188     tx1(i) = cmb2pa * (prs(i,ib+1) - prs(i,ib))/(grav*dt)
1189    
1190     c ratio of detraining cloud mass to mass in detrainment layer
1191     c -----------------------------------------------------------
1192     tx1(i) = rhfrac(i)*rknob * cmass(i,ib) / tx1(i)
1193     if(cmass(i,K).gt.0.) clf(i,ib) = clf(i,ib) + tx1(i)
1194     if( clf(i,ib).gt.1.) clf(i,ib) = 1.
1195     enddo
1196    
1197     c Compute Total Cloud Mass Flux
1198     c *****************************
1199     do L=ib,k
1200     do i=1,lenc
1201     cmass(i,L) = rhfrac(i)*cmass(i,L) * dt
1202     enddo
1203     enddo
1204    
1205     do L=ib,k
1206     do i=1,lenc
1207     cldmas(i,L) = cldmas(i,L) + cmass(i,L)
1208     enddo
1209     enddo
1210    
1211     do i=1,lenc
1212     detrain(i,ib) = detrain(i,ib) + cmass(i,ib)
1213     enddo
1214    
1215     DO L=IB,K
1216     DO I=1,LENC
1217     POI(I,L) = POI(I,L) + TCU(I,L) * DT * rhfrac(i)
1218     QOI(I,L) = QOI(I,L) + QCU(I,L) * DT * rhfrac(i)
1219     ENDDO
1220     ENDDO
1221     DO NT=1,Ntracer
1222     DO L=IB,K
1223     DO I=1,LENC
1224     UOI(I,L+nltop-1,NT)=UOI(I,L+nltop-1,NT)+UCU(I,L,NT)*DT*rhfrac(i)
1225     ENDDO
1226     ENDDO
1227     ENDDO
1228     DO I=1,LENC
1229     rains(I,ib) = rains(I,ib) + PCU(I)*dt * rhfrac(i)
1230     ENDDO
1231    
1232     100 CONTINUE
1233    
1234     c Fill Convective Cloud Fractions based on 3-D Rain Amounts
1235     c ---------------------------------------------------------
1236     do L=k-1,1,-1
1237     do i=1,lenc
1238     tx1(i) = 100*(prs(i,L+1)-prs(i,L))/grav
1239     cln(i,L) = min(1600*rains(i,L)/tx1(i),rasmax )
1240     enddo
1241     enddo
1242    
1243     RETURN
1244     END
1245     subroutine rndcloud (iras,nrnd,rnd,myid)
1246     implicit none
1247     integer n,iras,nrnd,myid
1248 molod 1.13 _RL random_numbx
1249     _RL rnd(nrnd)
1250 molod 1.1 integer irm
1251     parameter (irm = 1000)
1252 molod 1.13 _RL random(irm)
1253 molod 1.1 integer i,mcheck,numrand,iseed,index
1254     logical first
1255     data first /.true./
1256     integer iras0
1257     data iras0 /0/
1258     save random, iras0
1259    
1260     if(nrnd.eq.0.)then
1261     do i = 1,nrnd
1262     rnd(i) = 0
1263     enddo
1264 molod 1.15 if(first .and. myid.eq.1) print *,' NO RANDOM CLOUDS IN RAS '
1265 molod 1.1 go to 100
1266     endif
1267    
1268     mcheck = mod(iras-1,irm/nrnd)
1269    
1270     c First Time In From a Continuing RESTART (IRAS.GT.1) or Reading a New RESTART
1271     c ----------------------------------------------------------------------------
1272     if( first.and.(iras.gt.1) .or. iras.ne.iras0+1 )then
1273 molod 1.15 if( myid.eq.1 ) print *, 'Recreating Rand Numb Array in RNDCLOUD'
1274     if( myid.eq.1 ) print *, 'IRAS: ',iras,' IRAS0: ',iras0
1275 molod 1.1 numrand = mod(iras,irm/nrnd) * nrnd
1276     iseed = iras * nrnd - numrand
1277     call random_seedx(iseed)
1278     do i = 1,irm
1279 molod 1.12 random(i) = random_numbx(iseed)
1280 molod 1.1 enddo
1281     index = (iras-1)*nrnd
1282    
1283     c Multiple Time In But have Used Up all 1000 numbers (MCHECK.EQ.0)
1284     c ----------------------------------------------------------------
1285     else if (mcheck.eq.0) then
1286     iseed = (iras-1)*nrnd
1287     call random_seedx(iseed)
1288     do i = 1,irm
1289 molod 1.12 random(i) = random_numbx(iseed)
1290 molod 1.1 enddo
1291     index = iseed
1292    
1293     c Multiple Time In But have NOT Used Up all 1000 numbers (MCHECK.NE.0)
1294     c --------------------------------------------------------------------
1295     else
1296     index = (iras-1)*nrnd
1297     endif
1298    
1299     index = mod(index,irm)
1300     if( index+nrnd.gt.1000 ) index=1000-nrnd
1301    
1302     do n = 1,nrnd
1303     rnd(n) = random(index+n)
1304     enddo
1305    
1306     100 continue
1307     first = .false.
1308     iras0 = iras
1309     return
1310     end
1311 molod 1.12 function random_numbx(iseed)
1312 molod 1.1 implicit none
1313 molod 1.12 integer iseed
1314     real *8 seed,port_rand
1315 molod 1.13 _RL random_numbx
1316 molod 1.11 random_numbx = 0
1317 molod 1.10 #ifdef CRAY
1318 molod 1.13 _RL ranf
1319 molod 1.1 random_numbx = ranf()
1320 molod 1.12 #else
1321 molod 1.10 #ifdef SGI
1322 molod 1.13 _RL rand
1323 molod 1.1 random_numbx = rand()
1324     #endif
1325 molod 1.12 random_numbx = port_rand(seed)
1326     #endif
1327 molod 1.1 return
1328     end
1329     subroutine random_seedx (iseed)
1330     implicit none
1331     integer iseed
1332 molod 1.10 #ifdef CRAY
1333 molod 1.1 call ranset (iseed)
1334     #endif
1335 molod 1.10 #ifdef SGI
1336 molod 1.1 integer*4 seed
1337     seed = iseed
1338     call srand (seed)
1339     #endif
1340     return
1341     end
1342 molod 1.9 SUBROUTINE CLOUD(nn,LEN, LENC, K, NLTOP, nlayr, IC, RASALF
1343 molod 1.1 *, SETRAS, FRAC
1344     *, CP, ALHL, RKAP, GRAV, P00, CRTMSF
1345     *, POI, QOI, UOI, Ntracer, PRS, PRJ
1346     *, PCU, CLN, TCU, QCU, UCU, CMASS
1347     *, ALF, BET, GAM, PRH, PRI, HOL, ETA
1348     *, HST, QOL, GMH
1349     *, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, ALM
1350     *, WFN, AKM, QS1, CLF, UHT, WLQ
1351     *, IA, I1, I2,rhfrac)
1352     C
1353     C*********************************************************************
1354     C******************** Relaxed Arakawa-Schubert ***********************
1355     C********************* Plug Compatible Version **********************
1356     C************************ SUBROUTINE CLOUD ***************************
1357     C************************* 23 JULY 1992 ***************************
1358     C*********************************************************************
1359     C*********************************************************************
1360     C*********************************************************************
1361     C************************** Developed By *****************************
1362     C************************** *****************************
1363     C************************ Shrinivas Moorthi **************************
1364     C************************ and **************************
1365     C************************ Max J. Suarez *****************************
1366     C************************ *****************************
1367     C******************** Laboratory for Atmospheres *********************
1368     C****************** NASA/GSFC, Greenbelt, MD 20771 *******************
1369     C*********************************************************************
1370     C*********************************************************************
1371     C
1372     C The calculations of Moorthi and Suarez (1992, MWR) are
1373     C contained in the CLOUD routine.
1374     C It is probably advisable, at least initially, to treat CLOUD
1375     C as a black box that computes the single cloud adjustments. RAS,
1376     C on the other hand, can be tailored to each GCMs configuration
1377     C (ie, number and placement of levels, nature of boundary layer,
1378     C time step and frequency with which RAS is called).
1379     C
1380     C
1381     C Input:
1382     C ------
1383     C
1384     C LEN : The inner dimension of update and input arrays.
1385     C
1386     C LENC : The run: the number of soundings processes in a single call.
1387     C RAS works on the first LENC of the LEN soundings
1388     C passed. This allows working on pieces of the world
1389     C say for multitasking, without declaring temporary arrays
1390     C and copying the data to and from them. This is an f77
1391     C version. An F90 version would have to allow more
1392     C flexibility in the argument declarations. Obviously
1393     C (LENC<=LEN).
1394     C
1395     C K : Number of vertical layers (increasing downwards).
1396     C Need not be the same as the number of layers in the
1397     C GCM, since it is the outer dimension. The bottom layer
1398     C (K) is the subcloud layer.
1399     C
1400     C IC : Detrainment level to check for presence of convection
1401     C
1402     C RASALF : Relaxation parameter (< 1.) for present cloud-type
1403     C
1404     C SETRAS : Logical parameter to control re-calculation of
1405     C saturation specific humidity and mid level P**kappa
1406     C
1407     C FRAC : Fraction of the PBL (layer K) mass allowed to be used
1408     C by a cloud-type in time DT
1409     C
1410     C CP : Specific heat at constant pressure
1411     C
1412     C ALHL : Latent Heat of condensation
1413     C
1414     C RKAP : R/Cp, where R is the gas constant
1415     C
1416     C GRAV : Acceleration due to gravity
1417     C
1418     C P00 : A reference pressure in hPa, useually 1000 hPa
1419     C
1420     C CRTMSF : Critical value of mass flux above which cloudiness at
1421     C the detrainment layer of that cloud-type is assumed.
1422     C Affects only cloudiness calculation.
1423     C
1424     C POI : 2D array of dimension (LEN,K) containing potential
1425     C temperature. Updated but not initialized by RAS.
1426     C
1427     C QOI : 2D array of dimension (LEN,K) containing specific
1428     C humidity. Updated but not initialized by RAS.
1429     C
1430     C UOI : 3D array of dimension (LEN,K,NTRACER) containing tracers
1431     C Updated but not initialized by RAS.
1432     C
1433     C PRS : 2D array of dimension (LEN,K+1) containing pressure
1434     C in hPa at the interfaces of K-layers from top of the
1435     C atmosphere to the bottom. Not modified.
1436     C
1437     C PRJ : 2D array of dimension (LEN,K+1) containing (PRS/P00) **
1438     C RKAP. i.e. Exner function at layer edges. Not modified.
1439     C
1440     C rhfrac : 1D array of dimension (LEN) containing a rel.hum. scaling
1441     C fraction. Not modified.
1442     C
1443     C Output:
1444     C -------
1445     C
1446     C PCU : 1D array of length LEN containing accumulated
1447     C precipitation in mm/sec.
1448     C
1449     C CLN : 2D array of dimension (LEN,K) containing cloudiness
1450     C Note: CLN is bumped but NOT initialized
1451     C
1452     C TCU : 2D array of dimension (LEN,K) containing accumulated
1453     C convective heating (K/sec).
1454     C
1455     C QCU : 2D array of dimension (LEN,K) containing accumulated
1456     C convective drying (kg/kg/sec).
1457     C
1458     C CMASS : 2D array of dimension (LEN,K) containing the
1459     C cloud mass flux (kg/sec). Filled from cloud top
1460     C to base.
1461     C
1462     C Temporaries:
1463     C
1464     C ALF, BET, GAM, ETA, PRH, PRI, HOI, HST, QOL, GMH are temporary
1465     C 2D real arrays of dimension of at least (LENC,K) where LENC is
1466     C the horizontal dimension over which convection is invoked.
1467     C
1468     C
1469     C TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9, AKM, QS1, CLF, UHT
1470     C VHT, WLQ WFN are temporary real arrays of length at least LENC
1471     C
1472     C IA, I1, and I2 are temporary integer arrays of length LENC
1473     C
1474     C
1475     C************************************************************************
1476 molod 1.9 implicit none
1477     C Argument List declarations
1478     integer nn,LEN,LENC,K,NLTOP,nlayr,ic,ntracer
1479 molod 1.13 _RL rasalf
1480 molod 1.9 LOGICAL SETRAS
1481 molod 1.13 _RL frac, cp, alhl, rkap, grav, p00, crtmsf
1482     _RL POI(LEN,K),QOI(LEN,K),PRS(LEN,K+1),PRJ(LEN,K+1)
1483     _RL uoi(len,nlayr,ntracer)
1484     _RL PCU(LENC), CLN(LEN)
1485     _RL TCU(LEN,K), QCU(LEN,K), ucu(len,k,ntracer), CMASS(LEN,K)
1486     _RL ALF(LEN,K), BET(LEN,K), GAM(LEN,K), PRH(LEN,K), PRI(LEN,K)
1487     _RL HOL(LENC,K), ETA(LENC,K), HST(LENC,K), QOL(LENC,K)
1488     _RL GMH(LENC,K)
1489     _RL TX1(LENC), TX2(LENC), TX3(LENC), TX4(LENC)
1490     _RL TX5(LENC), TX6(LENC), TX7(LENC), TX8(LENC)
1491     _RL ALM(LENC), WFN(LENC), AKM(LENC), QS1(LENC)
1492     _RL WLQ(LENC), CLF(LENC)
1493     _RL uht(len,ntracer)
1494 molod 1.9 integer IA(LENC), I1(LENC),I2(LENC)
1495 molod 1.13 _RL rhfrac(len)
1496 molod 1.1
1497 molod 1.9 C Local Variables
1498 molod 1.13 _RL daylen,half,one,zero,cmb2pa,rhmax
1499 molod 1.1 PARAMETER (DAYLEN=86400.0, HALF=0.5, ONE=1.0, ZERO=0.0)
1500     PARAMETER (CMB2PA=100.0)
1501     PARAMETER (RHMAX=0.9999)
1502 molod 1.13 _RL rkapp1,onebcp,albcp,onebg,cpbg,twobal
1503 molod 1.1 C
1504 molod 1.9 integer nt,km1,ic1,i,L,len1,len2,isav,len11,ii
1505     integer lena,lena1,lenb,tem,tem1
1506 molod 1.1
1507     c Explicit Inline Directives
1508     c --------------------------
1509 molod 1.10 #ifdef CRAY
1510     #ifdef f77
1511 molod 1.1 cfpp$ expand (qsat)
1512     #endif
1513     #endif
1514    
1515     RKAPP1 = 1.0 + RKAP
1516     ONEBCP = 1.0 / CP
1517     ALBCP = ALHL * ONEBCP
1518     ONEBG = 1.0 / GRAV
1519     CPBG = CP * ONEBG
1520     TWOBAL = 2.0 / ALHL
1521     C
1522     KM1 = K - 1
1523     IC1 = IC + 1
1524     C
1525 molod 1.9 C SETTING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T.
1526 molod 1.1 C
1527    
1528     IF (SETRAS) THEN
1529    
1530     DO 2050 L=1,K
1531     DO 2030 I=1,LENC
1532     PRH(I,L) = (PRJ(I,L+1)*PRS(I,L+1) - PRJ(I,L)*PRS(I,L))
1533     * / ((PRS(I,L+1)-PRS(I,L)) * RKAPP1)
1534     2030 CONTINUE
1535     2050 CONTINUE
1536    
1537     DO 2070 L=1,K
1538     DO 2060 I=1,LENC
1539     TX5(I) = POI(I,L) * PRH(I,L)
1540     TX1(I) = (PRS(I,L) + PRS(I,L+1)) * 0.5
1541     TX3(I) = TX5(I)
1542     CALL QSAT(TX3(I), TX1(I), TX2(I), TX4(I), .TRUE.)
1543     ALF(I,L) = TX2(I) - TX4(I) * TX5(I)
1544     BET(I,L) = TX4(I) * PRH(I,L)
1545     GAM(I,L) = 1.0 / ((1.0 + TX4(I)*ALBCP) * PRH(I,L))
1546     PRI(I,L) = (CP/CMB2PA) / (PRS(I,L+1) - PRS(I,L))
1547     2060 CONTINUE
1548     2070 CONTINUE
1549    
1550     ENDIF
1551     C
1552     C
1553     DO 10 L=1,K
1554     DO 10 I=1,LEN
1555     TCU(I,L) = 0.0
1556     QCU(I,L) = 0.0
1557     CMASS(I,L) = 0.0
1558     10 CONTINUE
1559    
1560     do nt = 1,ntracer
1561     do L=1,K
1562     do I=1,LENC
1563     ucu(I,L,nt) = 0.0
1564     enddo
1565     enddo
1566     enddo
1567     C
1568     DO 30 I=1,LENC
1569     TX1(I) = PRJ(I,K+1) * POI(I,K)
1570     QS1(I) = ALF(I,K) + BET(I,K)*POI(I,K)
1571     QOL(I,K) = MIN(QS1(I)*RHMAX,QOI(I,K))
1572    
1573     HOL(I,K) = TX1(I)*CP + QOL(I,K)*ALHL
1574     ETA(I,K) = ZERO
1575     TX2(I) = (PRJ(I,K+1) - PRJ(I,K)) * POI(I,K) * CP
1576     30 CONTINUE
1577     C
1578     IF (IC .LT. KM1) THEN
1579     DO 3703 L=KM1,IC1,-1
1580     DO 50 I=1,LENC
1581     QS1(I) = ALF(I,L) + BET(I,L)*POI(I,L)
1582     QOL(I,L) = MIN(QS1(I)*RHMAX,QOI(I,L))
1583     C
1584     TEM1 = TX2(I) + PRJ(I,L+1) * POI(I,L) * CP
1585     HOL(I,L) = TEM1 + QOL(I,L )* ALHL
1586     HST(I,L) = TEM1 + QS1(I) * ALHL
1587    
1588     TX1(I) = (PRJ(I,L+1) - PRJ(I,L)) * POI(I,L)
1589     ETA(I,L) = ETA(I,L+1) + TX1(I)*CPBG
1590     TX2(I) = TX2(I) + TX1(I)*CP
1591     50 CONTINUE
1592     C
1593     3703 CONTINUE
1594     ENDIF
1595    
1596    
1597     DO 70 I=1,LENC
1598     HOL(I,IC) = TX2(I)
1599     QS1(I) = ALF(I,IC) + BET(I,IC)*POI(I,IC)
1600     QOL(I,IC) = MIN(QS1(I)*RHMAX,QOI(I,IC))
1601     c
1602     TEM1 = TX2(I) + PRJ(I,IC1) * POI(I,IC) * CP
1603     HOL(I,IC) = TEM1 + QOL(I,IC) * ALHL
1604     HST(I,IC) = TEM1 + QS1(I) * ALHL
1605     C
1606     TX3(I ) = (PRJ(I,IC1) - PRH(I,IC)) * POI(I,IC)
1607     ETA(I,IC) = ETA(I,IC1) + CPBG * TX3(I)
1608     70 CONTINUE
1609     C
1610     DO 130 I=1,LENC
1611     TX2(I) = HOL(I,K) - HST(I,IC)
1612     TX1(I) = ZERO
1613    
1614     130 CONTINUE
1615     C
1616     C ENTRAINMENT PARAMETER
1617     C
1618     DO 160 L=IC,KM1
1619     DO 160 I=1,LENC
1620     TX1(I) = TX1(I) + (HST(I,IC) - HOL(I,L)) * (ETA(I,L) - ETA(I,L+1))
1621     160 CONTINUE
1622     C
1623     LEN1 = 0
1624     LEN2 = 0
1625     ISAV = 0
1626     DO 195 I=1,LENC
1627     IF (TX1(I) .GT. ZERO .AND. TX2(I) .GT. ZERO
1628     . .AND. rhfrac(i).ne.0.0 ) THEN
1629     LEN1 = LEN1 + 1
1630     IA(LEN1) = I
1631     ALM(LEN1) = TX2(I) / TX1(I)
1632     ENDIF
1633     195 CONTINUE
1634     C
1635     LEN2 = LEN1
1636     if (IC1 .lt. K) then
1637     DO 196 I=1,LENC
1638     IF (TX2(I) .LE. 0.0 .AND. (HOL(I,K) .GT. HST(I,IC1))
1639     . .AND. rhfrac(i).ne.0.0 ) THEN
1640     LEN2 = LEN2 + 1
1641     IA(LEN2) = I
1642     ALM(LEN2) = 0.0
1643     ENDIF
1644     196 CONTINUE
1645     endif
1646     C
1647     IF (LEN2 .EQ. 0) THEN
1648     DO 5010 I=1,LENC*K
1649     HST(I,1) = 0.0
1650     QOL(I,1) = 0.0
1651     5010 CONTINUE
1652     DO 5020 I=1,LENC
1653     PCU(I) = 0.0
1654     5020 CONTINUE
1655     RETURN
1656     ENDIF
1657     LEN11 = LEN1 + 1
1658     C
1659     C NORMALIZED MASSFLUX
1660     C
1661     DO 250 I=1,LEN2
1662     ETA(I,K) = 1.0
1663     II = IA(I)
1664     TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
1665     TX4(I) = PRS(II,K)
1666     250 CONTINUE
1667     C
1668     DO 252 I=LEN11,LEN2
1669     WFN(I) = 0.0
1670     II = IA(I)
1671     IF (HST(II,IC1) .LT. HST(II,IC)) THEN
1672     TX6(I) = (HST(II,IC1)-HOL(II,K))/(HST(II,IC1)-HST(II,IC))
1673     ELSE
1674     TX6(I) = 0.0
1675     ENDIF
1676     TX2(I) = 0.5 * (PRS(II,IC1)+PRS(II,IC1+1)) * (1.0-TX6(I))
1677     * + TX2(I) * TX6(I)
1678     252 CONTINUE
1679     C
1680     CALL ACRITN(LEN2, TX2, TX4, TX3)
1681     C
1682     DO 260 L=KM1,IC,-1
1683     DO 255 I=1,LEN2
1684     TX1(I) = ETA(IA(I),L)
1685     255 CONTINUE
1686     DO 260 I=1,LEN2
1687     ETA(I,L) = 1.0 + ALM(I) * TX1(I)
1688     260 CONTINUE
1689     C
1690     C CLOUD WORKFUNCTION
1691     C
1692     IF (LEN1 .GT. 0) THEN
1693     DO 270 I=1,LEN1
1694     II = IA(I)
1695     WFN(I) = - GAM(II,IC) * (PRJ(II,IC1) - PRH(II,IC))
1696     * * HST(II,IC) * ETA(I,IC1)
1697     270 CONTINUE
1698     ENDIF
1699     C
1700     DO 290 I=1,LEN2
1701     II = IA(I)
1702     TX1(I) = HOL(II,K)
1703     290 CONTINUE
1704     C
1705     IF (IC1 .LE. KM1) THEN
1706    
1707     DO 380 L=KM1,IC1,-1
1708     DO 380 I=1,LEN2
1709     II = IA(I)
1710     TEM = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * HOL(II,L)
1711     C
1712     PCU(I) = PRJ(II,L+1) - PRH(II,L)
1713     TEM1 = ETA(I,L+1) * PCU(I)
1714     TX1(I) = TX1(I)*PCU(I)
1715     C
1716     PCU(I) = PRH(II,L) - PRJ(II,L)
1717     TEM1 = (TEM1 + ETA(I,L) * PCU(I)) * HST(II,L)
1718     TX1(I) = TX1(I) + TEM*PCU(I)
1719     C
1720     WFN(I) = WFN(I) + (TX1(I) - TEM1) * GAM(II,L)
1721     TX1(I) = TEM
1722     380 CONTINUE
1723     ENDIF
1724     C
1725     LENA = 0
1726     IF (LEN1 .GT. 0) THEN
1727     DO 512 I=1,LEN1
1728     II = IA(I)
1729     WFN(I) = WFN(I) + TX1(I) * GAM(II,IC)*(PRJ(II,IC1)-PRH(II,IC))
1730     * - TX3(I)
1731     IF (WFN(I) .GT. 0.0) THEN
1732     LENA = LENA + 1
1733     I1(LENA) = IA(I)
1734     I2(LENA) = I
1735     TX1(LENA) = WFN(I)
1736     TX2(LENA) = QS1(IA(I))
1737     TX6(LENA) = 1.0
1738     ENDIF
1739     512 CONTINUE
1740     ENDIF
1741     LENB = LENA
1742     DO 515 I=LEN11,LEN2
1743     WFN(I) = WFN(I) - TX3(I)
1744     IF (WFN(I) .GT. 0.0 .AND. TX6(I) .GT. 0.0) THEN
1745     LENB = LENB + 1
1746     I1(LENB) = IA(I)
1747     I2(LENB) = I
1748     TX1(LENB) = WFN(I)
1749     TX2(LENB) = QS1(IA(I))
1750     TX4(LENB) = TX6(I)
1751     ENDIF
1752     515 CONTINUE
1753     C
1754     IF (LENB .LE. 0) THEN
1755     DO 5030 I=1,LENC*K
1756     HST(I,1) = 0.0
1757     QOL(I,1) = 0.0
1758     5030 CONTINUE
1759     DO 5040 I=1,LENC
1760     PCU(I) = 0.0
1761     5040 CONTINUE
1762     RETURN
1763     ENDIF
1764    
1765     C
1766     DO 516 I=1,LENB
1767     WFN(I) = TX1(I)
1768     QS1(I) = TX2(I)
1769     516 CONTINUE
1770     C
1771     DO 520 L=IC,K
1772     DO 517 I=1,LENB
1773     TX1(I) = ETA(I2(I),L)
1774     517 CONTINUE
1775     DO 520 I=1,LENB
1776     ETA(I,L) = TX1(I)
1777     520 CONTINUE
1778     C
1779     LENA1 = LENA + 1
1780     C
1781     DO 510 I=1,LENA
1782     II = I1(I)
1783     TX8(I) = HST(II,IC) - HOL(II,IC)
1784     510 CONTINUE
1785     DO 530 I=LENA1,LENB
1786     II = I1(I)
1787     TX6(I) = TX4(I)
1788     TEM = TX6(I) * (HOL(II,IC)-HOL(II,IC1)) + HOL(II,IC1)
1789     TX8(I) = HOL(II,K) - TEM
1790    
1791     TEM1 = TX6(I) * (QOL(II,IC)-QOL(II,IC1)) + QOL(II,IC1)
1792     TX5(I) = TEM - TEM1 * ALHL
1793     QS1(I) = TEM1 + TX8(I)*(ONE/ALHL)
1794     TX3(I) = HOL(II,IC)
1795     530 CONTINUE
1796     C
1797     C
1798     DO 620 I=1,LENB
1799     II = I1(I)
1800     WLQ(I) = QOL(II,K) - QS1(I) * ETA(I,IC)
1801     TX7(I) = HOL(II,K)
1802     620 CONTINUE
1803     DO NT=1,Ntracer
1804     DO 621 I=1,LENB
1805     II = I1(I)
1806     UHT(I,NT) = UOI(II,K+nltop-1,NT)-UOI(II,IC+nltop-1,NT) * ETA(I,IC)
1807     621 CONTINUE
1808     ENDDO
1809     C
1810     DO 635 L=KM1,IC,-1
1811     DO 630 I=1,LENB
1812     II = I1(I)
1813     TEM = ETA(I,L) - ETA(I,L+1)
1814     WLQ(I) = WLQ(I) + TEM * QOL(II,L)
1815     630 CONTINUE
1816     635 CONTINUE
1817     DO NT=1,Ntracer
1818     DO L=KM1,IC,-1
1819     DO I=1,LENB
1820     II = I1(I)
1821     TEM = ETA(I,L) - ETA(I,L+1)
1822     UHT(I,NT) = UHT(I,NT) + TEM * UOI(II,L+nltop-1,NT)
1823     ENDDO
1824     ENDDO
1825     ENDDO
1826     C
1827     C CALCULATE GS AND PART OF AKM (THAT REQUIRES ETA)
1828     C
1829     DO 690 I=1,LENB
1830     II = I1(I)
1831     c TX7(I) = HOL(II,K)
1832     TEM = (POI(II,KM1) - POI(II,K)) / (PRH(II,K) - PRH(II,KM1))
1833     HOL(I,K) = TEM * (PRJ(II,K)-PRH(II,KM1))*PRH(II,K)*PRI(II,K)
1834     HOL(I,KM1) = TEM * (PRH(II,K)-PRJ(II,K))*PRH(II,KM1)*PRI(II,KM1)
1835     AKM(I) = ZERO
1836     TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
1837     690 CONTINUE
1838    
1839     IF (IC1 .LE. KM1) THEN
1840     DO 750 L=KM1,IC1,-1
1841     DO 750 I=1,LENB
1842     II = I1(I)
1843     TEM = (POI(II,L-1) - POI(II,L)) * ETA(I,L)
1844     * / (PRH(II,L) - PRH(II,L-1))
1845     C
1846     HOL(I,L) = TEM * (PRJ(II,L)-PRH(II,L-1)) * PRH(II,L)
1847     * * PRI(II,L) + HOL(I,L)
1848     HOL(I,L-1) = TEM * (PRH(II,L)-PRJ(II,L)) * PRH(II,L-1)
1849     * * PRI(II,L-1)
1850     C
1851     AKM(I) = AKM(I) - HOL(I,L)
1852     * * (ETA(I,L) * (PRH(II,L)-PRJ(II,L)) +
1853     * ETA(I,L+1) * (PRJ(II,L+1)-PRH(II,L))) / PRH(II,L)
1854     750 CONTINUE
1855     ENDIF
1856     C
1857     C
1858     CALL RNCL(LENB, TX2, TX1, CLF)
1859    
1860     DO 770 I=1,LENB
1861     TX2(I) = (ONE - TX1(I)) * WLQ(I)
1862     WLQ(I) = TX1(I) * WLQ(I)
1863     C
1864     TX1(I) = HOL(I,IC)
1865     770 CONTINUE
1866     DO 790 I=LENA1, LENB
1867     II = I1(I)
1868     TX1(I) = TX1(I) + (TX5(I)-TX3(I)+QOL(II,IC)*ALHL)*(PRI(II,IC)/CP)
1869     790 CONTINUE
1870    
1871     DO 800 I=1,LENB
1872     HOL(I,IC) = TX1(I) - TX2(I) * ALBCP * PRI(I1(I),IC)
1873     800 CONTINUE
1874    
1875     IF (LENA .GT. 0) THEN
1876     DO 810 I=1,LENA
1877     II = I1(I)
1878     AKM(I) = AKM(I) - ETA(I,IC1) * (PRJ(II,IC1) - PRH(II,IC))
1879     * * TX1(I) / PRH(II,IC)
1880     810 CONTINUE
1881     ENDIF
1882     c
1883     C CALCULATE GH
1884     C
1885     DO 830 I=1,LENB
1886     II = I1(I)
1887     TX3(I) = QOL(II,KM1) - QOL(II,K)
1888     GMH(I,K) = HOL(I,K) + TX3(I) * PRI(II,K) * (ALBCP)
1889    
1890     AKM(I) = AKM(I) + GAM(II,KM1)*(PRJ(II,K)-PRH(II,KM1))
1891     * * GMH(I,K)
1892     TX3(I) = zero
1893     830 CONTINUE
1894     C
1895     IF (IC1 .LE. KM1) THEN
1896     DO 840 L=KM1,IC1,-1
1897     DO 840 I=1,LENB
1898     II = I1(I)
1899     TX2(I) = TX3(I)
1900     TX3(I) = (QOL(II,L-1) - QOL(II,L)) * ETA(I,L)
1901     TX2(I) = TX2(I) + TX3(I)
1902     C
1903     GMH(I,L) = HOL(I,L) + TX2(I) * PRI(II,L) * (ALBCP*HALF)
1904     840 CONTINUE
1905     C
1906     C
1907     ENDIF
1908     DO 850 I=LENA1,LENB
1909     TX3(I) = TX3(I) + TWOBAL
1910     * * (TX7(I) - TX8(I) - TX5(I) - QOL(I1(I),IC)*ALHL)
1911     850 CONTINUE
1912     DO 860 I=1,LENB
1913     GMH(I,IC) = TX1(I) + PRI(I1(I),IC) * ONEBCP
1914     * * (TX3(I)*(ALHL*HALF) + ETA(I,IC) * TX8(I))
1915     860 CONTINUE
1916     C
1917     C CALCULATE HC PART OF AKM
1918     C
1919     IF (IC1 .LE. KM1) THEN
1920     DO 870 I=1,LENB
1921     TX1(I) = GMH(I,K)
1922     870 CONTINUE
1923     DO 3725 L=KM1,IC1,-1
1924     DO 880 I=1,LENB
1925     II = I1(I)
1926     TX1(I) = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * GMH(I,L)
1927     TX2(I) = GAM(II,L-1) * (PRJ(II,L) - PRH(II,L-1))
1928     880 CONTINUE
1929     C
1930     IF (L .EQ. IC1) THEN
1931     DO 890 I=LENA1,LENB
1932     TX2(I) = ZERO
1933     890 CONTINUE
1934     ENDIF
1935     DO 900 I=1,LENB
1936     II = I1(I)
1937     AKM(I) = AKM(I) + TX1(I) *
1938     * (TX2(I) + GAM(II,L)*(PRH(II,L)-PRJ(II,L)))
1939     900 CONTINUE
1940     3725 CONTINUE
1941     ENDIF
1942     C
1943     DO 920 I=LENA1,LENB
1944     II = I1(I)
1945     TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
1946     * + 0.5*(PRS(II,IC+2) - PRS(II,IC)) * (ONE-TX6(I))
1947     c
1948     TX1(I) = PRS(II,IC1)
1949     TX5(I) = 0.5 * (PRS(II,IC1) + PRS(II,IC+2))
1950     C
1951     IF ((TX2(I) .GE. TX1(I)) .AND. (TX2(I) .LT. TX5(I))) THEN
1952     TX6(I) = ONE - (TX2(I) - TX1(I)) / (TX5(I) - TX1(I))
1953     C
1954     TEM = PRI(II,IC1) / PRI(II,IC)
1955     HOL(I,IC1) = HOL(I,IC1) + HOL(I,IC) * TEM
1956     HOL(I,IC) = ZERO
1957     C
1958     GMH(I,IC1) = GMH(I,IC1) + GMH(I,IC) * TEM
1959     GMH(I,IC) = ZERO
1960     ELSEIF (TX2(I) .LT. TX1(I)) THEN
1961     TX6(I) = 1.0
1962     ELSE
1963     TX6(I) = 0.0
1964     ENDIF
1965     920 CONTINUE
1966     C
1967     C
1968     DO I=1,LENC
1969     PCU(I) = 0.0
1970     ENDDO
1971    
1972     DO 970 I=1,LENB
1973     II = I1(I)
1974     IF (AKM(I) .LT. ZERO .AND. WLQ(I) .GE. 0.0) THEN
1975     WFN(I) = - TX6(I) * WFN(I) * RASALF / AKM(I)
1976     ELSE
1977     WFN(I) = ZERO
1978     ENDIF
1979     TEM = (PRS(II,K+1)-PRS(II,K))*(CMB2PA*FRAC)
1980     WFN(I) = MIN(WFN(I), TEM)
1981     C
1982     C compute cloud amount
1983     C
1984     CC TX1(I) = CLN(II)
1985     CC IF (WFN(I) .GT. CRTMSF) TX1(I) = TX1(I) + CLF(I)
1986     CC IF (TX1(I) .GT. ONE) TX1(I) = ONE
1987     C
1988     C PRECIPITATION
1989     C
1990     PCU(II) = WLQ(I) * WFN(I) * ONEBG
1991     C
1992     C CUMULUS FRICTION AT THE BOTTOM LAYER
1993     C
1994     TX4(I) = WFN(I) * (1.0/ALHL)
1995     TX5(I) = WFN(I) * ONEBCP
1996     970 CONTINUE
1997     C
1998     C compute cloud mass flux for diagnostic output
1999     C
2000     DO L = IC,K
2001     DO I=1,LENB
2002     II = I1(I)
2003     if(L.lt.K)then
2004     CMASS(II,L) = ETA(I,L+1) * WFN(I) * ONEBG
2005     else
2006     CMASS(II,L) = WFN(I) * ONEBG
2007     endif
2008     ENDDO
2009     ENDDO
2010     C
2011     CC DO 975 I=1,LENB
2012     CC II = I1(I)
2013     CC CLN(II) = TX1(I)
2014     CC975 CONTINUE
2015     C
2016     C THETA AND Q CHANGE DUE TO CLOUD TYPE IC
2017     C
2018    
2019     c TEMA = 0.0
2020     c TEMB = 0.0
2021     DO 990 L=IC,K
2022     DO 980 I=1,LENB
2023     II = I1(I)
2024     TEM = (GMH(I,L) - HOL(I,L)) * TX4(I)
2025     TEM1 = HOL(I,L) * TX5(I)
2026     C
2027     TCU(II,L) = TEM1 / PRH(II,L)
2028     QCU(II,L) = TEM
2029     980 CONTINUE
2030    
2031     c I = I1(IP1)
2032     c
2033     c TEM = (PRS(I,L+1)-PRS(I,L)) * (ONEBG*100.0)
2034     c TEMA = TEMA + TCU(I,L) * PRH(I,L) * TEM * (CP/ALHL)
2035     c TEMB = TEMB + QCU(I,L) * TEM
2036     C
2037     990 CONTINUE
2038     C
2039     c Compute Tracer Tendencies
2040     c -------------------------
2041     do nt = 1,ntracer
2042    
2043     c Tracer Tendency at the Bottom Layer
2044     c -----------------------------------
2045     DO 995 I=1,LENB
2046     II = I1(I)
2047     TEM = half*TX5(I) * PRI(II,K)
2048     TX1(I) = (UOI(II,KM1+nltop-1,nt) - UOI(II,K+nltop-1,nt))
2049     ucu(II,K,nt) = TEM * TX1(I)
2050     995 CONTINUE
2051    
2052     c Tracer Tendency at all other Levels
2053     c -----------------------------------
2054     DO 1020 L=KM1,IC1,-1
2055     DO 1010 I=1,LENB
2056     II = I1(I)
2057     TEM = half*TX5(I) * PRI(II,L)
2058     TEM1 = TX1(I)
2059     TX1(I) = (UOI(II,L-1+nltop-1,nt)-UOI(II,L+nltop-1,nt)) * ETA(I,L)
2060     TX3(I) = (TX1(I) + TEM1) * TEM
2061     1010 CONTINUE
2062     DO 1020 I=1,LENB
2063     II = I1(I)
2064     ucu(II,L,nt) = TX3(I)
2065     1020 CONTINUE
2066    
2067     DO 1030 I=1,LENB
2068     II = I1(I)
2069     IF (TX6(I) .GE. 1.0) THEN
2070     TEM = half*TX5(I) * PRI(II,IC)
2071     ELSE
2072     TEM = 0.0
2073     ENDIF
2074     TX1(I) = (TX1(I) + UHT(I,nt) + UHT(I,nt)) * TEM
2075     1030 CONTINUE
2076     DO 1040 I=1,LENB
2077     II = I1(I)
2078     ucu(II,IC,nt) = TX1(I)
2079     1040 CONTINUE
2080    
2081     enddo
2082     C
2083     C PENETRATIVE CONVECTION CALCULATION OVER
2084     C
2085    
2086     RETURN
2087     END
2088     SUBROUTINE RNCL(LEN, PL, RNO, CLF)
2089     C
2090     C*********************************************************************
2091     C********************** Relaxed Arakawa-Schubert *********************
2092     C************************ SUBROUTINE RNCL ************************
2093     C**************************** 23 July 1992 ***************************
2094     C*********************************************************************
2095 molod 1.9 implicit none
2096     C Argument List declarations
2097     integer len
2098 molod 1.13 _RL PL(LEN), RNO(LEN), CLF(LEN)
2099 molod 1.1
2100 molod 1.9 C Local Variables
2101 molod 1.13 _RL p5,p8,pt8,pt2,pfac,p4,p6,p7,p9,cucld,cfac
2102 molod 1.1 PARAMETER (P5=500.0, P8=800.0, PT8=0.8, PT2=0.2)
2103     PARAMETER (PFAC=PT2/(P8-P5))
2104     PARAMETER (P4=400.0, P6=401.0)
2105     PARAMETER (P7=700.0, P9=900.0)
2106     PARAMETER (CUCLD=0.5,CFAC=CUCLD/(P6-P4))
2107 molod 1.9
2108     integer i
2109 molod 1.1 C
2110     DO 10 I=1,LEN
2111     rno(i) = 1.0
2112     ccc if( pl(i).le.400.0 ) rno(i) = max( 0.75, 1.0-0.0025*(400.0-pl(i)) )
2113    
2114     ccc IF ( PL(I).GE.P7 .AND. PL(I).LE.P9 ) THEN
2115     ccc RNO(I) = ((P9-PL(I))/(P9-P7)) **2
2116     ccc ELSE IF (PL(I).GT.P9) THEN
2117     ccc RNO(I) = 0.
2118     ccc ENDIF
2119    
2120     CLF(I) = CUCLD
2121     C
2122     CARIESIF (PL(I) .GE. P5 .AND. PL(I) .LE. P8) THEN
2123     CARIES RNO(I) = (P8-PL(I))*PFAC + PT8
2124     CARIESELSEIF (PL(I) .GT. P8 ) THEN
2125     CARIES RNO(I) = PT8
2126     CARIESENDIF
2127     CARIES
2128     IF (PL(I) .GE. P4 .AND. PL(I) .LE. P6) THEN
2129     CLF(I) = (P6-PL(I))*CFAC
2130     ELSEIF (PL(I) .GT. P6 ) THEN
2131     CLF(I) = 0.0
2132     ENDIF
2133     10 CONTINUE
2134     C
2135     RETURN
2136     END
2137     SUBROUTINE ACRITN ( LEN,PL,PLB,ACR )
2138    
2139     C*********************************************************************
2140     C********************** Relaxed Arakawa-Schubert *********************
2141     C************************** SUBROUTINE ACRIT *********************
2142     C****************** modified August 28, 1996 L.Takacs ************
2143     C**** *****
2144     C**** Note: Data obtained from January Mean After-Analysis *****
2145     C**** from 4x5 46-layer GEOS Assimilation *****
2146     C**** *****
2147     C*********************************************************************
2148 molod 1.9 implicit none
2149     C Argument List declarations
2150     integer len
2151 molod 1.13 _RL PL(LEN), PLB(LEN), ACR(LEN)
2152 molod 1.1
2153 molod 1.9 C Local variables
2154     integer lma
2155 molod 1.1 parameter (lma=18)
2156 molod 1.13 _RL p(lma)
2157     _RL a(lma)
2158 molod 1.9 integer i,L
2159 molod 1.13 _RL temp
2160 molod 1.1
2161     data p / 93.81, 111.65, 133.46, 157.80, 186.51,
2162     . 219.88, 257.40, 301.21, 352.49, 409.76,
2163     . 471.59, 535.04, 603.33, 672.79, 741.12,
2164     . 812.52, 875.31, 930.20/
2165    
2166     data a / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364,
2167     . 1.01971, .65846, .45867, .38687, .31002,
2168     . .25574, .20347, .17254, .15260, .16756,
2169     . .09916, .10360, .05880/
2170    
2171    
2172     do L=1,lma-1
2173     do i=1,len
2174     if( pl(i).ge.p(L) .and.
2175     . pl(i).le.p(L+1)) then
2176     temp = ( pl(i)-p(L) )/( p(L+1)-p(L) )
2177     acr(i) = a(L+1)*temp + a(L)*(1-temp)
2178     endif
2179     enddo
2180     enddo
2181    
2182     do i=1,len
2183     if( pl(i).lt.p(1) ) acr(i) = a(1)
2184     if( pl(i).gt.p(lma) ) acr(i) = a(lma)
2185     enddo
2186    
2187     do i=1,len
2188     acr(i) = acr(i) * (plb(i)-pl(i))
2189     enddo
2190    
2191     RETURN
2192     END
2193 molod 1.6 SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DP,PLKE,
2194 molod 1.1 1 PLK,TH,TEMP1,TEMP2,TEMP3,ITMP1,ITMP2,RCON,RLAR,CLSBTH,tmscl,
2195     2 tmfrc,cp,gravity,alhl,gamfac,cldlz,RHCRIT,offset,alpha)
2196    
2197 molod 1.9 implicit none
2198     C Argument List declarations
2199     integer nn,irun,nlay
2200 molod 1.13 _RL TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY),
2201 molod 1.9 . PL(IRUN,NLAY),CLFRAC(IRUN,NLAY),SP(IRUN),TEMP1(IRUN,NLAY),
2202     . TEMP2(IRUN,NLAY),PLKE(IRUN,NLAY+1),
2203     . RCON(IRUN),RLAR(IRUN),DP(IRUN,NLAY),PLK(IRUN,NLAY),TH(IRUN,NLAY),
2204     . TEMP3(IRUN,NLAY)
2205     integer ITMP1(IRUN,NLAY),ITMP2(IRUN,NLAY)
2206 molod 1.13 _RL CLSBTH(IRUN,NLAY)
2207     _RL tmscl,tmfrc,cp,gravity,alhl,gamfac,offset,alpha
2208     _RL cldlz(irun,nlay)
2209     _RL rhcrit(irun,nlay)
2210 molod 1.9 C
2211     C Local Variables
2212 molod 1.13 _RL zm1p04,zero,two89,zp44,zp01,half,zp578,one,thousand,z3600
2213     _RL zp1,zp001
2214 molod 1.1 PARAMETER (ZM1P04 = -1.04E-4 )
2215     PARAMETER (ZERO = 0.)
2216     PARAMETER (TWO89= 2.89E-5)
2217     PARAMETER ( ZP44= 0.44)
2218     PARAMETER ( ZP01= 0.01)
2219     PARAMETER ( ZP1 = 0.1 )
2220     PARAMETER ( ZP001= 0.001)
2221     PARAMETER ( HALF= 0.5)
2222     PARAMETER ( ZP578 = 0.578 )
2223     PARAMETER ( ONE = 1.)
2224     PARAMETER ( THOUSAND = 1000.)
2225     PARAMETER ( Z3600 = 3600.)
2226     C
2227 molod 1.13 _RL EVP9(IRUN,NLAY)
2228     _RL water(irun),crystal(irun)
2229     _RL watevap(irun),iceevap(irun)
2230     _RL fracwat,fracice, tice,rh,fact,dum
2231     _RL rainmax(irun)
2232     _RL getcon,rphf,elocp,cpog,relax
2233     _RL exparg,arearat,rpow
2234 molod 1.9
2235     integer i,L,n,nlaym1,irnlay,irnlm1
2236 molod 1.1
2237     c Explicit Inline Directives
2238     c --------------------------
2239 molod 1.10 #ifdef CRAY
2240     #ifdef f77
2241 molod 1.1 cfpp$ expand (qsat)
2242     #endif
2243     #endif
2244    
2245     tice = getcon('FREEZING-POINT')
2246    
2247     fracwat = 0.70
2248     fracice = 0.01
2249    
2250     NLAYM1 = NLAY - 1
2251     IRNLAY = IRUN*NLAY
2252     IRNLM1 = IRUN*(NLAY-1)
2253    
2254     RPHF = Z3600/tmscl
2255    
2256     ELOCP = alhl/cp
2257     CPOG = cp/gravity
2258    
2259     DO I = 1,IRUN
2260     RLAR(I) = 0.
2261     water(i) = 0.
2262     crystal(i) = 0.
2263     ENDDO
2264    
2265     do L = 1,nlay
2266     do i = 1,irun
2267     EVP9(i,L) = 0.
2268     TEMP1(i,L) = 0.
2269     TEMP2(i,L) = 0.
2270     TEMP3(i,L) = 0.
2271     CLSBTH(i,L) = 0.
2272     cldlz(i,L) = 0.
2273     enddo
2274     enddo
2275    
2276     C RHO(ZERO) / RHO FOR TERMINAL VELOCITY APPROX.
2277     c ---------------------------------------------
2278     DO L = 1,NLAY
2279     DO I = 1,IRUN
2280     TEMP2(I,L) = PL(I,L)*ZP001
2281     TEMP2(I,L) = SQRT(TEMP2(I,L))
2282     ENDDO
2283     ENDDO
2284    
2285     C INVERSE OF MASS IN EACH LAYER
2286     c -----------------------------
2287     DO L = 1,NLAY
2288     DO I = 1,IRUN
2289 molod 1.6 TEMP3(I,L) = GRAVITY*ZP01 / DP(I,L)
2290 molod 1.1 ENDDO
2291     ENDDO
2292    
2293     C DO LOOP FOR MOISTURE EVAPORATION ABILITY AND CONVEC EVAPORATION.
2294     c ----------------------------------------------------------------
2295     DO 100 L=1,NLAY
2296    
2297     DO I = 1,IRUN
2298     TEMP1(I,3) = TL(I,L)
2299     TEMP1(I,4) = QL(I,L)
2300     ENDDO
2301    
2302     DO 50 N=1,2
2303     IF(N.EQ.1)RELAX=HALF
2304     IF(N.GT.1)RELAX=ONE
2305    
2306     DO I = 1,IRUN
2307     call qsat ( temp1(i,3),pl(i,L),temp1(i,2),temp1(i,6),.true. )
2308     TEMP1(I,5)=TEMP1(I,2)-TEMP1(I,4)
2309     TEMP1(I,6)=TEMP1(I,6)*ELOCP
2310     TEMP1(I,5)=TEMP1(I,5)/(ONE+TEMP1(I,6))
2311     TEMP1(I,4)=TEMP1(I,4)+TEMP1(I,5)*RELAX
2312     TEMP1(I,3)=TEMP1(I,3)-TEMP1(I,5)*ELOCP*RELAX
2313     ENDDO
2314     50 CONTINUE
2315    
2316     DO I = 1,IRUN
2317     EVP9(I,L) = (TEMP1(I,4) - QL(I,L))/TEMP3(I,L)
2318     C convective detrained water
2319     cldlz(i,L) = rain(i,L)*temp3(i,L)
2320     if( tl(i,L).gt.tice-20.) then
2321     water(i) = water(i) + rain(i,L)
2322     else
2323     crystal(i) = crystal(i) + rain(i,L)
2324     endif
2325     ENDDO
2326    
2327     C**********************************************************************
2328     C FOR CONVECTIVE PRECIP, FIND THE "EVAPORATION EFFICIENCY" USING *
2329     C KESSLERS PARAMETERIZATION *
2330     C**********************************************************************
2331    
2332     DO 20 I=1,IRUN
2333    
2334     iceevap(i) = 0.
2335     watevap(i) = 0.
2336    
2337     if( (evp9(i,L).gt.0.) .and. (crystal(i).gt.0.) ) then
2338     iceevap(I) = EVP9(I,L)*fracice
2339     IF(iceevap(i).GE.crystal(i)) iceevap(i) = crystal(i)
2340     EVP9(I,L)=EVP9(I,L)-iceevap(I)
2341     crystal(i) = crystal(i) - iceevap(i)
2342     endif
2343    
2344     C and now warm precipitate
2345     if( (evp9(i,L).gt.0.) .and. (water(i).gt.0.) ) then
2346     exparg = ZM1P04*tmscl*((water(i)*RPHF*TEMP2(I,L))**ZP578)
2347     AREARAT = ONE-(EXP(EXPARG))
2348     watevap(I) = EVP9(I,L)*AREARAT*fracwat
2349     IF(watevap(I).GE.water(i)) watevap(I) = water(i)
2350     EVP9(I,L)=EVP9(I,L)-watevap(I)
2351     water(i) = water(i) - watevap(i)
2352     endif
2353    
2354     QL(I,L) = QL(I,L)+(iceevap(i)+watevap(i))*TEMP3(I,L)
2355     TL(I,L) = TL(I,L)-(iceevap(i)+watevap(i))*TEMP3(I,L)*ELOCP
2356    
2357     20 CONTINUE
2358    
2359     100 CONTINUE
2360    
2361     do i = 1,irun
2362     rcon(i) = water(i) + crystal(i)
2363     enddo
2364    
2365     C**********************************************************************
2366     C Large Scale Precip
2367     C**********************************************************************
2368    
2369     DO 200 L=1,NLAY
2370     DO I = 1,IRUN
2371     rainmax(i) = rhcrit(i,L)*evp9(i,L) +
2372     . ql(i,L)*(rhcrit(i,L)-1.)/temp3(i,L)
2373    
2374     if (rainmax(i).LE.0.0) then
2375     call qsat( tl(i,L),pl(i,L),rh,dum,.false.)
2376     rh = ql(i,L)/rh
2377    
2378     if( rhcrit(i,L).eq.1.0 ) then
2379     fact = 1.0
2380     else
2381     fact = min( 1.0, alpha + (1.0-alpha)*( rh-rhcrit(i,L)) /
2382     1 (1.0-rhcrit(i,L)) )
2383     endif
2384    
2385     C Do not allow clouds above 10 mb
2386     if( pl(i,L).ge.10.0 ) CLSBTH(I,L) = fact
2387     RLAR(I) = RLAR(I)-rainmax(I)
2388     QL(I,L) = QL(I,L)+rainmax(I)*TEMP3(I,L)
2389     TL(I,L) = TL(I,L)-rainmax(I)*TEMP3(I,L)*ELOCP
2390     C Large-scale water
2391     cldlz(i,L) = cldlz(i,L) - rainmax(i)*temp3(i,L)
2392     ENDIF
2393     ENDDO
2394    
2395     DO I=1,IRUN
2396     IF((RLAR(I).GT.0.0).AND.(rainmax(I).GT.0.0))THEN
2397     RPOW=(RLAR(I)*RPHF*TEMP2(I,L))**ZP578
2398     EXPARG = ZM1P04*tmscl*RPOW
2399     AREARAT = ONE-(EXP(EXPARG))
2400     TEMP1(I,7) = rainmax(I)*AREARAT
2401     IF(TEMP1(I,7).GE.RLAR(I)) TEMP1(I,7) = RLAR(I)
2402     RLAR(I) = RLAR(I)-TEMP1(I,7)
2403     QL(I,L) = QL(I,L)+TEMP1(I,7)*TEMP3(I,L)
2404     TL(I,L) = TL(I,L)-TEMP1(I,7)*TEMP3(I,L)*ELOCP
2405     ENDIF
2406     ENDDO
2407    
2408     200 CONTINUE
2409    
2410     RETURN
2411     END
2412     subroutine srclouds (th,q,plk,pl,plke,cloud,cldwat,irun,irise,
2413     1 rhc,offset,alpha)
2414     C***********************************************************************
2415     C
2416     C PURPOSE:
2417     C ========
2418     C Compute non-precipitating cloud fractions
2419     C based on Slingo and Ritter (1985).
2420     C Remove cloudiness where conditionally unstable.
2421     C
2422     C INPUT:
2423     C ======
2424     C th ......... Potential Temperature (irun,irise)
2425     C q .......... Specific Humidity (irun,irise)
2426     C plk ........ P**Kappa at mid-layer (irun,irise)
2427     C pl ......... Pressure at mid-layer (irun,irise)
2428     C plke ....... P**Kappa at edge (irun,irise+1)
2429     C irun ....... Horizontal dimension
2430     C irise ...... Vertical dimension
2431     C
2432     C OUTPUT:
2433     C =======
2434     C cloud ...... Cloud Fraction (irun,irise)
2435     C
2436     C***********************************************************************
2437    
2438     implicit none
2439     integer irun,irise
2440    
2441 molod 1.13 _RL th(irun,irise)
2442     _RL q(irun,irise)
2443     _RL plk(irun,irise)
2444     _RL pl(irun,irise)
2445     _RL plke(irun,irise+1)
2446    
2447     _RL cloud(irun,irise)
2448     _RL cldwat(irun,irise)
2449     _RL qs(irun,irise)
2450    
2451     _RL cp, alhl, getcon, akap
2452     _RL ratio, temp, elocp
2453     _RL rhcrit,rh,dum
2454 molod 1.9 integer i,L
2455 molod 1.1
2456 molod 1.13 _RL rhc(irun,irise)
2457     _RL offset,alpha
2458 molod 1.1
2459     c Explicit Inline Directives
2460     c --------------------------
2461 molod 1.10 #ifdef CRAY
2462     #ifdef f77
2463 molod 1.1 cfpp$ expand (qsat)
2464     #endif
2465     #endif
2466    
2467     cp = getcon('CP')
2468     alhl = getcon('LATENT HEAT COND')
2469     elocp = alhl/cp
2470     akap = getcon('KAPPA')
2471    
2472     do L = 1,irise
2473     do i = 1,irun
2474     temp = th(i,L)*plk(i,L)
2475     call qsat ( temp,pl(i,L),qs(i,L),dum,.false. )
2476     enddo
2477     enddo
2478    
2479     do L = 2,irise
2480     do i = 1,irun
2481     rh = q(i,L)/qs(i,L)
2482    
2483     rhcrit = rhc(i,L) - offset
2484     ratio = alpha*(rh-rhcrit)/offset
2485    
2486     if(cloud(i,L).eq. 0.0 .and. ratio.gt.0.0 ) then
2487     cloud(i,L) = min( ratio,1.0 )
2488     endif
2489    
2490     enddo
2491     enddo
2492    
2493     c Reduce clouds from conditionally unstable layer
2494     c -----------------------------------------------
2495     call ctei ( th,q,cloud,cldwat,pl,plk,plke,irun,irise )
2496    
2497     return
2498     end
2499    
2500     subroutine ctei ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm )
2501     implicit none
2502     integer im,lm
2503 molod 1.13 _RL th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm)
2504     _RL plk(im,lm),pl(im,lm),cldfrc(im,lm)
2505 molod 1.1 integer i,L
2506 molod 1.13 _RL getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq
2507     _RL k,krd,kmm,f
2508 molod 1.1
2509     cp = getcon('CP')
2510     alhl = getcon('LATENT HEAT COND')
2511     cpoel = cp/alhl
2512     elocp = alhl/cp
2513    
2514     do L=lm,2,-1
2515     do i=1,im
2516     dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1)
2517     if( dq.eq.0.0 ) dq = 1.0e-20
2518     k = 1.0 + cpoel*plke(i,L)*( th(i,L)-th(i,L-1) ) / dq
2519    
2520     t = th(i,L)*plk(i,L)
2521     p = pl(i,L)
2522     call qsat ( t,p,qs,dqsdt,.true. )
2523    
2524     krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t )
2525    
2526     kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t )
2527     . / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt )
2528    
2529     s = ( (k-krd)/(kmm-krd) )
2530     f = 1.0 - min( 1.0, max(0.0,1.0-exp(-s)) )
2531    
2532     cldfrc(i,L) = cldfrc(i,L)*f
2533     cldwat(i,L) = cldwat(i,L)*f
2534    
2535     enddo
2536     enddo
2537    
2538     return
2539     end
2540    
2541     subroutine back2grd(gathered,indeces,scattered,irun)
2542     implicit none
2543     integer i,irun,indeces(irun)
2544 molod 1.13 _RL gathered(irun),scattered(irun)
2545     _RL temp(irun)
2546 molod 1.1 do i = 1,irun
2547     temp(indeces(i)) = gathered(i)
2548     enddo
2549     do i = 1,irun
2550     scattered(i) = temp(i)
2551     enddo
2552     return
2553     end
  ViewVC Help
Powered by ViewVC 1.1.22