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