224 module module_radiation_clouds
229 &
ivflip, kind_phys, kind_io4
230 use physcons
, only : con_fvirt,
con_ttp, con_rocp, &
233 use module_microphysics
, only : rsipath2
234 use module_iounitdef
, only : nicltun
241 character(40),
parameter :: &
242 & VTAGCLD=
'NCEP-Radiation_clouds v5.1 Nov 2012 '
246 real (kind=kind_phys),
parameter :: gfac=1.0e5/
con_g
257 data ptopc / 1050., 650., 400., 0.0, 1050., 750., 500., 0.0 /
260 real (kind=kind_phys),
parameter :: climit = 0.001, climit2=0.05
261 real (kind=kind_phys),
parameter :: ovcst = 1.0 - 1.0e-8
295 real (kind=kind_phys),
parameter ::
xlim=5.0
297 data xlabdy / 30.0, 0.0, -30.0 /,
xlobdy / 0.0, 180., 360. /
300 real (kind=kind_phys),
parameter ::
vvcld1= 0.0003e0
302 real (kind=kind_phys),
parameter ::
vvcld2=-0.0005e0
373 integer,
intent(in) :: NLAY, me
375 real (kind=kind_phys),
intent(in) :: si(:)
380 integer :: k, kl, ier
389 if (me == 0) print *, vtagcld
392 if (me == 0) print *,
' - Using Diagnostic Cloud Method'
400 99
format(3x,
' *** Error in finding tuned RH table ***'
401 ' STOP at calling subroutine RHTABLE !!'/)
406 print *,
' - Using Prognostic Cloud Method'
408 print *,
' --- Zhao/Carr/Sundqvist microphysics'
410 print *,
' --- Ferrier cloud microphysics'
412 print *,
' --- zhao/carr/sundqvist + pdf cloud'
414 print *,
' !!! ERROR in cloud microphysc specification!!!',
426 lab_do_k0 :
do k = nlay, 2, -1
428 if (si(k) < 0.9e0)
exit lab_do_k0
433 lab_do_k1 :
do k = 2, nlay
435 if (si(k) < 0.9e0)
exit lab_do_k1
483 & xlat,xlon,slmsk, ix, nlay, nlp1,
561 integer,
intent(in) :: IX, NLAY, NLP1
563 logical,
intent(in) :: shoc_cld, lmfshal, lmfdeep2
565 real (kind=kind_phys),
dimension(:,:),
intent(in) :: plvl, plyr, &
566 & tlyr, tvly, qlyr, qstl, rhly, clw, cldcov
568 real (kind=kind_phys),
dimension(:),
intent(in) :: xlat, xlon, &
572 real (kind=kind_phys),
dimension(:,:,:),
intent(out) :: clouds
574 real (kind=kind_phys),
dimension(:,:),
intent(out) :: clds
576 integer,
dimension(:,:),
intent(out) :: mtop,mbot
579 real (kind=kind_phys),
dimension(IX,NLAY) :: cldtot, cldcnv, &
580 & cwp, cip, crp, csp, rew, rei, res, rer, delp, tem2d, clwf
582 real (kind=kind_phys) :: ptop1(ix,
nk_clds+1)
584 real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh,
value, &
587 integer :: i, k, id, nf
591 real (kind=kind_phys),
parameter :: xrc3 = 100.
617 tem2d(i,k) = min( 1.0, max( 0.0, (
con_ttp-tlyr(i,k))*0.05 ) )
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)
629 clwf(i,k) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1)
651 ptop1(i,id) =
ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 )
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)
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)
680 if (nint(slmsk(i)) == 1)
then
682 rew(i,k) = 5.0 + 5.0 * tem2d(i,k)
689 cldtot(i,k) = cldcov(i,k)
700 if (.not. lmfshal)
then
703 clwt = 1.0e-6 * (plyr(i,k)*0.001)
706 if (clwf(i,k) > clwt)
then
708 onemrh= max( 1.e-10, 1.0-rhly(i,k) )
709 clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
711 tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
715 value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 )
716 tem2 = sqrt( sqrt(rhly(i,k)) )
718 cldtot(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
725 clwt = 1.0e-6 * (plyr(i,k)*0.001)
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 )
732 tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0)
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 )
750 if (.not. lmfshal)
then
753 clwt = 1.0e-6 * (plyr(i,k)*0.001)
756 if (clwf(i,k) > clwt)
then
758 onemrh= max( 1.e-10, 1.0-rhly(i,k) )
759 clwm = clwmin / max( 0.01, plyr(i,k)*0.001 )
761 tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
766 value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 )
767 tem2 = sqrt( sqrt(rhly(i,k)) )
769 cldtot(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
776 clwt = 1.0e-6 * (plyr(i,k)*0.001)
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 )
783 tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0)
790 value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 )
791 tem2 = sqrt( sqrt(rhly(i,k)) )
793 cldtot(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
804 if (cldtot(i,k) < climit)
then
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
835 if (cip(i,k) > 0.0)
then
836 tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k))
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
845 rei(i,k) = (1250.0/9.387) * tem3 ** 0.031
849 rei(i,k) = max(10.0, min(rei(i,k), 150.0))
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)
864 clouds(i,k,7) = rer(i,k)
866 clouds(i,k,9) = rei(i,k)
882 & ( plyr, ptop1, cldtot, cldcnv,
935 & xlat,xlon,slmsk, f_ice,f_rain,r_rime,flgmin,
1022 integer,
intent(in) :: IX, NLAY, NLP1
1024 logical,
intent(in) :: lmfshal, lmfdeep2
1026 real (kind=kind_phys),
dimension(:,:),
intent(in) :: plvl, plyr, &
1027 & tlyr, tvly, qlyr, qstl, rhly, clw, f_ice, f_rain, r_rime
1029 real (kind=kind_phys),
dimension(:),
intent(in) :: xlat, xlon, &
1031 real (kind=kind_phys),
dimension(:),
intent(in) :: flgmin
1034 real (kind=kind_phys),
dimension(:,:,:),
intent(out) :: clouds
1036 real (kind=kind_phys),
dimension(:,:),
intent(out) :: clds
1038 integer,
dimension(:,:),
intent(out) :: mtop,mbot
1041 real (kind=kind_phys),
dimension(IX,NLAY) :: cldtot, cldcnv, &
1042 & cwp, cip, crp, csp, rew, rei, res, rer, tem2d, clw2,
1045 real (kind=kind_phys) :: ptop1(ix,
nk_clds+1)
1047 real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh,
value, &
1054 real (kind=kind_phys),
parameter :: xrc3 = 100.
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
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)
1087 clwf(i,k) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1)
1093 clwf(i,k) = clw(i,k)
1109 ptop1(i,id) =
ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 )
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)
1126 qcice(i,k) = clwf(i,k)
1129 clw2(i,k) = clwf(i,k)
1140 & ( plyr, plvl, tlyr, qlyr, qcwat, qcice, qrain, rrime,
1143 & cwp, cip, crp, csp, rew, rer, res, rsden
1150 tem2d(i,k) = (
con_g * plyr(i,k))
1157 tem2d(i,k) = (
con_g * plyr(i,k))
1168 if (.not. lmfshal)
then
1173 clwt = 2.0e-6 * (plyr(i,k)*0.001)
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 )
1184 tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
1185 tem1 = 2000.0 / tem1
1193 value = max( min( tem1*(clw2(i,k)-clwm), 50.0 ), 0.0 )
1194 tem2 = sqrt( sqrt(rhly(i,k)) )
1196 cldtot(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
1204 clwt = 2.0e-6 * (plyr(i,k)*0.001)
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 )
1210 tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0)
1217 value = max( min( tem1*(clw2(i,k)-clwm), 50.0 ), 0.0 )
1218 tem2 = sqrt( sqrt(rhly(i,k)) )
1220 cldtot(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
1229 if (.not. lmfshal)
then
1234 clwt = 2.0e-6 * (plyr(i,k)*0.001)
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 )
1245 tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0)
1246 tem1 = 2000.0 / tem1
1254 value = max( min( tem1*(clw2(i,k)-clwm), 50.0 ), 0.0 )
1255 tem2 = sqrt( sqrt(rhly(i,k)) )
1257 cldtot(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
1265 clwt = 2.0e-6 * (plyr(i,k)*0.001)
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 )
1271 tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0)
1278 value = max( min( tem1*(clw2(i,k)-clwm), 50.0 ), 0.0 )
1279 tem2 = sqrt( sqrt(rhly(i,k)) )
1281 cldtot(i,k) = max( tem2*(1.0-exp(-
value)), 0.0 )
1291 if (cldtot(i,k) < climit)
then
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
1332 if (tem2 > 0.0)
then
1333 tem3 = tem2d(i,k) * tem2 / tvly(i,k)
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
1342 rei(i,k) = (1250.0/9.387) * tem3 ** 0.031
1348 rei(i,k) = max(10.0, min(rei(i,k), 300.0))
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)
1367 clouds(i,k,8) = csp(i,k) * rsden(i,k)
1368 clouds(i,k,9) = rei(i,k)
1382 & ( plyr, ptop1, cldtot, cldcnv,
1516 integer,
intent(in) :: ix, nlay, nlp1,kdt
1518 real (kind=kind_phys),
dimension(:,:),
intent(in) :: plvl, plyr,
1519 & tlyr, tvly, qlyr, qstl, rhly, clw
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
1527 real (kind=kind_phys),
dimension(:),
intent(in) :: xlat, xlon,
1532 real (kind=kind_phys),
dimension(:,:,:),
intent(out) :: clouds
1534 real (kind=kind_phys),
dimension(:,:),
intent(out) :: clds
1536 integer,
dimension(:,:),
intent(out) :: mtop,mbot
1539 real (kind=kind_phys),
dimension(ix,nlay) :: cldtot, cldcnv, &
1540 & cwp, cip, crp, csp, rew, rei, res, rer, delp, tem2d, clwf
1542 real (kind=kind_phys) :: ptop1(ix,
nk_clds+1)
1544 real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh,
value, &
1547 integer :: i, k, id, nf
1555 clouds(i,k,nf) = 0.0
1573 tem2d(i,k) = min( 1.0, max( 0.0, (
con_ttp-tlyr(i,k))*0.05 ) )
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)
1585 clwf(i,k) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1)
1591 clwf(i,k) = clw(i,k)
1599 deltaq(i,k) = (1.-0.95)*qstl(i,k)
1615 ptop1(i,id) =
ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 )
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)
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)
1644 if (nint(slmsk(i)) == 1)
then
1646 rew(i,k) = 5.0 + 5.0 * tem2d(i,k)
1656 tem1 = tlyr(i,k) - 273.16
1658 qsc = sup * qstl(i,k)
1664 if(rhly(i,k) >= rhs)
then
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
1672 elseif(qtmp >= deltaq(i,k))
then
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)
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.)
1695 tem1 = tlyr(i,k) - 273.16
1697 qsc = sup * qstl(i,k)
1703 if(rhly(i,k) >= rhs)
then
1706 qtmp = qlyr(i,k) + clwf(i,k) - qsc
1707 if(deltaq(i,k) > epsq)
then
1709 if(qtmp <= -deltaq(i,k))
then
1711 elseif(qtmp >= deltaq(i,k))
then
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)
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.)
1736 if (cldtot(i,k) < climit)
then
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
1767 if (cip(i,k) > 0.0)
then
1769 tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k))
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
1778 rei(i,k) = (1250.0/9.387) * tem3 ** 0.031
1782 rei(i,k) = max(10.0, min(rei(i,k), 150.0))
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)
1797 clouds(i,k,7) = rer(i,k)
1799 clouds(i,k,9) = rei(i,k)
1814 & ( plyr, ptop1, cldtot, cldcnv,
1921 integer,
intent(in) :: IX, NLAY, NLP1
1923 real (kind=kind_phys),
dimension(:,:),
intent(in) :: plvl, plyr, &
1926 real (kind=kind_phys),
dimension(:),
intent(in) :: xlat, xlon, &
1927 & slmsk, cv, cvt, cvb
1930 real (kind=kind_phys),
dimension(:,:,:),
intent(out) :: clouds
1932 real (kind=kind_phys),
dimension(:,:),
intent(out) :: clds
1934 integer,
dimension(:,:),
intent(out) :: mtop,mbot
1937 real (kind=kind_phys),
dimension(IX,NLAY) :: cldtot, cldcnv, &
1938 & cldtau, taufac, dthdp, tem2d
1940 real (kind=kind_phys) :: ptop1(ix,
nk_clds+1)
1942 real (kind=kind_phys) :: cr1, cr2, crk, pval, cval, omeg,
value, &
1945 integer,
dimension(IX):: idom, kcut
1948 real (kind=kind_phys) :: xlatdg, xlondg, xlnn, xlss, xrgt, xlft, &
1949 & rhcla(NBIN,NLON,MCLD,NSEAL), rhcld(IX,NBIN,MCLD)
1951 integer :: ireg, ib, ic, id, id1, il, is, nhalf
1953 integer :: i, j, k, klowt
1965 lab_do_i_ix :
do i = 1, ix
1967 xlatdg = xlat(i) * tem1
1970 xlondg = xlon(i) * tem1
1971 if (xlondg < 0.0) xlondg = xlondg + 360.0
1977 lab_do_j :
do j = 1, 3
1978 if (xlatdg >
xlabdy(j))
then
1988 rhcla(ib,il,ic,is) =
rhcl(ib,il,ireg,ic,is)
1999 if (xlatdg < xlnn .and. xlatdg > xlss)
then
2004 rhcla(ib,il,ic,is) =
rhcl(ib,il,j+1,ic,is)
2018 if (slmsk(i) < 1.0)
then
2026 if (xlondg > 180.e0)
then
2034 rhcld(i,ib,ic) = rhcla(ib,il,ic,is)
2037 lab_do_k :
do k = 1, 3
2038 tem2 = abs(xlondg -
xlobdy(k))
2040 if (tem2 <
xlim)
then
2043 if (id1 >
nlon) id1 = 1
2049 rhcld(i,ib,ic) = rhcla(ib,id1,ic,is)
2070 ptop1(i,j) =
ptopc(j,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 )
2079 if (tem1 <= -10.0)
then
2080 cldtau(i,k) = max( 0.1e-3, 2.0e-6*(tem1+82.5)**2 )
2082 cldtau(i,k) = min( 0.08, 6.949e-3*tem1+0.08 )
2093 tem1 = (plyr(i,k)*0.001) ** (-con_rocp)
2094 tem2d(i,k) = tem1 * tlyr(i,k)
2100 dthdp(i,k) = (tem2d(i,k+1)-tem2d(i,k))/(plyr(i,k+1)-plyr(i,k))
2117 if (plvl(i,k) < ptop1(i,2)) klowt = max(klowt,k)
2118 taufac(i,k) = plvl(i,k+1) - plvl(i,k)
2129 lab_do_kcut0 :
do k = klowt-1, 2, -1
2130 if (plyr(i,k) <= ptop1(i,3) .and.
2139 if (cv(i) >= climit .and. cvt(i) < cvb(i))
then
2143 lab_do_k_cvt0 :
do k = 2, nlay
2144 if (cvt(i) <= plyr(i,k))
then
2150 lab_do_k_cvb0 :
do k = nlay-1, 1, -1
2151 if (cvb(i) >= plyr(i,k))
then
2157 tem1 = plyr(i,id1) - plyr(i,id)
2160 taufac(i,k) = taufac(i,k) * max( 0.25, 1.0-0.125*tem1 )
2176 do k = nlay-1, 2, -1
2184 lab_do_ic0 :
do ic = 2, 4
2185 if(plyr(i,k) >= ptop1(i,ic))
then
2199 nhalf = (
nbin + 1) / 2
2201 if (id <= 0 .or. k < kcut(i))
then
2203 elseif (rhly(i,k) <= rhcld(i,1,id))
then
2205 elseif (rhly(i,k) >= rhcld(i,
nbin,id))
then
2212 do while ( notstop )
2213 nhalf = (nhalf + 1) / 2
2214 cr1 = rhcld(i,ib, id)
2215 cr2 = rhcld(i,ib+1,id)
2217 if (crk <= cr1)
then
2218 ib = max( ib-nhalf, 1 )
2219 elseif (crk > cr2)
then
2220 ib = min( ib+nhalf,
nbin-1 )
2222 cldtot(i,k) = 0.01 * (ib + (crk - cr1)/(cr2 - cr1))
2234 do k = klowt,
llyr+1
2243 if (cval >= climit .and. pval >= ptop1(i,2))
then
2246 elseif (omeg >
vvcld2)
then
2247 tem1 = (
vvcld1 - omeg) /
value
2248 cldtot(i,k) = cldtot(i,k) * sqrt(tem1)
2265 if (plvl(i,k) < ptop1(i,2)) klowt = min(klowt,k)
2266 taufac(i,k) = plvl(i,k) - plvl(i,k+1)
2277 lab_do_kcut1 :
do k = klowt+1, nlay-1
2278 if (plyr(i,k) <= ptop1(i,3) .and.
2287 if (cv(i) >= climit .and. cvt(i) < cvb(i))
then
2291 lab_do_k_cvt :
do k = nlay-1, 1, -1
2292 if (cvt(i) <= plyr(i,k))
then
2298 lab_do_k_cvb :
do k = 2, nlay
2299 if (cvb(i) >= plyr(i,k))
then
2305 tem1 = plyr(i,id1) - plyr(i,id)
2308 taufac(i,k) = taufac(i,k) * max( 0.25, 1.0-0.125*tem1 )
2332 lab_do_ic1 :
do ic = 2, 4
2333 if(plyr(i,k) >= ptop1(i,ic))
then
2347 nhalf = (
nbin + 1) / 2
2349 if (id <= 0 .or. k > kcut(i))
then
2351 elseif (rhly(i,k) <= rhcld(i,1,id))
then
2353 elseif (rhly(i,k) >= rhcld(i,
nbin,id))
then
2360 do while ( notstop )
2361 nhalf = (nhalf + 1) / 2
2362 cr1 = rhcld(i,ib, id)
2363 cr2 = rhcld(i,ib+1,id)
2365 if (crk <= cr1)
then
2366 ib = max( ib-nhalf, 1 )
2367 elseif (crk > cr2)
then
2368 ib = min( ib+nhalf,
nbin-1 )
2370 cldtot(i,k) = 0.01 * (ib + (crk - cr1)/(cr2 - cr1))
2382 do k =
llyr-1, klowt
2391 if (cval >= climit .and. pval >= ptop1(i,2))
then
2394 elseif (omeg >
vvcld2)
then
2395 tem1 = (
vvcld1 - omeg) /
value
2396 cldtot(i,k) = cldtot(i,k) * sqrt(tem1)
2409 where (cldtot < climit)
2412 where (cldcnv < climit)
2416 where (cldtot < climit .and. cldcnv < climit)
2422 clouds(i,k,1) = max(cldtot(i,k), cldcnv(i,k))
2423 clouds(i,k,2) = cldtau(i,k) * taufac(i,k)
2438 & ( plyr, ptop1, cldtot, cldcnv,
2526 integer,
intent(in) :: IX, NLAY
2528 real (kind=kind_phys),
dimension(:,:),
intent(in) :: plyr, ptop1, &
2532 real (kind=kind_phys),
dimension(:,:),
intent(out) :: clds
2534 integer,
dimension(:,:),
intent(out) :: mtop, mbot
2537 real (kind=kind_phys) :: cl1(ix), cl2(ix)
2538 real (kind=kind_phys) :: pcur, pnxt, ccur, cnxt
2540 integer,
dimension(IX):: idom, kbt1, kth1, kbt2, kth2
2541 integer :: i, k, id, id1, kstr, kend, kinc
2569 if (
iovr == 0 )
then
2571 do k = kstr, kend, kinc
2573 ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
2574 if (ccur >= climit) cl1(i) = cl1(i) * (1.0 - ccur)
2579 clds(i,5) = 1.0 - cl1(i)
2585 clds(i,4) = 1.0 - cl1(i)
2590 do k = kstr, kend, kinc
2592 ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
2593 if (ccur >= climit)
then
2594 cl2(i) = min( cl2(i), (1.0 - ccur) )
2596 cl1(i) = cl1(i) * cl2(i)
2603 clds(i,5) = 1.0 - cl1(i) * cl2(i)
2609 clds(i,4) = 1.0 - cl1(i) * cl2(i)
2633 mbot(i,2) = nlay - 1
2634 mtop(i,2) = nlay - 1
2635 mbot(i,3) = nlay - 1
2636 mtop(i,3) = nlay - 1
2646 ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
2650 cnxt = min( ovcst, max( cldtot(i,k-1), cldcnv(i,k-1) ))
2656 if (pcur < ptop1(i,id1))
then
2662 if (ccur >= climit)
then
2663 if (kth2(i) == 0) kbt2(i) = k
2664 kth2(i) = kth2(i) + 1
2666 if (
iovr == 0 )
then
2667 cl2(i) = cl2(i) + ccur - cl2(i)*ccur
2669 cl2(i) = max( cl2(i), ccur )
2672 if (cnxt < climit .or. pnxt < ptop1(i,id1))
then
2673 kbt1(i) = nint( (cl1(i)*kbt1(i) + cl2(i)*kbt2(i) )
2685 if (pnxt < ptop1(i,id1))
then
2687 mtop(i,id) = min( kbt1(i), kbt1(i)-kth1(i)+1 )
2688 mbot(i,id) = kbt1(i)
2695 mbot(i,id1) = kbt1(i)
2696 mtop(i,id1) = kbt1(i)
2728 ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
2732 cnxt = min( ovcst, max( cldtot(i,k+1), cldcnv(i,k+1) ))
2738 if (pcur < ptop1(i,id1))
then
2744 if (ccur >= climit)
then
2745 if (kth2(i) == 0) kbt2(i) = k
2746 kth2(i) = kth2(i) + 1
2748 if (
iovr == 0 )
then
2749 cl2(i) = cl2(i) + ccur - cl2(i)*ccur
2751 cl2(i) = max( cl2(i), ccur )
2754 if (cnxt < climit .or. pnxt < ptop1(i,id1))
then
2755 kbt1(i) = nint( (cl1(i)*kbt1(i) + cl2(i)*kbt2(i))
2767 if (pnxt < ptop1(i,id1))
then
2769 mtop(i,id) = max( kbt1(i), kbt1(i)+kth1(i)-1 )
2770 mbot(i,id) = kbt1(i)
2773 kbt1(i) = min(k+1, nlay)
2777 mbot(i,id1) = kbt1(i)
2778 mtop(i,id1) = kbt1(i)
2817 integer,
intent(in) :: me
2820 integer,
intent(out) :: ier
2823 real (kind=kind_phys),
dimension(NBIN,NLON,NLAT,MCLD,NSEAL) :: &
2824 & rhfd, rtnfd, rhcf, rtncf, rhcla
2826 real (kind=kind_io4),
dimension(NBIN,NLON,NLAT,MCLD,NSEAL) :: &
2829 real(kind=kind_io4) :: fhour
2831 real(kind=kind_phys) :: binscl, cfrac, clsat, rhsat, cstem
2833 integer,
dimension(NLON,NLAT,MCLD,NSEAL) :: kpts, kkpts
2835 integer :: icdays(15), idate(4), nbdayi, isat
2837 integer :: i, i1, j, k, l, m, id, im, iy
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
2869 read (nicltun,err=998,end=999) nbdayi, icdays
2871 if (me == 0) print 11, nbdayi
2872 11
format(
' from rhtable DAYS ON FILE =',i5)
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)
2882 read (nicltun,err=998,end=999) fhour,idate
2885 read (nicltun) rhfd4
2888 read (nicltun) rtnfd4
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)
2906 kkpts = kkpts + kpts
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)
2922 if (kkpts(j,k,l,m) > 0)
then
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))
2932 rhcf(
nbin,j,k,l,m) = 1.0
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
2947 rhcf(i,j,k,l,m) = -0.1
2948 rtncf(i,j,k,l,m) = -0.1
2953 210
format(
' NO CRIT RH FOR LAT=',i3,
' AND LON BAND=',i3,
2954 ' LAND(=1) SEA=',i3//
' MODEL RH ',
' OBS RTCLD')
2956 print 203, rhcf(i,j,k,l,m), rtncf(i,j,k,l,m)
2975 cfrac = binscl * (i - 1)
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
2984 if (rtncf(i,j,k,l,m) >= 1.0)
then
2987 rhsat = rhcla(i,j,k,l,m)
2991 rhcla(i,j,k,l,m) = rhsat + (1.0 - rhsat)
2996 rhcla(
nbin,j,k,l,m) = 1.0
3015 lab_do_i :
do i = 1,
nbin - 2
3018 if (rhcla(i,j,k,l,m) >= 0.98)
then
3022 rhcl(i1,j,k,l,m) = rhcla(i,j,k,l,m)
3036 print *,
'completed rhtable for cloud tuninig tables'
3041 988
format(
' from rhtable ERROR READING TABLES')
3046 989
format(
' from rhtable E.O.F READING TABLES')
3057 end module module_radiation_clouds
real(kind=kind_phys), dimension(3) xlabdy
lat bndry between tuning regions
integer, save iovrsw
cloud overlapping control flag for SW
integer, parameter nlon
=1,2 for eastern and western hemispheres
real(kind=kind_phys), parameter con_pi
pi
integer, parameter, public nk_clds
number of cloud vertical domains
real(kind=kind_phys), dimension(3) xlobdy
lon bndry between tuning regions
real(kind=kind_phys), parameter vvcld2
low cloud vertical velocity adjustment boundaries in mb/sec
real(kind=kind_phys), parameter con_g
gravity ( )
subroutine, public diagcld1
This subroutine computes cloud fractions for radiation calculations.
real(kind=kind_phys), parameter xlim
+/- xlim for transition
logical, save lcnorm
in-cld condensate control flag
logical, save lcrick
eliminating CRICK control flag
subroutine, public progcld2
This subroutine computes cloud related quantities using ferrier's prognostic cloud microphysics schem...
real(kind=kind_phys), dimension(nk_clds+1, 2), save ptopc
pressure limits of cloud domain interfaces (low,mid,high) in mb (0.1kPa)
real(kind=kind_phys), parameter con_thgni
temperature the H.G.Nuc. ice starts
subroutine gethml
This subroutine computes high, mid, low, total, and boundary cloud fractions and cloud top/bottom lay...
integer iovr
maximum-random cloud overlapping method
integer, parameter nlat
=1,4 for 60n-30n,30n-equ,equ-30s,30s-60s
real(kind=kind_phys), parameter cldssa_def
default cld single scat albedo
integer, save icmphys
cloud micorphysics scheme control flag
real(kind=kind_phys), parameter con_t0c
temp at 0C (K)
integer, save icldflg
cloud optical property scheme control flag
real(kind=kind_phys), parameter vvcld1
low cloud vertical velocity adjustment boundaries in mb/sec
integer, parameter nbin
rh in one percent interval
real(kind=kind_phys), parameter cldasy_def
default cld asymmetry factor
logical, save lnoprec
precip effect on radiation flag (Ferrier microphysics)
subroutine rhtable
cld-rh relations obtained from mitchell-hahn procedure.
subroutine, public progcld1
This subroutine computes cloud related quantities using zhao/moorthi's prognostic cloud microphysics ...
integer, parameter nseal
=1,2 for land,sea
real(kind=kind_phys), parameter con_ttp
temp at H2O 3pt (K)
real(kind=kind_phys), parameter rsnow_def
default snow radius to 250 micron
real(kind=kind_phys), parameter con_rd
gas constant air ( )
subroutine, public progcld3
This subroutine computes cloud related quantities using zhao/moorthi's prognostic cloud microphysics ...
subroutine, public cld_init
This subroutine is an initialization program for cloud-radiation calculations and sets up boundary la...
real(kind=kind_phys), parameter reice_def
default ice radius to 50 micron
real(kind=kind_phys), parameter reliq_def
default liq radius to 10 micron
integer, parameter, public nf_clds
number of fields in cloud array
integer, parameter mcld
=1,4 for bl,low,mid,hi cld type
integer, save ivflip
vertical profile indexing flag
integer, save iovrlw
cloud overlapping control flag for LW
real(kind=kind_phys), parameter rrain_def
default rain radius to 1000 micron
integer llyr
upper limit of boundary layer clouds
real(kind=kind_phys), dimension(nbin, nlon, nlat, mcld, nseal) rhcl
tuned relative humidity relation table for diagnostic cloud scheme