diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/gfdl_cloud_microphys.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/gfdl_cloud_microphys.F90 index bcb4df19d..9e82c4564 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/gfdl_cloud_microphys.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/gfdl_cloud_microphys.F90 @@ -52,7 +52,7 @@ module gfdl2_cloud_microphys_mod public cloud_diagnosis public ICE_LSC_VFALL_PARAM, ICE_CNV_VFALL_PARAM - integer :: ICE_LSC_VFALL_PARAM = 1 + integer :: ICE_LSC_VFALL_PARAM = 1 integer :: ICE_CNV_VFALL_PARAM = 2 real :: missing_value = - 1.e10 @@ -118,7 +118,7 @@ module gfdl2_cloud_microphys_mod real, parameter :: sfcrho = 1.2 !< surface air density real, parameter :: rhor = 1.e3 !< density of rain water, lin83 - + real, parameter :: rc = (4. / 3.) * pi * rhor real :: cracs, csacr, cgacr, cgacs, csacw, craci, csaci, cgacw, cgaci, cracw !< constants for accretions @@ -206,7 +206,7 @@ module gfdl2_cloud_microphys_mod real :: rthreshu = 1.0e-6 !< critical cloud drop radius (micro m) real :: rthreshs = 10.0e-6 !< critical cloud drop radius (micro m) - + real :: sat_adj0 = 0.90 !< adjustment factor (0: no, 1: full) during fast_sat_adj real :: qc_crt = 5.0e-8 !< mini condensate mixing ratio to allow partial cloudiness @@ -215,7 +215,7 @@ module gfdl2_cloud_microphys_mod real :: ql_mlt = 2.0e-3 !< max value of cloud water allowed from melted cloud ice real :: qs_mlt = 1.0e-6 !< max cloud water due to snow melt - + real :: qi_gen = 9.82679e-5 !< max cloud ice generation at -40 C ! cloud condensate upper bounds: "safety valves" for ql & qi @@ -261,9 +261,9 @@ module gfdl2_cloud_microphys_mod real :: vi_max = 1.0 !< max fall speed for ice real :: vs_max = 3.0 !< max fall speed for snow - real :: vg_max = 6.0 !< max fall speed for graupel + real :: vg_max = 6.0 !< max fall speed for graupel real :: vr_max = 9.0 !< max fall speed for rain - real :: vh_max = 19.0 !< max fall speed for hail + real :: vh_max = 19.0 !< max fall speed for hail ! cloud microphysics switchers @@ -524,12 +524,12 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & real, intent (in), dimension (is:) :: eis real, intent (in), dimension (is:, js:, ks:) :: rhcrit - + real, intent (in) :: anv_icefall, lsc_icefall real, intent (in), dimension (is:, js:, ks:) :: uin, vin, delp, pt, dz real, intent (in), dimension (is:, js:, ks:) :: qv, qi, ql, qr, qs, qg, qa, qn - + real, intent (inout), dimension (is:, js:, ks:) :: u_dt, v_dt, w, pt_dt, qa_dt real, intent (inout), dimension (is:, js:, ks:) :: qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt real, intent ( out), dimension (is:, js:, ks:) :: revap, isubl @@ -541,8 +541,8 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & real, intent (out), dimension (is:, js:, ks:) :: vt_r, vt_s, vt_g, vt_i, qn2 real, intent (out), dimension (is:, js:, ks:) :: m2_rain, m2_sol - - real, dimension (ktop:kbot) :: h_var1d + + real, dimension (ktop:kbot) :: h_var1d real, dimension (ktop:kbot) :: qvz, qlz, qrz, qiz, qsz, qgz, qaz real, dimension (ktop:kbot) :: vtiz, vtsz, vtgz, vtrz real, dimension (ktop:kbot) :: dp0, dp1, dz0, dz1 @@ -551,7 +551,7 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & real, dimension (ktop:kbot) :: ccn, c_praut, m1_rain, m1_sol, m1, evap1, subl1 real, dimension (ktop:kbot) :: u0, v0, u1, v1, w1 - real :: onemsig, fac_eis + real :: onemsig, fac_eis real :: cpaut, rh_adj, rh_rain real :: r1, s1, i1, g1, rdt, ccn0 real :: dts @@ -576,7 +576,7 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & tz (k) = t0 (k) dp1 (k) = delp (i, j, k) dp0 (k) = dp1 (k) ! moist air mass * grav - + ! ----------------------------------------------------------------------- ! import horizontal subgrid variability with pressure dependence ! total water subgrid deviation in horizontal direction @@ -587,7 +587,7 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & ! ----------------------------------------------------------------------- ! convert moist mixing ratios to dry mixing ratios ! ----------------------------------------------------------------------- - + qvz (k) = qv (i, j, k) qlz (k) = ql (i, j, k) qiz (k) = qi (i, j, k) @@ -606,7 +606,7 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & qiz (k) = qiz (k) * omq qsz (k) = qsz (k) * omq qgz (k) = qgz (k) * omq - + qa0 (k) = qa (i, j, k) qaz (k) = qa (i, j, k) dz0 (k) = dz (i, j, k) @@ -650,14 +650,14 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & ! calculate cloud condensation nuclei (ccn) ! the following is based on klein eq. 15 ! ----------------------------------------------------------------------- - if (srf_type(i) < 2.0) then ! exclude snow/ice covered regions + if (srf_type(i) < 2.0) then ! exclude snow/ice covered regions fac_eis = min(1.0,eis(i)/10.0)**2 ! Estimated inversion strength determine stable regime cpaut = c_paut * (0.75*fac_eis + (1.0-fac_eis)) ! scaling autoconversion for stable->unstable else fac_eis = 0.0 cpaut = c_paut endif - ! ccn needs units #/m^3 + ! ccn needs units #/m^3 do k = ktop, kbot ! qn has units # / m^3 ccn (k) = qn (i, j, k) @@ -722,7 +722,7 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & ! ----------------------------------------------------------------------- ! warm rain processes ! ----------------------------------------------------------------------- - + call warm_rain (dts, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, & qgz, qaz, fac_eis, onemsig, den, denfac, ccn, c_praut, vtrz, & r1, evap1, m1_rain, w1, h_var1d) @@ -774,9 +774,9 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & ! ----------------------------------------------------------------------- ! fix all negative water species ! ----------------------------------------------------------------------- - + if (fix_negative) & - call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz) + call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz) ! ----------------------------------------------------------------------- ! update moist air mass (actually hydrostatic pressure) @@ -926,7 +926,7 @@ subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, qa, & real, intent (in) :: onemsig real, intent (in) :: fac_eis !< estimated inversion strength - + real, intent (inout), dimension (ktop:kbot) :: tz, vtr real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qi, qs, qg, qa real, intent (inout), dimension (ktop:kbot) :: evap1, m1_rain, w1 @@ -963,7 +963,7 @@ subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, qa, & m1_rain (:) = 0. call check_column (ktop, kbot, qr, no_fall) - + ! ----------------------------------------------------------------------- ! auto - conversion ! assuming linear subgrid vertical distribution of cloud water @@ -986,7 +986,7 @@ subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, qa, & ! ----------------------------------------------------------------------- ! no subgrid varaibility ! ----------------------------------------------------------------------- - + do k = ktop, kbot if (tz (k) > t_wfr) then qc = fac_rc * ccn (k) / den (k) @@ -1002,7 +1002,7 @@ subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, qa, & endif endif enddo - + else ! ----------------------------------------------------------------------- @@ -1042,11 +1042,11 @@ subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, qa, & ! Revert In-Cloud condensate ql = ql*qadum qi = qi*qadum - + ! ----------------------------------------------------------------------- ! fall speed of rain ! ----------------------------------------------------------------------- - + if (no_fall) then vtr (:) = vf_min elseif (const_vr) then @@ -1085,7 +1085,7 @@ subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, qa, & ! ----------------------------------------------------------------------- ! mass flux induced by falling rain ! ----------------------------------------------------------------------- - + if (no_fall) then r1 = 0.0 elseif (use_ppm) then @@ -1142,7 +1142,7 @@ subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, qa, revap, de integer, intent (in) :: ktop, kbot real, intent (in) :: dt ! time step (s) - + real, intent (in), dimension (ktop:kbot) :: h_var real, intent (in), dimension (ktop:kbot) :: den, denfac @@ -1154,17 +1154,20 @@ subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, qa, revap, de real :: dqv, qsat, dqsdt, evap, t2, qden, q_plus, q_minus, sink real :: qpz, dq, dqh, tin - real :: fac_revp + real :: fac_revp real :: TOT_PREC_LS, AREA_LS_PRC, AREA_LS_PRC_K integer :: k revap(:) = 0. + TOT_PREC_LS = 0.0 + AREA_LS_PRC = 0.0 + do k = ktop, kbot TOT_PREC_LS = TOT_PREC_LS + ( ( qr (k) + qs (k) + qg (k) ) * den (k) ) AREA_LS_PRC = AREA_LS_PRC + ( qa (k) * ( qr (k) + qs (k) + qg (k) ) * den (k) ) - + if (tz (k) > t_wfr .and. qr (k) > qpmin) then ! area and timescale efficiency on revap @@ -1226,7 +1229,7 @@ subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, qa, revap, de ! ----------------------------------------------------------------------- ! accretion: pracc ! ----------------------------------------------------------------------- - + if (qr (k) > qpmin .and. ql (k) > qcmin .and. qsat < q_minus) then sink = dt * denfac (k) * cracw * exp (0.95 * log (qr (k) * den (k))) sink = sink / (1. + sink) * ql (k) @@ -1314,13 +1317,13 @@ end subroutine linear_prof subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, & den, denfac, vts, vtg, vtr, qak, dts, subl1, h_var, ccn, cnv_fraction, srf_type, onemsig) - + implicit none integer, intent (in) :: ktop, kbot real, intent (in), dimension (ktop:kbot) :: p1, dp1, den, denfac, vts, vtg, vtr - + real, intent (inout), dimension (ktop:kbot) :: tzk, qvk, qlk, qrk, qik, qsk, qgk, qak real, intent (out), dimension (ktop:kbot) :: subl1 @@ -1328,7 +1331,7 @@ subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, & real, intent (in) :: dts, cnv_fraction, srf_type, onemsig real, intent (in), dimension (ktop:kbot) :: h_var, ccn - + real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, di, lhl, lhi real, dimension (ktop:kbot) :: cvm, q_liq, q_sol @@ -1343,9 +1346,9 @@ subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, & real :: qadum integer :: k, it - + rdts = 1. / dts - + ! ----------------------------------------------------------------------- ! define conversion scalar / factor ! ----------------------------------------------------------------------- @@ -1353,7 +1356,7 @@ subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, & fac_i2s = 1. - exp (- dts / tau_i2s) fac_imlt = 1. - exp (- dts / tau_imlt) fac_frz = 1. - exp (- dts / tau_frz) - + ! ----------------------------------------------------------------------- ! define heat capacity and latend heat coefficient ! ----------------------------------------------------------------------- @@ -1378,7 +1381,7 @@ subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, & ! Use In-Cloud condensates if (in_cloud) then qadum = max(qak (k),max(qcmin,onemsig)) - else + else qadum = 1.0 endif @@ -1388,7 +1391,7 @@ subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, & newice = max(0.0,qi + new_ice_condensate(tzk (k), ql, qi, cnv_fraction, srf_type)) newliq = max(0.0,ql + qi - newice) - melt = max(0.0,newliq - ql) + melt = max(0.0,newliq - ql) frez = max(0.0,newice - qi) if (melt > 0.0 .and. tzk (k) > tice .and. qi > qcmin) then @@ -1824,7 +1827,7 @@ end subroutine icloud subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & ql, qr, qi, qs, qg, qa, subl1, h_var, ccn, cnv_fraction, srf_type, onemsig) - + implicit none integer, intent (in) :: ktop, kbot @@ -1858,13 +1861,13 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & real :: pssub, pgsub, tsq, qden, fac_g2v, fac_v2g, fac_s2v, fac_v2s real :: ifrac, newqi, fac_frz real :: rh_adj, rh_rain - + integer :: k - + ! ----------------------------------------------------------------------- ! define conversion scalar / factor ! ----------------------------------------------------------------------- - + fac_l2v = 1. - exp (- dts / tau_l2v) fac_i2v = 1. - exp (- dts / tau_i2v) fac_s2v = 1. - exp (- dts / tau_s2v) @@ -1890,7 +1893,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & enddo do k = ktop, kbot - + rh_adj = 1. - h_var(k) - rh_inc rh_rain = max (0.35, 1. - h_var(k) - rh_inr) @@ -1923,7 +1926,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & icpk (k) = lhi (k) / cvm (k) tcpk (k) = lcpk (k) + icpk (k) tcp3 (k) = lcpk (k) + icpk (k) * min (1., dim (tice, tz (k)) / (tice - t_wfr)) - + ! ----------------------------------------------------------------------- ! cloud water < -- > vapor adjustment: LS evaporation ! ----------------------------------------------------------------------- @@ -1949,7 +1952,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & endif endif evap = evap*onemsig ! resolution dependent evap 0:1 coarse:fine - ! new total condensate / old condensate + ! new total condensate / old condensate qa(k) = max(0.0,min(1.,qa(k) * max(qi(k)+ql(k)-evap,0.0 ) / & max(qi(k)+ql(k) ,qcmin) ) ) qv (k) = qv (k) + evap @@ -2002,7 +2005,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice tz (k) = tz (k) + sink * lhi (k) / cvm (k) endif ! significant ql existed - + ! ----------------------------------------------------------------------- ! update capacity heat and latend heat coefficient ! ----------------------------------------------------------------------- @@ -2046,7 +2049,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & endif sink = sink*onemsig ! resolution dependent subl 0:1 coarse:fine endif - ! new total condensate / old condensate + ! new total condensate / old condensate qa(k) = max(0.0,min(1.,qa(k) * max(qi(k)+ql(k)+sink,0.0 ) / & max(qi(k)+ql(k) ,qcmin) ) ) qv (k) = qv (k) - sink @@ -2070,7 +2073,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & ! sublimation / deposition of snow ! this process happens for all temp rage ! ----------------------------------------------------------------------- - + if (qs (k) > qpmin) then qsi = iqs2 (tz (k), den (k), dqsdt) qden = qs (k) * den (k) @@ -2180,7 +2183,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & q_cond (k) = q_liq (k) + q_sol (k) qpz = qv (k) + q_cond (k) ! qpz is conserved - + ! ----------------------------------------------------------------------- ! use the "liquid - frozen water temperature" (tin) to compute saturated specific humidity ! ----------------------------------------------------------------------- @@ -2188,7 +2191,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & tin = tz (k) - (lcpk (k) * q_cond (k) + icpk (k) * q_sol (k)) ! minimum temperature ! tin = tz (k) - ((lv00 + d0_vap * tz (k)) * q_cond (k) + & ! (li00 + dc_ice * tz (k)) * q_sol (k)) / (c_air + qpz * c_vap) - + ! ----------------------------------------------------------------------- ! determine saturated specific humidity ! ----------------------------------------------------------------------- @@ -2206,7 +2209,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & if (q_cond (k) > 3.e-6) then rqi = q_sol (k) / q_cond (k) else - ! WMP impose CALIPSO ice polynomial from 0 C to -40 C + ! WMP impose CALIPSO ice polynomial from 0 C to -40 C rqi = ice_fraction(tin,cnv_fraction,srf_type) endif qstar = rqi * qsi + (1. - rqi) * qsw @@ -2224,7 +2227,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & if (icloud_f == 3) then ! triangular if(q_plus.le.qstar) then - ! little/no cloud cover + ! little/no cloud cover elseif ( (qpz.le.qstar).and.(qstar.lt.q_plus) ) then ! partial cloud cover qa (k) = max(qcmin, min(1., qa (k) + (q_plus-qstar)*(q_plus-qstar) / ( (q_plus-q_minus)*(q_plus-qpz) ))) elseif ( (q_minus.le.qstar).and.(qstar.lt.qpz) ) then ! partial cloud cover @@ -2235,7 +2238,7 @@ subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, tz, qv, & else ! top-hat if(q_plus.le.qstar) then - ! little/no cloud cover + ! little/no cloud cover elseif (qstar < q_plus .and. q_cond (k) > qc_crt) then qa (k) = max(qcmin, min(1., qa (k) + (q_plus - qstar) / (dq + dq) )) ! partial cloud cover elseif (qstar .le. q_minus) then @@ -2269,7 +2272,7 @@ subroutine terminal_fall (dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, & real, intent (out) :: r1, g1, s1, i1 real, dimension (ktop:kbot + 1) :: ze, zt - + real :: qsat, dqsdt, evap, dtime real :: factor, frac real :: tmp, precip, tc, sink @@ -2283,9 +2286,9 @@ subroutine terminal_fall (dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, & integer :: k, k0, m logical :: no_fall - + fac_imlt = 1. - exp (- dtm / tau_imlt) - + ! ----------------------------------------------------------------------- ! define heat capacity and latend heat coefficient ! ----------------------------------------------------------------------- @@ -3047,7 +3050,7 @@ subroutine fall_speed (ktop, kbot, pl, cnv_fraction, anv_icefall, lsc_icefall, & ! https://doi.org/10.1029/2008GL035054 ! ----------------------------------------------------------------------- viLSC = lsc_icefall*10.0**(log10(IWC) * (tc * (aaL * tc + bbL) + ccL) + ddL * tc + eeL) - else + else ! ----------------------------------------------------------------------- ! use Mishra et al (2014, JGR) 'Parameterization of ice fall speeds in ! ice clouds: Results from SPartICus' @@ -3178,7 +3181,7 @@ subroutine setupm pisq = pie * pie scm3 = (visk / vdifu) ** (1. / 3.) - c_paut = c_paut * 0.104 * grav / 1.717e-5 + c_paut = c_paut * 0.104 * grav / 1.717e-5 cracs = pisq * rnzr * rnzs * rhos csacr = pisq * rnzr * rnzs * rhor @@ -3216,7 +3219,7 @@ subroutine setupm ! decreasing gcon will reduce accretion of graupel from cloud ice/water cgacw = pie * rnzg * gcon * gam350 / (4. * act (6) ** 0.875) - cgaci = c_pgaci * cgacw + cgaci = c_pgaci * cgacw ! subl and revp: five constants for three separate processes