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