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

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

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

Revision 1.20 - (show 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 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.19 2004/09/30 19:45:24 molod Exp $
2 C $Name: $
3
4 #include "FIZHI_OPTIONS.h"
5 subroutine moistio (ndmoist,istrip,npcs,
6 . lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup,
7 . pz,plz,plze,dpres,pkht,pkl,tz,qz,bi,bj,ntracer,ptracer,
8 . qqz,dumoist,dvmoist,dtmoist,dqmoist,
9 . im,jm,lm,ptop,
10 . 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 . lpnt,myid)
14
15 implicit none
16
17 #ifdef ALLOW_DIAGNOSTICS
18 #include "SIZE.h"
19 #include "diagnostics_SIZE.h"
20 #include "diagnostics.h"
21 #endif
22
23 c Input Variables
24 c ---------------
25 integer im,jm,lm
26 integer ndmoist,istrip,npcs
27 integer bi,bj,ntracer,ptracer
28 integer lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup
29 _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 integer iras
37 _RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
38 integer nswcld,nswlz
39 _RL cldlsp_sw(im,jm,lm),cldras_sw(im,jm,lm)
40 _RL cldtot_sw(im,jm,lm),swlz(im,jm,lm)
41 integer nlwcld,nlwlz
42 _RL cldlsp_lw(im,jm,lm),cldras_lw(im,jm,lm)
43 _RL cldtot_lw(im,jm,lm),lwlz(im,jm,lm)
44 logical lpnt
45 integer myid
46
47 c Local Variables
48 c ---------------
49 integer ncrnd,nsecf
50
51 _RL fracqq, dum
52 integer snowcrit
53 parameter (fracqq = 0.1)
54 _RL one
55 parameter (one=1.)
56
57 _RL cldsr(im,jm,lm)
58 _RL srcld(istrip,lm)
59
60 _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
73 _RL qqcolmin(im,jm)
74 _RL qqcolmax(im,jm)
75 integer levpbl(im,jm)
76
77 c Gathered Arrays for Variable Cloud Base
78 c ---------------------------------------
79 _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
96 integer nindeces(lm)
97 integer pblindex(im*jm)
98 integer levgather(im*jm)
99
100 c Stripped Arrays
101 c ---------------
102 _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 integer ITMP1(ISTRIP,lm), ITMP2(ISTRIP,lm)
122
123 _RL PRECIP(ISTRIP), PCNET(ISTRIP)
124 _RL SP(ISTRIP), PREP(ISTRIP)
125 _RL PCPEN (ISTRIP,lm)
126 integer pbl(istrip),depths(lm)
127
128 _RL cldlz(istrip,lm), cldwater(im,jm,lm)
129 _RL rhfrac(istrip), rhmin, pup, ppbl, rhcrit(istrip,lm)
130 _RL offset, alpha, rasmax
131
132 logical first
133 logical lras
134 _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 DATA FIRST /.TRUE./
141
142 integer imstp,nsubcl,nlras
143 integer i,j,iloop,indx,indgath,l,nn,num,numdeps,nt
144 _RL tmstp,tminv,sday,grav,alhl,cp,elocp,gamfac
145 _RL rkappa,p0kappa,p0kinv,ptopkap,pcheck
146 _RL tice,getcon,pi
147 integer ntracedim
148
149 C **********************************************************************
150 C **** INITIALIZATION ****
151 C **********************************************************************
152
153 ntracedim= max(ntracer-ptracer,1)
154 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 C Difference in fraction between SR and LS Threshold
161 offset = 0.10
162 C Large-Scale Relative Humidity Threshold in PBL
163 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 c Determine Total number of Random Clouds to Check
188 c ---------------------------------------------
189 C ncrnd = (lm-nltop+1)/2
190 ncrnd = 0
191
192 if(first .and. myid.eq.1 .and. bi.eq.1 ) then
193 print *
194 print *,'Top Level Allowed for Convection : ',nltop
195 print *,' Highest Sub-Cloud Level: ',nsubmax
196 print *,' Lowest Sub-Cloud Level: ',nsubmin
197 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 indx = 0
237 do L = nsubmin,nltop,-1
238 do j = 1,jm
239 do i = 1,im
240 if(levpbl(i,j).eq.L) then
241 indx = indx + 1
242 pblindex(indx) = (j-1)*im + i
243 endif
244 enddo
245 enddo
246 enddo
247
248 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 enddo
254
255 do L = 1,lm
256 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 enddo
265 enddo
266 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
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 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 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 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
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 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
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 TMP5(I,L) = PLKE(I,L)*P0KINV
414 ENDDO
415 ENDDO
416 DO I=num,num+nindeces(nsubcl)-1
417 TMP4(I,lm+1) = PLE (I,lm+1)
418 TMP5(I,lm+1) = PLKE(I,lm+1)*P0KINV
419 ENDDO
420 DO 113 I=num,num+nindeces(nsubcl)-1
421 TMP4(I,NSUBCL+1) = PLE (I,lm+1)
422 TMP5(I,NSUBCL+1) = PLKE(I,lm+1)*P0KINV
423 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 do L = nsubcl, lm
456 rhcrit(i,L) = 1.
457 enddo
458 do L = 1, nsubcl-1
459 pcheck = pl(i,L)
460 if (pcheck .le. pup) then
461 rhcrit(i,L) = rhmin
462 else
463 ppbl = pl(i,nsubcl)
464 rhcrit(i,L) = rhmin + (1.-rhmin)/(19.) *
465 . ((atan( (2.*(pcheck-pup)/(ppbl-pup)-1.) *
466 . tan(20.*pi/21.-0.5*pi) )
467 . + 0.5*pi) * 21./pi - 1.)
468 endif
469 enddo
470 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 1, UL(num,1,1),ntracedim,TH(num,NLTOP),SHL(num,NLTOP)
485 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 dum = dp(i,L)/(ple(i,lm+1)-ple(i,nsubcl))
495 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 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
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 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
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 CALL RNEVP (NN,ISTRIP,lm,TL,SHL,PCPEN,PL,CLFRAC,SP,DP,PLKE,
609 . PLK,TH,TMP1,TMP2,TMP3,ITMP1,ITMP2,PCNET,PRECIP,
610 . CLSBTH,TMSTP,one,cp,grav,alhl,gamfac,cldlz,rhcrit,offset,alpha)
611
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 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
777
778 C **********************************************************************
779 C BUMP DIAGNOSTICS
780 C **********************************************************************
781
782
783 c Sub-Cloud Layer
784 c -------------------------
785 if( ipsubcld.ne.0 ) then
786 do j = 1,jm
787 do i = 1,im
788 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 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 qdiag(i,j,icldnp+L-1,bi,bj) = qdiag(i,j,icldnp+L-1,bi,bj) +
803 . cldsr(i,j,L)
804 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 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 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 qdiag(i,j,imoistq+L-1,bi,bj) = qdiag(i,j,imoistq+L-1,bi,bj) +
831 . (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 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 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 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 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 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 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 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 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 qdiag(i,j,ipreacc,bi,bj) = qdiag(i,j,ipreacc,bi,bj)
899 . + ( 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 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 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 plev = pl(i,L)
956
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 qdiag(i,j,icldtmp,bi,bj) = qdiag(i,j,icldtmp,bi,bj) +
1037 . cldtmp(i,j)*totcld(i,j)/tmpimjm(i,j)
1038 qdiag(i,j,icttcnt,bi,bj) = qdiag(i,j,icttcnt,bi,bj) +
1039 . totcld(i,j)
1040 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 qdiag(i,j,icldprs,bi,bj) = qdiag(i,j,icldprs,bi,bj) +
1050 . cldprs(i,j)*totcld(i,j)/tmpimjm(i,j)
1051 qdiag(i,j,ictpcnt,bi,bj) = qdiag(i,j,ictpcnt,bi,bj) +
1052 . totcld(i,j)
1053 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 #ifdef ALLOW_DIAGNOSTICS
1069 if( (bi.eq.1) .and. (bj.eq.1) ) then
1070 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 endif
1081 #endif
1082
1083 RETURN
1084 END
1085 SUBROUTINE RAS( NN, LNG, LENC, K, NLTOP, nlayr, DT
1086 *, 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 implicit none
1096
1097 C Argument List
1098 integer nn,lng,lenc,k,nltop,nlayr
1099 integer ntracer
1100 integer ncrnd
1101 _RL dt
1102 _RL UOI(lng,nlayr,ntracer), POI(lng,K)
1103 _RL QOI(lng,K), PRS(lng,K+1), PRJ(lng,K+1)
1104 _RL rnd(ncrnd)
1105 _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
1109 C Local Variables
1110 _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
1125 integer krmin,icm
1126 _RL rknob, cmb2pa
1127 PARAMETER (KRMIN=01)
1128 PARAMETER (ICM=1000)
1129 PARAMETER (CMB2PA=100.0)
1130 PARAMETER (rknob = 10.)
1131
1132 integer IC(ICM), IRND(icm)
1133 _RL cmass(lng,K)
1134 LOGICAL SETRAS
1135
1136 integer i,L,nc,ib,nt
1137 integer km1,kp1,kprv,kcr,kfx,ncmx
1138 _RL p00, crtmsf, frac, rasblf
1139
1140 do L = 1,k
1141 do I = 1,LENC
1142 rains(i,l) = 0.
1143 enddo
1144 enddo
1145
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
1192 c Initialize Cloud Fraction Array
1193 c -------------------------------
1194 do i = 1,lenc
1195 cloudn(i) = 0.0
1196 enddo
1197
1198 CALL CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IB, RASBLF,SETRAS,FRAC
1199 *, 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 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 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 _RL random_numbx
1275 _RL rnd(nrnd)
1276 integer irm
1277 parameter (irm = 1000)
1278 _RL random(irm)
1279 integer i,mcheck,numrand,iseed,indx
1280 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 if(first .and. myid.eq.1) print *,' NO RANDOM CLOUDS IN RAS '
1291 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 print *,' first ',first,' iras ',iras,' iras0 ',iras0
1300 if( myid.eq.1 ) print *, 'Recreating Rand Numb Array in RNDCLOUD'
1301 if( myid.eq.1 ) print *, 'IRAS: ',iras,' IRAS0: ',iras0
1302 numrand = mod(iras,irm/nrnd) * nrnd
1303 iseed = iras * nrnd - numrand
1304 call random_seedx(iseed)
1305 do i = 1,irm
1306 random(i) = random_numbx(iseed)
1307 enddo
1308 indx = (iras-1)*nrnd
1309
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 random(i) = random_numbx(iseed)
1317 enddo
1318 indx = iseed
1319
1320 c Multiple Time In But have NOT Used Up all 1000 numbers (MCHECK.NE.0)
1321 c --------------------------------------------------------------------
1322 else
1323 indx = (iras-1)*nrnd
1324 endif
1325
1326 indx = mod(indx,irm)
1327 if( indx+nrnd.gt.1000 ) indx=1000-nrnd
1328
1329 do n = 1,nrnd
1330 rnd(n) = random(indx+n)
1331 enddo
1332
1333 100 continue
1334 first = .false.
1335 iras0 = iras
1336 return
1337 end
1338 function random_numbx(iseed)
1339 implicit none
1340 integer iseed
1341 real *8 seed,port_rand
1342 _RL random_numbx
1343 random_numbx = 0
1344 #ifdef CRAY
1345 _RL ranf
1346 random_numbx = ranf()
1347 #else
1348 #ifdef SGI
1349 _RL rand
1350 random_numbx = rand()
1351 #endif
1352 random_numbx = port_rand(seed)
1353 #endif
1354 return
1355 end
1356 subroutine random_seedx (iseed)
1357 implicit none
1358 integer iseed
1359 #ifdef CRAY
1360 call ranset (iseed)
1361 #endif
1362 #ifdef SGI
1363 integer*4 seed
1364 seed = iseed
1365 call srand (seed)
1366 #endif
1367 return
1368 end
1369 SUBROUTINE CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IC, RASALF
1370 *, 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 C lng : The inner dimension of update and input arrays.
1412 C
1413 C LENC : The run: the number of soundings processes in a single call.
1414 C RAS works on the first LENC of the lng soundings
1415 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 C (LENC<=lng).
1421 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 C POI : 2D array of dimension (lng,K) containing potential
1452 C temperature. Updated but not initialized by RAS.
1453 C
1454 C QOI : 2D array of dimension (lng,K) containing specific
1455 C humidity. Updated but not initialized by RAS.
1456 C
1457 C UOI : 3D array of dimension (lng,K,NTRACER) containing tracers
1458 C Updated but not initialized by RAS.
1459 C
1460 C PRS : 2D array of dimension (lng,K+1) containing pressure
1461 C in hPa at the interfaces of K-layers from top of the
1462 C atmosphere to the bottom. Not modified.
1463 C
1464 C PRJ : 2D array of dimension (lng,K+1) containing (PRS/P00) **
1465 C RKAP. i.e. Exner function at layer edges. Not modified.
1466 C
1467 C rhfrac : 1D array of dimension (lng) containing a rel.hum. scaling
1468 C fraction. Not modified.
1469 C
1470 C Output:
1471 C -------
1472 C
1473 C PCU : 1D array of length lng containing accumulated
1474 C precipitation in mm/sec.
1475 C
1476 C CLN : 2D array of dimension (lng,K) containing cloudiness
1477 C Note: CLN is bumped but NOT initialized
1478 C
1479 C TCU : 2D array of dimension (lng,K) containing accumulated
1480 C convective heating (K/sec).
1481 C
1482 C QCU : 2D array of dimension (lng,K) containing accumulated
1483 C convective drying (kg/kg/sec).
1484 C
1485 C CMASS : 2D array of dimension (lng,K) containing the
1486 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 implicit none
1504 C Argument List declarations
1505 integer nn,lng,LENC,K,NLTOP,nlayr,ic,ntracer
1506 _RL rasalf
1507 LOGICAL SETRAS
1508 _RL frac, cp, alhl, rkap, grav, p00, crtmsf
1509 _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 _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 _RL uht(lng,ntracer)
1521 integer IA(LENC), I1(LENC),I2(LENC)
1522 _RL rhfrac(lng)
1523
1524 C Local Variables
1525 _RL daylen,half,one,zero,cmb2pa,rhmax
1526 PARAMETER (DAYLEN=86400.0, HALF=0.5, ONE=1.0, ZERO=0.0)
1527 PARAMETER (CMB2PA=100.0)
1528 PARAMETER (RHMAX=0.9999)
1529 _RL rkapp1,onebcp,albcp,onebg,cpbg,twobal
1530 C
1531 integer nt,km1,ic1,i,L,len1,len2,isav,len11,ii
1532 integer lena,lena1,lenb
1533 _RL tem,tem1
1534
1535 c Explicit Inline Directives
1536 c --------------------------
1537 #ifdef CRAY
1538 #ifdef f77
1539 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 C SETTING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T.
1554 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 DO 10 I=1,lng
1583 TCU(I,L) = 0.0
1584 QCU(I,L) = 0.0
1585 CMASS(I,L) = 0.0
1586 10 CONTINUE
1587
1588 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 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 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 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 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 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 c do nt = 1,ntracer
2070 c
2071 c Tracer Tendency at the Bottom Layer
2072 c -----------------------------------
2073 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 c Tracer Tendency at all other Levels
2081 c -----------------------------------
2082 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 C
2111 C PENETRATIVE CONVECTION CALCULATION OVER
2112 C
2113
2114 RETURN
2115 END
2116 SUBROUTINE RNCL(lng, PL, RNO, CLF)
2117 C
2118 C*********************************************************************
2119 C********************** Relaxed Arakawa-Schubert *********************
2120 C************************ SUBROUTINE RNCL ************************
2121 C**************************** 23 July 1992 ***************************
2122 C*********************************************************************
2123 implicit none
2124 C Argument List declarations
2125 integer lng
2126 _RL PL(lng), RNO(lng), CLF(lng)
2127
2128 C Local Variables
2129 _RL p5,p8,pt8,pt2,pfac,p4,p6,p7,p9,cucld,cfac
2130 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
2136 integer i
2137 C
2138 DO 10 I=1,lng
2139 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 SUBROUTINE ACRITN ( lng,PL,PLB,ACR )
2166
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 implicit none
2177 C Argument List declarations
2178 integer lng
2179 _RL PL(lng), PLB(lng), ACR(lng)
2180
2181 C Local variables
2182 integer lma
2183 parameter (lma=18)
2184 _RL p(lma)
2185 _RL a(lma)
2186 integer i,L
2187 _RL temp
2188
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 do i=1,lng
2202 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 do i=1,lng
2211 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 do i=1,lng
2216 acr(i) = acr(i) * (plb(i)-pl(i))
2217 enddo
2218
2219 RETURN
2220 END
2221 SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DP,PLKE,
2222 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 implicit none
2226 C Argument List declarations
2227 integer nn,irun,nlay
2228 _RL TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY),
2229 . 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 _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 C
2239 C Local Variables
2240 _RL zm1p04,zero,two89,zp44,zp01,half,zp578,one,thousand,z3600
2241 _RL zp1,zp001
2242 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 _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
2263 integer i,L,n,nlaym1,irnlay,irnlm1
2264
2265 c Explicit Inline Directives
2266 c --------------------------
2267 #ifdef CRAY
2268 #ifdef f77
2269 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 TEMP3(I,L) = GRAVITY*ZP01 / DP(I,L)
2318 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 _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 integer i,L
2483
2484 _RL rhc(irun,irise)
2485 _RL offset,alpha
2486
2487 c Explicit Inline Directives
2488 c --------------------------
2489 #ifdef CRAY
2490 #ifdef f77
2491 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 _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 integer i,L
2534 _RL getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq
2535 _RL k,krd,kmm,f
2536
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 _RL gathered(irun),scattered(irun)
2573 _RL temp(irun)
2574 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