Radiation Scheme in CCPP
radiation_clouds.f
Go to the documentation of this file.
1 !!!!! module_radiation_clouds description !!!!!
2 !!!!! ========================================================== !!!!!
3 ! !
4 ! the 'radiation_clouds.f' contains: !
5 ! !
6 ! 'module_radiation_clouds' --- compute cloud related quantities!
7 ! for radiation computations !
8 ! !
9 ! the following are the externally accessable subroutines: !
10 ! !
11 ! 'cld_init' --- initialization routine !
12 ! inputs: !
13 ! ( si, NLAY, me ) !
14 ! outputs: !
15 ! ( none ) !
16 ! !
17 ! 'progcld1' --- zhao/moorthi prognostic cloud scheme !
18 ! inputs: !
19 ! (plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, !
20 ! xlat,xlon,slmsk, !
21 ! IX, NLAY, NLP1, !
22 ! outputs: !
23 ! clouds,clds,mtop,mbot) !
24 ! !
25 ! 'progcld2' --- ferrier prognostic cloud microphysics !
26 ! inputs: !
27 ! (plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, !
28 ! xlat,xlon,slmsk, f_ice,f_rain,r_rime,flgmin, !
29 ! IX, NLAY, NLP1, !
30 ! outputs: !
31 ! clouds,clds,mtop,mbot) !
32 ! !
33 ! 'progcld3' --- zhao/moorthi prognostic cloud + pdfcld!
34 ! inputs: !
35 ! (plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, cnvw,cnvc, !
36 ! xlat,xlon,slmsk, !
37 ! ix, nlay, nlp1, !
38 ! deltaq,sup,kdt,me, !
39 ! outputs: !
40 ! clouds,clds,mtop,mbot)
41 ! !
42 ! 'diagcld1' --- diagnostic cloud calc routine !
43 ! inputs: !
44 ! (plyr,plvl,tlyr,rhly,vvel,cv,cvt,cvb, !
45 ! xlat,xlon,slmsk, !
46 ! IX, NLAY, NLP1, !
47 ! outputs: !
48 ! clouds,clds,mtop,mbot) !
49 ! !
50 ! internal accessable only subroutines: !
51 ! 'gethml' --- get diagnostic hi, mid, low clouds !
52 ! !
53 ! 'rhtable' --- rh lookup table for diag cloud scheme !
54 ! !
55 ! !
56 ! cloud array description: !
57 ! --- for prognostic cloud: icldflg=1 --- !
58 ! clouds(:,:,1) - layer total cloud fraction !
59 ! clouds(:,:,2) - layer cloud liq water path !
60 ! clouds(:,:,3) - mean effective radius for liquid cloud !
61 ! clouds(:,:,4) - layer cloud ice water path !
62 ! clouds(:,:,5) - mean effective radius for ice cloud !
63 ! clouds(:,:,6) - layer rain drop water path !
64 ! clouds(:,:,7) - mean effective radius for rain drop !
65 ! ** clouds(:,:,8) - layer snow flake water path !
66 ! clouds(:,:,9) - mean effective radius for snow flake !
67 ! ** fu's scheme need to be normalized by snow density (g/m**3/1.0e6)!
68 ! --- for diagnostic cloud: icldflg=0 --- !
69 ! clouds(:,:,1) - layer total cloud fraction !
70 ! clouds(:,:,2) - layer cloud optical depth !
71 ! clouds(:,:,3) - layer cloud single scattering albedo !
72 ! clouds(:,:,4) - layer cloud asymmetry factor !
73 ! !
74 ! external modules referenced: !
75 ! !
76 ! 'module physparam' in 'physparam.f' !
77 ! 'module physcons' in 'physcons.f' !
78 ! 'module module_microphysics' in 'module_bfmicrophysics.f' !
79 ! !
80 ! !
81 ! program history log: !
82 ! nov 1992, y.h., k.a.c, a.k. - cloud parameterization !
83 ! 'cldjms' patterned after slingo and slingo's work (jgr, !
84 ! 1992), stratiform clouds are allowed in any layer except !
85 ! the surface and upper stratosphere. the relative humidity !
86 ! criterion may cery in different model layers. !
87 ! mar 1993, kenneth campana - created original crhtab tables !
88 ! for critical rh look up references.
89 ! feb 1994, kenneth campana - modified to use only one table !
90 ! for all forecast hours. !
91 ! oct 1995, kenneth campana - tuned cloud rh curves !
92 ! rh-cld relation from tables created using mitchell-hahn !
93 ! tuning technique on airforce rtneph observations. !
94 ! nov 1995, kenneth campana - the bl relationships used !
95 ! below llyr, except in marine stratus regions. !
96 ! apr 1996, kenneth campana - save bl cld amt in cld(,5) !
97 ! aug 1997, kenneth campana - smooth out last bunch of bins !
98 ! of the rh look up tables !
99 ! dec 1998, s. moorthi - added prognostic cloud method !
100 ! apr 2003, yu-tai hou - rewritten in frotran 90 !
101 ! modulized form 'module_rad_clouds' from combining the original!
102 ! subroutines 'cldjms', 'cldprp', and 'gcljms'. and seperated !
103 ! prognostic and diagnostic methods into two packages. !
104 ! --- 2003, s. moorthi - adapted b. ferrier's prognostic !
105 ! cloud scheme to ncep gfs model as additional option. !
106 ! apr 2004, yu-tai hou - separated calculation of the !
107 ! averaged h,m,l,bl cloud amounts from each of the cld schemes !
108 ! to become an shared individule subprogram 'gethml'. !
109 ! may 2004, yu-tai hou - rewritten ferrier's scheme as a !
110 ! separated program 'progcld2' in the cloud module. !
111 ! apr 2005, yu-tai hou - modified cloud array and module !
112 ! structures. !
113 ! dec 2008, yu-tai hou - changed low-cld calculation, !
114 ! now cantains clds from sfc layer and upward to the low/mid !
115 ! boundary (include bl-cld). h,m,l clds domain boundaries are !
116 ! adjusted for better agreement with observations. !
117 ! jan 2011, yu-tai hou - changed virtual temperature !
118 ! as input variable instead of originally computed inside the !
119 ! two prognostic cld schemes 'progcld1' and 'progcld2'. !
120 ! aug 2012, yu-tai hou - modified subroutine cld_init !
121 ! to pass all fixed control variables at the start. and set !
122 ! their correponding internal module variables to be used by !
123 ! module subroutines. !
124 ! nov 2012, yu-tai hou - modified control parameters !
125 ! thru module 'physparam'. !
126 ! apr 2013, h.lin/y.hou - corrected error in calculating !
127 ! llyr for different vertical indexing directions. !
128 ! jul 2013, r. sun/h. pan - modified to use pdf cloud and !
129 ! convective cloud cover and water for radiation !
130 ! !
131 ! jul 2014 s. moorthi - merging with gfs version !
132 ! !
133 !!!!! ========================================================== !!!!!
134 !!!!! end descriptions !!!!!
135 !!!!! ========================================================== !!!!!
136 
138 !========================================!
139  module module_radiation_clouds !
140 !........................................!
141 !
142  use physparam, only : icldflg, icmphys, iovrsw, iovrlw, &
143  & lcrick, lcnorm, lnoprec, lsashal, &
144  & ivflip, kind_phys, kind_io4
145  use physcons, only : con_fvirt, con_ttp, con_rocp, &
146  & con_t0c, con_pi, con_g, con_rd, &
147  & con_thgni
148  use module_microphysics, only : rsipath2
149  use module_iounitdef, only : nicltun
150 !
151  implicit none
152 !
153  private
154 
155 ! --- version tag and last revision date
156  character(40), parameter :: &
157  & VTAGCLD='NCEP-Radiation_clouds v5.1 Nov 2012 '
158 ! & VTAGCLD='NCEP-Radiation_clouds v5.0 Aug 2012 '
159 
160 ! --- set constant parameters
161  real (kind=kind_phys), parameter :: gfac=1.0e5/con_g &
162  &, gord=con_g/con_rd
164  integer, parameter, public :: nf_clds = 9 ! number of fields in cloud array
166  integer, parameter, public :: nk_clds = 3 ! number of cloud vertical domains
167 
168 ! --- pressure limits of cloud domain interfaces (low,mid,high) in mb (0.1kPa)
169  real (kind=kind_phys), save :: ptopc(nk_clds+1,2)
170 
171 !org data ptopc / 1050., 642., 350., 0.0, 1050., 750., 500., 0.0 /
172  data ptopc / 1050., 650., 400., 0.0, 1050., 750., 500., 0.0 /
173 
174 ! real (kind=kind_phys), parameter :: climit = 0.01
175  real (kind=kind_phys), parameter :: climit = 0.001, climit2=0.05
176  real (kind=kind_phys), parameter :: ovcst = 1.0 - 1.0e-8
177 
178 ! --- set default quantities as parameters (for prognostic cloud)
179  real (kind=kind_phys), parameter :: reliq_def = 10.0 ! default liq radius to 10 micron
180  real (kind=kind_phys), parameter :: reice_def = 50.0 ! default ice radius to 50 micron
181  real (kind=kind_phys), parameter :: rrain_def = 1000.0 ! default rain radius to 1000 micron
182  real (kind=kind_phys), parameter :: rsnow_def = 250.0 ! default snow radius to 250 micron
183 
184 ! --- set look-up table dimensions and other parameters (for diagnostic cloud)
185  integer, parameter :: nbin=100 ! rh in one percent interval
186  integer, parameter :: nlon=2 ! =1,2 for eastern and western hemispheres
187  integer, parameter :: nlat=4 ! =1,4 for 60n-30n,30n-equ,equ-30s,30s-60s
188  integer, parameter :: mcld=4 ! =1,4 for bl,low,mid,hi cld type
189  integer, parameter :: nseal=2 ! =1,2 for land,sea
190 
191  real (kind=kind_phys), parameter :: cldssa_def = 0.99 ! default cld single scat albedo
192  real (kind=kind_phys), parameter :: cldasy_def = 0.84 ! default cld asymmetry factor
193 
194 ! --- xlabdy: lat bndry between tuning regions, +/- xlim for transition
195 ! xlobdy: lon bndry between tuning regions
196  real (kind=kind_phys), parameter :: xlim=5.0
197  real (kind=kind_phys) :: xlabdy(3), xlobdy(3)
198 
199  data xlabdy / 30.0, 0.0, -30.0 /, xlobdy / 0.0, 180., 360. /
200 
201 ! --- low cloud vertical velocity adjustment boundaries in mb/sec
202  real (kind=kind_phys), parameter :: vvcld1= 0.0003e0
203  real (kind=kind_phys), parameter :: vvcld2=-0.0005e0
204 
205 ! --- those data will be set up by "cld_init"
206 ! rhcl : tuned rh relation table for diagnostic cloud scheme
207 
208  real (kind=kind_phys) :: rhcl(nbin,nlon,nlat,mcld,nseal)
209 
210  integer :: llyr = 2 ! upper limit of boundary layer clouds
211  integer :: iovr = 1 ! cloud over lapping method for diagnostic 3-domain
212  ! output calc (see iovrsw/iovrlw description)
213 
214  public progcld1, progcld2, progcld3, diagcld1, cld_init
215 
216 
217 ! =================
218  contains
219 ! =================
220 
228 !-----------------------------------
229  subroutine cld_init &
230 !...................................
231 
232 ! --- inputs:
233  & ( si, nlay, me )
234 ! --- outputs:
235 ! ( none )
236 
237 ! =================================================================== !
238 ! !
239 ! abstract: cld_init is an initialization program for cloud-radiation !
240 ! calculations. it sets up boundary layer cloud top. !
241 ! !
242 ! !
243 ! inputs: !
244 ! si (L+1) : model vertical sigma layer interface !
245 ! NLAY : vertical layer number !
246 ! me : print control flag !
247 ! !
248 ! outputs: (none) !
249 ! to module variables !
250 ! !
251 ! external module variables: (in physparam) !
252 ! icldflg : cloud optical property scheme control flag !
253 ! =0: model use diagnostic cloud method !
254 ! =1: model use prognostic cloud method !
255 ! icmphys : cloud microphysics scheme control flag !
256 ! =1: zhao/carr/sundqvist microphysics cloud !
257 ! =2: ferrier microphysics cloud scheme !
258 ! =3: zhao/carr/sundqvist microphysics cloud +pdfcld!
259 ! iovrsw/iovrlw : sw/lw control flag for cloud overlapping scheme !
260 ! =0: random overlapping clouds !
261 ! =1: max/ran overlapping clouds !
262 ! ivflip : control flag for direction of vertical index !
263 ! =0: index from toa to surface !
264 ! =1: index from surface to toa !
265 ! usage: call cld_init !
266 ! !
267 ! subroutines called: rhtable !
268 ! !
269 ! =================================================================== !
270 !
271  implicit none
272 
273 ! --- inputs:
274  integer, intent(in) :: NLAY, me
275 
276  real (kind=kind_phys), intent(in) :: si(:)
277 
278 ! --- outputs: (none)
279 
280 ! --- locals:
281  integer :: k, kl, ier
282 
283 !
284 !===> ... begin here
285 !
286 ! --- set up module variables
287 
288  iovr = max( iovrsw, iovrlw ) !cld ovlp used for diag HML cld output
289 
290  if (me == 0) print *, vtagcld !print out version tag
291 
292  if ( icldflg == 0 ) then
293  if (me == 0) print *,' - Using Diagnostic Cloud Method'
294 
297 ! --- set up tuned rh table
298 
299  call rhtable( me, ier )
300 
301  if (ier < 0) then
302  write(6,99) ier
303  99 format(3x,' *** Error in finding tuned RH table ***' &
304  &, /3x,' STOP at calling subroutine RHTABLE !!'/)
305  stop 99
306  endif
307  else
308  if (me == 0) then
309  print *,' - Using Prognostic Cloud Method'
310  if (icmphys == 1) then
311  print *,' --- Zhao/Carr/Sundqvist microphysics'
312  elseif (icmphys == 2) then
313  print *,' --- Ferrier cloud microphysics'
314  elseif (icmphys == 3) then
315  print *,' --- zhao/carr/sundqvist + pdf cloud'
316  else
317  print *,' !!! ERROR in cloud microphysc specification!!!', &
318  & ' icmphys (NP3D) =',icmphys
319  stop
320  endif
321  endif
322  endif
323 
325 ! --- compute llyr - the top of bl cld and is topmost non cld(low) layer
326 ! for stratiform (at or above lowest 0.1 of the atmosphere)
327 
328  if ( ivflip == 0 ) then ! data from toa to sfc
329  lab_do_k0 : do k = nlay, 2, -1
330  kl = k
331  if (si(k) < 0.9e0) exit lab_do_k0
332  enddo lab_do_k0
333 
334  llyr = kl
335  else ! data from sfc to top
336  lab_do_k1 : do k = 2, nlay
337  kl = k
338  if (si(k) < 0.9e0) exit lab_do_k1
339  enddo lab_do_k1
340 
341  llyr = kl - 1
342  endif ! end_if_ivflip
343 
344 !
345  return
346 !...................................
347  end subroutine cld_init
348 !-----------------------------------
350 
381 !-----------------------------------
382  subroutine progcld1 &
383 !...................................
384 
385 ! --- inputs:
386  & ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, &
387  & xlat,xlon,slmsk, &
388  & ix, nlay, nlp1, shoc_cld, cldcov, &
389 ! --- outputs:
390  & clouds,clds,mtop,mbot &
391  & )
392 
393 ! ================= subprogram documentation block ================ !
394 ! !
395 ! subprogram: progcld1 computes cloud related quantities using !
396 ! zhao/moorthi's prognostic cloud microphysics scheme. !
397 ! !
398 ! abstract: this program computes cloud fractions from cloud !
399 ! condensates, calculates liquid/ice cloud droplet effective radius, !
400 ! and computes the low, mid, high, total and boundary layer cloud !
401 ! fractions and the vertical indices of low, mid, and high cloud !
402 ! top and base. the three vertical cloud domains are set up in the !
403 ! initial subroutine "cld_init". !
404 ! !
405 ! usage: call progcld1 !
406 ! !
407 ! subprograms called: gethml !
408 ! !
409 ! attributes: !
410 ! language: fortran 90 !
411 ! machine: ibm-sp, sgi !
412 ! !
413 ! !
414 ! ==================== definition of variables ==================== !
415 ! !
416 ! input variables: !
417 ! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) !
418 ! plvl (IX,NLP1) : model level pressure in mb (100Pa) !
419 ! tlyr (IX,NLAY) : model layer mean temperature in k !
420 ! tvly (IX,NLAY) : model layer virtual temperature in k !
421 ! qlyr (IX,NLAY) : layer specific humidity in gm/gm !
422 ! qstl (IX,NLAY) : layer saturate humidity in gm/gm !
423 ! rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) !
424 ! clw (IX,NLAY) : layer cloud condensate amount !
425 ! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2!
426 ! range, otherwise see in-line comment !
427 ! xlon (IX) : grid longitude in radians (not used) !
428 ! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) !
429 ! IX : horizontal dimention !
430 ! NLAY,NLP1 : vertical layer/level dimensions !
431 ! !
432 ! output variables: !
433 ! clouds(IX,NLAY,NF_CLDS) : cloud profiles !
434 ! clouds(:,:,1) - layer total cloud fraction !
435 ! clouds(:,:,2) - layer cloud liq water path (g/m**2) !
436 ! clouds(:,:,3) - mean eff radius for liq cloud (micron) !
437 ! clouds(:,:,4) - layer cloud ice water path (g/m**2) !
438 ! clouds(:,:,5) - mean eff radius for ice cloud (micron) !
439 ! clouds(:,:,6) - layer rain drop water path not assigned !
440 ! clouds(:,:,7) - mean eff radius for rain drop (micron) !
441 ! *** clouds(:,:,8) - layer snow flake water path not assigned !
442 ! clouds(:,:,9) - mean eff radius for snow flake (micron) !
443 ! *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) !
444 ! clds (IX,5) : fraction of clouds for low, mid, hi, tot, bl !
445 ! mtop (IX,3) : vertical indices for low, mid, hi cloud tops !
446 ! mbot (IX,3) : vertical indices for low, mid, hi cloud bases !
447 ! !
448 ! module variables: !
449 ! ivflip : control flag of vertical index direction !
450 ! =0: index from toa to surface !
451 ! =1: index from surface to toa !
452 ! lsashal : control flag for shallow convection !
453 ! lcrick : control flag for eliminating CRICK !
454 ! =t: apply layer smoothing to eliminate CRICK !
455 ! =f: do not apply layer smoothing !
456 ! lcnorm : control flag for in-cld condensate !
457 ! =t: normalize cloud condensate !
458 ! =f: not normalize cloud condensate !
459 ! !
460 ! ==================== end of description ===================== !
461 !
462  implicit none
463 
464 ! --- inputs
465  integer, intent(in) :: IX, NLAY, NLP1
466 
467  logical, intent(in) :: shoc_cld
468  real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, &
469  & tlyr, tvly, qlyr, qstl, rhly, clw, cldcov
470 
471  real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, &
472  & slmsk
473 
474 ! --- outputs
475  real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds
476 
477  real (kind=kind_phys), dimension(:,:), intent(out) :: clds
478 
479  integer, dimension(:,:), intent(out) :: mtop,mbot
480 
481 ! --- local variables:
482  real (kind=kind_phys), dimension(IX,NLAY) :: cldtot, cldcnv, &
483  & cwp, cip, crp, csp, rew, rei, res, rer, delp, tem2d, clwf
484 
485  real (kind=kind_phys) :: ptop1(ix,nk_clds+1)
486 
487  real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, &
488  & tem1, tem2, tem3
489 
490  integer :: i, k, id, nf
491 
492 !
493 !===> ... begin here
494 !
495  do nf=1,nf_clds
496  do k=1,nlay
497  do i=1,ix
498  clouds(i,k,nf) = 0.0
499  enddo
500  enddo
501  enddo
502 ! clouds(:,:,:) = 0.0
503 
504  do k = 1, nlay
505  do i = 1, ix
506  cldtot(i,k) = 0.0
507  cldcnv(i,k) = 0.0
508  cwp (i,k) = 0.0
509  cip (i,k) = 0.0
510  crp (i,k) = 0.0
511  csp (i,k) = 0.0
512  rew (i,k) = reliq_def ! default liq radius to 10 micron
513  rei (i,k) = reice_def ! default ice radius to 50 micron
514  rer (i,k) = rrain_def ! default rain radius to 1000 micron
515  res (i,k) = rsnow_def ! default snow radius to 250 micron
516  tem2d(i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) )
517  clwf(i,k) = 0.0
518  enddo
519  enddo
520 !
521  if ( lcrick ) then
522  do i = 1, ix
523  clwf(i,1) = 0.75*clw(i,1) + 0.25*clw(i,2)
524  clwf(i,nlay) = 0.75*clw(i,nlay) + 0.25*clw(i,nlay-1)
525  enddo
526  do k = 2, nlay-1
527  do i = 1, ix
528  clwf(i,k) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1)
529  enddo
530  enddo
531  else
532  do k = 1, nlay
533  do i = 1, ix
534  clwf(i,k) = clw(i,k)
535  enddo
536  enddo
537  endif
538 
540 ! --- find top pressure for each cloud domain for given latitude
541 ! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h;
542 ! --- i=1,2 are low-lat (<45 degree) and pole regions)
543 
544  do id = 1, 4
545  tem1 = ptopc(id,2) - ptopc(id,1)
546 
547  do i =1, ix
548  tem2 = xlat(i) / con_pi ! if xlat in pi/2 -> -pi/2 range
549 ! tem2 = 0.5 - xlat(i)/con_pi ! if xlat in 0 -> pi range
550 
551  ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 )
552  enddo
553  enddo
554 
556 ! --- compute liquid/ice condensate path in g/m**2
557 
558  if ( ivflip == 0 ) then ! input data from toa to sfc
559  do k = 1, nlay
560  do i = 1, ix
561  delp(i,k) = plvl(i,k+1) - plvl(i,k)
562  clwt = max(0.0, clwf(i,k)) * gfac * delp(i,k)
563  cip(i,k) = clwt * tem2d(i,k)
564  cwp(i,k) = clwt - cip(i,k)
565  enddo
566  enddo
567  else ! input data from sfc to toa
568  do k = 1, nlay
569  do i = 1, ix
570  delp(i,k) = plvl(i,k) - plvl(i,k+1)
571  clwt = max(0.0, clwf(i,k)) * gfac * delp(i,k)
572  cip(i,k) = clwt * tem2d(i,k)
573  cwp(i,k) = clwt - cip(i,k)
574  enddo
575  enddo
576  endif ! end_if_ivflip
577 
579 ! --- effective liquid cloud droplet radius over land
580 
581  do i = 1, ix
582  if (nint(slmsk(i)) == 1) then
583  do k = 1, nlay
584  rew(i,k) = 5.0 + 5.0 * tem2d(i,k)
585  enddo
586  endif
587  enddo
588  if (shoc_cld) then ! use shoc generated sgs clouds
589  do k = 1, nlay
590  do i = 1, ix
591  cldtot(i,k) = cldcov(i,k)
592  enddo
593  enddo
594 
595  else
596 ! --- layer cloud fraction
597 
598  if ( ivflip == 0 ) then ! input data from toa to sfc
599 
600  clwmin = 0.0
601  if (.not. lsashal) then
602  do k = nlay, 1, -1
603  do i = 1, ix
604  clwt = 1.0e-6 * (plyr(i,k)*0.001)
605 ! clwt = 2.0e-6 * (plyr(i,k)*0.001)
606 
607  if (clwf(i,k) > clwt) then
608 
609  onemrh= max( 1.e-10, 1.0-rhly(i,k) )
610  clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
611 
612  tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
613  tem1 = 2000.0 / tem1
614 ! tem1 = 1000.0 / tem1
615 
616  value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 )
617  tem2 = sqrt( sqrt(rhly(i,k)) )
618 
619  cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
620  endif
621  enddo
622  enddo
623  else
624  do k = nlay, 1, -1
625  do i = 1, ix
626  clwt = 1.0e-6 * (plyr(i,k)*0.001)
627 ! clwt = 2.0e-6 * (plyr(i,k)*0.001)
628 
629  if (clwf(i,k) > clwt) then
630  onemrh= max( 1.e-10, 1.0-rhly(i,k) )
631  clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
632 
633 ! tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
634 ! tem1 = 2000.0 / tem1
635 
636  tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan
637  tem1 = 100.0 / tem1
638 !
639 ! tem1 = 2000.0 / tem1
640 ! tem1 = 1000.0 / tem1
641 !
642 
643  value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 )
644  tem2 = sqrt( sqrt(rhly(i,k)) )
645  cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
646  endif
647  enddo
648  enddo
649  endif
650 
651  else ! input data from sfc to toa
652 
653  clwmin = 0.0
654  if (.not. lsashal) then
655  do k = 1, nlay
656  do i = 1, ix
657  clwt = 1.0e-6 * (plyr(i,k)*0.001)
658 ! clwt = 2.0e-6 * (plyr(i,k)*0.001)
659 
660  if (clwf(i,k) > clwt) then
661 
662  onemrh= max( 1.e-10, 1.0-rhly(i,k) )
663  clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
664 
665  tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
666  tem1 = 2000.0 / tem1
667 
668 ! tem1 = 1000.0 / tem1
669 
670  value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 )
671  tem2 = sqrt( sqrt(rhly(i,k)) )
672 
673  cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
674  endif
675  enddo
676  enddo
677  else
678  do k = 1, nlay
679  do i = 1, ix
680  clwt = 1.0e-6 * (plyr(i,k)*0.001)
681 ! clwt = 2.0e-6 * (plyr(i,k)*0.001)
682 
683  if (clwf(i,k) > clwt) then
684  onemrh= max( 1.e-10, 1.0-rhly(i,k) )
685  clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
686 
687 ! tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
688 ! tem1 = 2000.0 / tem1
689 
690  tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan
691  tem1 = 100.0 / tem1
692 !
693 ! tem1 = 2000.0 / tem1
694 ! tem1 = 1000.0 / tem1
695 !
696 
697  value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 )
698  tem2 = sqrt( sqrt(rhly(i,k)) )
699 
700  cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
701  endif
702  enddo
703  enddo
704  endif
705 
706  endif ! end_if_flip
707  endif ! if (shoc_cld) then
708 
709  do k = 1, nlay
710  do i = 1, ix
711  if (cldtot(i,k) < climit) then
712  cldtot(i,k) = 0.0
713  cwp(i,k) = 0.0
714  cip(i,k) = 0.0
715  crp(i,k) = 0.0
716  csp(i,k) = 0.0
717  endif
718  enddo
719  enddo
720 
721  if ( lcnorm ) then
722  do k = 1, nlay
723  do i = 1, ix
724  if (cldtot(i,k) >= climit) then
725  tem1 = 1.0 / max(climit2, cldtot(i,k))
726  cwp(i,k) = cwp(i,k) * tem1
727  cip(i,k) = cip(i,k) * tem1
728  crp(i,k) = crp(i,k) * tem1
729  csp(i,k) = csp(i,k) * tem1
730  endif
731  enddo
732  enddo
733  endif
734 
742 ! --- effective ice cloud droplet radius
743 
744  do k = 1, nlay
745  do i = 1, ix
746  tem2 = tlyr(i,k) - con_ttp
747 
748  if (cip(i,k) > 0.0) then
749  tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k))
750 
751  if (tem2 < -50.0) then
752  rei(i,k) = (1250.0/9.917) * tem3 ** 0.109
753  elseif (tem2 < -40.0) then
754  rei(i,k) = (1250.0/9.337) * tem3 ** 0.08
755  elseif (tem2 < -30.0) then
756  rei(i,k) = (1250.0/9.208) * tem3 ** 0.055
757  else
758  rei(i,k) = (1250.0/9.387) * tem3 ** 0.031
759  endif
760 ! rei(i,k) = max(20.0, min(rei(i,k), 300.0))
761 ! rei(i,k) = max(10.0, min(rei(i,k), 100.0))
762  rei(i,k) = max(10.0, min(rei(i,k), 150.0))
763 ! rei(i,k) = max(5.0, min(rei(i,k), 130.0))
764  endif
765  enddo
766  enddo
767 
768 !
769  do k = 1, nlay
770  do i = 1, ix
771  clouds(i,k,1) = cldtot(i,k)
772  clouds(i,k,2) = cwp(i,k)
773  clouds(i,k,3) = rew(i,k)
774  clouds(i,k,4) = cip(i,k)
775  clouds(i,k,5) = rei(i,k)
776 ! clouds(i,k,6) = 0.0
777  clouds(i,k,7) = rer(i,k)
778 ! clouds(i,k,8) = 0.0
779  clouds(i,k,9) = rei(i,k)
780  enddo
781  enddo
782 
784 ! --- compute low, mid, high, total, and boundary layer cloud fractions
785 ! and clouds top/bottom layer indices for low, mid, and high clouds.
786 ! The three cloud domain boundaries are defined by ptopc. The cloud
787 ! overlapping method is defined by control flag 'iovr', which may
788 ! be different for lw and sw radiation programs.
789 
790  call gethml &
791 ! --- inputs:
792  & ( plyr, ptop1, cldtot, cldcnv, &
793  & ix,nlay, &
794 ! --- outputs:
795  & clds, mtop, mbot &
796  & )
797 
798 
799 !
800  return
801 !...................................
802  end subroutine progcld1
803 !-----------------------------------
805 
838 
841 !-----------------------------------
842  subroutine progcld2 &
843 !...................................
844 
845 ! --- inputs:
846  & ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, &
847  & xlat,xlon,slmsk, f_ice,f_rain,r_rime,flgmin, &
848  & ix, nlay, nlp1, &
849 ! --- outputs:
850  & clouds,clds,mtop,mbot &
851  & )
852 
853 ! ================= subprogram documentation block ================ !
854 ! !
855 ! subprogram: progcld2 computes cloud related quantities using !
856 ! ferrier's prognostic cloud microphysics scheme. !
857 ! !
858 ! abstract: this program computes cloud fractions from cloud !
859 ! condensates, calculates liquid/ice cloud droplet effective radius, !
860 ! and computes the low, mid, high, total and boundary layer cloud !
861 ! fractions and the vertical indices of low, mid, and high cloud !
862 ! top and base. the three vertical cloud domains are set up in the !
863 ! initial subroutine "cld_init". !
864 ! !
865 ! usage: call progcld2 !
866 ! !
867 ! subprograms called: gethml !
868 ! !
869 ! attributes: !
870 ! language: fortran 90 !
871 ! machine: ibm-sp, sgi !
872 ! !
873 ! !
874 ! ==================== definition of variables ==================== !
875 ! !
876 ! input variables: !
877 ! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) !
878 ! plvl (IX,NLP1) : model level pressure in mb (100Pa) !
879 ! tlyr (IX,NLAY) : model layer mean temperature in k !
880 ! tvly (IX,NLAY) : model layer virtual temperature in k !
881 ! qlyr (IX,NLAY) : layer specific humidity in gm/gm !
882 ! qstl (IX,NLAY) : layer saturate humidity in gm/gm !
883 ! rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) !
884 ! clw (IX,NLAY) : layer cloud condensate amount !
885 ! f_ice (IX,NLAY) : fraction of layer cloud ice (ferrier micro-phys) !
886 ! f_rain(IX,NLAY) : fraction of layer rain water (ferrier micro-phys) !
887 ! r_rime(IX,NLAY) : mass ratio of total ice to unrimed ice (>=1) !
888 ! flgmin(IX) : minimim large ice fraction !
889 ! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2!
890 ! range, otherwise see in-line comment !
891 ! xlon (IX) : grid longitude in radians (not used) !
892 ! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) !
893 ! IX : horizontal dimention !
894 ! NLAY,NLP1 : vertical layer/level dimensions !
895 ! !
896 ! output variables: !
897 ! clouds(IX,NLAY,NF_CLDS) : cloud profiles !
898 ! clouds(:,:,1) - layer total cloud fraction !
899 ! clouds(:,:,2) - layer cloud liq water path (g/m**2) !
900 ! clouds(:,:,3) - mean eff radius for liq cloud (micron) !
901 ! clouds(:,:,4) - layer cloud ice water path (g/m**2) !
902 ! clouds(:,:,5) - mean eff radius for ice cloud (micron) !
903 ! clouds(:,:,6) - layer rain drop water path (g/m**2) !
904 ! clouds(:,:,7) - mean eff radius for rain drop (micron) !
905 ! *** clouds(:,:,8) - layer snow flake water path (g/m**2) !
906 ! clouds(:,:,9) - mean eff radius for snow flake (micron) !
907 ! *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) !
908 ! clds (IX,5) : fraction of clouds for low, mid, hi, tot, bl !
909 ! mtop (IX,3) : vertical indices for low, mid, hi cloud tops !
910 ! mbot (IX,3) : vertical indices for low, mid, hi cloud bases !
911 ! !
912 ! external module variables: !
913 ! ivflip : control flag of vertical index direction !
914 ! =0: index from toa to surface !
915 ! =1: index from surface to toa !
916 ! lsashal : control flag for shallow convection !
917 ! lcrick : control flag for eliminating CRICK !
918 ! =t: apply layer smoothing to eliminate CRICK !
919 ! =f: do not apply layer smoothing !
920 ! lcnorm : control flag for in-cld condensate !
921 ! =t: normalize cloud condensate !
922 ! =f: not normalize cloud condensate !
923 ! lnoprec : precip effect in radiation flag (ferrier scheme) !
924 ! =t: snow/rain has no impact on radiation !
925 ! =f: snow/rain has impact on radiation !
926 ! !
927 ! ==================== end of description ===================== !
928 !
929  implicit none
930 
931 ! --- constants
932 
933 ! --- inputs
934  integer, intent(in) :: IX, NLAY, NLP1
935 
936  real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, &
937  & tlyr, tvly, qlyr, qstl, rhly, clw, f_ice, f_rain, r_rime
938 
939  real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, &
940  & slmsk
941  real (kind=kind_phys), dimension(:), intent(in) :: flgmin
942 
943 ! --- outputs
944  real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds
945 
946  real (kind=kind_phys), dimension(:,:), intent(out) :: clds
947 
948  integer, dimension(:,:), intent(out) :: mtop,mbot
949 
950 ! --- local variables:
951  real (kind=kind_phys), dimension(IX,NLAY) :: cldtot, cldcnv, &
952  & cwp, cip, crp, csp, rew, rei, res, rer, tem2d, clw2, &
953  & qcwat, qcice, qrain, fcice, frain, rrime, rsden, clwf
954 
955  real (kind=kind_phys) :: ptop1(ix,nk_clds+1)
956 
957  real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, &
958  & tem1, tem2, tem3
959 
960  integer :: i, k, id
961 
962 !
963 !===> ... begin here
964 !
965 ! clouds(:,:,:) = 0.0
966 
967  do k = 1, nlay
968  do i = 1, ix
969  cldtot(i,k) = 0.0
970  cldcnv(i,k) = 0.0
971  cwp (i,k) = 0.0
972  cip (i,k) = 0.0
973  crp (i,k) = 0.0
974  csp (i,k) = 0.0
975  rew (i,k) = reliq_def ! default liq radius to 10 micron
976  rei (i,k) = reice_def ! default ice radius to 50 micron
977  rer (i,k) = rrain_def ! default rain radius to 1000 micron
978  res (i,k) = rsnow_def ! default snow radius to 250 micron
979  fcice(i,k) = max(0.0, min(1.0, f_ice(i,k)))
980  frain(i,k) = max(0.0, min(1.0, f_rain(i,k)))
981  rrime(i,k) = max(1.0, r_rime(i,k))
982  tem2d(i,k) = tlyr(i,k) - con_t0c
983  enddo
984  enddo
985 !
986  if ( lcrick ) then
987  do i = 1, ix
988  clwf(i,1) = 0.75*clw(i,1) + 0.25*clw(i,2)
989  clwf(i,nlay) = 0.75*clw(i,nlay) + 0.25*clw(i,nlay-1)
990  enddo
991  do k = 2, nlay-1
992  do i = 1, ix
993  clwf(i,k) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1)
994  enddo
995  enddo
996  else
997  do k = 1, nlay
998  do i = 1, ix
999  clwf(i,k) = clw(i,k)
1000  enddo
1001  enddo
1002  endif
1005 ! --- find top pressure for each cloud domain for given latitude
1006 ! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h;
1007 ! --- i=1,2 are low-lat (<45 degree) and pole regions)
1008 
1009  do id = 1, 4
1010  tem1 = ptopc(id,2) - ptopc(id,1)
1011 
1012  do i =1, ix
1013  tem2 = xlat(i) / con_pi ! if xlat in pi/2 -> -pi/2 range
1014 ! tem2 = 0.5 - xlat(i)/con_pi ! if xlat in 0 -> pi range
1015 
1016  ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 )
1017  enddo
1018  enddo
1019 
1022 ! --- separate cloud condensate into liquid, ice, and rain types, and
1023 ! save the liquid+ice condensate in array clw2 for later
1024 ! calculation of cloud fraction
1025 
1026  do k = 1, nlay
1027  do i = 1, ix
1028  if (tem2d(i,k) > -40.0) then
1029  qcice(i,k) = clwf(i,k) * fcice(i,k)
1030  tem1 = clwf(i,k) - qcice(i,k)
1031  qrain(i,k) = tem1 * frain(i,k)
1032  qcwat(i,k) = tem1 - qrain(i,k)
1033  clw2(i,k) = qcwat(i,k) + qcice(i,k)
1034  else
1035  qcice(i,k) = clwf(i,k)
1036  qrain(i,k) = 0.0
1037  qcwat(i,k) = 0.0
1038  clw2(i,k) = clwf(i,k)
1039  endif
1040  enddo
1041  enddo
1042 
1046  call rsipath2 &
1047 ! --- inputs:
1048  & ( plyr, plvl, tlyr, qlyr, qcwat, qcice, qrain, rrime, &
1049  & ix, nlay, ivflip, flgmin, &
1050 ! --- outputs:
1051  & cwp, cip, crp, csp, rew, rer, res, rsden &
1052  & )
1053 
1054 
1055  if ( ivflip == 0 ) then ! input data from toa to sfc
1056  do k = 1, nlay
1057  do i = 1, ix
1058  tem2d(i,k) = (con_g * plyr(i,k)) &
1059  & / (con_rd* (plvl(i,k+1) - plvl(i,k)))
1060  enddo
1061  enddo
1062  else ! input data from sfc to toa
1063  do k = 1, nlay
1064  do i = 1, ix
1065  tem2d(i,k) = (con_g * plyr(i,k)) &
1066  & / (con_rd* (plvl(i,k) - plvl(i,k+1)))
1067  enddo
1068  enddo
1069  endif ! end_if_ivflip
1071 ! --- layer cloud fraction
1072 
1073  if ( ivflip == 0 ) then ! input data from toa to sfc
1074 
1075  clwmin = 0.0
1076  if (.not. lsashal) then
1077  do k = nlay, 1, -1
1078  do i = 1, ix
1079 ! clwt = 1.0e-7 * (plyr(i,k)*0.001)
1080 ! clwt = 1.0e-6 * (plyr(i,k)*0.001)
1081  clwt = 2.0e-6 * (plyr(i,k)*0.001)
1082 ! clwt = 5.0e-6 * (plyr(i,k)*0.001)
1083 ! clwt = 5.0e-6
1084 
1085  if (clw2(i,k) > clwt) then
1086  onemrh= max( 1.e-10, 1.0-rhly(i,k) )
1087  clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
1088 
1089 ! tem1 = min(max(sqrt(onemrh*qstl(i,k)),0.0001),1.0)
1090 ! tem1 = 100.0 / tem1
1091 
1092  tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
1093  tem1 = 2000.0 / tem1
1094 ! tem1 = 2400.0 / tem1
1095 !cnt tem1 = 2500.0 / tem1
1096 ! tem1 = min(max(sqrt(onemrh*qstl(i,k)),0.0001),1.0)
1097 ! tem1 = 2000.0 / tem1
1098 ! tem1 = 1000.0 / tem1
1099 ! tem1 = 100.0 / tem1
1100 
1101  value = max( min( tem1*(clw2(i,k)-clwm), 50.0 ), 0.0 )
1102  tem2 = sqrt( sqrt(rhly(i,k)) )
1103 
1104  cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
1105  endif
1106  enddo
1107  enddo
1108  else
1109  do k = nlay, 1, -1
1110  do i = 1, ix
1111 ! clwt = 1.0e-6 * (plyr(i,k)*0.001)
1112  clwt = 2.0e-6 * (plyr(i,k)*0.001)
1113 
1114  if (clw2(i,k) > clwt) then
1115  onemrh= max( 1.e-10, 1.0-rhly(i,k) )
1116  clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
1117 
1118  tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan
1119  tem1 = 100.0 / tem1
1120 
1121 ! tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
1122 ! tem1 = 2000.0 / tem1
1123 !
1124 ! tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
1125 ! tem1 = 2200.0 / tem1
1126 ! tem1 = 2400.0 / tem1
1127 ! tem1 = 2500.0 / tem1
1128 ! tem1 = min(max(sqrt(onemrh*qstl(i,k)),0.0001),1.0)
1129 ! tem1 = 2000.0 / tem1
1130 ! tem1 = 1000.0 / tem1
1131 ! tem1 = 100.0 / tem1
1132 
1133  value = max( min( tem1*(clw2(i,k)-clwm), 50.0 ), 0.0 )
1134  tem2 = sqrt( sqrt(rhly(i,k)) )
1135 
1136  cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
1137  endif
1138  enddo
1139  enddo
1140  endif
1141 
1142  else ! input data from sfc to toa
1143 
1144  clwmin = 0.0e-6
1145  if (.not. lsashal) then
1146  do k = 1, nlay
1147  do i = 1, ix
1148 ! clwt = 1.0e-7 * (plyr(i,k)*0.001)
1149 ! clwt = 1.0e-6 * (plyr(i,k)*0.001)
1150  clwt = 2.0e-6 * (plyr(i,k)*0.001)
1151 ! clwt = 5.0e-6 * (plyr(i,k)*0.001)
1152 ! clwt = 5.0e-6
1153 
1154  if (clw2(i,k) > clwt) then
1155  onemrh= max( 1.e-10, 1.0-rhly(i,k) )
1156  clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
1157 
1158 ! tem1 = min(max(sqrt(onemrh*qstl(i,k)),0.0001),1.0)
1159 ! tem1 = 100.0 / tem1
1160 
1161  tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
1162  tem1 = 2000.0 / tem1
1163 ! tem1 = 2400.0 / tem1
1164 !cnt tem1 = 2500.0 / tem1
1165 ! tem1 = min(max(sqrt(onemrh*qstl(i,k)),0.0001),1.0)
1166 ! tem1 = 2000.0 / tem1
1167 ! tem1 = 1000.0 / tem1
1168 ! tem1 = 100.0 / tem1
1169 
1170  value = max( min( tem1*(clw2(i,k)-clwm), 50.0 ), 0.0 )
1171  tem2 = sqrt( sqrt(rhly(i,k)) )
1172 
1173  cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
1174  endif
1175  enddo
1176  enddo
1177  else
1178  do k = 1, nlay
1179  do i = 1, ix
1180 ! clwt = 1.0e-6 * (plyr(i,k)*0.001)
1181  clwt = 2.0e-6 * (plyr(i,k)*0.001)
1182 
1183  if (clw2(i,k) > clwt) then
1184  onemrh= max( 1.e-10, 1.0-rhly(i,k) )
1185  clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
1186 
1187  tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan
1188  tem1 = 100.0 / tem1
1189 
1190 ! tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
1191 ! tem1 = 2000.0 / tem1
1192 !
1193 ! tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
1194 ! tem1 = 2200.0 / tem1
1195 ! tem1 = 2400.0 / tem1
1196 ! tem1 = 2500.0 / tem1
1197 ! tem1 = min(max(sqrt(onemrh*qstl(i,k)),0.0001),1.0)
1198 ! tem1 = 2000.0 / tem1
1199 ! tem1 = 1000.0 / tem1
1200 ! tem1 = 100.0 / tem1
1201 
1202  value = max( min( tem1*(clw2(i,k)-clwm), 50.0 ), 0.0 )
1203  tem2 = sqrt( sqrt(rhly(i,k)) )
1204 
1205  cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
1206  endif
1207  enddo
1208  enddo
1209  endif
1210 
1211  endif ! end_if_flip
1212 
1213  do k = 1, nlay
1214  do i = 1, ix
1215  if (cldtot(i,k) < climit) then
1216  cldtot(i,k) = 0.0
1217  cwp(i,k) = 0.0
1218  cip(i,k) = 0.0
1219  crp(i,k) = 0.0
1220  csp(i,k) = 0.0
1221  endif
1222  enddo
1223  enddo
1224 
1225 ! When lnoprec = .true. snow/rain has no impact on radiation
1226  if ( lnoprec ) then
1227  do k = 1, nlay
1228  do i = 1, ix
1229  crp(i,k) = 0.0
1230  csp(i,k) = 0.0
1231  enddo
1232  enddo
1233  endif
1234 !
1235  if ( lcnorm ) then
1236  do k = 1, nlay
1237  do i = 1, ix
1238  if (cldtot(i,k) >= climit) then
1239  tem1 = 1.0 / max(climit2, cldtot(i,k))
1240  cwp(i,k) = cwp(i,k) * tem1
1241  cip(i,k) = cip(i,k) * tem1
1242  crp(i,k) = crp(i,k) * tem1
1243  csp(i,k) = csp(i,k) * tem1
1244  endif
1245  enddo
1246  enddo
1247  endif
1248 
1250 ! --- effective ice cloud droplet radius
1251 
1252  do k = 1, nlay
1253  do i = 1, ix
1254  tem1 = tlyr(i,k) - con_ttp
1255  tem2 = cip(i,k)
1256 
1257  if (tem2 > 0.0) then
1258  tem3 = tem2d(i,k) * tem2 / tvly(i,k)
1259 
1260  if (tem1 < -50.0) then
1261  rei(i,k) = (1250.0/9.917) * tem3 ** 0.109
1262  elseif (tem1 < -40.0) then
1263  rei(i,k) = (1250.0/9.337) * tem3 ** 0.08
1264  elseif (tem1 < -30.0) then
1265  rei(i,k) = (1250.0/9.208) * tem3 ** 0.055
1266  else
1267  rei(i,k) = (1250.0/9.387) * tem3 ** 0.031
1268  endif
1269 
1270 ! if (lprnt .and. k == l) print *,' reiL=',rei(i,k),' icec=', &
1271 ! & icec,' cip=',cip(i,k),' tem=',tem,' delt=',delt
1272 
1273  rei(i,k) = max(10.0, min(rei(i,k), 300.0))
1274 ! rei(i,k) = max(20.0, min(rei(i,k), 300.0))
1275 !!!! rei(i,k) = max(30.0, min(rei(i,k), 300.0))
1276 ! rei(i,k) = max(50.0, min(rei(i,k), 300.0))
1277 ! rei(i,k) = max(100.0, min(rei(i,k), 300.0))
1278  endif
1279  enddo
1280  enddo
1281 !
1282  do k = 1, nlay
1283  do i = 1, ix
1284  clouds(i,k,1) = cldtot(i,k)
1285  clouds(i,k,2) = cwp(i,k)
1286  clouds(i,k,3) = rew(i,k)
1287  clouds(i,k,4) = cip(i,k)
1288  clouds(i,k,5) = rei(i,k)
1289  clouds(i,k,6) = crp(i,k)
1290  clouds(i,k,7) = rer(i,k)
1291 ! clouds(i,k,8) = csp(i,k) !ncar scheme
1292  clouds(i,k,8) = csp(i,k) * rsden(i,k) !fu's scheme
1293  clouds(i,k,9) = rei(i,k)
1294  enddo
1295  enddo
1296 
1299 ! --- compute low, mid, high, total, and boundary layer cloud fractions
1300 ! and clouds top/bottom layer indices for low, mid, and high clouds.
1301 ! The three cloud domain boundaries are defined by ptopc. The cloud
1302 ! overlapping method is defined by control flag 'iovr', which may
1303 ! be different for lw and sw radiation programs.
1304 
1305  call gethml &
1306 ! --- inputs:
1307  & ( plyr, ptop1, cldtot, cldcnv, &
1308  & ix,nlay, &
1309 ! --- outputs:
1310  & clds, mtop, mbot &
1311  & )
1312 
1313 
1314 !
1315  return
1316 !...................................
1317  end subroutine progcld2
1318 !-----------------------------------
1320 
1354 !-----------------------------------
1355  subroutine progcld3 &
1356 !...................................
1357 
1358 ! --- inputs:
1359  & ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw,cnvw,cnvc, &
1360  & xlat,xlon,slmsk, &
1361  & ix, nlay, nlp1, &
1362  & deltaq,sup,kdt,me, &
1363 ! --- outputs:
1364  & clouds,clds,mtop,mbot &
1365  & )
1366 
1367 ! ================= subprogram documentation block ================ !
1368 ! !
1369 ! subprogram: progcld3 computes cloud related quantities using !
1370 ! zhao/moorthi's prognostic cloud microphysics scheme. !
1371 ! !
1372 ! abstract: this program computes cloud fractions from cloud !
1373 ! condensates, calculates liquid/ice cloud droplet effective radius, !
1374 ! and computes the low, mid, high, total and boundary layer cloud !
1375 ! fractions and the vertical indices of low, mid, and high cloud !
1376 ! top and base. the three vertical cloud domains are set up in the !
1377 ! initial subroutine "cld_init". !
1378 ! !
1379 ! usage: call progcld3 !
1380 ! !
1381 ! subprograms called: gethml !
1382 ! !
1383 ! attributes: !
1384 ! language: fortran 90 !
1385 ! machine: ibm-sp, sgi !
1386 ! !
1387 ! !
1388 ! ==================== defination of variables ==================== !
1389 ! !
1390 ! input variables: !
1391 ! plyr (ix,nlay) : model layer mean pressure in mb (100pa) !
1392 ! plvl (ix,nlp1) : model level pressure in mb (100pa) !
1393 ! tlyr (ix,nlay) : model layer mean temperature in k !
1394 ! tvly (ix,nlay) : model layer virtual temperature in k !
1395 ! qlyr (ix,nlay) : layer specific humidity in gm/gm !
1396 ! qstl (ix,nlay) : layer saturate humidity in gm/gm !
1397 ! rhly (ix,nlay) : layer relative humidity (=qlyr/qstl) !
1398 ! clw (ix,nlay) : layer cloud condensate amount !
1399 ! xlat (ix) : grid latitude in radians, default to pi/2 -> -pi/2!
1400 ! range, otherwise see in-line comment !
1401 ! xlon (ix) : grid longitude in radians (not used) !
1402 ! slmsk (ix) : sea/land mask array (sea:0,land:1,sea-ice:2) !
1403 ! ix : horizontal dimention !
1404 ! nlay,nlp1 : vertical layer/level dimensions !
1405 ! cnvw (ix,nlay) : layer convective cloud condensate !
1406 ! cnvc (ix,nlay) : layer convective cloud cover !
1407 ! deltaq(ix,nlay) : half total water distribution width !
1408 ! sup : supersaturation !
1409 
1410 ! !
1411 ! output variables: !
1412 ! clouds(ix,nlay,nf_clds) : cloud profiles !
1413 ! clouds(:,:,1) - layer total cloud fraction !
1414 ! clouds(:,:,2) - layer cloud liq water path (g/m**2) !
1415 ! clouds(:,:,3) - mean eff radius for liq cloud (micron) !
1416 ! clouds(:,:,4) - layer cloud ice water path (g/m**2) !
1417 ! clouds(:,:,5) - mean eff radius for ice cloud (micron) !
1418 ! clouds(:,:,6) - layer rain drop water path not assigned !
1419 ! clouds(:,:,7) - mean eff radius for rain drop (micron) !
1420 ! *** clouds(:,:,8) - layer snow flake water path not assigned !
1421 ! clouds(:,:,9) - mean eff radius for snow flake (micron) !
1422 ! *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) !
1423 ! clds (ix,5) : fraction of clouds for low, mid, hi, tot, bl !
1424 ! mtop (ix,3) : vertical indices for low, mid, hi cloud tops !
1425 ! mbot (ix,3) : vertical indices for low, mid, hi cloud bases !
1426 ! !
1427 ! module variables: !
1428 ! ivflip : control flag of vertical index direction !
1429 ! =0: index from toa to surface !
1430 ! =1: index from surface to toa !
1431 ! lcrick : control flag for eliminating crick !
1432 ! =t: apply layer smoothing to eliminate crick !
1433 ! =f: do not apply layer smoothing !
1434 ! lcnorm : control flag for in-cld condensate !
1435 ! =t: normalize cloud condensate !
1436 ! =f: not normalize cloud condensate !
1437 ! !
1438 ! ==================== end of description ===================== !
1439 !
1440  implicit none
1441 
1442 ! --- inputs
1443  integer, intent(in) :: ix, nlay, nlp1,kdt
1444 
1445  real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, &
1446  & tlyr, tvly, qlyr, qstl, rhly, clw
1447 ! & tlyr, tvly, qlyr, qstl, rhly, clw, cnvw, cnvc
1448 ! real (kind=kind_phys), dimension(:,:), intent(in) :: deltaq
1449  real (kind=kind_phys), dimension(:,:) :: deltaq, cnvw, cnvc
1450  real (kind=kind_phys) qtmp,qsc,rhs
1451  real (kind=kind_phys), intent(in) :: sup
1452  real (kind=kind_phys), parameter :: epsq = 1.0e-12
1453 
1454  real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, &
1455  & slmsk
1456  integer :: me
1457 
1458 ! --- outputs
1459  real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds
1460 
1461  real (kind=kind_phys), dimension(:,:), intent(out) :: clds
1462 
1463  integer, dimension(:,:), intent(out) :: mtop,mbot
1464 
1465 ! --- local variables:
1466  real (kind=kind_phys), dimension(ix,nlay) :: cldtot, cldcnv, &
1467  & cwp, cip, crp, csp, rew, rei, res, rer, delp, tem2d, clwf
1468 
1469  real (kind=kind_phys) :: ptop1(ix,nk_clds+1)
1470 
1471  real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, &
1472  & tem1, tem2, tem3
1473 
1474  integer :: i, k, id, nf
1475 
1476 !
1477 !===> ... begin here
1478 !
1479  do nf=1,nf_clds
1480  do k=1,nlay
1481  do i=1,ix
1482  clouds(i,k,nf) = 0.0
1483  enddo
1484  enddo
1485  enddo
1486 ! clouds(:,:,:) = 0.0
1487 
1488  do k = 1, nlay
1489  do i = 1, ix
1490  cldtot(i,k) = 0.0
1491  cldcnv(i,k) = 0.0
1492  cwp (i,k) = 0.0
1493  cip (i,k) = 0.0
1494  crp (i,k) = 0.0
1495  csp (i,k) = 0.0
1496  rew (i,k) = reliq_def ! default liq radius to 10 micron
1497  rei (i,k) = reice_def ! default ice radius to 50 micron
1498  rer (i,k) = rrain_def ! default rain radius to 1000 micron
1499  res (i,k) = rsnow_def ! default snow radius to 250 micron
1500  tem2d(i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) )
1501  clwf(i,k) = 0.0
1502  enddo
1503  enddo
1504 !
1505  if ( lcrick ) then
1506  do i = 1, ix
1507  clwf(i,1) = 0.75*clw(i,1) + 0.25*clw(i,2)
1508  clwf(i,nlay) = 0.75*clw(i,nlay) + 0.25*clw(i,nlay-1)
1509  enddo
1510  do k = 2, nlay-1
1511  do i = 1, ix
1512  clwf(i,k) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1)
1513  enddo
1514  enddo
1515  else
1516  do k = 1, nlay
1517  do i = 1, ix
1518  clwf(i,k) = clw(i,k)
1519  enddo
1520  enddo
1521  endif
1522 
1523  if(kdt==1) then
1524  do k = 1, nlay
1525  do i = 1, ix
1526  deltaq(i,k) = (1.-0.95)*qstl(i,k)
1527  enddo
1528  enddo
1529  endif
1530 
1532 ! --- find top pressure for each cloud domain for given latitude
1533 ! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,l,m,h;
1534 ! --- i=1,2 are low-lat (<45 degree) and pole regions)
1535 
1536  do id = 1, 4
1537  tem1 = ptopc(id,2) - ptopc(id,1)
1538 
1539  do i =1, ix
1540  tem2 = xlat(i) / con_pi ! if xlat in pi/2 -> -pi/2 range
1541 ! tem2 = 0.5 - xlat(i)/con_pi ! if xlat in 0 -> pi range
1542 
1543  ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 )
1544  enddo
1545  enddo
1546 
1548 ! --- compute liquid/ice condensate path in g/m**2
1549 
1550  if ( ivflip == 0 ) then ! input data from toa to sfc
1551  do k = 1, nlay
1552  do i = 1, ix
1553  delp(i,k) = plvl(i,k+1) - plvl(i,k)
1554  clwt = max(0.0,(clwf(i,k)+cnvw(i,k))) * gfac * delp(i,k)
1555  cip(i,k) = clwt * tem2d(i,k)
1556  cwp(i,k) = clwt - cip(i,k)
1557  enddo
1558  enddo
1559  else ! input data from sfc to toa
1560  do k = 1, nlay
1561  do i = 1, ix
1562  delp(i,k) = plvl(i,k) - plvl(i,k+1)
1563  clwt = max(0.0,(clwf(i,k)+cnvw(i,k))) * gfac * delp(i,k)
1564  cip(i,k) = clwt * tem2d(i,k)
1565  cwp(i,k) = clwt - cip(i,k)
1566  enddo
1567  enddo
1568  endif ! end_if_ivflip
1569 
1571 ! --- effective liquid cloud droplet radius over land
1572 
1573  do i = 1, ix
1574  if (nint(slmsk(i)) == 1) then
1575  do k = 1, nlay
1576  rew(i,k) = 5.0 + 5.0 * tem2d(i,k)
1577  enddo
1578  endif
1579  enddo
1580 
1581 ! --- layer cloud fraction
1582 
1583  if ( ivflip == 0 ) then ! input data from toa to sfc
1584  do k = nlay, 1, -1
1585  do i = 1, ix
1586  tem1 = tlyr(i,k) - 273.16
1587  if(tem1 < con_thgni) then ! for pure ice, has to be consistent with gscond
1588  qsc = sup * qstl(i,k)
1589  rhs = sup
1590  else
1591  qsc = qstl(i,k)
1592  rhs = 1.0
1593  endif
1594  if(rhly(i,k) >= rhs) then
1595  cldtot(i,k) = 1.0
1596  else
1597  qtmp = qlyr(i,k) + clwf(i,k) - qsc
1598  if(deltaq(i,k) > epsq) then
1599  if(qtmp < -deltaq(i,k) .or. clwf(i,k) < epsq) then
1600 ! if(qtmp < -deltaq(i,k)) then
1601  cldtot(i,k) = 0.0
1602  elseif(qtmp >= deltaq(i,k)) then
1603  cldtot(i,k) = 1.0
1604  else
1605  cldtot(i,k) = 0.5*qtmp/deltaq(i,k) + 0.5
1606  cldtot(i,k) = max(cldtot(i,k),0.0)
1607  cldtot(i,k) = min(cldtot(i,k),1.0)
1608  endif
1609  else
1610  if(qtmp.gt.0) then
1611  cldtot(i,k) = 1.0
1612  else
1613  cldtot(i,k) = 0.0
1614  endif
1615  endif
1616  endif
1617  cldtot(i,k) = cnvc(i,k)+(1-cnvc(i,k))*cldtot(i,k)
1618  cldtot(i,k) = max(cldtot(i,k),0.)
1619  cldtot(i,k) = min(cldtot(i,k),1.)
1620  enddo
1621  enddo
1622  else ! input data from sfc to toa
1623  do k = 1, nlay
1624  do i = 1, ix
1625  tem1 = tlyr(i,k) - 273.16
1626  if(tem1 < con_thgni) then ! for pure ice, has to be consistent with gscond
1627  qsc = sup * qstl(i,k)
1628  rhs = sup
1629  else
1630  qsc = qstl(i,k)
1631  rhs = 1.0
1632  endif
1633  if(rhly(i,k) >= rhs) then
1634  cldtot(i,k) = 1.0
1635  else
1636  qtmp = qlyr(i,k) + clwf(i,k) - qsc
1637  if(deltaq(i,k) > epsq) then
1638 ! if(qtmp <= -deltaq(i,k) .or. cwmik < epsq) then
1639  if(qtmp <= -deltaq(i,k)) then
1640  cldtot(i,k) = 0.0
1641  elseif(qtmp >= deltaq(i,k)) then
1642  cldtot(i,k) = 1.0
1643  else
1644  cldtot(i,k) = 0.5*qtmp/deltaq(i,k) + 0.5
1645  cldtot(i,k) = max(cldtot(i,k),0.0)
1646  cldtot(i,k) = min(cldtot(i,k),1.0)
1647  endif
1648  else
1649  if(qtmp > 0.) then
1650  cldtot(i,k) = 1.0
1651  else
1652  cldtot(i,k) = 0.0
1653  endif
1654  endif
1655  endif
1656  cldtot(i,k) = cnvc(i,k) + (1-cnvc(i,k))*cldtot(i,k)
1657  cldtot(i,k) = max(cldtot(i,k),0.)
1658  cldtot(i,k) = min(cldtot(i,k),1.)
1659 
1660  enddo
1661  enddo
1662  endif ! end_if_flip
1663 
1664  do k = 1, nlay
1665  do i = 1, ix
1666  if (cldtot(i,k) < climit) then
1667  cldtot(i,k) = 0.0
1668  cwp(i,k) = 0.0
1669  cip(i,k) = 0.0
1670  crp(i,k) = 0.0
1671  csp(i,k) = 0.0
1672  endif
1673  enddo
1674  enddo
1675 
1676  if ( lcnorm ) then
1677  do k = 1, nlay
1678  do i = 1, ix
1679  if (cldtot(i,k) >= climit) then
1680  tem1 = 1.0 / max(climit2, cldtot(i,k))
1681  cwp(i,k) = cwp(i,k) * tem1
1682  cip(i,k) = cip(i,k) * tem1
1683  crp(i,k) = crp(i,k) * tem1
1684  csp(i,k) = csp(i,k) * tem1
1685  endif
1686  enddo
1687  enddo
1688  endif
1689 
1697 ! --- effective ice cloud droplet radius
1698 
1699  do k = 1, nlay
1700  do i = 1, ix
1701  tem2 = tlyr(i,k) - con_ttp
1702 
1703  if (cip(i,k) > 0.0) then
1704 ! tem3 = gord * cip(i,k) * (plyr(i,k)/delp(i,k)) / tvly(i,k)
1705  tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k))
1706 
1707  if (tem2 < -50.0) then
1708  rei(i,k) = (1250.0/9.917) * tem3 ** 0.109
1709  elseif (tem2 < -40.0) then
1710  rei(i,k) = (1250.0/9.337) * tem3 ** 0.08
1711  elseif (tem2 < -30.0) then
1712  rei(i,k) = (1250.0/9.208) * tem3 ** 0.055
1713  else
1714  rei(i,k) = (1250.0/9.387) * tem3 ** 0.031
1715  endif
1716 ! rei(i,k) = max(20.0, min(rei(i,k), 300.0))
1717 ! rei(i,k) = max(10.0, min(rei(i,k), 100.0))
1718  rei(i,k) = max(10.0, min(rei(i,k), 150.0))
1719 ! rei(i,k) = max(5.0, min(rei(i,k), 130.0))
1720  endif
1721  enddo
1722  enddo
1723 
1724 !
1725  do k = 1, nlay
1726  do i = 1, ix
1727  clouds(i,k,1) = cldtot(i,k)
1728  clouds(i,k,2) = cwp(i,k)
1729  clouds(i,k,3) = rew(i,k)
1730  clouds(i,k,4) = cip(i,k)
1731  clouds(i,k,5) = rei(i,k)
1732 ! clouds(i,k,6) = 0.0
1733  clouds(i,k,7) = rer(i,k)
1734 ! clouds(i,k,8) = 0.0
1735  clouds(i,k,9) = rei(i,k)
1736  enddo
1737  enddo
1738 
1740 ! --- compute low, mid, high, total, and boundary layer cloud fractions
1741 ! and clouds top/bottom layer indices for low, mid, and high clouds.
1742 ! the three cloud domain boundaries are defined by ptopc. the cloud
1743 ! overlapping method is defined by control flag 'iovr', which may
1744 ! be different for lw and sw radiation programs.
1745 
1746 
1747  call gethml &
1748 ! --- inputs:
1749  & ( plyr, ptop1, cldtot, cldcnv, &
1750  & ix,nlay, &
1751 ! --- outputs:
1752  & clds, mtop, mbot &
1753  & )
1754 
1755 
1756 !
1757  return
1758 !...................................
1759  end subroutine progcld3
1760 !-----------------------------------
1762 
1786 !-----------------------------------
1787  subroutine diagcld1 &
1788 
1789 !...................................
1790 
1791 ! --- inputs:
1792  & ( plyr,plvl,tlyr,rhly,vvel,cv,cvt,cvb, &
1793  & xlat,xlon,slmsk, &
1794  & ix, nlay, nlp1, &
1795 ! --- outputs:
1796  & clouds,clds,mtop,mbot &
1797  & )
1798 
1799 ! ================= subprogram documentation block ================ !
1800 ! !
1801 ! subprogram: diagcld1 computes cloud fractions for radiation !
1802 ! calculations. !
1803 ! !
1804 ! abstract: clouds are diagnosed from layer relative humidity, and !
1805 ! estimate cloud optical depth from temperature and layer thickness. !
1806 ! then computes the low, mid, high, total and boundary layer cloud !
1807 ! fractions and the vertical indices of low, mid, and high cloud top !
1808 ! and base. the three vertical cloud domains are set up in the !
1809 ! initial subroutine "cld_init". !
1810 ! !
1811 ! usage: call diagcld1 !
1812 ! !
1813 ! subprograms called: gethml !
1814 ! !
1815 ! attributes: !
1816 ! language: fortran 90 !
1817 ! machine: ibm-sp, sgi !
1818 ! !
1819 ! !
1820 ! ==================== definition of variables ==================== !
1821 ! !
1822 ! input variables: !
1823 ! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) !
1824 ! plvl (IX,NLP1) : model level pressure in mb (100Pa) !
1825 ! tlyr (IX,NLAY) : model layer mean temperature in k !
1826 ! rhly (IX,NLAY) : layer relative humidity !
1827 ! vvel (IX,NLAY) : layer mean vertical velocity in mb/sec !
1828 ! clw (IX,NLAY) : layer cloud condensate amount (not used) !
1829 ! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2!
1830 ! range, otherwise see in-line comment !
1831 ! xlon (IX) : grid longitude in radians, ok for both 0->2pi or !
1832 ! -pi -> +pi ranges !
1833 ! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) !
1834 ! cv (IX) : fraction of convective cloud !
1835 ! cvt, cvb (IX) : conv cloud top/bottom pressure in mb !
1836 ! IX : horizontal dimention !
1837 ! NLAY,NLP1 : vertical layer/level dimensions !
1838 ! !
1839 ! output variables: !
1840 ! clouds(IX,NLAY,NF_CLDS) : cloud profiles !
1841 ! clouds(:,:,1) - layer total cloud fraction !
1842 ! clouds(:,:,2) - layer cloud optical depth !
1843 ! clouds(:,:,3) - layer cloud single scattering albedo !
1844 ! clouds(:,:,4) - layer cloud asymmetry factor !
1845 ! clds (IX,5) : fraction of clouds for low, mid, hi, tot, bl !
1846 ! mtop (IX,3) : vertical indices for low, mid, hi cloud tops !
1847 ! mbot (IX,3) : vertical indices for low, mid, hi cloud bases !
1848 ! !
1849 ! external module variables: !
1850 ! ivflip : control flag of vertical index direction !
1851 ! =0: index from toa to surface !
1852 ! =1: index from surface to toa !
1853 ! !
1854 ! ==================== end of description ===================== !
1855 !
1856  implicit none
1857 
1858 ! --- inputs
1859  integer, intent(in) :: IX, NLAY, NLP1
1860 
1861  real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, &
1862  & tlyr, rhly, vvel
1863 
1864  real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, &
1865  & slmsk, cv, cvt, cvb
1866 
1867 ! --- outputs
1868  real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds
1869 
1870  real (kind=kind_phys), dimension(:,:), intent(out) :: clds
1871 
1872  integer, dimension(:,:), intent(out) :: mtop,mbot
1873 
1874 ! --- local variables:
1875  real (kind=kind_phys), dimension(IX,NLAY) :: cldtot, cldcnv, &
1876  & cldtau, taufac, dthdp, tem2d
1877 
1878  real (kind=kind_phys) :: ptop1(ix,nk_clds+1)
1879 
1880  real (kind=kind_phys) :: cr1, cr2, crk, pval, cval, omeg, value, &
1881  & tem1, tem2
1882 
1883  integer, dimension(IX):: idom, kcut
1884 
1885 ! --- for rh-cl calculation
1886  real (kind=kind_phys) :: xlatdg, xlondg, xlnn, xlss, xrgt, xlft, &
1887  & rhcla(NBIN,NLON,MCLD,NSEAL), rhcld(IX,NBIN,MCLD)
1888 
1889  integer :: ireg, ib, ic, id, id1, il, is, nhalf
1890 
1891  integer :: i, j, k, klowt
1892 ! integer :: klowb
1893 
1894  logical :: notstop
1895 
1896 !
1897 !===> ... begin here
1898 !
1899  clouds(:,:,:) = 0.0
1900 
1901  tem1 = 180.0 / con_pi
1902 
1903  lab_do_i_ix : do i = 1, ix
1904 
1905  xlatdg = xlat(i) * tem1 ! if xlat in pi/2 -> -pi/2 range
1906 ! xlatdg = 90.0 - xlat(i)*tem1 ! if xlat in 0 -> pi range
1907 
1908  xlondg = xlon(i) * tem1
1909  if (xlondg < 0.0) xlondg = xlondg + 360.0 ! if in -180->180, chg to 0->360 range
1910 
1911  ireg = 4
1912 
1914 ! --- get rh-cld relation for this lat
1915 
1916  lab_do_j : do j = 1, 3
1917  if (xlatdg > xlabdy(j)) then
1918  ireg = j
1919  exit lab_do_j
1920  endif
1921  enddo lab_do_j
1922 
1923  do is = 1, nseal
1924  do ic = 1, mcld
1925  do il = 1, nlon
1926  do ib = 1, nbin
1927  rhcla(ib,il,ic,is) = rhcl(ib,il,ireg,ic,is)
1928  enddo
1929  enddo
1930  enddo
1931  enddo
1932 
1934 ! --- linear transition between latitudinal regions...
1935  do j = 1, 3
1936  xlnn = xlabdy(j) + xlim
1937  xlss = xlabdy(j) - xlim
1938 
1939  if (xlatdg < xlnn .and. xlatdg > xlss) then
1940  do is = 1, nseal
1941  do ic = 1, mcld
1942  do il = 1, nlon
1943  do ib = 1, nbin
1944  rhcla(ib,il,ic,is) = rhcl(ib,il,j+1,ic,is) &
1945  & + (rhcl(ib,il,j,ic,is)-rhcl(ib,il,j+1,ic,is)) &
1946  & * (xlatdg-xlss) / (xlnn-xlss)
1947  enddo
1948  enddo
1949  enddo
1950  enddo
1951  endif
1952 
1953  enddo ! end_j_loop
1954 
1956 ! --- get rh-cld relationship for each grid point, interpolating
1957 ! longitudinally between regions if necessary..
1958 
1959  if (slmsk(i) < 1.0) then
1960  is = 2
1961  else
1962  is = 1
1963  endif
1964 
1965 ! --- which hemisphere (e,w)
1966 
1967  if (xlondg > 180.e0) then
1968  il = 2
1969  else
1970  il = 1
1971  endif
1972 
1973  do ic = 1, mcld
1974  do ib = 1, nbin
1975  rhcld(i,ib,ic) = rhcla(ib,il,ic,is)
1976  enddo
1977 
1978  lab_do_k : do k = 1, 3
1979  tem2 = abs(xlondg - xlobdy(k))
1980 
1981  if (tem2 < xlim) then
1982  id = il
1983  id1= id + 1
1984  if (id1 > nlon) id1 = 1
1985 
1986  xlft = xlobdy(k) - xlim
1987  xrgt = xlobdy(k) + xlim
1988 
1989  do ib = 1, nbin
1990  rhcld(i,ib,ic) = rhcla(ib,id1,ic,is) &
1991  & + (rhcla(ib,id,ic,is) - rhcla(ib,id1,ic,is)) &
1992  & * (xlondg-xrgt)/(xlft-xrgt)
1993  enddo
1994  exit lab_do_k
1995  endif
1996 
1997  enddo lab_do_k
1998 
1999  enddo ! end_do_ic_loop
2000  enddo lab_do_i_ix
2001 
2003 ! --- find top pressure for each cloud domain
2004 
2005  do j = 1, 4
2006  tem1 = ptopc(j,2) - ptopc(j,1)
2007 
2008  do i = 1, ix
2009  tem2 = xlat(i) / con_pi ! if xlat in pi/2 -> -pi/2 range
2010 ! tem2 = 0.5 - xlat(i)/con_pi ! if xlat in 0 -> pi range
2011 
2012  ptop1(i,j) = ptopc(j,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 )
2013  enddo
2014  enddo
2015 
2017 ! --- stratiform cloud optical depth
2018 
2019  do k = 1, nlay
2020  do i = 1, ix
2021  tem1 = tlyr(i,k) - con_ttp
2022  if (tem1 <= -10.0) then
2023  cldtau(i,k) = max( 0.1e-3, 2.0e-6*(tem1+82.5)**2 )
2024  else
2025  cldtau(i,k) = min( 0.08, 6.949e-3*tem1+0.08 )
2026  endif
2027  enddo
2028  enddo
2029 
2031 ! --- potential temperature and its lapse rate
2032 
2033  do k = 1, nlay
2034  do i = 1, ix
2035  cldtot(i,k) = 0.0
2036  cldcnv(i,k) = 0.0
2037  tem1 = (plyr(i,k)*0.001) ** (-con_rocp)
2038  tem2d(i,k) = tem1 * tlyr(i,k)
2039  enddo
2040  enddo
2041 
2042  do k = 1, nlay-1
2043  do i = 1, ix
2044  dthdp(i,k) = (tem2d(i,k+1)-tem2d(i,k))/(plyr(i,k+1)-plyr(i,k))
2045  enddo
2046  enddo
2047 
2049 !===> ... diagnostic method to find cloud amount cldtot, cldcnv
2050 !
2051 
2052  if ( ivflip == 0 ) then ! input data from toa to sfc
2055 ! --- find the lowest low cloud top sigma level, computed for each lat cause
2056 ! domain definition changes with latitude...
2057 
2058 ! klowb = 1
2059  klowt = 1
2060  do k = 1, nlay
2061  do i = 1, ix
2062 ! if (plvl(i,k) < ptop1(i,2)) klowb = k
2063  if (plvl(i,k) < ptop1(i,2)) klowt = max(klowt,k)
2064  taufac(i,k) = plvl(i,k+1) - plvl(i,k)
2065  enddo
2066  enddo
2067 
2068  do i = 1, ix
2069 
2073 ! --- find the stratosphere cut off layer for high cloud (about 250mb).
2074 ! it is assumed to be above the layer with dthdp less than -0.25 in
2075 ! the high cloud domain
2076 
2077  kcut(i) = 2
2078  lab_do_kcut0 : do k = klowt-1, 2, -1
2079  if (plyr(i,k) <= ptop1(i,3) .and. &
2080  & dthdp(i,k) < -0.25e0) then
2081  kcut(i) = k
2082  exit lab_do_kcut0
2083  endif
2084  enddo lab_do_kcut0
2085 
2087 ! --- put convective cloud into 'cldcnv', no merge at this point..
2088 
2089  if (cv(i) >= climit .and. cvt(i) < cvb(i)) then
2090  id = nlay
2091  id1 = nlay
2092 
2093  lab_do_k_cvt0 : do k = 2, nlay
2094  if (cvt(i) <= plyr(i,k)) then
2095  id = k - 1
2096  exit lab_do_k_cvt0
2097  endif
2098  enddo lab_do_k_cvt0
2099 
2100  lab_do_k_cvb0 : do k = nlay-1, 1, -1
2101  if (cvb(i) >= plyr(i,k)) then
2102  id1 = k + 1
2103  exit lab_do_k_cvb0
2104  endif
2105  enddo lab_do_k_cvb0
2106 
2107  tem1 = plyr(i,id1) - plyr(i,id)
2108  do k = id, id1
2109  cldcnv(i,k) = cv(i)
2110  taufac(i,k) = taufac(i,k) * max( 0.25, 1.0-0.125*tem1 )
2111  cldtau(i,k) = 0.06
2112  enddo
2113  endif
2114 
2115  enddo ! end_do_i_loop
2116 
2119 ! --- calculate stratiform cloud and put into array 'cldtot' using
2120 ! the cloud-rel.humidity relationship from table look-up..where
2121 ! tables obtained using k.mitchell frequency distribution tuning
2122 !bl (observations are daily means from us af rtneph).....k.a.c.
2123 !bl tables created without lowest 10 percent of atmos.....k.a.c.
2124 ! (observations are synoptic using -6,+3 window from rtneph)
2125 ! tables are created with lowest 10-percent-of-atmos, and are
2126 ! --- now used.. 25 october 1995 ... kac.
2127 
2128  do k = nlay-1, 2, -1
2129 
2130  if (k < llyr) then
2131  do i = 1, ix
2132  idom(i) = 0
2133  enddo
2134 
2135  do i = 1, ix
2136  lab_do_ic0 : do ic = 2, 4
2137  if(plyr(i,k) >= ptop1(i,ic)) then
2138  idom(i) = ic
2139  exit lab_do_ic0
2140  endif
2141  enddo lab_do_ic0
2142  enddo
2143  else
2144  do i = 1, ix
2145  idom(i) = 1
2146  enddo
2147  endif
2148 
2149  do i = 1, ix
2150  id = idom(i)
2151  nhalf = (nbin + 1) / 2
2152 
2153  if (id <= 0 .or. k < kcut(i)) then
2154  cldtot(i,k) = 0.0
2155  elseif (rhly(i,k) <= rhcld(i,1,id)) then
2156  cldtot(i,k) = 0.0
2157  elseif (rhly(i,k) >= rhcld(i,nbin,id)) then
2158  cldtot(i,k) = 1.0
2159  else
2160  ib = nhalf
2161  crk = rhly(i,k)
2162 
2163  notstop = .true.
2164  do while ( notstop )
2165  nhalf = (nhalf + 1) / 2
2166  cr1 = rhcld(i,ib, id)
2167  cr2 = rhcld(i,ib+1,id)
2168 
2169  if (crk <= cr1) then
2170  ib = max( ib-nhalf, 1 )
2171  elseif (crk > cr2) then
2172  ib = min( ib+nhalf, nbin-1 )
2173  else
2174  cldtot(i,k) = 0.01 * (ib + (crk - cr1)/(cr2 - cr1))
2175  notstop = .false.
2176  endif
2177  enddo ! end_do_while
2178  endif
2179  enddo ! end_do_i_loop
2180 
2181  enddo ! end_do_k_loop
2182 
2184 ! --- vertical velocity adjustment on low clouds
2185 
2186  value = vvcld1 - vvcld2
2187  do k = klowt, llyr+1
2188  do i = 1, ix
2189 
2190  omeg = vvel(i,k)
2191  cval = cldtot(i,k)
2192  pval = plyr(i,k)
2193 
2194 ! --- vertical velocity adjustment on low clouds
2195 
2196  if (cval >= climit .and. pval >= ptop1(i,2)) then
2197  if (omeg >= vvcld1) then
2198  cldtot(i,k) = 0.0
2199  elseif (omeg > vvcld2) then
2200  tem1 = (vvcld1 - omeg) / value
2201  cldtot(i,k) = cldtot(i,k) * sqrt(tem1)
2202  endif
2203  endif
2204 
2205  enddo ! end_do_i_loop
2206  enddo ! end_do_k_loop
2207 
2208  else ! input data from sfc to toa
2209 
2210 ! --- find the lowest low cloud top sigma level, computed for each lat cause
2211 ! domain definition changes with latitude...
2212 
2213 ! klowb = NLAY
2214  klowt = nlay
2215  do k = nlay, 1, -1
2216  do i = 1, ix
2217 ! if (plvl(i,k) < ptop1(i,2)) klowb = k
2218  if (plvl(i,k) < ptop1(i,2)) klowt = min(klowt,k)
2219  taufac(i,k) = plvl(i,k) - plvl(i,k+1) ! dp for later cal cldtau use
2220  enddo
2221  enddo
2222 
2223  do i = 1, ix
2224 
2225 ! --- find the stratosphere cut off layer for high cloud (about 250mb).
2226 ! it is assumed to be above the layer with dthdp less than -0.25 in
2227 ! the high cloud domain
2228 
2229  kcut(i) = nlay - 1
2230  lab_do_kcut1 : do k = klowt+1, nlay-1
2231  if (plyr(i,k) <= ptop1(i,3) .and. &
2232  & dthdp(i,k) < -0.25e0) then
2233  kcut(i) = k
2234  exit lab_do_kcut1
2235  endif
2236  enddo lab_do_kcut1
2237 
2238 ! --- put convective cloud into 'cldcnv', no merge at this point..
2239 
2240  if (cv(i) >= climit .and. cvt(i) < cvb(i)) then
2241  id = 1
2242  id1 = 1
2243 
2244  lab_do_k_cvt : do k = nlay-1, 1, -1
2245  if (cvt(i) <= plyr(i,k)) then
2246  id = k + 1
2247  exit lab_do_k_cvt
2248  endif
2249  enddo lab_do_k_cvt
2250 
2251  lab_do_k_cvb : do k = 2, nlay
2252  if (cvb(i) >= plyr(i,k)) then
2253  id1 = k - 1
2254  exit lab_do_k_cvb
2255  endif
2256  enddo lab_do_k_cvb
2257 
2258  tem1 = plyr(i,id1) - plyr(i,id)
2259  do k = id1, id
2260  cldcnv(i,k) = cv(i)
2261  taufac(i,k) = taufac(i,k) * max( 0.25, 1.0-0.125*tem1 )
2262  cldtau(i,k) = 0.06
2263  enddo
2264  endif
2265 
2266  enddo ! end_do_i_loop
2267 
2268 ! --- calculate stratiform cloud and put into array 'cldtot' using
2269 ! the cloud-rel.humidity relationship from table look-up..where
2270 ! tables obtained using k.mitchell frequency distribution tuning
2271 !bl (observations are daily means from us af rtneph).....k.a.c.
2272 !bl tables created without lowest 10 percent of atmos.....k.a.c.
2273 ! (observations are synoptic using -6,+3 window from rtneph)
2274 ! tables are created with lowest 10-percent-of-atmos, and are
2275 ! --- now used.. 25 october 1995 ... kac.
2276 
2277  do k = 2, nlay-1
2278 
2279  if (k > llyr) then
2280  do i = 1, ix
2281  idom(i) = 0
2282  enddo
2283 
2284  do i = 1, ix
2285  lab_do_ic1 : do ic = 2, 4
2286  if(plyr(i,k) >= ptop1(i,ic)) then
2287  idom(i) = ic
2288  exit lab_do_ic1
2289  endif
2290  enddo lab_do_ic1
2291  enddo
2292  else
2293  do i = 1, ix
2294  idom(i) = 1
2295  enddo
2296  endif
2297 
2298  do i = 1, ix
2299  id = idom(i)
2300  nhalf = (nbin + 1) / 2
2301 
2302  if (id <= 0 .or. k > kcut(i)) then
2303  cldtot(i,k) = 0.0
2304  elseif (rhly(i,k) <= rhcld(i,1,id)) then
2305  cldtot(i,k) = 0.0
2306  elseif (rhly(i,k) >= rhcld(i,nbin,id)) then
2307  cldtot(i,k) = 1.0
2308  else
2309  ib = nhalf
2310  crk = rhly(i,k)
2311 
2312  notstop = .true.
2313  do while ( notstop )
2314  nhalf = (nhalf + 1) / 2
2315  cr1 = rhcld(i,ib, id)
2316  cr2 = rhcld(i,ib+1,id)
2317 
2318  if (crk <= cr1) then
2319  ib = max( ib-nhalf, 1 )
2320  elseif (crk > cr2) then
2321  ib = min( ib+nhalf, nbin-1 )
2322  else
2323  cldtot(i,k) = 0.01 * (ib + (crk - cr1)/(cr2 - cr1))
2324  notstop = .false.
2325  endif
2326  enddo ! end_do_while
2327  endif
2328  enddo ! end_do_i_loop
2329 
2330  enddo ! end_do_k_loop
2331 
2332 ! --- vertical velocity adjustment on low clouds
2333 
2334  value = vvcld1 - vvcld2
2335  do k = llyr-1, klowt
2336  do i = 1, ix
2337 
2338  omeg = vvel(i,k)
2339  cval = cldtot(i,k)
2340  pval = plyr(i,k)
2341 
2342 ! --- vertical velocity adjustment on low clouds
2343 
2344  if (cval >= climit .and. pval >= ptop1(i,2)) then
2345  if (omeg >= vvcld1) then
2346  cldtot(i,k) = 0.0
2347  elseif (omeg > vvcld2) then
2348  tem1 = (vvcld1 - omeg) / value
2349  cldtot(i,k) = cldtot(i,k) * sqrt(tem1)
2350  endif
2351  endif
2352 
2353  enddo ! end_do_i_loop
2354  enddo ! end_do_k_loop
2355 
2356  endif ! end_if_ivflip
2358 ! --- diagnostic cloud optical depth
2359 ! cldtau = cldtau * taufac
2360 
2361  where (cldtot < climit)
2362  cldtot = 0.0
2363  endwhere
2364  where (cldcnv < climit)
2365  cldcnv = 0.0
2366  endwhere
2367 
2368  where (cldtot < climit .and. cldcnv < climit)
2369  cldtau = 0.0
2370  endwhere
2371 
2372  do k = 1, nlay
2373  do i = 1, ix
2374  clouds(i,k,1) = max(cldtot(i,k), cldcnv(i,k))
2375  clouds(i,k,2) = cldtau(i,k) * taufac(i,k)
2376  clouds(i,k,3) = cldssa_def
2377  clouds(i,k,4) = cldasy_def
2378  enddo
2379  enddo
2380 
2383 !
2384 !===> ... compute low, mid, high, total, and boundary layer cloud fractions
2385 ! and clouds top/bottom layer indices for low, mid, and high clouds.
2386 ! the three cloud domain boundaries are defined by ptopc. the cloud
2387 ! overlapping method is defined by control flag 'iovr', which may
2388 ! be different for lw and sw radiation programs.
2389 
2390  call gethml &
2391 ! --- inputs:
2392  & ( plyr, ptop1, cldtot, cldcnv, &
2393  & ix, nlay, &
2394 ! --- outputs:
2395  & clds, mtop, mbot &
2396  & )
2397 
2398 !
2399  return
2400 !...................................
2401  end subroutine diagcld1
2402 !-----------------------------------
2404 
2407 !----------------------------------- !
2408  subroutine gethml &
2409 !................................... !
2410 
2411 ! --- inputs:
2412  & ( plyr, ptop1, cldtot, cldcnv, &
2413  & ix, nlay, &
2414 ! --- outputs:
2415  & clds, mtop, mbot &
2416  & )
2417 
2418 ! =================================================================== !
2419 ! !
2420 ! abstract: compute high, mid, low, total, and boundary cloud fractions !
2421 ! and cloud top/bottom layer indices for model diagnostic output. !
2422 ! the three cloud domain boundaries are defined by ptopc. the cloud !
2423 ! overlapping method is defined by control flag 'iovr', which is also !
2424 ! used by lw and sw radiation programs. !
2425 ! !
2426 ! usage: call gethml !
2427 ! !
2428 ! subprograms called: none !
2429 ! !
2430 ! attributes: !
2431 ! language: fortran 90 !
2432 ! machine: ibm-sp, sgi !
2433 ! !
2434 ! !
2435 ! ==================== definition of variables ==================== !
2436 ! !
2437 ! input variables: !
2438 ! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) !
2439 ! ptop1 (IX,4) : pressure limits of cloud domain interfaces !
2440 ! (sfc,low,mid,high) in mb (100Pa) !
2441 ! cldtot(IX,NLAY) : total or straiform cloud profile in fraction !
2442 ! cldcnv(IX,NLAY) : convective cloud (for diagnostic scheme only) !
2443 ! IX : horizontal dimention !
2444 ! NLAY : vertical layer dimensions !
2445 ! !
2446 ! output variables: !
2447 ! clds (IX,5) : fraction of clouds for low, mid, hi, tot, bl !
2448 ! mtop (IX,3) : vertical indices for low, mid, hi cloud tops !
2449 ! mbot (IX,3) : vertical indices for low, mid, hi cloud bases !
2450 ! !
2451 ! external module variables: (in physparam) !
2452 ! ivflip : control flag of vertical index direction !
2453 ! =0: index from toa to surface !
2454 ! =1: index from surface to toa !
2455 ! !
2456 ! internal module variables: !
2457 ! iovr : control flag for cloud overlap !
2458 ! =0 random overlapping clouds !
2459 ! =1 max/ran overlapping clouds !
2460 ! !
2461 ! !
2462 ! ==================== end of description ===================== !
2463 !
2464  implicit none!
2465 
2466 ! --- inputs:
2467  integer, intent(in) :: IX, NLAY
2468 
2469  real (kind=kind_phys), dimension(:,:), intent(in) :: plyr, ptop1, &
2470  & cldtot, cldcnv
2471 
2472 ! --- outputs
2473  real (kind=kind_phys), dimension(:,:), intent(out) :: clds
2474 
2475  integer, dimension(:,:), intent(out) :: mtop, mbot
2476 
2477 ! --- local variables:
2478  real (kind=kind_phys) :: cl1(ix), cl2(ix)
2479  real (kind=kind_phys) :: pcur, pnxt, ccur, cnxt
2480 
2481  integer, dimension(IX):: idom, kbt1, kth1, kbt2, kth2
2482  integer :: i, k, id, id1, kstr, kend, kinc
2483 
2484 !
2485 !===> ... begin here
2486 !
2487  clds(:,:) = 0.0
2488 
2489  do i = 1, ix
2490  cl1(i) = 1.0
2491  cl2(i) = 1.0
2492  enddo
2493 
2494 ! --- total and bl clouds, where cl1, cl2 are fractions of clear-sky view
2495 ! layer processed from surface and up
2496 
2497  if ( ivflip == 0 ) then ! input data from toa to sfc
2498  kstr = nlay
2499  kend = 1
2500  kinc = -1
2501  else ! input data from sfc to toa
2502  kstr = 1
2503  kend = nlay
2504  kinc = 1
2505  endif ! end_if_ivflip
2506 
2507  if ( iovr == 0 ) then ! random overlap
2508 
2509  do k = kstr, kend, kinc
2510  do i = 1, ix
2511  ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
2512  if (ccur >= climit) cl1(i) = cl1(i) * (1.0 - ccur)
2513  enddo
2514 
2515  if (k == llyr) then
2516  do i = 1, ix
2517  clds(i,5) = 1.0 - cl1(i) ! save bl cloud
2518  enddo
2519  endif
2520  enddo
2521 
2522  do i = 1, ix
2523  clds(i,4) = 1.0 - cl1(i) ! save total cloud
2524  enddo
2525 
2526  else ! max/ran overlap
2527 
2528  do k = kstr, kend, kinc
2529  do i = 1, ix
2530  ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
2531  if (ccur >= climit) then ! cloudy layer
2532  cl2(i) = min( cl2(i), (1.0 - ccur) )
2533  else ! clear layer
2534  cl1(i) = cl1(i) * cl2(i)
2535  cl2(i) = 1.0
2536  endif
2537  enddo
2538 
2539  if (k == llyr) then
2540  do i = 1, ix
2541  clds(i,5) = 1.0 - cl1(i) * cl2(i) ! save bl cloud
2542  enddo
2543  endif
2544  enddo
2545 
2546  do i = 1, ix
2547  clds(i,4) = 1.0 - cl1(i) * cl2(i) ! save total cloud
2548  enddo
2549 
2550  endif ! end_if_iovr
2551 
2552 ! --- high, mid, low clouds, where cl1, cl2 are cloud fractions
2553 ! layer processed from one layer below llyr and up
2554 ! --- change! layer processed from surface to top, so low clouds will
2555 ! contains both bl and low clouds.
2556 
2557  if ( ivflip == 0 ) then ! input data from toa to sfc
2558 
2559  do i = 1, ix
2560  cl1(i) = 0.0
2561  cl2(i) = 0.0
2562  kbt1(i) = nlay
2563  kbt2(i) = nlay
2564  kth1(i) = 0
2565  kth2(i) = 0
2566  idom(i) = 1
2567  mbot(i,1) = nlay
2568  mtop(i,1) = nlay
2569  mbot(i,2) = nlay - 1
2570  mtop(i,2) = nlay - 1
2571  mbot(i,3) = nlay - 1
2572  mtop(i,3) = nlay - 1
2573  enddo
2574 
2575 !org do k = llyr-1, 1, -1
2576  do k = nlay, 1, -1
2577  do i = 1, ix
2578  id = idom(i)
2579  id1= id + 1
2580 
2581  pcur = plyr(i,k)
2582  ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
2583 
2584  if (k > 1) then
2585  pnxt = plyr(i,k-1)
2586  cnxt = min( ovcst, max( cldtot(i,k-1), cldcnv(i,k-1) ))
2587  else
2588  pnxt = -1.0
2589  cnxt = 0.0
2590  endif
2591 
2592  if (pcur < ptop1(i,id1)) then
2593  id = id + 1
2594  id1= id1 + 1
2595  idom(i) = id
2596  endif
2597 
2598  if (ccur >= climit) then
2599  if (kth2(i) == 0) kbt2(i) = k
2600  kth2(i) = kth2(i) + 1
2601 
2602  if ( iovr == 0 ) then
2603  cl2(i) = cl2(i) + ccur - cl2(i)*ccur
2604  else
2605  cl2(i) = max( cl2(i), ccur )
2606  endif
2607 
2608  if (cnxt < climit .or. pnxt < ptop1(i,id1)) then
2609  kbt1(i) = nint( (cl1(i)*kbt1(i) + cl2(i)*kbt2(i) ) &
2610  & / (cl1(i) + cl2(i)) )
2611  kth1(i) = nint( (cl1(i)*kth1(i) + cl2(i)*kth2(i) ) &
2612  & / (cl1(i) + cl2(i)) )
2613  cl1(i) = cl1(i) + cl2(i) - cl1(i)*cl2(i)
2614 
2615  kbt2(i) = k - 1
2616  kth2(i) = 0
2617  cl2(i) = 0.0
2618  endif ! end_if_cnxt_or_pnxt
2619  endif ! end_if_ccur
2620 
2621  if (pnxt < ptop1(i,id1)) then
2622  clds(i,id) = cl1(i)
2623  mtop(i,id) = min( kbt1(i), kbt1(i)-kth1(i)+1 )
2624  mbot(i,id) = kbt1(i)
2625 
2626  cl1(i) = 0.0
2627  kbt1(i) = k - 1
2628  kth1(i) = 0
2629 
2630  if (id1 <= nk_clds) then
2631  mbot(i,id1) = kbt1(i)
2632  mtop(i,id1) = kbt1(i)
2633  endif
2634  endif ! end_if_pnxt
2635 
2636  enddo ! end_do_i_loop
2637  enddo ! end_do_k_loop
2638 
2639  else ! input data from sfc to toa
2640 
2641  do i = 1, ix
2642  cl1(i) = 0.0
2643  cl2(i) = 0.0
2644  kbt1(i) = 1
2645  kbt2(i) = 1
2646  kth1(i) = 0
2647  kth2(i) = 0
2648  idom(i) = 1
2649  mbot(i,1) = 1
2650  mtop(i,1) = 1
2651  mbot(i,2) = 2
2652  mtop(i,2) = 2
2653  mbot(i,3) = 2
2654  mtop(i,3) = 2
2655  enddo
2656 
2657 !org do k = llyr+1, NLAY
2658  do k = 1, nlay
2659  do i = 1, ix
2660  id = idom(i)
2661  id1= id + 1
2662 
2663  pcur = plyr(i,k)
2664  ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
2665 
2666  if (k < nlay) then
2667  pnxt = plyr(i,k+1)
2668  cnxt = min( ovcst, max( cldtot(i,k+1), cldcnv(i,k+1) ))
2669  else
2670  pnxt = -1.0
2671  cnxt = 0.0
2672  endif
2673 
2674  if (pcur < ptop1(i,id1)) then
2675  id = id + 1
2676  id1= id1 + 1
2677  idom(i) = id
2678  endif
2679 
2680  if (ccur >= climit) then
2681  if (kth2(i) == 0) kbt2(i) = k
2682  kth2(i) = kth2(i) + 1
2683 
2684  if ( iovr == 0 ) then
2685  cl2(i) = cl2(i) + ccur - cl2(i)*ccur
2686  else
2687  cl2(i) = max( cl2(i), ccur )
2688  endif
2689 
2690  if (cnxt < climit .or. pnxt < ptop1(i,id1)) then
2691  kbt1(i) = nint( (cl1(i)*kbt1(i) + cl2(i)*kbt2(i)) &
2692  & / (cl1(i) + cl2(i)) )
2693  kth1(i) = nint( (cl1(i)*kth1(i) + cl2(i)*kth2(i)) &
2694  & / (cl1(i) + cl2(i)) )
2695  cl1(i) = cl1(i) + cl2(i) - cl1(i)*cl2(i)
2696 
2697  kbt2(i) = k + 1
2698  kth2(i) = 0
2699  cl2(i) = 0.0
2700  endif ! end_if_cnxt_or_pnxt
2701  endif ! end_if_ccur
2702 
2703  if (pnxt < ptop1(i,id1)) then
2704  clds(i,id) = cl1(i)
2705  mtop(i,id) = max( kbt1(i), kbt1(i)+kth1(i)-1 )
2706  mbot(i,id) = kbt1(i)
2707 
2708  cl1(i) = 0.0
2709  kbt1(i) = min(k+1, nlay)
2710  kth1(i) = 0
2711 
2712  if (id1 <= nk_clds) then
2713  mbot(i,id1) = kbt1(i)
2714  mtop(i,id1) = kbt1(i)
2715  endif
2716  endif ! end_if_pnxt
2717 
2718  enddo ! end_do_i_loop
2719  enddo ! end_do_k_loop
2720 
2721  endif ! end_if_ivflip
2722 
2723 !
2724  return
2725 !...................................
2726  end subroutine gethml
2727 !-----------------------------------
2728 
2730 !----------------------------------- !
2731  subroutine rhtable &
2732 !................................... !
2733 
2734 ! --- inputs:
2735  & ( me &
2736 ! --- outputs:
2737  &, ier )
2738 
2739 ! =================================================================== !
2740 ! !
2741 ! abstract: cld-rh relations obtained from mitchell-hahn procedure, !
2742 ! here read cld/rh tuning tables for day 0,1,...,5 and merge into 1 !
2743 ! file. !
2744 ! !
2745 ! inputs: !
2746 ! me : check print control flag !
2747 ! !
2748 ! outputs: !
2749 ! ier : error flag !
2750 ! !
2751 ! =================================================================== !
2752 !
2753  implicit none!
2754 
2755 ! --- inputs:
2756  integer, intent(in) :: me
2757 
2758 ! --- output:
2759  integer, intent(out) :: ier
2760 
2761 ! --- locals:
2762  real (kind=kind_phys), dimension(NBIN,NLON,NLAT,MCLD,NSEAL) :: &
2763  & rhfd, rtnfd, rhcf, rtncf, rhcla
2764 
2765  real (kind=kind_io4), dimension(NBIN,NLON,NLAT,MCLD,NSEAL) :: &
2766  & rhfd4, rtnfd4
2767 
2768  real(kind=kind_io4) :: fhour
2769 
2770  real(kind=kind_phys) :: binscl, cfrac, clsat, rhsat, cstem
2771 
2772  integer, dimension(NLON,NLAT,MCLD,NSEAL) :: kpts, kkpts
2773 
2774  integer :: icdays(15), idate(4), nbdayi, isat
2775 
2776  integer :: i, i1, j, k, l, m, id, im, iy
2777 
2778 !
2779 !===> ... begin here
2780 !
2781 
2782  ier = 1
2783 
2784  rewind nicltun
2785 
2786  binscl = 1.0 / nbin
2787 
2788 ! --- array initializations
2789 
2790  do m=1,nseal
2791  do l=1,mcld
2792  do k=1,nlat
2793  do j=1,nlon
2794  do i=1,nbin
2795  rhcf(i,j,k,l,m) = 0.0
2796  rtncf(i,j,k,l,m) = 0.0
2797  rhcla(i,j,k,l,m) = -0.1
2798  enddo
2799  enddo
2800  enddo
2801  enddo
2802  enddo
2803 
2804  kkpts = 0
2805 
2806 ! --- read the data off the rotating file
2807 
2808  read (nicltun,err=998,end=999) nbdayi, icdays
2809 
2810  if (me == 0) print 11, nbdayi
2811  11 format(' from rhtable DAYS ON FILE =',i5)
2812 
2813  do i = 1, nbdayi
2814  id = icdays(i) / 10000
2815  im = (icdays(i)-id*10000) / 100
2816  iy = icdays(i)-id*10000-im*100
2817  if (me == 0) print 51, id,im,iy
2818  51 format(' from rhtable ARCHV DATA FROM DA,MO,YR=',3i4)
2819  enddo
2820 
2821  read (nicltun,err=998,end=999) fhour,idate
2822 
2823  do i1 = 1, nbdayi
2824  read (nicltun) rhfd4
2825  rhfd = rhfd4
2826 
2827  read (nicltun) rtnfd4
2828  rtnfd = rtnfd4
2829 
2830  read (nicltun) kpts
2831 
2832  do m=1,nseal
2833  do l=1,mcld
2834  do k=1,nlat
2835  do j=1,nlon
2836  do i=1,nbin
2837  rhcf(i,j,k,l,m) = rhcf(i,j,k,l,m) + rhfd(i,j,k,l,m)
2838  rtncf(i,j,k,l,m) = rtncf(i,j,k,l,m) + rtnfd(i,j,k,l,m)
2839  enddo
2840  enddo
2841  enddo
2842  enddo
2843  enddo
2844 
2845  kkpts = kkpts + kpts
2846 
2847  enddo ! end_do_i1_loop
2848 
2849  do m = 1, nseal
2850  do l = 1, mcld
2851  do k = 1, nlat
2852  do j = 1, nlon
2853 
2854 ! --- compute the cumulative frequency distribution
2855 
2856  do i = 2, nbin
2857  rhcf(i,j,k,l,m) = rhcf(i-1,j,k,l,m) + rhcf(i,j,k,l,m)
2858  rtncf(i,j,k,l,m) = rtncf(i-1,j,k,l,m) + rtncf(i,j,k,l,m)
2859  enddo ! end_do_i_loop
2860 
2861  if (kkpts(j,k,l,m) > 0) then
2862  do i = 1, nbin
2863  rhcf(i,j,k,l,m)= rhcf(i,j,k,l,m)/kkpts(j,k,l,m)
2864  rtncf(i,j,k,l,m)=min(1., rtncf(i,j,k,l,m)/kkpts(j,k,l,m))
2865  enddo
2866 
2867 ! --- cause we mix calculations of rh retune with cray and ibm words
2868 ! the last value of rhcf is close to but ne 1.0,
2869 ! --- so we reset it in order that the 360 loop gives complete tabl
2870 
2871  rhcf(nbin,j,k,l,m) = 1.0
2872 
2873  do i = 1, nbin
2874  lab_do_i1 : do i1 = 1, nbin
2875  if (rhcf(i1,j,k,l,m) >= rtncf(i,j,k,l,m)) then
2876  rhcla(i,j,k,l,m) = i1 * binscl
2877  exit lab_do_i1
2878  endif
2879  enddo lab_do_i1
2880  enddo
2881 
2882  else ! if_kkpts
2883 ! --- no critical rh
2884 
2885  do i = 1, nbin
2886  rhcf(i,j,k,l,m) = -0.1
2887  rtncf(i,j,k,l,m) = -0.1
2888  enddo
2889 
2890  if (me == 0) then
2891  print 210, k,j,m
2892  210 format(' NO CRIT RH FOR LAT=',i3,' AND LON BAND=',i3, &
2893  & ' LAND(=1) SEA=',i3//' MODEL RH ',' OBS RTCLD')
2894  do i = 1, nbin
2895  print 203, rhcf(i,j,k,l,m), rtncf(i,j,k,l,m)
2896  203 format(2f10.2)
2897  enddo
2898  endif
2899 
2900  endif ! if_kkpts
2901 
2902  enddo ! end_do_j_loop
2903  enddo ! end_do_k_loop
2904  enddo ! end_do_l_loop
2905  enddo ! end_do_m_loop
2906 
2907  do m = 1, nseal
2908  do l = 1, mcld
2909  do k = 1, nlat
2910  do j = 1, nlon
2911 
2912  isat = 0
2913  do i = 1, nbin-1
2914  cfrac = binscl * (i - 1)
2915 
2916  if (rhcla(i,j,k,l,m) < 0.0) then
2917  print 1941, i,m,l,k,j
2918  1941 format(' NEG RHCLA FOR IT,NSL,NC,LAT,LON=',5i4 &
2919  &, '...STOPPP..')
2920  stop
2921  endif
2922 
2923  if (rtncf(i,j,k,l,m) >= 1.0) then
2924  if (isat <= 0) then
2925  isat = i
2926  rhsat = rhcla(i,j,k,l,m)
2927  clsat = cfrac
2928  endif
2929 
2930  rhcla(i,j,k,l,m) = rhsat + (1.0 - rhsat) &
2931  & * (cfrac - clsat) / (1.0 - clsat)
2932  endif
2933  enddo
2934 
2935  rhcla(nbin,j,k,l,m) = 1.0
2936 
2937  enddo ! end_do_j_loop
2938  enddo ! end_do_k_loop
2939  enddo ! end_do_l_loop
2940  enddo ! end_do_m_loop
2941 
2942 ! --- smooth out the table as it reaches rh=1.0, via linear interpolation
2943 ! between location of rh ge .98 and the NBIN bin (where rh=1.0)
2944 ! previously rh=1.0 occurred for many of the latter bins in the
2945 ! --- table, thereby giving a cloud value of less then 1.0 for rh=1.0
2946 
2947  rhcl = rhcla
2948 
2949  do m = 1, nseal
2950  do l = 1, mcld
2951  do k = 1, nlat
2952  do j = 1, nlon
2953 
2954  lab_do_i : do i = 1, nbin - 2
2955  cfrac = binscl * i
2956 
2957  if (rhcla(i,j,k,l,m) >= 0.98) then
2958  do i1 = i, nbin
2959  cstem = binscl * i1
2960 
2961  rhcl(i1,j,k,l,m) = rhcla(i,j,k,l,m) &
2962  & + (rhcla(nbin,j,k,l,m) - rhcla(i,j,k,l,m)) &
2963  & * (cstem - cfrac) / (1.0 - cfrac)
2964  enddo
2965  exit lab_do_i
2966  endif
2967  enddo lab_do_i
2968 
2969  enddo ! end_do_j_loop
2970  enddo ! end_do_k_loop
2971  enddo ! end_do_l_loop
2972  enddo ! end_do_m_loop
2973 
2974  if (me == 0) then
2975  print *,'completed rhtable for cloud tuninig tables'
2976  endif
2977  return
2978 
2979  998 print 988
2980  988 format(' from rhtable ERROR READING TABLES')
2981  ier = -1
2982  return
2983 
2984  999 print 989
2985  989 format(' from rhtable E.O.F READING TABLES')
2986  ier = -1
2987  return
2988 
2989 !...................................
2990  end subroutine rhtable
2991 !-----------------------------------
2992 
2993 
2994 !
2995 !........................................!
2996  end module module_radiation_clouds !
2997 !========================================!
2998 
module_radiation_clouds::progcld1
subroutine, public progcld1
This subroutine computes cloud related quantities using zhao/moorthi&#39;s prognostic cloud microphysics ...
Definition: radiation_clouds.f:384
module_radiation_clouds::rhcl
real(kind=kind_phys), dimension(nbin, nlon, nlat, mcld, nseal) rhcl
Definition: radiation_clouds.f:208
module_radiation_clouds::nbin
integer, parameter nbin
Definition: radiation_clouds.f:185
module_radiation_clouds::ptopc
real(kind=kind_phys), dimension(nk_clds+1, 2), save ptopc
Definition: radiation_clouds.f:169
module_radiation_clouds::mcld
integer, parameter mcld
Definition: radiation_clouds.f:188
module_microphysics::rsipath2
subroutine rsipath2
This subroutine is a modified version of ferrier&#39;s original "rsipath" subprogram. It computes layer&#39;s...
Definition: module_bfmicrophysics.f:2927
physcons::con_t0c
real(kind=kind_phys), parameter con_t0c
temp at 0C (K)
Definition: physcons.f:89
module_radiation_clouds::climit
real(kind=kind_phys), parameter climit
Definition: radiation_clouds.f:175
module_radiation_clouds::iovr
integer iovr
Definition: radiation_clouds.f:211
module_radiation_clouds::rrain_def
real(kind=kind_phys), parameter rrain_def
Definition: radiation_clouds.f:181
module_radiation_clouds::nk_clds
integer, parameter, public nk_clds
number of cloud vertical domains
Definition: radiation_clouds.f:166
module_radiation_clouds
This module computes cloud related quantities for radiation computations.
Definition: radiation_clouds.f:139
physparam::iovrsw
integer, save iovrsw
cloud overlapping control flag for sw
Definition: physparam.f:197
module_radiation_clouds::rsnow_def
real(kind=kind_phys), parameter rsnow_def
Definition: radiation_clouds.f:182
physcons
This module contains some the most frequently used math and physics constants for gcm models...
Definition: physcons.f:29
module_radiation_clouds::xlim
real(kind=kind_phys), parameter xlim
Definition: radiation_clouds.f:196
module_radiation_clouds::reice_def
real(kind=kind_phys), parameter reice_def
Definition: radiation_clouds.f:180
module_radiation_clouds::cldasy_def
real(kind=kind_phys), parameter cldasy_def
Definition: radiation_clouds.f:192
physcons::con_ttp
real(kind=kind_phys), parameter con_ttp
temp at H2O 3pt (K)
Definition: physcons.f:91
physparam::iovrlw
integer, save iovrlw
cloud overlapping control flag for lw
Definition: physparam.f:199
physcons::con_thgni
real(kind=kind_phys), parameter con_thgni
temperature the H.G.Nuc. ice starts
Definition: physcons.f:145
module_radiation_clouds::vtagcld
character(40), parameter vtagcld
Definition: radiation_clouds.f:156
physparam
This module defines commonly used control variables/parameters in physics related programs...
Definition: physparam.f:27
module_radiation_clouds::xlobdy
real(kind=kind_phys), dimension(3) xlobdy
Definition: radiation_clouds.f:197
module_radiation_clouds::cldssa_def
real(kind=kind_phys), parameter cldssa_def
Definition: radiation_clouds.f:191
physcons::con_rd
real(kind=kind_phys), parameter con_rd
gas constant air (J/kg/K)
Definition: physcons.f:69
module_radiation_clouds::gord
real(kind=kind_phys), parameter gord
Definition: radiation_clouds.f:161
module_radiation_clouds::llyr
integer llyr
Definition: radiation_clouds.f:210
module_radiation_clouds::rhtable
subroutine rhtable
cld-rh relations obtained from mitchell-hahn procedure.
Definition: radiation_clouds.f:2733
physparam::icmphys
integer, save icmphys
cloud microphysics scheme control flag
Definition: physparam.f:195
module_radiation_clouds::nlon
integer, parameter nlon
Definition: radiation_clouds.f:186
physparam::lnoprec
logical, save lnoprec
precip effect on radiation flag (ferrier microphysics)
Definition: physparam.f:205
module_radiation_clouds::reliq_def
real(kind=kind_phys), parameter reliq_def
Definition: radiation_clouds.f:179
physcons::con_rocp
real(kind=kind_phys), parameter con_rocp
Definition: physcons.f:104
physparam::lcrick
logical, save lcrick
eliminating CRICK control flag
Definition: physparam.f:201
physparam::ivflip
integer, save ivflip
vertical profile indexing flag
Definition: physparam.f:224
module_radiation_clouds::vvcld1
real(kind=kind_phys), parameter vvcld1
Definition: radiation_clouds.f:202
module_radiation_clouds::vvcld2
real(kind=kind_phys), parameter vvcld2
Definition: radiation_clouds.f:203
module_radiation_clouds::climit2
real(kind=kind_phys), parameter climit2
Definition: radiation_clouds.f:175
module_radiation_clouds::xlabdy
real(kind=kind_phys), dimension(3) xlabdy
Definition: radiation_clouds.f:197
module_radiation_clouds::nseal
integer, parameter nseal
Definition: radiation_clouds.f:189
module_radiation_clouds::gethml
subroutine gethml
This subroutine computes high, mid, low, total, and boundary cloud fractions and cloud top/bottom lay...
Definition: radiation_clouds.f:2410
module_microphysics
Definition: module_bfmicrophysics.f:1
physcons::con_g
real(kind=kind_phys), parameter con_g
gravity (m/s2)
Definition: physcons.f:52
module_radiation_clouds::gfac
real(kind=kind_phys), parameter gfac
Definition: radiation_clouds.f:161
physcons::con_pi
real(kind=kind_phys), parameter con_pi
Definition: physcons.f:41
physparam::lsashal
logical, save lsashal
shallow convection flag
Definition: physparam.f:207
module_radiation_clouds::nf_clds
integer, parameter, public nf_clds
number of fields in cloud array
Definition: radiation_clouds.f:164
module_radiation_clouds::ovcst
real(kind=kind_phys), parameter ovcst
Definition: radiation_clouds.f:176
module_radiation_clouds::cld_init
subroutine, public cld_init
This subroutine is an initialization program for cloud-radiation calculations. It sets up boundary la...
Definition: radiation_clouds.f:231
module_radiation_clouds::diagcld1
subroutine, public diagcld1
This subroutine computes cloud fractions for radiation calculations.
Definition: radiation_clouds.f:1789
physcons::con_fvirt
real(kind=kind_phys), parameter con_fvirt
Definition: physcons.f:107
module_radiation_clouds::progcld2
subroutine, public progcld2
This subroutine computes cloud related quantities using ferrier&#39;s prognostic cloud microphysics schem...
Definition: radiation_clouds.f:844
module_radiation_clouds::progcld3
subroutine, public progcld3
This subroutine computes cloud related quantities using zhao/moorthi&#39;s prognostic cloud microphysics ...
Definition: radiation_clouds.f:1357
module_radiation_clouds::nlat
integer, parameter nlat
Definition: radiation_clouds.f:187
physparam::lcnorm
logical, save lcnorm
in-cld condensate control flag
Definition: physparam.f:203
physparam::icldflg
integer, save icldflg
cloud optical property scheme control flag
Definition: physparam.f:193