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