/[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.27 - (show annotations) (download)
Tue Mar 1 16:53:05 2005 UTC (19 years, 4 months ago) by jmc
Branch: MAIN
Changes since 1.26: +9 -3 lines 
updated after changing port_rand function.
1 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.26 2005/03/01 15:25:11 jmc 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 = -1.d0
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 #else
1470 #ifdef SGI
1471 integer*4 seed
1472 seed = iseed
1473 call srand (seed)
1474 #else
1475 real*8 tmpRdN
1476 real*8 seed
1477 seed = iseed
1478 tmpRdN = port_rand(seed)
1479 #endif
1480 #endif
1481 return
1482 end
1483 SUBROUTINE CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IC, RASALF
1484 *, SETRAS, FRAC
1485 *, CP, ALHL, RKAP, GRAV, P00, CRTMSF
1486 *, POI, QOI, UOI, ntracedim, Ntracer, PRS, PRJ
1487 *, PCU, CLN, TCU, QCU, UCU, CMASS
1488 *, ALF, BET, GAM, PRH, PRI, HOL, ETA
1489 *, HST, QOL, GMH
1490 *, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, ALM
1491 *, WFN, AKM, QS1, CLF, UHT, WLQ
1492 *, IA, I1, I2,rhfrac)
1493 C
1494 C*********************************************************************
1495 C******************** Relaxed Arakawa-Schubert ***********************
1496 C********************* Plug Compatible Version **********************
1497 C************************ SUBROUTINE CLOUD ***************************
1498 C************************* 23 JULY 1992 ***************************
1499 C*********************************************************************
1500 C*********************************************************************
1501 C*********************************************************************
1502 C************************** Developed By *****************************
1503 C************************** *****************************
1504 C************************ Shrinivas Moorthi **************************
1505 C************************ and **************************
1506 C************************ Max J. Suarez *****************************
1507 C************************ *****************************
1508 C******************** Laboratory for Atmospheres *********************
1509 C****************** NASA/GSFC, Greenbelt, MD 20771 *******************
1510 C*********************************************************************
1511 C*********************************************************************
1512 C
1513 C The calculations of Moorthi and Suarez (1992, MWR) are
1514 C contained in the CLOUD routine.
1515 C It is probably advisable, at least initially, to treat CLOUD
1516 C as a black box that computes the single cloud adjustments. RAS,
1517 C on the other hand, can be tailored to each GCMs configuration
1518 C (ie, number and placement of levels, nature of boundary layer,
1519 C time step and frequency with which RAS is called).
1520 C
1521 C
1522 C Input:
1523 C ------
1524 C
1525 C lng : The inner dimension of update and input arrays.
1526 C
1527 C LENC : The run: the number of soundings processes in a single call.
1528 C RAS works on the first LENC of the lng soundings
1529 C passed. This allows working on pieces of the world
1530 C say for multitasking, without declaring temporary arrays
1531 C and copying the data to and from them. This is an f77
1532 C version. An F90 version would have to allow more
1533 C flexibility in the argument declarations. Obviously
1534 C (LENC<=lng).
1535 C
1536 C K : Number of vertical layers (increasing downwards).
1537 C Need not be the same as the number of layers in the
1538 C GCM, since it is the outer dimension. The bottom layer
1539 C (K) is the subcloud layer.
1540 C
1541 C IC : Detrainment level to check for presence of convection
1542 C
1543 C RASALF : Relaxation parameter (< 1.) for present cloud-type
1544 C
1545 C SETRAS : Logical parameter to control re-calculation of
1546 C saturation specific humidity and mid level P**kappa
1547 C
1548 C FRAC : Fraction of the PBL (layer K) mass allowed to be used
1549 C by a cloud-type in time DT
1550 C
1551 C CP : Specific heat at constant pressure
1552 C
1553 C ALHL : Latent Heat of condensation
1554 C
1555 C RKAP : R/Cp, where R is the gas constant
1556 C
1557 C GRAV : Acceleration due to gravity
1558 C
1559 C P00 : A reference pressure in hPa, useually 1000 hPa
1560 C
1561 C CRTMSF : Critical value of mass flux above which cloudiness at
1562 C the detrainment layer of that cloud-type is assumed.
1563 C Affects only cloudiness calculation.
1564 C
1565 C POI : 2D array of dimension (lng,K) containing potential
1566 C temperature. Updated but not initialized by RAS.
1567 C
1568 C QOI : 2D array of dimension (lng,K) containing specific
1569 C humidity. Updated but not initialized by RAS.
1570 C
1571 C UOI : 3D array of dimension (lng,K,NTRACER) containing tracers
1572 C Updated but not initialized by RAS.
1573 C
1574 C PRS : 2D array of dimension (lng,K+1) containing pressure
1575 C in hPa at the interfaces of K-layers from top of the
1576 C atmosphere to the bottom. Not modified.
1577 C
1578 C PRJ : 2D array of dimension (lng,K+1) containing (PRS/P00) **
1579 C RKAP. i.e. Exner function at layer edges. Not modified.
1580 C
1581 C rhfrac : 1D array of dimension (lng) containing a rel.hum. scaling
1582 C fraction. Not modified.
1583 C
1584 C Output:
1585 C -------
1586 C
1587 C PCU : 1D array of length lng containing accumulated
1588 C precipitation in mm/sec.
1589 C
1590 C CLN : 2D array of dimension (lng,K) containing cloudiness
1591 C Note: CLN is bumped but NOT initialized
1592 C
1593 C TCU : 2D array of dimension (lng,K) containing accumulated
1594 C convective heating (K/sec).
1595 C
1596 C QCU : 2D array of dimension (lng,K) containing accumulated
1597 C convective drying (kg/kg/sec).
1598 C
1599 C CMASS : 2D array of dimension (lng,K) containing the
1600 C cloud mass flux (kg/sec). Filled from cloud top
1601 C to base.
1602 C
1603 C Temporaries:
1604 C
1605 C ALF, BET, GAM, ETA, PRH, PRI, HOI, HST, QOL, GMH are temporary
1606 C 2D real arrays of dimension of at least (LENC,K) where LENC is
1607 C the horizontal dimension over which convection is invoked.
1608 C
1609 C
1610 C TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9, AKM, QS1, CLF, UHT
1611 C VHT, WLQ WFN are temporary real arrays of length at least LENC
1612 C
1613 C IA, I1, and I2 are temporary integer arrays of length LENC
1614 C
1615 C
1616 C************************************************************************
1617 implicit none
1618 C Argument List declarations
1619 integer nn,lng,LENC,K,NLTOP,nlayr,ic,ntracedim, ntracer
1620 _RL rasalf
1621 LOGICAL SETRAS
1622 _RL frac, cp, alhl, rkap, grav, p00, crtmsf
1623 _RL POI(lng,K),QOI(lng,K),PRS(lng,K+1),PRJ(lng,K+1)
1624 _RL uoi(lng,nlayr,ntracedim)
1625 _RL PCU(LENC), CLN(lng)
1626 _RL TCU(lng,K),QCU(lng,K),ucu(lng,k,ntracedim),CMASS(lng,K)
1627 _RL ALF(lng,K), BET(lng,K), GAM(lng,K), PRH(lng,K), PRI(lng,K)
1628 _RL HOL(LENC,K), ETA(LENC,K), HST(LENC,K), QOL(LENC,K)
1629 _RL GMH(LENC,K)
1630 _RL TX1(LENC), TX2(LENC), TX3(LENC), TX4(LENC)
1631 _RL TX5(LENC), TX6(LENC), TX7(LENC), TX8(LENC)
1632 _RL ALM(LENC), WFN(LENC), AKM(LENC), QS1(LENC)
1633 _RL WLQ(LENC), CLF(LENC)
1634 _RL uht(lng,ntracedim)
1635 integer IA(LENC), I1(LENC),I2(LENC)
1636 _RL rhfrac(lng)
1637
1638 C Local Variables
1639 _RL daylen,half,one,zero,cmb2pa,rhmax
1640 PARAMETER (DAYLEN=86400.0, HALF=0.5, ONE=1.0, ZERO=0.0)
1641 PARAMETER (CMB2PA=100.0)
1642 PARAMETER (RHMAX=0.9999)
1643 _RL rkapp1,onebcp,albcp,onebg,cpbg,twobal
1644 C
1645 integer nt,km1,ic1,i,L,len1,len2,isav,len11,ii
1646 integer lena,lena1,lenb
1647 _RL tem,tem1
1648
1649 c Explicit Inline Directives
1650 c --------------------------
1651 #ifdef CRAY
1652 #ifdef f77
1653 cfpp$ expand (qsat)
1654 #endif
1655 #endif
1656
1657 RKAPP1 = 1.0 + RKAP
1658 ONEBCP = 1.0 / CP
1659 ALBCP = ALHL * ONEBCP
1660 ONEBG = 1.0 / GRAV
1661 CPBG = CP * ONEBG
1662 TWOBAL = 2.0 / ALHL
1663
1664 C
1665 KM1 = K - 1
1666 IC1 = IC + 1
1667 C
1668 C SETTING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T.
1669 C
1670
1671 IF (SETRAS) THEN
1672
1673 DO 2050 L=1,K
1674 DO 2030 I=1,LENC
1675 PRH(I,L) = (PRJ(I,L+1)*PRS(I,L+1) - PRJ(I,L)*PRS(I,L))
1676 * / ((PRS(I,L+1)-PRS(I,L)) * RKAPP1)
1677 2030 CONTINUE
1678 2050 CONTINUE
1679
1680 DO 2070 L=1,K
1681 DO 2060 I=1,LENC
1682 TX5(I) = POI(I,L) * PRH(I,L)
1683 TX1(I) = (PRS(I,L) + PRS(I,L+1)) * 0.5
1684 TX3(I) = TX5(I)
1685 CALL QSAT(TX3(I), TX1(I), TX2(I), TX4(I), .TRUE.)
1686 ALF(I,L) = TX2(I) - TX4(I) * TX5(I)
1687 BET(I,L) = TX4(I) * PRH(I,L)
1688 GAM(I,L) = 1.0 / ((1.0 + TX4(I)*ALBCP) * PRH(I,L))
1689 PRI(I,L) = (CP/CMB2PA) / (PRS(I,L+1) - PRS(I,L))
1690 2060 CONTINUE
1691 2070 CONTINUE
1692
1693 ENDIF
1694 C
1695 C
1696 DO 10 L=1,K
1697 DO 10 I=1,lng
1698 TCU(I,L) = 0.0
1699 QCU(I,L) = 0.0
1700 CMASS(I,L) = 0.0
1701 10 CONTINUE
1702
1703 do nt = 1,ntracer
1704 do L=1,K
1705 do I=1,LENC
1706 ucu(I,L,nt) = 0.0
1707 enddo
1708 enddo
1709 enddo
1710 C
1711 DO 30 I=1,LENC
1712 TX1(I) = PRJ(I,K+1) * POI(I,K)
1713 QS1(I) = ALF(I,K) + BET(I,K)*POI(I,K)
1714 QOL(I,K) = MIN(QS1(I)*RHMAX,QOI(I,K))
1715
1716 HOL(I,K) = TX1(I)*CP + QOL(I,K)*ALHL
1717 ETA(I,K) = ZERO
1718 TX2(I) = (PRJ(I,K+1) - PRJ(I,K)) * POI(I,K) * CP
1719 30 CONTINUE
1720 C
1721 IF (IC .LT. KM1) THEN
1722 DO 3703 L=KM1,IC1,-1
1723 DO 50 I=1,LENC
1724 QS1(I) = ALF(I,L) + BET(I,L)*POI(I,L)
1725 QOL(I,L) = MIN(QS1(I)*RHMAX,QOI(I,L))
1726 C
1727 TEM1 = TX2(I) + PRJ(I,L+1) * POI(I,L) * CP
1728 HOL(I,L) = TEM1 + QOL(I,L )* ALHL
1729 HST(I,L) = TEM1 + QS1(I) * ALHL
1730
1731 TX1(I) = (PRJ(I,L+1) - PRJ(I,L)) * POI(I,L)
1732 ETA(I,L) = ETA(I,L+1) + TX1(I)*CPBG
1733 TX2(I) = TX2(I) + TX1(I)*CP
1734 50 CONTINUE
1735 C
1736 3703 CONTINUE
1737 ENDIF
1738
1739
1740 DO 70 I=1,LENC
1741 HOL(I,IC) = TX2(I)
1742 QS1(I) = ALF(I,IC) + BET(I,IC)*POI(I,IC)
1743 QOL(I,IC) = MIN(QS1(I)*RHMAX,QOI(I,IC))
1744 c
1745 TEM1 = TX2(I) + PRJ(I,IC1) * POI(I,IC) * CP
1746 HOL(I,IC) = TEM1 + QOL(I,IC) * ALHL
1747 HST(I,IC) = TEM1 + QS1(I) * ALHL
1748 C
1749 TX3(I ) = (PRJ(I,IC1) - PRH(I,IC)) * POI(I,IC)
1750 ETA(I,IC) = ETA(I,IC1) + CPBG * TX3(I)
1751 70 CONTINUE
1752 C
1753 DO 130 I=1,LENC
1754 TX2(I) = HOL(I,K) - HST(I,IC)
1755 TX1(I) = ZERO
1756
1757 130 CONTINUE
1758 C
1759 C ENTRAINMENT PARAMETER
1760 C
1761 DO 160 L=IC,KM1
1762 DO 160 I=1,LENC
1763 TX1(I) = TX1(I) + (HST(I,IC) - HOL(I,L)) * (ETA(I,L) - ETA(I,L+1))
1764 160 CONTINUE
1765 C
1766 LEN1 = 0
1767 LEN2 = 0
1768 ISAV = 0
1769 DO 195 I=1,LENC
1770 IF (TX1(I) .GT. ZERO .AND. TX2(I) .GT. ZERO
1771 . .AND. rhfrac(i).ne.0.0 ) THEN
1772 LEN1 = LEN1 + 1
1773 IA(LEN1) = I
1774 ALM(LEN1) = TX2(I) / TX1(I)
1775 ENDIF
1776 195 CONTINUE
1777 C
1778 LEN2 = LEN1
1779 if (IC1 .lt. K) then
1780 DO 196 I=1,LENC
1781 IF (TX2(I) .LE. 0.0 .AND. (HOL(I,K) .GT. HST(I,IC1))
1782 . .AND. rhfrac(i).ne.0.0 ) THEN
1783 LEN2 = LEN2 + 1
1784 IA(LEN2) = I
1785 ALM(LEN2) = 0.0
1786 ENDIF
1787 196 CONTINUE
1788 endif
1789 C
1790 IF (LEN2 .EQ. 0) THEN
1791 DO 5010 I=1,LENC*K
1792 HST(I,1) = 0.0
1793 QOL(I,1) = 0.0
1794 5010 CONTINUE
1795 DO 5020 I=1,LENC
1796 PCU(I) = 0.0
1797 5020 CONTINUE
1798 RETURN
1799 ENDIF
1800 LEN11 = LEN1 + 1
1801 C
1802 C NORMALIZED MASSFLUX
1803 C
1804 DO 250 I=1,LEN2
1805 ETA(I,K) = 1.0
1806 II = IA(I)
1807 TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
1808 TX4(I) = PRS(II,K)
1809 250 CONTINUE
1810 C
1811 DO 252 I=LEN11,LEN2
1812 WFN(I) = 0.0
1813 II = IA(I)
1814 IF (HST(II,IC1) .LT. HST(II,IC)) THEN
1815 TX6(I) = (HST(II,IC1)-HOL(II,K))/(HST(II,IC1)-HST(II,IC))
1816 ELSE
1817 TX6(I) = 0.0
1818 ENDIF
1819 TX2(I) = 0.5 * (PRS(II,IC1)+PRS(II,IC1+1)) * (1.0-TX6(I))
1820 * + TX2(I) * TX6(I)
1821 252 CONTINUE
1822 C
1823 CALL ACRITN(LEN2, TX2, TX4, TX3)
1824 C
1825 DO 260 L=KM1,IC,-1
1826 DO 255 I=1,LEN2
1827 TX1(I) = ETA(IA(I),L)
1828 255 CONTINUE
1829 DO 260 I=1,LEN2
1830 ETA(I,L) = 1.0 + ALM(I) * TX1(I)
1831 260 CONTINUE
1832 C
1833 C CLOUD WORKFUNCTION
1834 C
1835 IF (LEN1 .GT. 0) THEN
1836 DO 270 I=1,LEN1
1837 II = IA(I)
1838 WFN(I) = - GAM(II,IC) * (PRJ(II,IC1) - PRH(II,IC))
1839 * * HST(II,IC) * ETA(I,IC1)
1840 270 CONTINUE
1841 ENDIF
1842 C
1843 DO 290 I=1,LEN2
1844 II = IA(I)
1845 TX1(I) = HOL(II,K)
1846 290 CONTINUE
1847 C
1848 IF (IC1 .LE. KM1) THEN
1849
1850 DO 380 L=KM1,IC1,-1
1851 DO 380 I=1,LEN2
1852 II = IA(I)
1853 TEM = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * HOL(II,L)
1854 C
1855 PCU(I) = PRJ(II,L+1) - PRH(II,L)
1856 TEM1 = ETA(I,L+1) * PCU(I)
1857 TX1(I) = TX1(I)*PCU(I)
1858 C
1859 PCU(I) = PRH(II,L) - PRJ(II,L)
1860 TEM1 = (TEM1 + ETA(I,L) * PCU(I)) * HST(II,L)
1861 TX1(I) = TX1(I) + TEM*PCU(I)
1862 C
1863 WFN(I) = WFN(I) + (TX1(I) - TEM1) * GAM(II,L)
1864 TX1(I) = TEM
1865 380 CONTINUE
1866 ENDIF
1867 C
1868 LENA = 0
1869 IF (LEN1 .GT. 0) THEN
1870 DO 512 I=1,LEN1
1871 II = IA(I)
1872 WFN(I) = WFN(I) + TX1(I) * GAM(II,IC)*(PRJ(II,IC1)-PRH(II,IC))
1873 * - TX3(I)
1874 IF (WFN(I) .GT. 0.0) THEN
1875 LENA = LENA + 1
1876 I1(LENA) = IA(I)
1877 I2(LENA) = I
1878 TX1(LENA) = WFN(I)
1879 TX2(LENA) = QS1(IA(I))
1880 TX6(LENA) = 1.0
1881 ENDIF
1882 512 CONTINUE
1883 ENDIF
1884 LENB = LENA
1885 DO 515 I=LEN11,LEN2
1886 WFN(I) = WFN(I) - TX3(I)
1887 IF (WFN(I) .GT. 0.0 .AND. TX6(I) .GT. 0.0) THEN
1888 LENB = LENB + 1
1889 I1(LENB) = IA(I)
1890 I2(LENB) = I
1891 TX1(LENB) = WFN(I)
1892 TX2(LENB) = QS1(IA(I))
1893 TX4(LENB) = TX6(I)
1894 ENDIF
1895 515 CONTINUE
1896 C
1897 IF (LENB .LE. 0) THEN
1898 DO 5030 I=1,LENC*K
1899 HST(I,1) = 0.0
1900 QOL(I,1) = 0.0
1901 5030 CONTINUE
1902 DO 5040 I=1,LENC
1903 PCU(I) = 0.0
1904 5040 CONTINUE
1905 RETURN
1906 ENDIF
1907
1908 C
1909 DO 516 I=1,LENB
1910 WFN(I) = TX1(I)
1911 QS1(I) = TX2(I)
1912 516 CONTINUE
1913 C
1914 DO 520 L=IC,K
1915 DO 517 I=1,LENB
1916 TX1(I) = ETA(I2(I),L)
1917 517 CONTINUE
1918 DO 520 I=1,LENB
1919 ETA(I,L) = TX1(I)
1920 520 CONTINUE
1921 C
1922 LENA1 = LENA + 1
1923 C
1924 DO 510 I=1,LENA
1925 II = I1(I)
1926 TX8(I) = HST(II,IC) - HOL(II,IC)
1927 510 CONTINUE
1928 DO 530 I=LENA1,LENB
1929 II = I1(I)
1930 TX6(I) = TX4(I)
1931 TEM = TX6(I) * (HOL(II,IC)-HOL(II,IC1)) + HOL(II,IC1)
1932 TX8(I) = HOL(II,K) - TEM
1933
1934 TEM1 = TX6(I) * (QOL(II,IC)-QOL(II,IC1)) + QOL(II,IC1)
1935 TX5(I) = TEM - TEM1 * ALHL
1936 QS1(I) = TEM1 + TX8(I)*(ONE/ALHL)
1937 TX3(I) = HOL(II,IC)
1938 530 CONTINUE
1939 C
1940 C
1941 DO 620 I=1,LENB
1942 II = I1(I)
1943 WLQ(I) = QOL(II,K) - QS1(I) * ETA(I,IC)
1944 TX7(I) = HOL(II,K)
1945 620 CONTINUE
1946 DO NT=1,Ntracer
1947 DO 621 I=1,LENB
1948 II = I1(I)
1949 UHT(I,NT) = UOI(II,K+nltop-1,NT)-UOI(II,IC+nltop-1,NT) * ETA(I,IC)
1950 621 CONTINUE
1951 ENDDO
1952 C
1953 DO 635 L=KM1,IC,-1
1954 DO 630 I=1,LENB
1955 II = I1(I)
1956 TEM = ETA(I,L) - ETA(I,L+1)
1957 WLQ(I) = WLQ(I) + TEM * QOL(II,L)
1958 630 CONTINUE
1959 635 CONTINUE
1960 DO NT=1,Ntracer
1961 DO L=KM1,IC,-1
1962 DO I=1,LENB
1963 II = I1(I)
1964 TEM = ETA(I,L) - ETA(I,L+1)
1965 UHT(I,NT) = UHT(I,NT) + TEM * UOI(II,L+nltop-1,NT)
1966 ENDDO
1967 ENDDO
1968 ENDDO
1969 C
1970 C CALCULATE GS AND PART OF AKM (THAT REQUIRES ETA)
1971 C
1972 DO 690 I=1,LENB
1973 II = I1(I)
1974 c TX7(I) = HOL(II,K)
1975 TEM = (POI(II,KM1) - POI(II,K)) / (PRH(II,K) - PRH(II,KM1))
1976 HOL(I,K) = TEM * (PRJ(II,K)-PRH(II,KM1))*PRH(II,K)*PRI(II,K)
1977 HOL(I,KM1) = TEM * (PRH(II,K)-PRJ(II,K))*PRH(II,KM1)*PRI(II,KM1)
1978 AKM(I) = ZERO
1979 TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
1980 690 CONTINUE
1981
1982 IF (IC1 .LE. KM1) THEN
1983 DO 750 L=KM1,IC1,-1
1984 DO 750 I=1,LENB
1985 II = I1(I)
1986 TEM = (POI(II,L-1) - POI(II,L)) * ETA(I,L)
1987 * / (PRH(II,L) - PRH(II,L-1))
1988 C
1989 HOL(I,L) = TEM * (PRJ(II,L)-PRH(II,L-1)) * PRH(II,L)
1990 * * PRI(II,L) + HOL(I,L)
1991 HOL(I,L-1) = TEM * (PRH(II,L)-PRJ(II,L)) * PRH(II,L-1)
1992 * * PRI(II,L-1)
1993 C
1994 AKM(I) = AKM(I) - HOL(I,L)
1995 * * (ETA(I,L) * (PRH(II,L)-PRJ(II,L)) +
1996 * ETA(I,L+1) * (PRJ(II,L+1)-PRH(II,L))) / PRH(II,L)
1997 750 CONTINUE
1998 ENDIF
1999 C
2000 C
2001 CALL RNCL(LENB, TX2, TX1, CLF)
2002
2003 DO 770 I=1,LENB
2004 TX2(I) = (ONE - TX1(I)) * WLQ(I)
2005 WLQ(I) = TX1(I) * WLQ(I)
2006 C
2007 TX1(I) = HOL(I,IC)
2008 770 CONTINUE
2009 DO 790 I=LENA1, LENB
2010 II = I1(I)
2011 TX1(I) = TX1(I) + (TX5(I)-TX3(I)+QOL(II,IC)*ALHL)*(PRI(II,IC)/CP)
2012 790 CONTINUE
2013
2014 DO 800 I=1,LENB
2015 HOL(I,IC) = TX1(I) - TX2(I) * ALBCP * PRI(I1(I),IC)
2016 800 CONTINUE
2017
2018 IF (LENA .GT. 0) THEN
2019 DO 810 I=1,LENA
2020 II = I1(I)
2021 AKM(I) = AKM(I) - ETA(I,IC1) * (PRJ(II,IC1) - PRH(II,IC))
2022 * * TX1(I) / PRH(II,IC)
2023 810 CONTINUE
2024 ENDIF
2025 c
2026 C CALCULATE GH
2027 C
2028 DO 830 I=1,LENB
2029 II = I1(I)
2030 TX3(I) = QOL(II,KM1) - QOL(II,K)
2031 GMH(I,K) = HOL(I,K) + TX3(I) * PRI(II,K) * (ALBCP)
2032
2033 AKM(I) = AKM(I) + GAM(II,KM1)*(PRJ(II,K)-PRH(II,KM1))
2034 * * GMH(I,K)
2035 TX3(I) = zero
2036 830 CONTINUE
2037 C
2038 IF (IC1 .LE. KM1) THEN
2039 DO 840 L=KM1,IC1,-1
2040 DO 840 I=1,LENB
2041 II = I1(I)
2042 TX2(I) = TX3(I)
2043 TX3(I) = (QOL(II,L-1) - QOL(II,L)) * ETA(I,L)
2044 TX2(I) = TX2(I) + TX3(I)
2045 C
2046 GMH(I,L) = HOL(I,L) + TX2(I) * PRI(II,L) * (ALBCP*HALF)
2047 840 CONTINUE
2048 C
2049 C
2050 ENDIF
2051 DO 850 I=LENA1,LENB
2052 TX3(I) = TX3(I) + TWOBAL
2053 * * (TX7(I) - TX8(I) - TX5(I) - QOL(I1(I),IC)*ALHL)
2054 850 CONTINUE
2055 DO 860 I=1,LENB
2056 GMH(I,IC) = TX1(I) + PRI(I1(I),IC) * ONEBCP
2057 * * (TX3(I)*(ALHL*HALF) + ETA(I,IC) * TX8(I))
2058 860 CONTINUE
2059 C
2060 C CALCULATE HC PART OF AKM
2061 C
2062 IF (IC1 .LE. KM1) THEN
2063 DO 870 I=1,LENB
2064 TX1(I) = GMH(I,K)
2065 870 CONTINUE
2066 DO 3725 L=KM1,IC1,-1
2067 DO 880 I=1,LENB
2068 II = I1(I)
2069 TX1(I) = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * GMH(I,L)
2070 TX2(I) = GAM(II,L-1) * (PRJ(II,L) - PRH(II,L-1))
2071 880 CONTINUE
2072 C
2073 IF (L .EQ. IC1) THEN
2074 DO 890 I=LENA1,LENB
2075 TX2(I) = ZERO
2076 890 CONTINUE
2077 ENDIF
2078 DO 900 I=1,LENB
2079 II = I1(I)
2080 AKM(I) = AKM(I) + TX1(I) *
2081 * (TX2(I) + GAM(II,L)*(PRH(II,L)-PRJ(II,L)))
2082 900 CONTINUE
2083 3725 CONTINUE
2084 ENDIF
2085 C
2086 DO 920 I=LENA1,LENB
2087 II = I1(I)
2088 TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
2089 * + 0.5*(PRS(II,IC+2) - PRS(II,IC)) * (ONE-TX6(I))
2090 c
2091 TX1(I) = PRS(II,IC1)
2092 TX5(I) = 0.5 * (PRS(II,IC1) + PRS(II,IC+2))
2093 C
2094 IF ((TX2(I) .GE. TX1(I)) .AND. (TX2(I) .LT. TX5(I))) THEN
2095 TX6(I) = ONE - (TX2(I) - TX1(I)) / (TX5(I) - TX1(I))
2096 C
2097 TEM = PRI(II,IC1) / PRI(II,IC)
2098 HOL(I,IC1) = HOL(I,IC1) + HOL(I,IC) * TEM
2099 HOL(I,IC) = ZERO
2100 C
2101 GMH(I,IC1) = GMH(I,IC1) + GMH(I,IC) * TEM
2102 GMH(I,IC) = ZERO
2103 ELSEIF (TX2(I) .LT. TX1(I)) THEN
2104 TX6(I) = 1.0
2105 ELSE
2106 TX6(I) = 0.0
2107 ENDIF
2108 920 CONTINUE
2109 C
2110 C
2111 DO I=1,LENC
2112 PCU(I) = 0.0
2113 ENDDO
2114
2115 DO 970 I=1,LENB
2116 II = I1(I)
2117 IF (AKM(I) .LT. ZERO .AND. WLQ(I) .GE. 0.0) THEN
2118 WFN(I) = - TX6(I) * WFN(I) * RASALF / AKM(I)
2119 ELSE
2120 WFN(I) = ZERO
2121 ENDIF
2122 TEM = (PRS(II,K+1)-PRS(II,K))*(CMB2PA*FRAC)
2123 WFN(I) = MIN(WFN(I), TEM)
2124 C
2125 C compute cloud amount
2126 C
2127 CC TX1(I) = CLN(II)
2128 CC IF (WFN(I) .GT. CRTMSF) TX1(I) = TX1(I) + CLF(I)
2129 CC IF (TX1(I) .GT. ONE) TX1(I) = ONE
2130 C
2131 C PRECIPITATION
2132 C
2133 PCU(II) = WLQ(I) * WFN(I) * ONEBG
2134 C
2135 C CUMULUS FRICTION AT THE BOTTOM LAYER
2136 C
2137 TX4(I) = WFN(I) * (1.0/ALHL)
2138 TX5(I) = WFN(I) * ONEBCP
2139 970 CONTINUE
2140 C
2141 C compute cloud mass flux for diagnostic output
2142 C
2143 DO L = IC,K
2144 DO I=1,LENB
2145 II = I1(I)
2146 if(L.lt.K)then
2147 CMASS(II,L) = ETA(I,L+1) * WFN(I) * ONEBG
2148 else
2149 CMASS(II,L) = WFN(I) * ONEBG
2150 endif
2151 ENDDO
2152 ENDDO
2153 C
2154 CC DO 975 I=1,LENB
2155 CC II = I1(I)
2156 CC CLN(II) = TX1(I)
2157 CC975 CONTINUE
2158 C
2159 C THETA AND Q CHANGE DUE TO CLOUD TYPE IC
2160 C
2161
2162 c TEMA = 0.0
2163 c TEMB = 0.0
2164 DO 990 L=IC,K
2165 DO 980 I=1,LENB
2166 II = I1(I)
2167 TEM = (GMH(I,L) - HOL(I,L)) * TX4(I)
2168 TEM1 = HOL(I,L) * TX5(I)
2169 C
2170 TCU(II,L) = TEM1 / PRH(II,L)
2171 QCU(II,L) = TEM
2172 980 CONTINUE
2173
2174 c I = I1(IP1)
2175 c
2176 c TEM = (PRS(I,L+1)-PRS(I,L)) * (ONEBG*100.0)
2177 c TEMA = TEMA + TCU(I,L) * PRH(I,L) * TEM * (CP/ALHL)
2178 c TEMB = TEMB + QCU(I,L) * TEM
2179 C
2180 990 CONTINUE
2181 C
2182 c Compute Tracer Tendencies
2183 c -------------------------
2184 do nt = 1,ntracer
2185 c
2186 c Tracer Tendency at the Bottom Layer
2187 c -----------------------------------
2188 DO 995 I=1,LENB
2189 II = I1(I)
2190 TEM = half*TX5(I) * PRI(II,K)
2191 TX1(I) = ( UOI(II,KM1+nltop-1,nt) - UOI(II,K+nltop-1,nt))
2192 ucu(II,K,nt) = TEM * TX1(I)
2193 995 CONTINUE
2194 c
2195 c Tracer Tendency at all other Levels
2196 c -----------------------------------
2197 DO 1020 L=KM1,IC1,-1
2198 DO 1010 I=1,LENB
2199 II = I1(I)
2200 TEM = half*TX5(I) * PRI(II,L)
2201 TEM1 = TX1(I)
2202 TX1(I) = (UOI(II,L-1+nltop-1,nt)-UOI(II,L+nltop-1,nt)) * ETA(I,L)
2203 TX3(I) = (TX1(I) + TEM1) * TEM
2204 1010 CONTINUE
2205 DO 1020 I=1,LENB
2206 II = I1(I)
2207 ucu(II,L,nt) = TX3(I)
2208 1020 CONTINUE
2209
2210 DO 1030 I=1,LENB
2211 II = I1(I)
2212 IF (TX6(I) .GE. 1.0) THEN
2213 TEM = half*TX5(I) * PRI(II,IC)
2214 ELSE
2215 TEM = 0.0
2216 ENDIF
2217 TX1(I) = (TX1(I) + UHT(I,nt) + UHT(I,nt)) * TEM
2218 1030 CONTINUE
2219 DO 1040 I=1,LENB
2220 II = I1(I)
2221 ucu(II,IC,nt) = TX1(I)
2222 1040 CONTINUE
2223
2224 enddo
2225 C
2226 C PENETRATIVE CONVECTION CALCULATION OVER
2227 C
2228
2229 RETURN
2230 END
2231 SUBROUTINE RNCL(lng, PL, RNO, CLF)
2232 C
2233 C*********************************************************************
2234 C********************** Relaxed Arakawa-Schubert *********************
2235 C************************ SUBROUTINE RNCL ************************
2236 C**************************** 23 July 1992 ***************************
2237 C*********************************************************************
2238 implicit none
2239 C Argument List declarations
2240 integer lng
2241 _RL PL(lng), RNO(lng), CLF(lng)
2242
2243 C Local Variables
2244 _RL p5,p8,pt8,pt2,pfac,p4,p6,p7,p9,cucld,cfac
2245 PARAMETER (P5=500.0, P8=800.0, PT8=0.8, PT2=0.2)
2246 PARAMETER (PFAC=PT2/(P8-P5))
2247 PARAMETER (P4=400.0, P6=401.0)
2248 PARAMETER (P7=700.0, P9=900.0)
2249 PARAMETER (CUCLD=0.5,CFAC=CUCLD/(P6-P4))
2250
2251 integer i
2252 C
2253 DO 10 I=1,lng
2254 rno(i) = 1.0
2255 ccc if( pl(i).le.400.0 ) rno(i) = max( 0.75, 1.0-0.0025*(400.0-pl(i)) )
2256
2257 ccc IF ( PL(I).GE.P7 .AND. PL(I).LE.P9 ) THEN
2258 ccc RNO(I) = ((P9-PL(I))/(P9-P7)) **2
2259 ccc ELSE IF (PL(I).GT.P9) THEN
2260 ccc RNO(I) = 0.
2261 ccc ENDIF
2262
2263 CLF(I) = CUCLD
2264 C
2265 CARIESIF (PL(I) .GE. P5 .AND. PL(I) .LE. P8) THEN
2266 CARIES RNO(I) = (P8-PL(I))*PFAC + PT8
2267 CARIESELSEIF (PL(I) .GT. P8 ) THEN
2268 CARIES RNO(I) = PT8
2269 CARIESENDIF
2270 CARIES
2271 IF (PL(I) .GE. P4 .AND. PL(I) .LE. P6) THEN
2272 CLF(I) = (P6-PL(I))*CFAC
2273 ELSEIF (PL(I) .GT. P6 ) THEN
2274 CLF(I) = 0.0
2275 ENDIF
2276 10 CONTINUE
2277 C
2278 RETURN
2279 END
2280 SUBROUTINE ACRITN ( lng,PL,PLB,ACR )
2281
2282 C*********************************************************************
2283 C********************** Relaxed Arakawa-Schubert *********************
2284 C************************** SUBROUTINE ACRIT *********************
2285 C****************** modified August 28, 1996 L.Takacs ************
2286 C**** *****
2287 C**** Note: Data obtained from January Mean After-Analysis *****
2288 C**** from 4x5 46-layer GEOS Assimilation *****
2289 C**** *****
2290 C*********************************************************************
2291 implicit none
2292 C Argument List declarations
2293 integer lng
2294 _RL PL(lng), PLB(lng), ACR(lng)
2295
2296 C Local variables
2297 integer lma
2298 parameter (lma=18)
2299 _RL p(lma)
2300 _RL a(lma)
2301 integer i,L
2302 _RL temp
2303
2304 data p / 93.81, 111.65, 133.46, 157.80, 186.51,
2305 . 219.88, 257.40, 301.21, 352.49, 409.76,
2306 . 471.59, 535.04, 603.33, 672.79, 741.12,
2307 . 812.52, 875.31, 930.20/
2308
2309 data a / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364,
2310 . 1.01971, .65846, .45867, .38687, .31002,
2311 . .25574, .20347, .17254, .15260, .16756,
2312 . .09916, .10360, .05880/
2313
2314
2315 do L=1,lma-1
2316 do i=1,lng
2317 if( pl(i).ge.p(L) .and.
2318 . pl(i).le.p(L+1)) then
2319 temp = ( pl(i)-p(L) )/( p(L+1)-p(L) )
2320 acr(i) = a(L+1)*temp + a(L)*(1-temp)
2321 endif
2322 enddo
2323 enddo
2324
2325 do i=1,lng
2326 if( pl(i).lt.p(1) ) acr(i) = a(1)
2327 if( pl(i).gt.p(lma) ) acr(i) = a(lma)
2328 enddo
2329
2330 do i=1,lng
2331 acr(i) = acr(i) * (plb(i)-pl(i))
2332 enddo
2333
2334 RETURN
2335 END
2336 SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DP,PLKE,
2337 1 PLK,TH,TEMP1,TEMP2,TEMP3,ITMP1,ITMP2,RCON,RLAR,CLSBTH,tmscl,
2338 2 tmfrc,cp,gravity,alhl,gamfac,cldlz,RHCRIT,offset,alpha)
2339
2340 implicit none
2341 C Argument List declarations
2342 integer nn,irun,nlay
2343 _RL TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY),
2344 . PL(IRUN,NLAY),CLFRAC(IRUN,NLAY),SP(IRUN),TEMP1(IRUN,NLAY),
2345 . TEMP2(IRUN,NLAY),PLKE(IRUN,NLAY+1),
2346 . RCON(IRUN),RLAR(IRUN),DP(IRUN,NLAY),PLK(IRUN,NLAY),TH(IRUN,NLAY),
2347 . TEMP3(IRUN,NLAY)
2348 integer ITMP1(IRUN,NLAY),ITMP2(IRUN,NLAY)
2349 _RL CLSBTH(IRUN,NLAY)
2350 _RL tmscl,tmfrc,cp,gravity,alhl,gamfac,offset,alpha
2351 _RL cldlz(irun,nlay)
2352 _RL rhcrit(irun,nlay)
2353 C
2354 C Local Variables
2355 _RL zm1p04,zero,two89,zp44,zp01,half,zp578,one,thousand,z3600
2356 _RL zp1,zp001
2357 PARAMETER (ZM1P04 = -1.04E-4 )
2358 PARAMETER (ZERO = 0.)
2359 PARAMETER (TWO89= 2.89E-5)
2360 PARAMETER ( ZP44= 0.44)
2361 PARAMETER ( ZP01= 0.01)
2362 PARAMETER ( ZP1 = 0.1 )
2363 PARAMETER ( ZP001= 0.001)
2364 PARAMETER ( HALF= 0.5)
2365 PARAMETER ( ZP578 = 0.578 )
2366 PARAMETER ( ONE = 1.)
2367 PARAMETER ( THOUSAND = 1000.)
2368 PARAMETER ( Z3600 = 3600.)
2369 C
2370 _RL EVP9(IRUN,NLAY)
2371 _RL water(irun),crystal(irun)
2372 _RL watevap(irun),iceevap(irun)
2373 _RL fracwat,fracice, tice,rh,fact,dum
2374 _RL rainmax(irun)
2375 _RL getcon,rphf,elocp,cpog,relax
2376 _RL exparg,arearat,rpow
2377
2378 integer i,L,n,nlaym1,irnlay,irnlm1
2379
2380 c Explicit Inline Directives
2381 c --------------------------
2382 #ifdef CRAY
2383 #ifdef f77
2384 cfpp$ expand (qsat)
2385 #endif
2386 #endif
2387
2388 tice = getcon('FREEZING-POINT')
2389
2390 fracwat = 0.70
2391 fracice = 0.01
2392
2393 NLAYM1 = NLAY - 1
2394 IRNLAY = IRUN*NLAY
2395 IRNLM1 = IRUN*(NLAY-1)
2396
2397 RPHF = Z3600/tmscl
2398
2399 ELOCP = alhl/cp
2400 CPOG = cp/gravity
2401
2402 DO I = 1,IRUN
2403 RLAR(I) = 0.
2404 water(i) = 0.
2405 crystal(i) = 0.
2406 ENDDO
2407
2408 do L = 1,nlay
2409 do i = 1,irun
2410 EVP9(i,L) = 0.
2411 TEMP1(i,L) = 0.
2412 TEMP2(i,L) = 0.
2413 TEMP3(i,L) = 0.
2414 CLSBTH(i,L) = 0.
2415 cldlz(i,L) = 0.
2416 enddo
2417 enddo
2418
2419 C RHO(ZERO) / RHO FOR TERMINAL VELOCITY APPROX.
2420 c ---------------------------------------------
2421 DO L = 1,NLAY
2422 DO I = 1,IRUN
2423 TEMP2(I,L) = PL(I,L)*ZP001
2424 TEMP2(I,L) = SQRT(TEMP2(I,L))
2425 ENDDO
2426 ENDDO
2427
2428 C INVERSE OF MASS IN EACH LAYER
2429 c -----------------------------
2430 DO L = 1,NLAY
2431 DO I = 1,IRUN
2432 TEMP3(I,L) = GRAVITY*ZP01 / DP(I,L)
2433 ENDDO
2434 ENDDO
2435
2436 C DO LOOP FOR MOISTURE EVAPORATION ABILITY AND CONVEC EVAPORATION.
2437 c ----------------------------------------------------------------
2438 DO 100 L=1,NLAY
2439
2440 DO I = 1,IRUN
2441 TEMP1(I,3) = TL(I,L)
2442 TEMP1(I,4) = QL(I,L)
2443 ENDDO
2444
2445 DO 50 N=1,2
2446 IF(N.EQ.1)RELAX=HALF
2447 IF(N.GT.1)RELAX=ONE
2448
2449 DO I = 1,IRUN
2450 call qsat ( temp1(i,3),pl(i,L),temp1(i,2),temp1(i,6),.true. )
2451 TEMP1(I,5)=TEMP1(I,2)-TEMP1(I,4)
2452 TEMP1(I,6)=TEMP1(I,6)*ELOCP
2453 TEMP1(I,5)=TEMP1(I,5)/(ONE+TEMP1(I,6))
2454 TEMP1(I,4)=TEMP1(I,4)+TEMP1(I,5)*RELAX
2455 TEMP1(I,3)=TEMP1(I,3)-TEMP1(I,5)*ELOCP*RELAX
2456 ENDDO
2457 50 CONTINUE
2458
2459 DO I = 1,IRUN
2460 EVP9(I,L) = (TEMP1(I,4) - QL(I,L))/TEMP3(I,L)
2461 C convective detrained water
2462 cldlz(i,L) = rain(i,L)*temp3(i,L)
2463 if( tl(i,L).gt.tice-20.) then
2464 water(i) = water(i) + rain(i,L)
2465 else
2466 crystal(i) = crystal(i) + rain(i,L)
2467 endif
2468 ENDDO
2469
2470 C**********************************************************************
2471 C FOR CONVECTIVE PRECIP, FIND THE "EVAPORATION EFFICIENCY" USING *
2472 C KESSLERS PARAMETERIZATION *
2473 C**********************************************************************
2474
2475 DO 20 I=1,IRUN
2476
2477 iceevap(i) = 0.
2478 watevap(i) = 0.
2479
2480 if( (evp9(i,L).gt.0.) .and. (crystal(i).gt.0.) ) then
2481 iceevap(I) = EVP9(I,L)*fracice
2482 IF(iceevap(i).GE.crystal(i)) iceevap(i) = crystal(i)
2483 EVP9(I,L)=EVP9(I,L)-iceevap(I)
2484 crystal(i) = crystal(i) - iceevap(i)
2485 endif
2486
2487 C and now warm precipitate
2488 if( (evp9(i,L).gt.0.) .and. (water(i).gt.0.) ) then
2489 exparg = ZM1P04*tmscl*((water(i)*RPHF*TEMP2(I,L))**ZP578)
2490 AREARAT = ONE-(EXP(EXPARG))
2491 watevap(I) = EVP9(I,L)*AREARAT*fracwat
2492 IF(watevap(I).GE.water(i)) watevap(I) = water(i)
2493 EVP9(I,L)=EVP9(I,L)-watevap(I)
2494 water(i) = water(i) - watevap(i)
2495 endif
2496
2497 QL(I,L) = QL(I,L)+(iceevap(i)+watevap(i))*TEMP3(I,L)
2498 TL(I,L) = TL(I,L)-(iceevap(i)+watevap(i))*TEMP3(I,L)*ELOCP
2499
2500 20 CONTINUE
2501
2502 100 CONTINUE
2503
2504 do i = 1,irun
2505 rcon(i) = water(i) + crystal(i)
2506 enddo
2507
2508 C**********************************************************************
2509 C Large Scale Precip
2510 C**********************************************************************
2511
2512 DO 200 L=1,NLAY
2513 DO I = 1,IRUN
2514 rainmax(i) = rhcrit(i,L)*evp9(i,L) +
2515 . ql(i,L)*(rhcrit(i,L)-1.)/temp3(i,L)
2516
2517 if (rainmax(i).LE.0.0) then
2518 call qsat( tl(i,L),pl(i,L),rh,dum,.false.)
2519 rh = ql(i,L)/rh
2520
2521 if( rhcrit(i,L).eq.1.0 ) then
2522 fact = 1.0
2523 else
2524 fact = min( 1.0, alpha + (1.0-alpha)*( rh-rhcrit(i,L)) /
2525 1 (1.0-rhcrit(i,L)) )
2526 endif
2527
2528 C Do not allow clouds above 10 mb
2529 if( pl(i,L).ge.10.0 ) CLSBTH(I,L) = fact
2530 RLAR(I) = RLAR(I)-rainmax(I)
2531 QL(I,L) = QL(I,L)+rainmax(I)*TEMP3(I,L)
2532 TL(I,L) = TL(I,L)-rainmax(I)*TEMP3(I,L)*ELOCP
2533 C Large-scale water
2534 cldlz(i,L) = cldlz(i,L) - rainmax(i)*temp3(i,L)
2535 ENDIF
2536 ENDDO
2537
2538 DO I=1,IRUN
2539 IF((RLAR(I).GT.0.0).AND.(rainmax(I).GT.0.0))THEN
2540 RPOW=(RLAR(I)*RPHF*TEMP2(I,L))**ZP578
2541 EXPARG = ZM1P04*tmscl*RPOW
2542 AREARAT = ONE-(EXP(EXPARG))
2543 TEMP1(I,7) = rainmax(I)*AREARAT
2544 IF(TEMP1(I,7).GE.RLAR(I)) TEMP1(I,7) = RLAR(I)
2545 RLAR(I) = RLAR(I)-TEMP1(I,7)
2546 QL(I,L) = QL(I,L)+TEMP1(I,7)*TEMP3(I,L)
2547 TL(I,L) = TL(I,L)-TEMP1(I,7)*TEMP3(I,L)*ELOCP
2548 ENDIF
2549 ENDDO
2550
2551 200 CONTINUE
2552
2553 RETURN
2554 END
2555 subroutine srclouds (th,q,plk,pl,plke,cloud,cldwat,irun,irise,
2556 1 rhc,offset,alpha)
2557 C***********************************************************************
2558 C
2559 C PURPOSE:
2560 C ========
2561 C Compute non-precipitating cloud fractions
2562 C based on Slingo and Ritter (1985).
2563 C Remove cloudiness where conditionally unstable.
2564 C
2565 C INPUT:
2566 C ======
2567 C th ......... Potential Temperature (irun,irise)
2568 C q .......... Specific Humidity (irun,irise)
2569 C plk ........ P**Kappa at mid-layer (irun,irise)
2570 C pl ......... Pressure at mid-layer (irun,irise)
2571 C plke ....... P**Kappa at edge (irun,irise+1)
2572 C irun ....... Horizontal dimension
2573 C irise ...... Vertical dimension
2574 C
2575 C OUTPUT:
2576 C =======
2577 C cloud ...... Cloud Fraction (irun,irise)
2578 C
2579 C***********************************************************************
2580
2581 implicit none
2582 integer irun,irise
2583
2584 _RL th(irun,irise)
2585 _RL q(irun,irise)
2586 _RL plk(irun,irise)
2587 _RL pl(irun,irise)
2588 _RL plke(irun,irise+1)
2589
2590 _RL cloud(irun,irise)
2591 _RL cldwat(irun,irise)
2592 _RL qs(irun,irise)
2593
2594 _RL cp, alhl, getcon, akap
2595 _RL ratio, temp, elocp
2596 _RL rhcrit,rh,dum
2597 integer i,L
2598
2599 _RL rhc(irun,irise)
2600 _RL offset,alpha
2601
2602 c Explicit Inline Directives
2603 c --------------------------
2604 #ifdef CRAY
2605 #ifdef f77
2606 cfpp$ expand (qsat)
2607 #endif
2608 #endif
2609
2610 cp = getcon('CP')
2611 alhl = getcon('LATENT HEAT COND')
2612 elocp = alhl/cp
2613 akap = getcon('KAPPA')
2614
2615 do L = 1,irise
2616 do i = 1,irun
2617 temp = th(i,L)*plk(i,L)
2618 call qsat ( temp,pl(i,L),qs(i,L),dum,.false. )
2619 enddo
2620 enddo
2621
2622 do L = 2,irise
2623 do i = 1,irun
2624 rh = q(i,L)/qs(i,L)
2625
2626 rhcrit = rhc(i,L) - offset
2627 ratio = alpha*(rh-rhcrit)/offset
2628
2629 if(cloud(i,L).eq. 0.0 .and. ratio.gt.0.0 ) then
2630 cloud(i,L) = min( ratio,1.0 )
2631 endif
2632
2633 enddo
2634 enddo
2635
2636 c Reduce clouds from conditionally unstable layer
2637 c -----------------------------------------------
2638 call ctei ( th,q,cloud,cldwat,pl,plk,plke,irun,irise )
2639
2640 return
2641 end
2642
2643 subroutine ctei ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm )
2644 implicit none
2645 integer im,lm
2646 _RL th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm)
2647 _RL plk(im,lm),pl(im,lm),cldfrc(im,lm)
2648 integer i,L
2649 _RL getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq
2650 _RL k,krd,kmm,f
2651
2652 cp = getcon('CP')
2653 alhl = getcon('LATENT HEAT COND')
2654 cpoel = cp/alhl
2655 elocp = alhl/cp
2656
2657 do L=lm,2,-1
2658 do i=1,im
2659 dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1)
2660 if( dq.eq.0.0 ) dq = 1.0e-20
2661 k = 1.0 + cpoel*plke(i,L)*( th(i,L)-th(i,L-1) ) / dq
2662
2663 t = th(i,L)*plk(i,L)
2664 p = pl(i,L)
2665 call qsat ( t,p,qs,dqsdt,.true. )
2666
2667 krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t )
2668
2669 kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t )
2670 . / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt )
2671
2672 s = ( (k-krd)/(kmm-krd) )
2673 f = 1.0 - min( 1.0, max(0.0,1.0-exp(-s)) )
2674
2675 cldfrc(i,L) = cldfrc(i,L)*f
2676 cldwat(i,L) = cldwat(i,L)*f
2677
2678 enddo
2679 enddo
2680
2681 return
2682 end
2683
2684 subroutine back2grd(gathered,indeces,scattered,irun)
2685 implicit none
2686 integer i,irun,indeces(irun)
2687 _RL gathered(irun),scattered(irun)
2688 _RL temp(irun)
2689 do i = 1,irun
2690 temp(indeces(i)) = gathered(i)
2691 enddo
2692 do i = 1,irun
2693 scattered(i) = temp(i)
2694 enddo
2695 return
2696 end
  ViewVC Help
Powered by ViewVC 1.1.22