diff --git a/schemes/cloud_fraction/compute_cloud_fraction.F90 b/schemes/cloud_fraction/compute_cloud_fraction.F90 new file mode 100644 index 00000000..29873730 --- /dev/null +++ b/schemes/cloud_fraction/compute_cloud_fraction.F90 @@ -0,0 +1,561 @@ +! Compute cloud fractions using relative humidity (RH) threshold and other methods. +! CCPP-ized: Haipeng Lin, February 2025 +module compute_cloud_fraction + use ccpp_kinds, only: kind_phys + private + save + + public :: compute_cloud_fraction_init + public :: compute_cloud_fraction_timestep_init + public :: compute_cloud_fraction_run + + ! Namelist variables + logical :: cldfrc_freeze_dry ! flag for Vavrus correction + logical :: cldfrc_ice ! flag to compute ice cloud fraction + integer :: iceopt ! option for ice cloud closure + ! 1=wang & sassen 2=schiller (iciwc) + ! 3=wood & field, 4=Wilson (based on smith) + real(kind_phys) :: rhminl ! minimum rh for low stable clouds + real(kind_phys) :: rhminl_adj_land ! rhminl adjustment for snowfree land + real(kind_phys) :: rhminh ! minimum rh for high stable clouds + real(kind_phys) :: premit ! top pressure bound for mid level cloud + real(kind_phys) :: premib ! bottom pressure bound for mid level cloud + real(kind_phys) :: icecrit ! critical RH for ice clouds in Wilson and Ballard closure (smaller = more ice clouds) + + ! flags + logical :: inversion_cld_off ! turn off stratification-based cloud fraction (inversion_cld) + + integer :: k700 ! model level nearest to 700 mb [index] + + ! tuning constants + real(kind_phys), parameter :: pretop = 1.0e2_kind_phys ! pressure bounding high cloud [Pa] + +contains + ! Initialize cloud_fraction from namelist parameters +!> \section arg_table_compute_cloud_fraction_init Argument Table +!! \htmlinclude compute_cloud_fraction_init.html + subroutine compute_cloud_fraction_init( & + amIRoot, iulog, & + pver, & + pref_mid, & + inversion_cld_off_in, & + cldfrc_freeze_dry_in, cldfrc_ice_in, iceopt_in, & + rhminl_in, rhminl_adj_land_in, & + rhminh_in, & + premit_in, premib_in, & + icecrit_in, & + errmsg, errflg) + + ! Input arguments + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog + integer, intent(in) :: pver + real(kind_phys), intent(in) :: pref_mid(:) ! reference_pressure [Pa] + logical, intent(in) :: inversion_cld_off_in ! flag to turn off inversion_cld? + logical, intent(in) :: cldfrc_freeze_dry_in ! flag for Vavrus correction + logical, intent(in) :: cldfrc_ice_in ! flag to compute ice cloud fraction + integer, intent(in) :: iceopt_in ! option for ice cloud closure + real(kind_phys), intent(in) :: rhminl_in ! minimum rh for low stable clouds + real(kind_phys), intent(in) :: rhminl_adj_land_in ! rhminl adjustment for snowfree land + real(kind_phys), intent(in) :: rhminh_in ! minimum rh for high stable clouds + real(kind_phys), intent(in) :: premit_in ! top pressure bound for mid level cloud + real(kind_phys), intent(in) :: premib_in ! bottom pressure bound for mid level cloud + real(kind_phys), intent(in) :: icecrit_in ! critical RH for ice clouds in Wilson and Ballard closure + + ! Output arguments + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + errflg = 0 + errmsg = '' + + ! Set module variables from input arguments + cldfrc_freeze_dry = cldfrc_freeze_dry_in + cldfrc_ice = cldfrc_ice_in + iceopt = iceopt_in + rhminl = rhminl_in + rhminl_adj_land = rhminl_adj_land_in + rhminh = rhminh_in + premit = premit_in + premib = premib_in + icecrit = icecrit_in + inversion_cld_off = inversion_cld_off_in + + if(amIRoot) then + write(iulog,*) 'tuning parameters compute_cloud_fraction_init: inversion_cld_off ', inversion_cld_off + write(iulog,*) 'tuning parameters compute_cloud_fraction_init: rhminl ',rhminl,' rhminl_adj_land ',rhminl_adj_land, & + ' rhminh ',rhminh,' premit ',premit,' premib ',premib + write(iulog,*) 'tuning parameters compute_cloud_fraction_init: iceopt ',iceopt,' icecrit ',icecrit + endif + + ! Find vertical level nearest 700 mb. + k700 = minloc(abs(pref_mid(:) - 7.e4_kind_phys), 1) + + if(amIRoot) then + write(iulog,*) 'compute_cloud_fraction: model level nearest 700 mb is ',k700,' which is ',pref_mid(k700),' pascals ' + endif + + end subroutine compute_cloud_fraction_init + + ! Timestep initialization for cloud fraction + ! Specify the perturbation to RH to none by default. + ! If a scheme has to perturb RH this quantity should be set by an interstitial + ! before calling cloud_fraction a second time. +!> \section arg_table_compute_cloud_fraction_timestep_init Argument Table +!! \htmlinclude compute_cloud_fraction_timestep_init.html + subroutine compute_cloud_fraction_timestep_init(rhpert_flag, errmsg, errflg) + logical, intent(out) :: rhpert_flag + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + rhpert_flag = .false. + end subroutine compute_cloud_fraction_timestep_init + + ! Compute cloud fraction. + ! + ! Calculates cloud fraction using a relative humidity threshold + ! The threshold depends upon pressure, and upon the presence or absence + ! of convection as defined by a reasonably large vertical mass flux + ! entering that layer from below. + ! + ! Original authors: Many, including Jim McCaa. +!> \section arg_table_compute_cloud_fraction_run Argument Table +!! \htmlinclude compute_cloud_fraction_run.html + subroutine compute_cloud_fraction_run( & + ncol, pver, & + cappa, gravit, rair, tmelt, pref, lapse_rate, & + top_lev_cloudphys, & + pmid, ps, temp, sst, & + q, cldice, & + phis, & + shallowcu, deepcu, concld, & ! convective cloud cover + landfrac, ocnfrac, snowh, & + rhpert_flag, & ! todo: decide what to do with this + cloud, rhcloud, & + cldst, & + rhu00, icecldf, liqcldf, relhum, & + errmsg, errflg) + + ! To-be-ccppized dependencies + use wv_saturation, only: qsat, qsat_water, svp_ice_vect + + ! Arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(kind_phys), intent(in) :: cappa + real(kind_phys), intent(in) :: gravit + real(kind_phys), intent(in) :: rair + real(kind_phys), intent(in) :: tmelt + real(kind_phys), intent(in) :: pref + real(kind_phys), intent(in) :: lapse_rate + integer, intent(in) :: top_lev_cloudphys ! vertical_layer_index_of_cloud_fraction_top [index] + real(kind_phys), intent(in) :: pmid(:, :) ! air_pressure [Pa] + real(kind_phys), intent(in) :: ps(:) ! surface_air_pressure [Pa] + real(kind_phys), intent(in) :: temp(:, :) ! air_temperature [K] + real(kind_phys), intent(in) :: sst(:) ! sea_surface_temperature [K] + real(kind_phys), intent(in) :: q(:, :) ! water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldice(:, :) ! cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: phis(:) ! surface_geopotential [m2 s-2] + real(kind_phys), intent(in) :: shallowcu(:, :) ! shallow convective cloud fraction [fraction] + real(kind_phys), intent(in) :: deepcu(:, :) ! deep convective cloud fraction [fraction] + real(kind_phys), intent(in) :: concld(:, :) ! convective_cloud_area_fraction [fraction] + real(kind_phys), intent(in) :: landfrac(:) ! land_area_fraction [fraction] + real(kind_phys), intent(in) :: ocnfrac(:) ! ocean_area_fraction [fraction] + real(kind_phys), intent(in) :: snowh(:) ! lwe_surface_snow_depth_over_land [m] + logical, intent(in) :: rhpert_flag ! 0 or 1 to perturb rh + + ! Output arguments + real(kind_phys), intent(out) :: cloud(:, :) ! cloud_area_fraction [fraction] + real(kind_phys), intent(out) :: rhcloud(:, :) ! cloud fraction + real(kind_phys), intent(out) :: cldst(:, :) ! stratiform_cloud_area_fraction [fraction] + real(kind_phys), intent(out) :: rhu00(:, :) ! RH threshold for cloud [fraction] + real(kind_phys), intent(out) :: relhum(:, :) ! RH for prognostic cldwat [fraction] + real(kind_phys), intent(out) :: icecldf(:, :) ! ice cloud fraction [fraction] + real(kind_phys), intent(out) :: liqcldf(:, :) ! liquid cloud fraction (combined into cloud) [fraction] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local variables: + real(kind_phys) :: cld ! intermediate scratch variable (low cld) + real(kind_phys) :: dthdpmn(ncol) ! most stable lapse rate below 750 mb + real(kind_phys) :: dthdp ! lapse rate (intermediate variable) + real(kind_phys) :: es(ncol, pver) ! saturation vapor pressure + real(kind_phys) :: qs(ncol, pver) ! saturation specific humidity + real(kind_phys) :: rhwght ! weighting function for rhlim transition + real(kind_phys) :: rh(ncol, pver) ! relative humidity + real(kind_phys) :: rhdif ! intermediate scratch variable + real(kind_phys) :: strat ! intermediate scratch variable + real(kind_phys) :: theta(ncol, pver) ! potential temperature + real(kind_phys) :: thetas(ncol) ! ocean surface potential temperature + real(kind_phys) :: rhlim ! local rel. humidity threshold estimate + real(kind_phys) :: clrsky(ncol) ! temporary used in random overlap calc + real(kind_phys) :: rpdeli(ncol, pver - 1) ! 1./(pmid(k+1)-pmid(k)) + real(kind_phys) :: rhpert ! the specified perturbation to rh + + logical :: cldbnd(ncol) ! region below high cloud boundary + + integer :: i, ierror, k ! column, level indices + integer :: kp1, ifld + integer :: kdthdp(ncol) + integer :: numkcld ! number of levels in which to allow clouds + + ! In Cloud Ice Content variables + real(kind_phys) :: a, b, c, as, bs, cs ! fit parameters + real(kind_phys) :: Kc ! constant for ice cloud calc (wood & field) + real(kind_phys) :: ttmp ! limited temperature + real(kind_phys) :: icicval ! empirical iwc value + real(kind_phys) :: rho ! local air density + real(kind_phys) :: esl(ncol, pver) ! liq sat vapor pressure + real(kind_phys) :: esi(ncol, pver) ! ice sat vapor pressure + real(kind_phys) :: ncf, phi ! Wilson and Ballard parameters + + ! Statement functions + logical land + land(i) = nint(landfrac(i)) == 1 + + errmsg = '' + errflg = 0 + + !================================================================================== + ! PHILOSOPHY OF PRESENT IMPLEMENTATION + ! Modification to philosophy for ice supersaturation + ! philosophy below is based on RH water only. This is 'liquid condensation' + ! or liquid cloud (even though it will freeze immediately to ice) + ! The idea is that the RH limits for condensation are strict only for + ! water saturation + ! + ! Ice clouds are formed by explicit parameterization of ice nucleation. + ! Closure for ice cloud fraction is done on available cloud ice, such that + ! the in-cloud ice content matches an empirical fit + ! thus, icecldf = min(cldice/icicval,1) where icicval = f(temp,cldice,numice) + ! for a first cut, icicval=f(temp) only. + ! Combined cloud fraction is maximum overlap cloud=max(1,max(icecldf,liqcldf)) + ! No dA/dt term for ice? + ! + ! There are three co-existing cloud types: convective, inversion related low-level + ! stratocumulus, and layered cloud (based on relative humidity). Layered and + ! stratocumulus clouds do not compete with convective cloud for which one creates + ! the most cloud. They contribute collectively to the total grid-box average cloud + ! amount. This is reflected in the way in which the total cloud amount is evaluated + ! (a sum as opposed to a logical "or" operation) + ! + !================================================================================== + ! set defaults for rhu00 - arbitrary number larger than 1.0 (100%) + rhu00(:, :) = 2.0_kind_phys + ! define rh perturbation in order to estimate rhdfda + rhpert = 0.01_kind_phys + + ! Parameters for ice cloud fraction methods + ! References for these parameters, and the corresponding methods, are in the compute section + ! below. + ! iceopt = 1. Wang and Sassen IWC + a = 26.87_kind_phys + b = 0.569_kind_phys + c = 0.002892_kind_phys + ! iceopt = 2. Schiller + as = -68.4202_kind_phys + bs = 0.983917_kind_phys + cs = 2.81795_kind_phys + ! iceopt = 3. Wood and Field + Kc = 75._kind_phys + + ! Evaluate potential temperature and relative humidity + if (cldfrc_ice) then + ! If computing ice cloud fraction (CAM5+ MG Morrison and Gettelman microphysics) then water RH + do k = top_lev_cloudphys, pver + call qsat_water(temp(1:ncol, k), pmid(1:ncol, k), esl(1:ncol, k), qs(1:ncol, k), ncol) + call svp_ice_vect(temp(1:ncol, k), esi(1:ncol, k), ncol) + end do + else + ! If not computing ice cloud fraction (CAM4 RK Rasch-Kristjansson microphysics) then hybrid RH, + do k = top_lev_cloudphys, pver + call qsat(temp(1:ncol, k), pmid(1:ncol, k), es(1:ncol, k), qs(1:ncol, k), ncol) + end do + end if + + cloud = 0._kind_phys + icecldf = 0._kind_phys + liqcldf = 0._kind_phys + rhcloud = 0._kind_phys + cldst = 0._kind_phys + + do k = top_lev_cloudphys, pver + theta(:ncol, k) = temp(:ncol, k)*(pref/pmid(:ncol, k))**cappa + + do i = 1, ncol + if(.not. rhpert_flag) then + ! default: no RH perturbation + rh(i, k) = q(i, k)/qs(i, k) + else + rh(i, k) = q(i, k)/qs(i, k)*(1.0_kind_phys + rhpert) + endif + ! record relhum, rh itself will later be modified related with concld + relhum(i, k) = rh(i, k) + end do + end do + + ! Initialize other temporary variables + ierror = 0 + do i = 1, ncol + ! Adjust thetas(i) in the presence of non-zero ocean heights. + ! This reduces the temperature for positive heights according to a standard lapse rate. + if (ocnfrac(i) > 0.01_kind_phys) then + thetas(i) = (sst(i) - lapse_rate*phis(i)/gravit)*(pref/ps(i))**cappa + endif + if (ocnfrac(i) > 0.01_kind_phys .and. sst(i) < 260._kind_phys) then + ierror = i + endif + end do + + if (ierror > 0) then + write (iulog, *) 'COLDSST: encountered in cldfrc:', ierror, ocnfrac(ierror), sst(ierror) + end if + + do k = top_lev_cloudphys, pver - 1 + rpdeli(:ncol, k) = 1._kind_phys/(pmid(:ncol, k + 1) - pmid(:ncol, k)) + end do + + ! shallow and deep convective cloudiness are calculated in the + ! convective_cloud_cover CCPP scheme. +#ifndef PERGRO + do k = top_lev_cloudphys, pver + do i = 1, ncol + rh(i, k) = (rh(i, k) - concld(i, k))/(1.0_kind_phys - concld(i, k)) + end do + end do +#endif + !================================================================================== + ! + ! ****** Compute layer cloudiness ****** + ! + !==================================================================== + ! Begin the evaluation of layered cloud amount based on (modified) RH + !==================================================================== + ! + numkcld = pver + do k = top_lev_cloudphys + 1, numkcld + kp1 = min(k + 1, pver) + do i = 1, ncol + ! This is now designed to apply FOR LIQUID CLOUDS (condensation > RH water) + cldbnd(i) = pmid(i, k) >= pretop + + if (pmid(i, k) >= premib) then + !============================================================== + ! This is the low cloud (below premib) block + !============================================================== + ! enhance low cloud activation over land with no snow cover + if (land(i) .and. (snowh(i) <= 0.000001_kind_phys)) then + rhlim = rhminl - rhminl_adj_land + else + rhlim = rhminl + end if + + rhdif = (rh(i, k) - rhlim)/(1.0_kind_phys - rhlim) + rhcloud(i, k) = min(0.999_kind_phys, (max(rhdif, 0.0_kind_phys))**2) + + ! SJV: decrease cloud amount if very low water vapor content + ! (thus very cold): "freeze dry" + if (cldfrc_freeze_dry) then + rhcloud(i, k) = rhcloud(i, k)*max(0.15_kind_phys, min(1.0_kind_phys, q(i, k)/0.0030_kind_phys)) + end if + + else if (pmid(i, k) < premit) then + !============================================================== + ! This is the high cloud (above premit) block + !============================================================== + ! + rhlim = rhminh + ! + rhdif = (rh(i, k) - rhlim)/(1.0_kind_phys - rhlim) + rhcloud(i, k) = min(0.999_kind_phys, (max(rhdif, 0.0_kind_phys))**2) + else + !============================================================== + ! This is the middle cloud block + !============================================================== + ! + ! linear rh threshold transition between thresholds for low & high cloud + ! + rhwght = (premib - (max(pmid(i, k), premit)))/(premib - premit) + + if (land(i) .and. (snowh(i) <= 0.000001_kind_phys)) then + rhlim = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_kind_phys - rhwght) + else + rhlim = rhminh*rhwght + rhminl*(1.0_kind_phys - rhwght) + end if + rhdif = (rh(i, k) - rhlim)/(1.0_kind_phys - rhlim) + rhcloud(i, k) = min(0.999_kind_phys, (max(rhdif, 0.0_kind_phys))**2) + end if + !================================================================================== + ! WE NEED TO DOCUMENT THE PURPOSE OF THIS TYPE OF CODE (ASSOCIATED WITH 2ND CALL) + !================================================================================== + ! ! + ! ! save rhlim to rhu00, it handles well by itself for low/high cloud + ! ! + rhu00(i, k) = rhlim + !================================================================================== + + if (cldfrc_ice) then + ! Evaluate ice cloud fraction based on in-cloud ice content + if (iceopt .eq. 1 .or. iceopt .eq. 2) then + if (iceopt .eq. 1) then + ! ICE CLOUD OPTION 1 - Wang & Sassen 2002 + ! Evaluate desired in-cloud water content + ! icicval = f(temp,cldice,numice) + ! Start with a function of temperature. + ! Wang & Sassen 2002 (JAS), based on ARM site MMCR (midlat cirrus) + ! https://doi.org/10.1175/1520-0469(2002)059<2291:CCMPRU>2.0.CO;2 + ! parameterization valid for 203-253K + ! icival > 0 for t>195K + ttmp = max(195._kind_phys, min(temp(i, k), 253._kind_phys)) - 273.16_kind_phys + icicval = a + b*ttmp + c*ttmp**2._kind_phys + !convert units + rho = pmid(i, k)/(rair*temp(i, k)) + icicval = icicval*1.e-6_kind_phys/rho + else + ! ICE CLOUD OPTION 2 - Schiller 2008 (JGR) + ! https://doi.org/10.1029/2008JD010342 + ! Use a curve based on FISH measurements in + ! tropics, mid-lats and arctic. Curve is for 180-250K (raise to 273K?) + ! use median all flights + ttmp = max(190._kind_phys, min(temp(i, k), 273.16_kind_phys)) + icicval = 10._kind_phys**(as*bs**ttmp + cs) + ! convert units from ppmv to kg kg-1 + icicval = icicval*1.e-6_kind_phys*18._kind_phys/28.97_kind_phys + end if + + ! Set ice cloud fraction for options 1 or 2 + icecldf(i, k) = max(0._kind_phys, min(cldice(i, k)/icicval, 1._kind_phys)) + else if (iceopt .eq. 3) then + ! ICE CLOUD OPTION 3 - Wood & Field 2000 (JAS) + ! https://doi.org/10.1175/1520-0469(2000)057<1888:RBTWCW>2.0.CO;2 + ! eq 6: cloud fraction = 1 - exp (-K * qc/qsati) + icecldf(i, k) = 1._kind_phys - exp(-Kc*cldice(i, k)/(qs(i, k)*(esi(i, k)/esl(i, k)))) + icecldf(i, k) = max(0._kind_phys, min(icecldf(i, k), 1._kind_phys)) + else + ! ICE CLOUD OPTION 4 - Wilson and Ballard 1999 + ! https://doi.org/10.1002/qj.49712555707 + ! inversion of smith.... + ! ncf = cldice / ((1-RHcrit)*qs) + ! then a function of ncf.... + ncf = cldice(i, k)/((1._kind_phys - icecrit)*qs(i, k)) + if (ncf <= 0._kind_phys) then + icecldf(i, k) = 0._kind_phys + else if (ncf > 0._kind_phys .and. ncf <= 1._kind_phys/6._kind_phys) then + icecldf(i, k) = 0.5_kind_phys*(6._kind_phys*ncf)**(2._kind_phys/3._kind_phys) + else if (ncf > 1._kind_phys/6._kind_phys .and. ncf < 1._kind_phys) then + phi = (acos(3._kind_phys*(1._kind_phys - ncf)/2._kind_phys**(3._kind_phys/2._kind_phys)) + 4._kind_phys*3.1415927_kind_phys)/3._kind_phys + icecldf(i, k) = (1._kind_phys - 4._kind_phys*cos(phi)*cos(phi)) + else + icecldf(i, k) = 1._kind_phys + end if + icecldf(i, k) = max(0._kind_phys, min(icecldf(i, k), 1._kind_phys)) + end if + + ! Test code: Combine ice and liquid cloud fraction assuming maximum overlap. + ! Combined cloud fraction is maximum overlap + ! cloud(i,k)=min(1._kind_phys,max(icecldf(i,k),rhcloud(i,k))) + + liqcldf(i, k) = (1._kind_phys - icecldf(i, k))*rhcloud(i, k) + cloud(i, k) = liqcldf(i, k) + icecldf(i, k) + else ! cldfrc_ice = .false. + ! For RK microphysics + cloud(i, k) = rhcloud(i, k) + end if + end do + end do + + ! + ! Add in the marine strat + ! MARINE STRATUS SHOULD BE A SPECIAL CASE OF LAYERED CLOUD + ! CLOUD CURRENTLY CONTAINS LAYERED CLOUD DETERMINED BY RH CRITERIA + ! TAKE THE MAXIMUM OF THE DIAGNOSED LAYERED CLOUD OR STRATOCUMULUS + ! + !=================================================================================== + ! + ! SOME OBSERVATIONS ABOUT THE FOLLOWING SECTION OF CODE (missed in earlier look) + ! K700 IS SET AS A CONSTANT BASED ON HYBRID COORDINATE: IT DOES NOT DEPEND ON + ! LOCAL PRESSURE; THERE IS NO PRESSURE RAMP => LOOKS LEVEL DEPENDENT AND + ! DISCONTINUOUS IN SPACE (I.E., STRATUS WILL END SUDDENLY WITH NO TRANSITION) + ! + ! IT APPEARS THAT STRAT IS EVALUATED ACCORDING TO KLEIN AND HARTMANN; HOWEVER, + ! THE ACTUAL STRATUS AMOUNT (CLDST) APPEARS TO DEPEND DIRECTLY ON THE RH BELOW + ! THE STRONGEST PART OF THE LOW LEVEL INVERSION. + !PJR answers: 1) the rh limitation is a physical/mathematical limitation + ! cant have more cloud than there is RH + ! allowed the cloud to exist two layers below the inversion + ! because the numerics frequently make 50% relative humidity + ! in level below the inversion which would allow no cloud + ! 2) since the cloud is only allowed over ocean, it should + ! be very insensitive to surface pressure (except due to + ! spectral ringing, which also causes so many other problems + ! I didnt worry about it. + ! + !================================================================================== + if (.not. inversion_cld_off) then + ! + ! Find most stable level below 750 mb for evaluating stratus regimes + ! + do i = 1, ncol + ! Nothing triggers unless a stability greater than this minimum threshold is found + dthdpmn(i) = -0.125_kind_phys + kdthdp(i) = 0 + end do + + do k = top_lev_cloudphys + 1, pver + do i = 1, ncol + if (pmid(i, k) >= premib .and. ocnfrac(i) > 0.01_kind_phys) then + ! I think this is done so that dtheta/dp is in units of dg/mb (JJH) + dthdp = 100.0_kind_phys*(theta(i, k) - theta(i, k - 1))*rpdeli(i, k - 1) + if (dthdp < dthdpmn(i)) then + dthdpmn(i) = dthdp + kdthdp(i) = k ! index of interface of max inversion + end if + end if + end do + end do + + ! Also check between the bottom layer and the surface + ! Only perform this check if the criteria were not met above + do i = 1, ncol + if (kdthdp(i) .eq. 0 .and. ocnfrac(i) > 0.01_kind_phys) then + dthdp = 100.0_kind_phys*(thetas(i) - theta(i, pver))/(ps(i) - pmid(i, pver)) + if (dthdp < dthdpmn(i)) then + dthdpmn(i) = dthdp + kdthdp(i) = pver ! index of interface of max inversion + end if + end if + end do + + do i = 1, ncol + if (kdthdp(i) /= 0) then + k = kdthdp(i) + kp1 = min(k + 1, pver) + ! Note: strat will be zero unless ocnfrac > 0.01 + strat = min(1._kind_phys, max(0._kind_phys, ocnfrac(i)*((theta(i, k700) - thetas(i))*.057_kind_phys - .5573_kind_phys))) + ! + ! assign the stratus to the layer just below max inversion + ! the relative humidity changes so rapidly across the inversion + ! that it is not safe to just look immediately below the inversion + ! so limit the stratus cloud by rh in both layers below the inversion + ! + cldst(i, k) = min(strat, max(rh(i, k), rh(i, kp1))) + end if + end do + end if ! .not.inversion_cld_off + + do k = top_lev_cloudphys, pver + do i = 1, ncol + ! which is greater; standard layered cloud amount or stratocumulus diagnosis + cloud(i, k) = max(rhcloud(i, k), cldst(i, k)) + + ! add in the contributions of convective cloud (determined separately and accounted + ! for by modifications to the large-scale relative humidity. + cloud(i, k) = min(cloud(i, k) + concld(i, k), 1.0_kind_phys) + end do + end do + + end subroutine compute_cloud_fraction_run +end module compute_cloud_fraction diff --git a/schemes/cloud_fraction/compute_cloud_fraction.meta b/schemes/cloud_fraction/compute_cloud_fraction.meta new file mode 100644 index 00000000..c314e2cf --- /dev/null +++ b/schemes/cloud_fraction/compute_cloud_fraction.meta @@ -0,0 +1,322 @@ +[ccpp-table-properties] + name = compute_cloud_fraction + type = scheme + dependencies = ../../to_be_ccppized/wv_saturation.F90 + +[ccpp-arg-table] + name = compute_cloud_fraction_init + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ pref_mid ] + standard_name = reference_pressure_in_atmosphere_layer + units = Pa + type = real | kind = kind_phys + dimensions = (vertical_layer_dimension) + intent = in +[ inversion_cld_off_in ] + standard_name = do_no_stratification_based_cloud_fraction + units = flag + type = logical + dimensions = () + intent = in +[ cldfrc_freeze_dry_in ] + standard_name = do_vavrus_freeze_dry_adjustment_for_cloud_fraction + units = flag + type = logical + dimensions = () + intent = in +[ cldfrc_ice_in ] + standard_name = do_ice_cloud_fraction_for_cloud_fraction + units = flag + type = logical + dimensions = () + intent = in +[ iceopt_in ] + standard_name = control_for_ice_cloud_fraction + units = 1 + type = integer + dimensions = () + intent = in +[ rhminl_in ] + standard_name = tunable_parameter_for_minimum_relative_humidity_for_low_stable_clouds_for_cloud_fraction + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rhminl_adj_land_in ] + standard_name = tunable_parameter_for_adjustment_to_minimum_relative_humidity_for_low_stable_clouds_for_land_without_snow_cover_for_cloud_fraction + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rhminh_in ] + standard_name = tunable_parameter_for_minimum_relative_humidity_for_high_stable_clouds_for_cloud_fraction + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ premit_in ] + standard_name = tunable_parameter_for_top_pressure_bound_for_mid_level_clouds_for_cloud_fraction + units = Pa + type = real | kind = kind_phys + dimensions = () + intent = in +[ premib_in ] + standard_name = tunable_parameter_for_bottom_pressure_bound_for_mid_level_liquid_stratus_for_cloud_fraction + units = Pa + type = real | kind = kind_phys + dimensions = () + intent = in +[ icecrit_in ] + standard_name = tunable_parameter_for_critical_relative_humidity_for_ice_clouds_for_cloud_fraction_using_wilson_and_ballard_scheme + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = compute_cloud_fraction_timestep_init + type = scheme +[ rhpert_flag ] + standard_name = do_relative_humidity_perturbation_for_cloud_fraction + units = flag + type = logical + dimensions = () + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = compute_cloud_fraction_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ cappa ] + standard_name = ratio_of_dry_air_gas_constant_to_specific_heat_of_dry_air_at_constant_pressure + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rair ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ tmelt ] + standard_name = freezing_point_of_water + units = K + type = real | kind = kind_phys + dimensions = () + intent = in +[ pref ] + standard_name = surface_reference_pressure + units = Pa + type = real | kind = kind_phys + dimensions = () + intent = in +[ lapse_rate ] + standard_name = reference_temperature_lapse_rate + units = K m-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ top_lev_cloudphys ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ pmid ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ ps ] + standard_name = surface_air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ temp ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ sst ] + standard_name = sea_surface_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ q ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldice ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ phis ] + standard_name = surface_geopotential + units = m2 s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ shallowcu ] + standard_name = shallow_convective_cloud_area_fraction_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ deepcu ] + standard_name = deep_convective_cloud_area_fraction_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ concld ] + standard_name = convective_cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ landfrac ] + standard_name = land_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ ocnfrac ] + standard_name = ocean_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snowh ] + standard_name = lwe_surface_snow_depth_over_land + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ rhpert_flag ] + standard_name = do_relative_humidity_perturbation_for_cloud_fraction + units = flag + type = logical + dimensions = () + intent = in +[ cloud ] + standard_name = cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ rhcloud ] + standard_name = cloud_area_fraction_from_relative_humidity_method_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ cldst ] + standard_name = stratiform_cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ rhu00 ] + standard_name = relative_humidity_threshold_for_prognostic_cloud_water_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ icecldf ] + standard_name = stratiform_cloud_ice_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ liqcldf ] + standard_name = stratiform_cloud_liquid_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ relhum ] + standard_name = relative_humidity_for_prognostic_cloud_water_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml b/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml new file mode 100644 index 00000000..18c511bd --- /dev/null +++ b/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml @@ -0,0 +1,229 @@ + + + + + + + + + logical + cldfrc + cldfrc_nl + do_vavrus_freeze_dry_adjustment_for_cloud_fraction + flag + + Turn on for Vavrus "freeze dry" adjustment in cloud fraction. + Default: TRUE + + + .true. + + + + logical + cldfrc + cldfrc_nl + do_ice_cloud_fraction_for_cloud_fraction + flag + + Turn on ice cloud fraction calculation + Default: TRUE (CAM5, 6, 7) or FALSE + + + .false. + .true. + .true. + + + + logical + cldfrc + cldfrc_nl + do_no_stratification_based_cloud_fraction + flag + + Turn off stratification-based cloud fraction + Default: FALSE + + + .false. + + + + + real + cldfrc + cldfrc_nl + tunable_parameter_for_minimum_relative_humidity_for_low_stable_clouds_for_cloud_fraction + 1 + + Minimum rh for low stable clouds. + + + 0.900D0 + 0.8975D0 + 0.8975D0 + + + + real + cldfrc + cldfrc_nl + tunable_parameter_for_adjustment_to_minimum_relative_humidity_for_low_stable_clouds_for_land_without_snow_cover_for_cloud_fraction + 1 + + Adjustment to rhminl for land without snow cover. + Default: 0.0 for CAM6, CAM7; 0.10 for all others + + + 0.1D0 + 0.0D0 + + + + real + cldfrc + cldfrc_nl + tunable_parameter_for_minimum_relative_humidity_for_high_stable_clouds_for_cloud_fraction + 1 + + Minimum rh for high stable clouds. + + + 0.800D0 + + + + real + cldfrc + cldfrc_nl + tunable_parameter_for_top_pressure_bound_for_mid_level_clouds_for_cloud_fraction + Pa + + Top pressure bound for mid level cloud. + + + 75000.0D0 + + + + real + cldfrc + cldfrc_nl + tunable_parameter_for_bottom_pressure_bound_for_mid_level_liquid_stratus_for_cloud_fraction + Pa + + Bottom pressure bound for mid-level liquid stratus fraction. + + + 750.0D2 + 700.0D2 + 700.0D2 + + + + integer + cldfrc + cldfrc_nl + control_for_ice_cloud_fraction + 1 + + Scheme for ice cloud fraction. + 1 = Wang and Sassen (https://doi.org/10.1175/1520-0469(2002)059<2291:CCMPRU>2.0.CO;2) + 2 = Schiller (iciwc) (https://doi.org/10.1029/2008JD010342) + 3 = Wood and Field (https://doi.org/10.1175/1520-0469(2000)057<1888:RBTWCW>2.0.CO;2) + 4 = Wilson (based on Smith) (https://doi.org/10.1002/qj.49712555707) + 5 = modified Slingo (ssat and empty cloud) + + + 1 + 5 + 5 + + + + real + cldfrc + cldfrc_nl + tunable_parameter_for_critical_relative_humidity_for_ice_clouds_for_cloud_fraction_using_wilson_and_ballard_scheme + 1 + + Critical RH for ice clouds using the Wilson and Ballard scheme + + + 0.95D0 + 0.93D0 + 0.93D0 + + + diff --git a/schemes/cloud_fraction/convective_cloud_cover.F90 b/schemes/cloud_fraction/convective_cloud_cover.F90 new file mode 100644 index 00000000..97730edd --- /dev/null +++ b/schemes/cloud_fraction/convective_cloud_cover.F90 @@ -0,0 +1,121 @@ +! Computes deep and shallow convective cloud fractions +! CCPP-ized: Haipeng Lin, February 2025 +module convective_cloud_cover + + use ccpp_kinds, only: kind_phys + + implicit none + private + save + + ! public CCPP-compliant subroutines + public :: convective_cloud_cover_init + public :: convective_cloud_cover_run + + ! Module variables for tuning parameters + real(kind_phys) :: sh1 ! Shallow convection tuning parameter 1 [fraction] + real(kind_phys) :: sh2 ! Shallow convection tuning parameter 2 [m2 s kg-1] + real(kind_phys) :: dp1 ! Deep convection tuning parameter 1 [fraction] + real(kind_phys) :: dp2 ! Deep convection tuning parameter 2 [m2 s kg-1] + +contains + +!> \section arg_table_convective_cloud_cover_init Argument Table +!! \htmlinclude convective_cloud_cover_init.html + subroutine convective_cloud_cover_init( & + amIRoot, iulog, & + sh1_in, sh2_in, dp1_in, dp2_in, errmsg, errflg) + + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog ! log output unit + real(kind_phys), intent(in) :: sh1_in ! Shallow convection parameter 1 + real(kind_phys), intent(in) :: sh2_in ! Shallow convection parameter 2 + real(kind_phys), intent(in) :: dp1_in ! Deep convection parameter 1 + real(kind_phys), intent(in) :: dp2_in ! Deep convection parameter 2 + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + ! Set module variables from namelist + sh1 = sh1_in + sh2 = sh2_in + dp1 = dp1_in + dp2 = dp2_in + + if(amIRoot) then + write(iulog,*) 'tuning parameters convective_cloud_cover: dp1 ',dp1,' dp2 ',dp2,' sh1 ',sh1,' sh2 ',sh2 + endif + + end subroutine convective_cloud_cover_init + + ! Compute convective cloud cover (deep and shallow) + ! Should produce typical numbers of 20% + ! Shallow and deep convective cloudiness are evaluated separately and summed + ! because the processes are evaluated separately. +!> \section arg_table_convective_cloud_cover_run Argument Table +!! \htmlinclude convective_cloud_cover_run.html + subroutine convective_cloud_cover_run( & + ncol, pver, & + top_lev_cloudphys, & + use_shfrc, shfrc, & + cmfmc_total, cmfmc_sh, & + shallowcu, deepcu, concld, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: top_lev_cloudphys ! Top vertical level for cloud physics [index] + + logical, intent(in) :: use_shfrc ! Use cloud area fraction provided by shallow convection? [flag] + real(kind_phys), intent(in) :: shfrc(:, :) ! Input shallow cloud fraction [fraction] + + real(kind_phys), intent(in) :: cmfmc_total(:, :) ! atmosphere_convective_mass_flux_due_to_all_convection [kg m-2 s-1] + real(kind_phys), intent(in) :: cmfmc_sh(:, :) ! atmosphere_convective_mass_flux_due_to_shallow_convection [kg m-2 s-1] + + ! Output arguments + real(kind_phys), intent(out) :: shallowcu(:, :) ! Shallow convective cloud fraction [fraction] + real(kind_phys), intent(out) :: deepcu(:, :) ! Deep convective cloud fraction [fraction] + real(kind_phys), intent(out) :: concld(:, :) ! Convective cloud cover [fraction] + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables + integer :: i, k + + errmsg = '' + errflg = 0 + + concld(:ncol, :) = 0.0_kind_phys + do k = top_lev_cloudphys, pver + do i = 1, ncol + if (.not. use_shfrc) then + ! Compute the shallow convection cloud cover using the shallow convective mass flux. + shallowcu(i, k) = max(0.0_kind_phys, & + min(sh1*log(1.0_kind_phys + sh2*cmfmc_sh(i, k + 1)), 0.30_kind_phys)) + else + ! If use shallow convective calculated clouds, then just assign + shallowcu(i, k) = shfrc(i, k) + end if + + ! REMOVECAM: This could be changed to use deep convective mass flux + ! since it is independently available in CCPP, once CAM is retired + deepcu(i, k) = max(0.0_kind_phys, & + min(dp1 * & + log(1.0_kind_phys + & + dp2 * (cmfmc_total(i, k + 1) - cmfmc_sh(i, k + 1)) & + ), & + 0.60_kind_phys)) + + ! Estimate of local convective cloud cover based on convective mass flux + ! Modify local large-scale relative humidity to account for presence of + ! convective cloud when evaluating relative humidity based layered cloud amount + concld(i, k) = min(shallowcu(i, k) + deepcu(i, k), 0.80_kind_phys) + end do + end do + + end subroutine convective_cloud_cover_run + +end module convective_cloud_cover diff --git a/schemes/cloud_fraction/convective_cloud_cover.meta b/schemes/cloud_fraction/convective_cloud_cover.meta new file mode 100644 index 00000000..7340b8f6 --- /dev/null +++ b/schemes/cloud_fraction/convective_cloud_cover.meta @@ -0,0 +1,131 @@ +[ccpp-table-properties] + name = convective_cloud_cover + type = scheme + +[ccpp-arg-table] + name = convective_cloud_cover_init + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ sh1_in ] + standard_name = tunable_parameter_for_shallow_convection_1_for_cloud_fraction + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ sh2_in ] + standard_name = tunable_parameter_for_shallow_convection_2_for_cloud_fraction + units = s kg-1 m2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ dp1_in ] + standard_name = tunable_parameter_for_deep_convection_1_for_cloud_fraction + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ dp2_in ] + standard_name = tunable_parameter_for_deep_convection_2_for_cloud_fraction + units = s kg-1 m2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = convective_cloud_cover_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ top_lev_cloudphys ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ use_shfrc ] + standard_name = flag_for_cloud_area_fraction_to_use_shallow_convection_calculated_cloud_area_fraction + units = flag + type = logical + dimensions = () + intent = in +[ shfrc ] + standard_name = shallow_convective_cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cmfmc_total ] + standard_name = atmosphere_convective_mass_flux_due_to_all_convection + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ cmfmc_sh ] + standard_name = atmosphere_convective_mass_flux_due_to_shallow_convection + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ shallowcu ] + standard_name = shallow_convective_cloud_area_fraction_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ deepcu ] + standard_name = deep_convective_cloud_area_fraction_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ concld ] + standard_name = convective_cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/cloud_fraction/convective_cloud_cover_namelist.xml b/schemes/cloud_fraction/convective_cloud_cover_namelist.xml new file mode 100644 index 00000000..5efac053 --- /dev/null +++ b/schemes/cloud_fraction/convective_cloud_cover_namelist.xml @@ -0,0 +1,138 @@ + + + + + + + + + real + kind_phys + cldfrc + convective_cloud_cover_nl + tunable_parameter_for_shallow_convection_1_for_cloud_fraction + 1 + + Tunable parameter 1 for shallow convection cloud fraction. + + + 0.07D0 + 0.04D0 + + + + real + kind_phys + cldfrc + convective_cloud_cover_nl + tunable_parameter_for_shallow_convection_2_for_cloud_fraction + s kg-1 m2 + + Tunable parameter 2 for shallow convection cloud fraction. + + + 500.0D0 + + + + real + kind_phys + cldfrc + convective_cloud_cover_nl + tunable_parameter_for_deep_convection_1_for_cloud_fraction + 1 + + Tunable parameter 1 for deep convection cloud fraction. + + + 0.14D0 + 0.10D0 + 0.10D0 + 0.10D0 + + + + real + kind_phys + cldfrc + convective_cloud_cover_nl + tunable_parameter_for_deep_convection_2_for_cloud_fraction + s kg-1 m2 + + Tunable parameter 2 for deep convection cloud fraction. + + + 500.0D0 + + + diff --git a/schemes/cloud_fraction/set_cloud_fraction_top.F90 b/schemes/cloud_fraction/set_cloud_fraction_top.F90 index f30e1045..a13bce30 100644 --- a/schemes/cloud_fraction/set_cloud_fraction_top.F90 +++ b/schemes/cloud_fraction/set_cloud_fraction_top.F90 @@ -1,6 +1,3 @@ -! Copyright (C) 2024 National Science Foundation-National Center for Atmospheric Research -! SPDX-License-Identifier: Apache-2.0 -! ! Stub scheme to set top of cloud physics to below top cloud level. ! Used for all macrophysical schemes except RK. module set_cloud_fraction_top diff --git a/schemes/dry_adiabatic_adjust/dadadj_namelist.xml b/schemes/dry_adiabatic_adjust/dadadj_namelist.xml index 27504b75..020aeaa2 100644 --- a/schemes/dry_adiabatic_adjust/dadadj_namelist.xml +++ b/schemes/dry_adiabatic_adjust/dadadj_namelist.xml @@ -67,7 +67,7 @@ to this variable when the namelist generator is run. If the tag has no attributes, it is a default value. In general, the namelist generator attempts to find a value with the - maximium number of attribute matches (and no non-matches). + maximum number of attribute matches (and no non-matches). standard_name This is the CCPP Standard Name of the variable diff --git a/schemes/hack_shallow/convect_shallow_sum_to_deep.F90 b/schemes/hack_shallow/convect_shallow_sum_to_deep.F90 index 988c6707..1ca4618b 100644 --- a/schemes/hack_shallow/convect_shallow_sum_to_deep.F90 +++ b/schemes/hack_shallow/convect_shallow_sum_to_deep.F90 @@ -1,6 +1,3 @@ -! Copyright (C) 2024-2025 National Science Foundation-National Center for Atmospheric Research -! SPDX-License-Identifier: Apache-2.0 -! ! Note: There is an implicit dependency in the shallow convective code ! that deep convection runs *before* shallow convection. ! Efforts have been made to make this dependency clear in the form that diff --git a/schemes/hack_shallow/hack_convect_shallow.F90 b/schemes/hack_shallow/hack_convect_shallow.F90 index 9ad58058..802a42e8 100644 --- a/schemes/hack_shallow/hack_convect_shallow.F90 +++ b/schemes/hack_shallow/hack_convect_shallow.F90 @@ -1,6 +1,3 @@ -! Copyright (C) 2024-2025 National Science Foundation-National Center for Atmospheric Research -! SPDX-License-Identifier: Apache-2.0 -! ! Hack shallow convective scheme. ! The main subroutine was formerly named "cmfmca", and its initialization "mfinti". ! diff --git a/schemes/hack_shallow/set_general_conv_fluxes_to_shallow.F90 b/schemes/hack_shallow/set_general_conv_fluxes_to_shallow.F90 index e945f878..028cf8cc 100644 --- a/schemes/hack_shallow/set_general_conv_fluxes_to_shallow.F90 +++ b/schemes/hack_shallow/set_general_conv_fluxes_to_shallow.F90 @@ -1,5 +1,3 @@ -! Copyright (C) 2025 National Science Foundation-National Center for Atmospheric Research -! SPDX-License-Identifier: Apache-2.0 module set_general_conv_fluxes_to_shallow use ccpp_kinds, only: kind_phys diff --git a/schemes/hack_shallow/set_shallow_conv_fluxes_to_general.F90 b/schemes/hack_shallow/set_shallow_conv_fluxes_to_general.F90 index fc1747c1..a18ac019 100644 --- a/schemes/hack_shallow/set_shallow_conv_fluxes_to_general.F90 +++ b/schemes/hack_shallow/set_shallow_conv_fluxes_to_general.F90 @@ -1,5 +1,3 @@ -! Copyright (C) 2025 National Science Foundation-National Center for Atmospheric Research -! SPDX-License-Identifier: Apache-2.0 module set_shallow_conv_fluxes_to_general use ccpp_kinds, only: kind_phys diff --git a/schemes/rasch_kristjansson/cloud_particle_sedimentation.F90 b/schemes/rasch_kristjansson/cloud_particle_sedimentation.F90 new file mode 100644 index 00000000..5732856a --- /dev/null +++ b/schemes/rasch_kristjansson/cloud_particle_sedimentation.F90 @@ -0,0 +1,671 @@ +! Compute tendencies from sedimentation of cloud liquid and ice particles +! Original authors: Byron Boville, September 2002 from code by P.J. Rasch +! CCPP-ized: Haipeng Lin, January 2025 +module cloud_particle_sedimentation + + use ccpp_kinds, only: kind_phys + + implicit none + private + save + + ! public CCPP-compliant subroutines + public :: cloud_particle_sedimentation_init + public :: cloud_particle_sedimentation_run + + ! tuning parameters for cloud liquid and ice particle sedimentation + real(kind_phys), parameter :: vland = 1.5_kind_phys ! liquid fall velocity over land [cm s-1] + real(kind_phys), parameter :: vocean = 2.8_kind_phys ! liquid fall velocity over ocean [cm s-1] + real(kind_phys), parameter :: mxsedfac = 0.99_kind_phys ! maximum sedimentation flux factor + + ! use Stokes velocity method + ! or McFarquhar and Heymsfield (https://doi.org/10.1175/1520-0469(1997)054<2187:POTCIC>2.0.CO;2) + logical, parameter :: stokes = .true. ! use Stokes velocity instead of McFarquhar and Heymsfield + + ! tuning parameters for Stokes terminal velocity + real(kind_phys) :: cldsed_ice_stokes_fac ! factor applied to ice fall vel + ! from Stokes terminal velocity + real(kind_phys), parameter :: eta = 1.7e-5_kind_phys ! viscosity of air [kg m s-1] + real(kind_phys), parameter :: r40 = 40._kind_phys ! 40 micron radius + real(kind_phys), parameter :: r400 = 400._kind_phys ! 400 micron radius + real(kind_phys), parameter :: v400 = 1.00_kind_phys ! fall velocity of 400 micron sphere [m s-1] + real(kind_phys) :: v40 ! Stokes fall velocity of 40 micron sphere [m s-1] + real(kind_phys) :: vslope ! linear slope for large particles [m s-1 micron-1] + + ! tuning parameters for McFarquhar and Heymsfield terminal velocity + real(kind_phys), parameter :: vice_small = 1._kind_phys ! ice fall velocity for small concentration [cm s-1] + real(kind_phys) :: lbound ! lower bound for iciwc + + ! misc + real(kind_phys), parameter :: um_to_m = 1.e-12_kind_phys + +contains + +!> \section arg_table_cloud_particle_sedimentation_init Argument Table +!! \htmlinclude arg_table_cloud_particle_sedimentation_init.html + subroutine cloud_particle_sedimentation_init(& + amIRoot, iulog, & + cldsed_ice_stokes_fac_in, & + rhoh2o, gravit, & + errmsg, errflg) + + ! Input arguments + logical, intent(in) :: amIRoot + integer, intent(in) :: iulog ! log output unit + real(kind_phys), intent(in) :: cldsed_ice_stokes_fac_in + real(kind_phys), intent(in) :: gravit + real(kind_phys), intent(in) :: rhoh2o + + ! Output arguments + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local variables + real(kind_phys) :: ac, bc, cc ! constants for McFarquhar and Heymsfield method + + errmsg = '' + errflg = 0 + cldsed_ice_stokes_fac = cldsed_ice_stokes_fac_in + + ! linear ramp variables for fall velocity + v40 = (2._kind_phys/9._kind_phys)*rhoh2o*gravit/eta*r40**2*1.e-12_kind_phys + vslope = (v400 - v40)/(r400 - r40) + + ! McFarquhar and Heymsfield lower bound for iciwc + cc = 128.64_kind_phys + bc = 53.242_kind_phys + ac = 5.4795_kind_phys + lbound = (-bc + sqrt(bc*bc - 4*ac*cc))/(2*ac) + lbound = 10._kind_phys**lbound + + if(amIRoot) then + write(iulog,*) 'cloud_particle_sedimentation_init: cldsed_ice_stokes_fac = ', cldsed_ice_stokes_fac + endif + + end subroutine cloud_particle_sedimentation_init + + ! Compute gravitational sedimentation velocities for cloud liquid water + ! and ice, based on Lawrence and Crutzen (1998). + ! https://doi.org/10.3402/tellusb.v50i3.16129 + ! + ! and apply cloud particle gravitational sedimentation to condensate. + ! + ! Original authors: mgl, March 1998 (MATCH-MPIC version 2.0) + ! adapted by P.J. Rasch + ! B.A. Boville, September 2002 + ! P.J. Rasch, May 2003 +!> \section arg_table_cloud_particle_sedimentation_run Argument Table +!! \htmlinclude arg_table_cloud_particle_sedimentation_run.html + subroutine cloud_particle_sedimentation_run( & + ncol, pver, pverp, dtime, & + tmelt, gravit, latvap, latice, rair, rhoh2o, & + icritc, & ! from prognostic_cloud_water namelist + pint, pmid, pdel, t, cloud, & + icefrac, landfrac, ocnfrac, & + cldliq, cldice, snowh, landm, & + ! below output for sedimentation velocity + pvliq, pvice, & + ! below output for sedimentation to condensate + liqtend, icetend, wvtend, htend, sfliq, sfice, & + errmsg, errflg) + ! To-be-CCPPized dependencies + use cloud_optical_properties, only: reltab, reitab + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: pverp + real(kind_phys), intent(in) :: dtime + real(kind_phys), intent(in) :: tmelt + real(kind_phys), intent(in) :: gravit + real(kind_phys), intent(in) :: latvap + real(kind_phys), intent(in) :: latice + real(kind_phys), intent(in) :: rair + real(kind_phys), intent(in) :: rhoh2o + real(kind_phys), intent(in) :: icritc ! tunable_parameter_for_autoconversion_of_cold_ice_for_rk_microphysics [kg kg-1] + real(kind_phys), intent(in) :: pint(:,:) ! air_pressure_at_interface [Pa] + real(kind_phys), intent(in) :: pmid(:,:) ! air_pressure [Pa] + real(kind_phys), intent(in) :: pdel(:,:) ! air_pressure_thickness [Pa] + real(kind_phys), intent(in) :: t(:,:) ! air_temperature [K] + real(kind_phys), intent(in) :: cloud(:,:) ! cloud_area_fraction [fraction] + real(kind_phys), intent(in) :: icefrac(:) ! sea_ice_area_fraction [fraction] + real(kind_phys), intent(in) :: landfrac(:) ! land_area_fraction [fraction] + real(kind_phys), intent(in) :: ocnfrac(:) ! ocean_area_fraction [fraction] + real(kind_phys), intent(in) :: cldliq(:,:) ! adv: cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldice(:,:) ! adv: cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: snowh(:) ! lwe_surface_snow_depth_over_land [m] + real(kind_phys), intent(in) :: landm(:) ! smoothed_land_area_fraction [fraction] + + ! Output arguments + ! note: pvliq, pvice are at the interfaces (loss from cell is based on pvel(k+1)) + real(kind_phys), intent(out) :: pvliq(:,:) ! vertical velocity of cloud liquid drops [Pa s-1] + real(kind_phys), intent(out) :: pvice(:,:) ! vertical velocity of cloud ice particles [Pa s-1] + + real(kind_phys), intent(out) :: liqtend(:,:) ! liquid condensate tendency [kg kg-1 s-1] -- to apply cldliq tendency + real(kind_phys), intent(out) :: icetend(:,:) ! ice condensate tendency [kg kg-1 s-1] -- to apply cldice tendency + real(kind_phys), intent(out) :: wvtend(:,:) ! water vapor tendency [kg kg-1 s-1] -- to apply wv tendency + real(kind_phys), intent(out) :: htend(:,:) ! heating rate [J kg-1 s-1] -- to apply s tendency + real(kind_phys), intent(out) :: sfliq(:) ! surface flux of liquid (rain) [kg m-2 s-1] + real(kind_phys), intent(out) :: sfice(:) ! lwe_cloud_ice_sedimentation_rate_at_surface_due_to_microphysics [m s-1] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local variables + real(kind_phys) :: rho(ncol, pver) ! air density [kg m-3] + real(kind_phys) :: vfall ! settling velocity of cloud particles [m s-1] + real(kind_phys) :: icice ! in cloud ice water content [kg kg-1] + real(kind_phys) :: iciwc ! in cloud ice water content [g m-3] + real(kind_phys) :: icefac + real(kind_phys) :: logiwc + real(kind_phys) :: rei(ncol, pver) ! effective radius of ice particles [um] + real(kind_phys) :: rel(ncol, pver) ! effective radius of liquid particles [um] + real(kind_phys) :: fxliq(ncol, pverp) ! fluxes at interface, liquid (positive is down) [Pa] + real(kind_phys) :: fxice(ncol, pverp) ! fluxes at interface, ice (positive is down) [Pa] + real(kind_phys) :: cldab(ncol) ! cloud in layer above [fraction] + real(kind_phys) :: evapliq ! evaporation of cloud liquid into environment [kg kg-1 s-1] + real(kind_phys) :: evapice ! evaporation of cloud ice into environment [kg kg-1 s-1] + real(kind_phys) :: cldovrl ! cloud overlap factor + integer :: i, k + + errmsg = '' + errflg = 0 + + !-------------------------------------------------- + ! COMPUTE GRAVITATIONAL SEDIMENTATION VELOCITIES + !-------------------------------------------------- + ! NEED TO BE CAREFUL - VELOCITIES SHOULD BE AT THE *LOWER* INTERFACE + ! (THAT IS, FOR K+1), FLUXES ARE ALSO AT THE LOWER INTERFACE (K+1), + ! BUT MIXING RATIOS ARE AT THE MIDPOINTS (K)... + ! + ! NOTE THAT PVEL IS ON INTERFACES, WHEREAS VFALL IS FOR THE CELL + ! AVERAGES (I.E., MIDPOINTS); ASSUME THE FALL VELOCITY APPLICABLE TO THE + ! LOWER INTERFACE (K+1) IS THE SAME AS THAT APPLICABLE FOR THE CELL (V(K)) + + ! LIQUID: + ! The fall velocities assume that droplets have a gamma distribution + ! with effective radii for land and ocean as assessed by Han et al.; + ! see Lawrence and Crutzen (1998) for a derivation. + rho(:ncol,:) = pmid(:ncol,:)/(rair*t(:ncol,:)) + pvliq(:ncol,:) = 0._kind_phys + + ! get effective radius of liquid drop (rel) + call reltab(ncol=ncol, pver=pver, tmelt=tmelt, t=t(:ncol,:), landfrac=landfrac(:ncol), landm=landm(:ncol), icefrac=icefrac(:ncol), snowh=snowh(:ncol), rel=rel(:ncol,:)) + + do k = 1, pver + do i = 1, ncol + if (cloud(i, k) > 0._kind_phys .and. cldliq(i, k) > 0._kind_phys) then + if (rel(i, k) < 40._kind_phys) then + vfall = 2._kind_phys/9._kind_phys*rhoh2o*gravit*rel(i, k)**2/eta*um_to_m ! um^2 -> m^2 + else + vfall = v40 + vslope*(rel(i, k) - r40) ! linear above 40 microns + end if + ! convert the fall speed to pressure units, but do not apply the traditional + ! negative convention for pvel. + pvliq(i, k + 1) = vfall*rho(i, k)*gravit ! m s-1 to Pa s-1 + end if + end do + end do + + ! ICE: + ! The fall velocities are based on data from McFarquhar and Heymsfield + ! or on Stokes terminal velocity for spheres and the effective radius. + pvice(:ncol,:) = 0._kind_phys + + if(stokes) then + ! stokes terminal velocity < 40 microns + + ! get effective radius (rei) + call reitab(ncol=ncol, pver=pver, t=t(:ncol,:), re=rei(:ncol,:)) + + do k = 1, pver + do i = 1, ncol + if (cloud(i, k) > 0._kind_phys .and. cldice(i, k) > 0._kind_phys) then + if (rei(i, k) < 40._kind_phys) then + vfall = 2._kind_phys/9._kind_phys*rhoh2o*gravit*rei(i, k)**2/eta*um_to_m ! microns^2 -> m^2 + vfall = vfall*cldsed_ice_stokes_fac + else + vfall = v40 + vslope*(rei(i, k) - r40) ! linear above 40 microns + end if + + ! convert the fall speed to pressure units, but do not apply the traditional + ! negative convention for pvel. + pvice(i, k + 1) = vfall*rho(i, k)*gravit ! m s-1 to Pa s-1 + end if + end do + end do + + else + ! McFarquhar and Heymsfield > icritc + do k = 1, pver + do i = 1, ncol + if (cloud(i, k) > 0._kind_phys .and. cldice(i, k) > 0._kind_phys) then + + ! compute the in-cloud ice concentration (kg/kg) + icice = cldice(i, k)/cloud(i, k) + + ! compute the ice water content in g/m3 + iciwc = icice*rho(i, k)*1.e3_kind_phys + + ! set the fall velocity (cm/s) to depend on the ice water content in g/m3, + if (iciwc > lbound) then ! need this because of log10 + logiwc = log10(iciwc) + ! Median - + vfall = 128.64_kind_phys + 53.242_kind_phys*logiwc + 5.4795_kind_phys*logiwc**2 + ! Average - + ! vfall = 122.63 + 44.111*logiwc + 4.2144*logiwc**2 + else + vfall = 0._kind_phys + end if + + ! set ice velocity to 1 cm/s if ice mixing ratio < icritc, ramp to value + ! calculated above at 2*icritc + if (icice <= icritc) then + vfall = vice_small + else if (icice < 2*icritc) then + icefac = (icice - icritc)/icritc + vfall = vice_small*(1._kind_phys - icefac) + vfall*icefac + end if + + ! bound the terminal velocity of ice particles at high concentration + vfall = min(100.0_kind_phys, vfall) + + ! convert the fall speed to pressure units, but do not apply the traditional + ! negative convention for pvel. + pvice(i, k + 1) = vfall & + *0.01_kind_phys & ! cm to meters + *rho(i, k)*gravit ! m s-1 to Pa s-1 + end if + end do + end do + endif + + !-------------------------------------------------- + ! APPLY CLOUD PARTICLE GRAVITATIONAL SEDIMENTATION TO CONDENSATE + !-------------------------------------------------- + ! Initialize variables + fxliq(:ncol, :) = 0._kind_phys ! flux at interfaces (liquid) + fxice(:ncol, :) = 0._kind_phys ! flux at interfaces (ice) + liqtend(:ncol, :) = 0._kind_phys ! condensate tend (liquid) + icetend(:ncol, :) = 0._kind_phys ! condensate tend (ice) + wvtend(:ncol, :) = 0._kind_phys ! environmental moistening + htend(:ncol, :) = 0._kind_phys ! evaporative cooling + sfliq(:ncol) = 0._kind_phys ! condensate sedimentation flux out bot of column (liquid) + sfice(:ncol) = 0._kind_phys ! condensate sedimentation flux out bot of column (ice) + + ! Get fluxes at interior points + call getflx(ncol, pver, pverp, pint, cldliq, pvliq, dtime, fxliq, errmsg, errflg) + if(errflg /= 0) then + return + endif + call getflx(ncol, pver, pverp, pint, cldice, pvice, dtime, fxice, errmsg, errflg) + if(errflg /= 0) then + return + endif + + ! Calculate fluxes at boundaries + do i = 1, ncol + fxliq(i, 1) = 0._kind_phys + fxice(i, 1) = 0._kind_phys + ! surface flux by upstream scheme + fxliq(i, pverp) = cldliq(i, pver)*pvliq(i, pverp)*dtime + fxice(i, pverp) = cldice(i, pver)*pvice(i, pverp)*dtime + end do + + ! filter out any negative fluxes from the getflx routine + ! (typical fluxes are of order > 1e-3 when clouds are present) + do k = 2, pver + fxliq(:ncol, k) = max(0._kind_phys, fxliq(:ncol, k)) + fxice(:ncol, k) = max(0._kind_phys, fxice(:ncol, k)) + end do + + ! Limit the flux out of the bottom of each cell to the water content in each phase. + ! Apply mxsedfac to prevent generating very small negative cloud water/ice + ! NOTE, REMOVED CLOUD FACTOR FROM AVAILABLE WATER. ALL CLOUD WATER IS IN CLOUDS. + ! ***Should we include the flux in the top, to allow for thin surface layers? + ! ***Requires simple treatment of cloud overlap, already included below. + do k = 1, pver + do i = 1, ncol + fxliq(i, k + 1) = min(fxliq(i, k + 1), mxsedfac*cldliq(i, k)*pdel(i, k)) + fxice(i, k + 1) = min(fxice(i, k + 1), mxsedfac*cldice(i, k)*pdel(i, k)) + ! fxliq(i,k+1) = min( fxliq(i,k+1), cldliq(i,k) * pdel(i,k) + fxliq(i,k)) + ! fxice(i,k+1) = min( fxice(i,k+1), cldice(i,k) * pdel(i,k) + fxice(i,k)) + ! fxliq(i,k+1) = min( fxliq(i,k+1), cloud(i,k) * cldliq(i,k) * pdel(i,k) ) + ! fxice(i,k+1) = min( fxice(i,k+1), cloud(i,k) * cldice(i,k) * pdel(i,k) ) + end do + end do + + ! Now calculate the tendencies assuming that condensate evaporates when + ! it falls into environment, and does not when it falls into cloud. + ! All flux out of the layer comes from the cloudy part. + ! Assume maximum overlap for stratiform clouds + ! if cloud above < cloud, all water falls into cloud below + ! if cloud above > cloud, water split between cloud and environment + do k = 1, pver + cldab(:ncol) = 0._kind_phys + do i = 1, ncol + ! cloud overlap cloud factor + cldovrl = min(cloud(i, k)/(cldab(i) + .0001_kind_phys), 1._kind_phys) + cldab(i) = cloud(i, k) + ! evaporation into environment cause moistening and cooling + evapliq = fxliq(i, k)*(1._kind_phys - cldovrl)/(dtime*pdel(i, k)) ! into env [kg kg-1 s-1] + evapice = fxice(i, k)*(1._kind_phys - cldovrl)/(dtime*pdel(i, k)) ! into env [kg kg-1 s-1] + wvtend(i, k) = evapliq + evapice ! evaporation into environment [kg kg-1 s-1] + htend(i, k) = -latvap*evapliq - (latvap + latice)*evapice ! evaporation [W kg-1] + ! net flux into cloud changes cloud liquid/ice (all flux is out of cloud) + liqtend(i, k) = (fxliq(i, k)*cldovrl - fxliq(i, k + 1))/(dtime*pdel(i, k)) + icetend(i, k) = (fxice(i, k)*cldovrl - fxice(i, k + 1))/(dtime*pdel(i, k)) + end do + end do + + ! convert flux out the bottom to mass units Pa -> kg/m2/s + sfliq(:ncol) = fxliq(:ncol, pverp)/(dtime*gravit) + sfice(:ncol) = fxice(:ncol, pverp)/(dtime*gravit) + + ! Convert lwe_cloud_ice_sedimentation_rate_at_surface_due_to_microphysics from kg m-2 s-1 to precip units m s-1 + sfice(:ncol) = sfice(:ncol)/1000._kind_phys + + end subroutine cloud_particle_sedimentation_run + + ! Compute fluxes at interior points + subroutine getflx(ncol, pver, pverp, & + xw, phi, vel, deltat, flux, & + errmsg, errflg) + !.....xw1.......xw2.......xw3.......xw4.......xw5.......xw6 + !....psiw1.....psiw2.....psiw3.....psiw4.....psiw5.....psiw6 + !....velw1.....velw2.....velw3.....velw4.....velw5.....velw6 + !.........phi1......phi2.......phi3.....phi4.......phi5....... + + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: pverp + real(kind_phys), intent(in) :: xw(ncol, pverp) + real(kind_phys), intent(in) :: phi(ncol, pver) + real(kind_phys), intent(in) :: vel(ncol, pverp) + real(kind_phys), intent(in) :: deltat + + real(kind_phys), intent(out) :: flux(ncol, pverp) + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + integer :: i, k + real(kind_phys) :: psi(ncol, pverp) + real(kind_phys) :: fdot(ncol, pverp) + real(kind_phys) :: xx(ncol) + real(kind_phys) :: fxdot(ncol) + real(kind_phys) :: fxdd(ncol) + real(kind_phys) :: psistar(ncol) + real(kind_phys) :: xxk(ncol, pver) + + do i = 1, ncol + ! integral of phi + psi(i, 1) = 0._kind_phys + ! fluxes at boundaries + flux(i, 1) = 0._kind_phys + flux(i, pverp) = 0._kind_phys + end do + + ! integral function + do k = 2, pverp + do i = 1, ncol + psi(i, k) = phi(i, k - 1)*(xw(i, k) - xw(i, k - 1)) + psi(i, k - 1) + end do + end do + + ! Calculate the derivatives for the interpolating polynomial + call cfdotmc(ncol, pver, pverp, xw, psi, fdot) + + ! NEW WAY + ! calculate fluxes at interior pts + do k = 2, pver + do i = 1, ncol + xxk(i, k) = xw(i, k) - vel(i, k)*deltat + end do + end do + do k = 2, pver + call cfint2(ncol, pverp, xw, psi, fdot, xxk(1, k), fxdot, fxdd, psistar, errmsg, errflg) + if(errflg /= 0) then + return + endif + do i = 1, ncol + flux(i, k) = (psi(i, k) - psistar(i)) + end do + end do + end subroutine getflx + + ! Vertical level interpolation + subroutine cfint2(ncol, pverp, & + x, f, fdot, xin, fxdot, fxdd, psistar, & + errmsg, errflg) + + integer, intent(in) :: ncol + integer, intent(in) :: pverp + real(kind_phys), intent(in) :: x(ncol, pverp) + real(kind_phys), intent(in) :: f(ncol, pverp) + real(kind_phys), intent(in) :: fdot(ncol, pverp) + real(kind_phys), intent(in) :: xin(ncol) + + real(kind_phys), intent(out) :: fxdot(ncol) + real(kind_phys), intent(out) :: fxdd(ncol) + real(kind_phys), intent(out) :: psistar(ncol) + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + integer :: i, k + integer :: intz(ncol) + real(kind_phys) :: dx + real(kind_phys) :: s + real(kind_phys) :: c2 + real(kind_phys) :: c3 + real(kind_phys) :: xx + real(kind_phys) :: xinf + real(kind_phys) :: psi1, psi2, psi3, psim + real(kind_phys) :: cfint + real(kind_phys) :: cfnew + real(kind_phys) :: xins(ncol) + logical :: found_error + + do i = 1, ncol + xins(i) = medan(x(i, 1), xin(i), x(i, pverp)) + intz(i) = 0 + end do + + ! first find the interval + do k = 1, pverp - 1 + do i = 1, ncol + if ((xins(i) - x(i, k))*(x(i, k + 1) - xins(i)) .ge. 0) then + intz(i) = k + end if + end do + end do + + found_error = .false. + do i = 1, ncol + if (intz(i) .eq. 0._kind_phys) found_error = .true. + end do + if (found_error) then + do i = 1, ncol + if (intz(i) .eq. 0._kind_phys) then + write(errmsg,'(a,i0)') 'cloud_particle_sedimentation/cfint2: interval was not found for column ', i + errflg = 1 + end if + end do + end if + + ! now interpolate + do i = 1, ncol + k = intz(i) + dx = (x(i, k + 1) - x(i, k)) + s = (f(i, k + 1) - f(i, k))/dx + c2 = (3*s - 2*fdot(i, k) - fdot(i, k + 1))/dx + c3 = (fdot(i, k) + fdot(i, k + 1) - 2*s)/dx**2 + xx = (xins(i) - x(i, k)) + fxdot(i) = (3*c3*xx + 2*c2)*xx + fdot(i, k) + fxdd(i) = 6*c3*xx + 2*c2 + cfint = ((c3*xx + c2)*xx + fdot(i, k))*xx + f(i, k) + + ! limit the interpolant + psi1 = f(i, k) + (f(i, k + 1) - f(i, k))*xx/dx + if (k .eq. 1) then + psi2 = f(i, 1) + else + psi2 = f(i, k) + (f(i, k) - f(i, k - 1))*xx/(x(i, k) - x(i, k - 1)) + end if + if (k + 1 .eq. pverp) then + psi3 = f(i, pverp) + else + psi3 = f(i, k + 1) - (f(i, k + 2) - f(i, k + 1))*(dx - xx)/(x(i, k + 2) - x(i, k + 1)) + end if + psim = medan(psi1, psi2, psi3) + cfnew = medan(cfint, psi1, psim) + !if (abs(cfnew - cfint)/(abs(cfnew) + abs(cfint) + 1.e-36_kind_phys) .gt. .03_kind_phys) then + ! CHANGE THIS BACK LATER!!! + ! $ .gt..1) then + ! UNCOMMENT THIS LATER!!! + ! write(iulog,*) ' cfint2 limiting important ', cfint, cfnew + !end if + psistar(i) = cfnew + end do + end subroutine cfint2 + + ! Calculate the derivative for the interpolating polynomial, multi-column version + subroutine cfdotmc(ncol, pver, pverp, x, f, fdot) + ! assumed variable distribution + ! x1.......x2.......x3.......x4.......x5.......x6 1,pverp points + ! f1.......f2.......f3.......f4.......f5.......f6 1,pverp points + ! ...sh1.......sh2......sh3......sh4......sh5.... 1,pver points + ! .........d2.......d3.......d4.......d5......... 2,pver points + ! .........s2.......s3.......s4.......s5......... 2,pver points + ! .............dh2......dh3......dh4............. 2,pver-1 points + ! .............eh2......eh3......eh4............. 2,pver-1 points + ! ..................e3.......e4.................. 3,pver-1 points + ! .................ppl3......ppl4................ 3,pver-1 points + ! .................ppr3......ppr4................ 3,pver-1 points + ! .................t3........t4.................. 3,pver-1 points + ! ................fdot3.....fdot4................ 3,pver-1 points + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: pverp + real(kind_phys), intent(in) :: x(ncol, pverp) + real(kind_phys), intent(in) :: f(ncol, pverp) + + real(kind_phys), intent(out) :: fdot(ncol, pverp) ! derivative at nodes + + integer :: i, k + real(kind_phys) :: a, b, c ! work var + real(kind_phys) :: s(ncol, pverp) ! first divided differences at nodes + real(kind_phys) :: sh(ncol, pverp) ! first divided differences between nodes + real(kind_phys) :: d(ncol, pverp) ! second divided differences at nodes + real(kind_phys) :: dh(ncol, pverp) ! second divided differences between nodes + real(kind_phys) :: e(ncol, pverp) ! third divided differences at nodes + real(kind_phys) :: eh(ncol, pverp) ! third divided differences between nodes + real(kind_phys) :: pp ! p prime + real(kind_phys) :: ppl(ncol, pverp) ! p prime on left + real(kind_phys) :: ppr(ncol, pverp) ! p prime on right + real(kind_phys) :: qpl + real(kind_phys) :: qpr + real(kind_phys) :: ttt + real(kind_phys) :: t + real(kind_phys) :: tmin + real(kind_phys) :: tmax + real(kind_phys) :: delxh(ncol, pverp) + + do k = 1, pver + ! first divided differences between nodes + do i = 1, ncol + delxh(i, k) = (x(i, k + 1) - x(i, k)) + sh(i, k) = (f(i, k + 1) - f(i, k))/delxh(i, k) + end do + + ! first and second divided differences at nodes + if (k .ge. 2) then + do i = 1, ncol + d(i, k) = (sh(i, k) - sh(i, k - 1))/(x(i, k + 1) - x(i, k - 1)) + s(i, k) = minmod(sh(i, k), sh(i, k - 1)) + end do + end if + end do + + ! second and third divided diffs between nodes + do k = 2, pver - 1 + do i = 1, ncol + eh(i, k) = (d(i, k + 1) - d(i, k))/(x(i, k + 2) - x(i, k - 1)) + dh(i, k) = minmod(d(i, k), d(i, k + 1)) + end do + end do + + ! treat the boundaries + do i = 1, ncol + e(i, 2) = eh(i, 2) + e(i, pver) = eh(i, pver - 1) + ! outside level + fdot(i, 1) = sh(i, 1) - d(i, 2)*delxh(i, 1) & + - eh(i, 2)*delxh(i, 1)*(x(i, 1) - x(i, 3)) + fdot(i, 1) = minmod(fdot(i, 1), 3*sh(i, 1)) + fdot(i, pverp) = sh(i, pver) + d(i, pver)*delxh(i, pver) & + + eh(i, pver - 1)*delxh(i, pver)*(x(i, pverp) - x(i, pver - 1)) + fdot(i, pverp) = minmod(fdot(i, pverp), 3*sh(i, pver)) + ! one in from boundary + fdot(i, 2) = sh(i, 1) + d(i, 2)*delxh(i, 1) - eh(i, 2)*delxh(i, 1)*delxh(i, 2) + fdot(i, 2) = minmod(fdot(i, 2), 3*s(i, 2)) + fdot(i, pver) = sh(i, pver) - d(i, pver)*delxh(i, pver) & + - eh(i, pver - 1)*delxh(i, pver)*delxh(i, pver - 1) + fdot(i, pver) = minmod(fdot(i, pver), 3*s(i, pver)) + end do + + do k = 3, pver - 1 + do i = 1, ncol + e(i, k) = minmod(eh(i, k), eh(i, k - 1)) + end do + end do + + do k = 3, pver - 1 + do i = 1, ncol + ! p prime at k-0.5 + ppl(i, k) = sh(i, k - 1) + dh(i, k - 1)*delxh(i, k - 1) + ! p prime at k+0.5 + ppr(i, k) = sh(i, k) - dh(i, k)*delxh(i, k) + + t = minmod(ppl(i, k), ppr(i, k)) + + ! derivate from parabola thru f(i,k-1), f(i,k), and f(i,k+1) + pp = sh(i, k - 1) + d(i, k)*delxh(i, k - 1) + + ! quartic estimate of fdot + fdot(i, k) = pp & + - delxh(i, k - 1)*delxh(i, k) & + *(eh(i, k - 1)*(x(i, k + 2) - x(i, k)) & + + eh(i, k)*(x(i, k) - x(i, k - 2)) & + )/(x(i, k + 2) - x(i, k - 2)) + + ! now limit it + qpl = sh(i, k - 1) & + + delxh(i, k - 1)*minmod(d(i, k - 1) + e(i, k - 1)*(x(i, k) - x(i, k - 2)), & + d(i, k) - e(i, k)*delxh(i, k)) + qpr = sh(i, k) & + + delxh(i, k)*minmod(d(i, k) + e(i, k)*delxh(i, k - 1), & + d(i, k + 1) + e(i, k + 1)*(x(i, k) - x(i, k + 2))) + + fdot(i, k) = medan(fdot(i, k), qpl, qpr) + + ttt = minmod(qpl, qpr) + tmin = min(0._kind_phys, 3*s(i, k), 1.5_kind_phys*t, ttt) + tmax = max(0._kind_phys, 3*s(i, k), 1.5_kind_phys*t, ttt) + + fdot(i, k) = fdot(i, k) + minmod(tmin - fdot(i, k), tmax - fdot(i, k)) + end do + end do + end subroutine cfdotmc + + pure function minmod(a, b) result(res) + real(kind_phys), intent(in) :: a, b + real(kind_phys) :: res + res = 0.5_kind_phys*(sign(1._kind_phys, a) + sign(1._kind_phys, b))*min(abs(a), abs(b)) + end function minmod + + pure function medan(a, b, c) result(res) + real(kind_phys), intent(in) :: a, b, c + real(kind_phys) :: res + res = a + minmod(b - a, c - a) + end function medan + +end module cloud_particle_sedimentation diff --git a/schemes/rasch_kristjansson/cloud_particle_sedimentation.meta b/schemes/rasch_kristjansson/cloud_particle_sedimentation.meta new file mode 100644 index 00000000..4b318031 --- /dev/null +++ b/schemes/rasch_kristjansson/cloud_particle_sedimentation.meta @@ -0,0 +1,257 @@ +[ccpp-table-properties] + name = cloud_particle_sedimentation + type = scheme + dependencies = ../../to_be_ccppized/cloud_optical_properties.F90 + +[ccpp-arg-table] + name = cloud_particle_sedimentation_init + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ cldsed_ice_stokes_fac_in ] + standard_name = tunable_parameter_for_ice_fall_velocity_for_rk_microphysics + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rhoh2o ] + standard_name = fresh_liquid_water_density_at_0c + units = kg m-3 + type = real | kind = kind_phys + dimensions = () + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = cloud_particle_sedimentation_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ pverp ] + standard_name = vertical_interface_dimension + units = count + type = integer + dimensions = () + intent = in +[ dtime ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ tmelt ] + standard_name = freezing_point_of_water + units = K + type = real | kind = kind_phys + dimensions = () + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ latvap ] + standard_name = latent_heat_of_vaporization_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ latice ] + standard_name = latent_heat_of_fusion_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rair ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rhoh2o ] + standard_name = fresh_liquid_water_density_at_0c + units = kg m-3 + type = real | kind = kind_phys + dimensions = () + intent = in +[ icritc ] + standard_name = tunable_parameter_for_autoconversion_of_cold_ice_for_rk_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ pint ] + standard_name = air_pressure_at_interface + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ pmid ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ pdel ] + standard_name = air_pressure_thickness + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cloud ] + standard_name = cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ icefrac ] + standard_name = sea_ice_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ landfrac ] + standard_name = land_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ ocnfrac ] + standard_name = ocean_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ cldliq ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + advected = True + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldice ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + advected = True + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ snowh ] + standard_name = lwe_surface_snow_depth_over_land + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ landm ] + standard_name = smoothed_land_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ pvliq ] + standard_name = vertical_velocity_of_cloud_liquid_water_due_to_sedimentation_tbd + units = Pa s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = out +[ pvice ] + standard_name = vertical_velocity_of_cloud_ice_due_to_sedimentation_tbd + units = Pa s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = out +[ liqtend ] + standard_name = tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = True +[ icetend ] + standard_name = tendency_of_cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = True +[ wvtend ] + standard_name = tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = True +[ htend ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ sfliq ] + standard_name = stratiform_rain_flux_at_surface_due_to_sedimentation + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ sfice ] + standard_name = lwe_cloud_ice_sedimentation_rate_at_surface_due_to_microphysics + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/rasch_kristjansson/cloud_particle_sedimentation_namelist.xml b/schemes/rasch_kristjansson/cloud_particle_sedimentation_namelist.xml new file mode 100644 index 00000000..ad09c34a --- /dev/null +++ b/schemes/rasch_kristjansson/cloud_particle_sedimentation_namelist.xml @@ -0,0 +1,91 @@ + + + + + + + + + real + cldsed + cldsed_nl + tunable_parameter_for_ice_fall_velocity_for_rk_microphysics + 1 + + Factor applied to the ice fall velocity computed from Stokes terminal velocity. + + + 1.0 + + + diff --git a/schemes/rasch_kristjansson/prognostic_cloud_water.F90 b/schemes/rasch_kristjansson/prognostic_cloud_water.F90 new file mode 100644 index 00000000..b2484d6b --- /dev/null +++ b/schemes/rasch_kristjansson/prognostic_cloud_water.F90 @@ -0,0 +1,1189 @@ +! Prognostic cloud water data and methods (cldwat) +! Original authors as marked in subroutines +! CCPP-ized: Haipeng Lin, January 2025 +module prognostic_cloud_water + use ccpp_kinds, only: kind_phys + + implicit none + private + save + + ! public CCPP-compliant subroutines + public :: prognostic_cloud_water_init + public :: prognostic_cloud_water_run + + ! tuning parameters used by prognostic cloud water (RK stratiform) + real(kind_phys), public :: icritc ! threshold for autoconversion of cold ice [kg kg-1] + ! REMOVECAM: icritc does not have to be public after CAM is retired (namelist option; CCPP framework can provide it) + real(kind_phys) :: icritw ! threshold for autoconversion of warm ice [kg kg-1] + real(kind_phys) :: conke ! tunable constant for evaporation of precipitation [kg-0.5 m s-0.5] + real(kind_phys) :: r3lcrit ! critical radius where liq conversion begins [m] + + logical :: do_psrhmin ! set rh in stratosphere poleward of 50 degrees [flag] + real(kind_phys) :: psrhmin ! rh set in stratosphere poleward of 50 degrees [%] + + ! module private variables + real(kind_phys) :: rhonot ! air density at surface [g cm-3] + real(kind_phys) :: rhos ! assumed snow density [g cm-3] + real(kind_phys) :: rhow ! water density [g cm-3] + real(kind_phys) :: rhoi ! ice density [g cm-3] + real(kind_phys) :: esi ! Collection efficiency for ice by snow [1] + real(kind_phys) :: esw ! Collection efficiency for water by snow [1] + real(kind_phys) :: t0 ! Approx. freezing temperature [K] + real(kind_phys) :: cldmin ! Assumed minimum cloud amount [1] + real(kind_phys) :: small ! Small number compared to unity [1] + real(kind_phys) :: c ! Constant for graupel-like snow [cm^(1-d)/s] + + real(kind_phys) :: d ! Constant for graupel-like snow [1] + real(kind_phys) :: thrpd ! 3+d + real(kind_phys) :: gam3pd ! gamma(3+d) + real(kind_phys) :: gam4pd ! gamma(4+d) + + real(kind_phys) :: nos ! Particles snow concentration [cm^-4] + real(kind_phys) :: prhonos ! pi * rhos * nos + + ! cloud microphysics constants + real(kind_phys) :: mcon01 + real(kind_phys) :: mcon02 + real(kind_phys) :: mcon03 + real(kind_phys) :: mcon04 + real(kind_phys) :: mcon05 + real(kind_phys) :: mcon06 + real(kind_phys) :: mcon07 + real(kind_phys) :: mcon08 + + ! Parameters for findmcnew + real(kind_phys) :: capnw ! Warm continental cloud particle number concentration [cm-3] + real(kind_phys) :: capnc ! Cold/oceanic cloud particle number concentration [cm-3] + real(kind_phys) :: capnsi ! Sea ice cloud particle number concentration [cm-3] + real(kind_phys) :: kconst ! Terminal velocity constant (Stokes regime) [1] + real(kind_phys) :: effc ! Autoconv collection efficiency [1] + real(kind_phys) :: alpha ! Ratio of 3rd moment radius to 2nd [1] + real(kind_phys) :: capc ! Autoconversion constant [1] + real(kind_phys) :: convfw ! Fall velocity calculation constant [1] + real(kind_phys) :: cracw ! Rain accreting water constant [1] + real(kind_phys) :: critpr ! Critical precip rate for collection efficiency [mm day-1] + real(kind_phys) :: ciautb ! Ice autoconversion coefficient [s-1] + +contains + + ! Initialize prognostic cloud water module and calculate microphysical constants +!> \section arg_table_prognostic_cloud_water_init Argument Table +!! \htmlinclude arg_table_prognostic_cloud_water_init.html + subroutine prognostic_cloud_water_init( & + amIRoot, iulog, & + tmelt, rhodair, pi, & + icritc_in, icritw_in, & + conke_in, r3lcrit_in, & + do_psrhmin_in, psrhmin_in, & + errmsg, errflg) + + ! Input arguments + logical, intent(in) :: amIRoot ! are we on the MPI root task? + integer, intent(in) :: iulog ! log output unit + + real(kind_phys), intent(in) :: tmelt ! freezing_point_of_water [K] + real(kind_phys), intent(in) :: rhodair ! density_of_dry_air_at_stp [kg m-3] + real(kind_phys), intent(in) :: pi + + real(kind_phys), intent(in) :: icritc_in + real(kind_phys), intent(in) :: icritw_in + real(kind_phys), intent(in) :: conke_in + real(kind_phys), intent(in) :: r3lcrit_in + logical, intent(in) :: do_psrhmin_in + real(kind_phys), intent(in) :: psrhmin_in + + ! Output arguments + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + errmsg = '' + errflg = 0 + + ! First populate tuning parameters in-module + icritc = icritc_in + icritw = icritw_in + conke = conke_in + r3lcrit = r3lcrit_in + do_psrhmin = do_psrhmin_in + psrhmin = psrhmin_in + + if(amIRoot) then + write(iulog,*) 'tuning parameters prognostic_cloud_water_init: icritw ',icritw,' icritc ',icritc,' conke ',conke,' r3lcrit ',r3lcrit + write(iulog,*) 'prognostic_cloud_water_init: do_psrhmin = ', do_psrhmin + endif + + !-------------------------------------------------- + ! Initialize constants used for prognostic condensate (inimc) + ! Original author: P. Rasch, April 1997 + !-------------------------------------------------- + + rhonot = rhodair/1000.0_kind_phys ! convert from kg m-3 to g cm-3 + + ! assumed densities of snow, water, ice [g cm-3] + rhos = 0.1_kind_phys + rhow = 1._kind_phys + rhoi = 1._kind_phys + + esi = 1._kind_phys + esw = 0.1_kind_phys + t0 = tmelt + cldmin = 0.02_kind_phys + small = 1.e-22_kind_phys + c = 152.93_kind_phys + d = 0.25_kind_phys ! constant for graupel like snow + ! values other than 0.25 are not supported. + if(d == 0.25_kind_phys) then + gam3pd = 2.549256966718531_kind_phys + gam4pd = 8.285085141835282_kind_phys + + ! on cray machines this can be calculated using gamma function: + ! call gamma(3._kind_phys+d, signgam, gam3pd) + ! gam3pd = sign(exp(gam3pd),signgam) + ! call gamma(4._kind_phys+d, signgam, gam4pd) + ! gam4pd = sign(exp(gam4pd),signgam) + ! write(iulog,*) ' d, gamma(3+d), gamma(4+d) =', gam3pd, gam4pd + else + errflg = 1 + errmsg = 'prognostic_cloud_water_init: cannot use d /= 0.25' + endif + + nos = 3.e-2_kind_phys + + ! note: pi was originally calculated here as 4._kind_phys * atan(1.0_kind_phys) + ! changed in ccpp-ization to use physconst value. + prhonos = pi*rhos*nos + thrpd = 3._kind_phys + d + + mcon01 = pi*nos*c*gam3pd/4._kind_phys + mcon02 = 1._kind_phys/(c*gam4pd*sqrt(rhonot)/(6*prhonos**(d/4._kind_phys))) + mcon03 = -(0.5_kind_phys+d/4._kind_phys) + mcon04 = 4._kind_phys/(4._kind_phys+d) + mcon05 = (3+d)/(4+d) + mcon06 = (3+d)/4._kind_phys + mcon07 = mcon01*sqrt(rhonot)*mcon02**mcon05/prhonos**mcon06 + mcon08 = -0.5_kind_phys/(4._kind_phys+d) + + ! Initialize parameters used by findmcnew + ! Cloud particle densities [cm-3] for ... + capnw = 400._kind_phys ! ...warm continental + capnc = 150._kind_phys ! ...cold and oceanic + capnsi = 75._kind_phys ! ...sea ice + + kconst = 1.18e6_kind_phys + + ! Autoconv collection efficiencies. + ! Tripoli and Cotton (default) = 0.55 + ! Boucher 96 = 1. + ! Baker 93 = 0.55*0.05 + ! Turn off = 0 + effc = 0.55_kind_phys + + alpha = 1.1_kind_phys ** 4 + ! constant for autoconversion + capc = pi**(-.333_kind_phys)*kconst*effc *(0.75_kind_phys)**(1.333_kind_phys)*alpha + + ! critical precip rate at which we assume the collector drops can change the + ! drop size enough to enhance the auto-conversion process (mm/day) + critpr = 0.5_kind_phys + convfw = 1.94_kind_phys*2.13_kind_phys*sqrt(rhow*1000._kind_phys*9.81_kind_phys*2.7e-4_kind_phys) + + ! liquid microphysics - rain accreting water constant + ! Tripoli and Cotton = default + ! Beheng = 6. + cracw = .884_kind_phys*sqrt(9.81_kind_phys/(rhow*1000._kind_phys*2.7e-4_kind_phys)) ! tripoli and cotton + + ! ice microphysics autoconversion coefficient + ciautb = 5.e-4_kind_phys + + if(amIRoot) then + write(iulog,*) 'tuning parameters prognostic_cloud_water_init: capnw ',capnw,' capnc ',capnc,' capnsi ',capnsi,' kconst ',kconst + write(iulog,*) 'tuning parameters prognostic_cloud_water_init: effc ',effc,' alpha ',alpha,' capc ',capc + write(iulog,*) 'tuning parameters prognostic_cloud_water_init: critpr ',critpr,' convfw ',convfw,' cracw ',cracw,' ciautb ',ciautb + endif + + end subroutine prognostic_cloud_water_init + + ! Calculate prognostic condensate. + ! Cloud water parameterization, returns tendencies to water vapor, temperature, + ! and cloud water variables. + ! + ! Rasch, P. J., and J. E. Kristjánsson, 1998: A Comparison of the CCM3 Model Climate Using Diagnosed and Predicted Condensate Parameterizations. J. Climate, 11, 1587–1614 + ! https://doi.org/10.1175/1520-0442(1998)011<1587:ACOTCM>2.0.CO;2 + ! + ! with modifications to improve the method of determining condensation/evaporation + ! Zhang, M., W. Lin, C. Bretherton, J. Hack, and P. J. Rasch, 2003, A modified formulation of fractional stratiform condensation rate in the NCAR Community Atmospheric Model (CAM2), J. Geophys. Res., 108(D1), 4035 + ! https://doi.org/10.1029/2002JD002523 + ! + ! Original authors: M. Zhang, W. Lin, P. Rasch and J.E. Kristjansson + ! B. A. Boville (latent heat of fusion) +!> \section arg_table_prognostic_cloud_water_run Argument Table +!! \htmlinclude arg_table_prognostic_cloud_water_run.html + subroutine prognostic_cloud_water_run( & + ncol, pver, top_lev, deltat, & + iulog, & + pi, gravit, rh2o, epsilo, latvap, latice, cpair, & + dlat, & + pmid, pdel, & + zi, & + troplev, & + ttend, tn, & + qtend, qn, & + ltend, & + cldice, cldliq, & + omega, & + cldn, & + fice, fsnow, & + rhdfda, rhu00, & + landm, seaicef, snowh, & + qme, & + prodprec, prodsnow, & + evapprec, evapsnow, & + evapheat, prfzheat, meltheat, & + precip, snowab, & + ice2pr, liq2pr, liq2snow, & + lsflxprc, & + lsflxsnw, & + pracwo, psacwo, psacio, & + fwaut, fsaut, fracw, fsacw, fsaci, & + errmsg, errflg) + ! Dependencies to-be-ccppized. + use wv_saturation, only: qsat, estblf, svp_to_qsat, findsp_vc + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: top_lev ! vertical_layer_index_of_troposphere_cloud_physics_top [index] + real(kind_phys), intent(in) :: deltat ! timestep [s] + integer, intent(in) :: iulog ! log output unit [1] + real(kind_phys), intent(in) :: pi ! pi_constant [1] + real(kind_phys), intent(in) :: gravit ! gravitational acceleration [m s-2] + real(kind_phys), intent(in) :: rh2o ! gas_constant_of_water_vapor [J K-1 kg-1] + real(kind_phys), intent(in) :: epsilo ! ratio_of_water_vapor_to_dry_air_molecular_weights [1] + real(kind_phys), intent(in) :: latvap ! latent_heat_of_vaporization_of_water_at_0c [J kg-1] + real(kind_phys), intent(in) :: latice ! latent_heat_of_fusion_of_water_at_0c [J kg-1] + real(kind_phys), intent(in) :: cpair ! specific_heat_of_dry_air_at_constant_pressure [J K-1 kg-1] + real(kind_phys), intent(in) :: dlat(:) ! latitude_degrees_north [degrees] + real(kind_phys), intent(in) :: pmid(:,:) ! air_pressure [Pa] + real(kind_phys), intent(in) :: pdel(:,:) ! air_pressure_thickness [Pa] + real(kind_phys), intent(in) :: zi(:,:) ! geopotential_height_wrt_surface_at_interface [m] + integer, intent(in) :: tropLev(:) ! tropopause_vertical_layer_index [index] + + real(kind_phys), intent(in) :: ttend(:,:) ! Temperature tendency [K s-1] -- from non-micro/macrophysics + real(kind_phys), intent(in) :: tn(:,:) ! New Temperature [K] + real(kind_phys), intent(in) :: qtend(:,:) ! Water vapor tendency [kg kg-1 s-1] -- from non-micro/macrophysics + real(kind_phys), intent(in) :: qn(:,:) ! New Water vapor mixing ratio [kg kg-1] + real(kind_phys), intent(in) :: ltend(:,:) ! Cloud liquid water tendency [kg kg-1 s-1] -- from non-micro/macrophysics + real(kind_phys), intent(in) :: cldice(:,:) ! adv: cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldliq(:,:) ! adv: cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + + real(kind_phys), intent(in) :: omega(:,:) ! lagrangian_tendency_of_air_pressure [Pa s-1] + real(kind_phys), intent(in) :: cldn(:,:) ! New Cloud fraction [fraction] + real(kind_phys), intent(in) :: fice(:,:) ! Ice fraction within cwat [fraction] + real(kind_phys), intent(in) :: fsnow(:,:) ! Fraction of rain that freezes to snow [fraction] + + real(kind_phys), intent(in) :: rhdfda(:,:) ! dG(a)/da, rh=G(a), when rh>u00 [??] + real(kind_phys), intent(in) :: rhu00(:,:) ! Rhlim for cloud [percent] + real(kind_phys), intent(in) :: landm(:) ! smoothed_land_area_fraction [fraction] + real(kind_phys), intent(in) :: seaicef(:) ! Sea ice fraction [1 (fraction)] + real(kind_phys), intent(in) :: snowh(:) ! Snow depth over land, water equivalent [m] + + ! Output arguments + real(kind_phys), intent(out) :: qme(:,:) ! Rate of condensation-evaporation of condensate (net_condensation_rate_due_to_microphysics) [kg kg-1 s-1] + real(kind_phys), intent(out) :: prodprec(:,:) ! Conversion rate of condensate to precip (precipitation_production_due_to_microphysics) [kg kg-1 s-1] + real(kind_phys), intent(out) :: prodsnow(:,:) ! Snow production rate (ignored in RK?) [kg kg-1 s-1] + real(kind_phys), intent(out) :: evapprec(:,:) ! Falling precipitation evaporation rate (precipitation_evaporation_due_to_microphysics) [kg kg-1 s-1] -- & combined to apply q(wv) tendency + real(kind_phys), intent(out) :: evapsnow(:,:) ! Falling snow evaporation rate [kg kg-1 s-1] + real(kind_phys), intent(out) :: evapheat(:,:) ! heating rate due to evaporation of precipitation [J kg-1 s-1] + real(kind_phys), intent(out) :: prfzheat(:,:) ! heating rate due to freezing of precipitation [J kg-1 s-1] + real(kind_phys), intent(out) :: meltheat(:,:) ! heating rate due to snow melt [J kg-1 s-1] + ! note -- these are precip units for atmosphere-surface exchange. + real(kind_phys), intent(out) :: precip(:) ! Precipitation rate (lwe_stratiform_precipitation_rate_at_surface) [m s-1] + real(kind_phys), intent(out) :: snowab(:) ! Snow rate (lwe_snow_precipitation_rate_at_surface_due_to_microphysics) [m s-1] + ! / note + real(kind_phys), intent(out) :: ice2pr(:,:) ! Conversion rate of ice to precip [kg kg-1 s-1] -- apply q(cldice) tendency + real(kind_phys), intent(out) :: liq2pr(:,:) ! Conversion rate of liquid to precip [kg kg-1 s-1] -- apply q(cldliq) tendency + real(kind_phys), intent(out) :: liq2snow(:,:) ! Conversion rate of liquid to snow [kg kg-1 s-1] -- discard?? + real(kind_phys), intent(out) :: lsflxprc(:,:) ! Grid-box mean, flux large scale cloud rain+snow at interfaces (stratiform_rain_and_snow_flux_at_interface) [kg m-2 s-1] + real(kind_phys), intent(out) :: lsflxsnw(:,:) ! Grid-box mean, flux large scale cloud snow at interfaces (stratiform_snow_flux_at_interface) [kg m-2 s-1] + real(kind_phys), intent(out) :: pracwo(:,:) ! accretion of cloud water by rain [s-1] + real(kind_phys), intent(out) :: psacwo(:,:) ! accretion of cloud water by snow [s-1] + real(kind_phys), intent(out) :: psacio(:,:) ! accretion of cloud ice by snow [s-1] + real(kind_phys), intent(out) :: fwaut(ncol,pver) ! Relative importance of liquid autoconversion [fraction] + real(kind_phys), intent(out) :: fsaut(ncol,pver) ! Relative importance of ice autoconversion [fraction] + real(kind_phys), intent(out) :: fracw(ncol,pver) ! Relative importance of liquid collection by rain [fraction] + real(kind_phys), intent(out) :: fsacw(ncol,pver) ! Relative importance of liquid collection by snow [fraction] + real(kind_phys), intent(out) :: fsaci(ncol,pver) ! Relative importance of ice collection by snow [fraction] + + ! Output arguments + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local variables + integer :: i, k, l ! Iteration index [1] + integer :: iter ! # of iterations for precipitation calculation [1] + logical :: error_found ! Flag for error detection [flag] + + ! Total cloud water (from cldice+cldliq) + real(kind_phys) :: cwat(ncol,pver) ! Cloud water mixing ratio [kg kg-1] + + ! Precipitation and conversion rates + real(kind_phys) :: nice2pr ! Rate of conversion from ice to snow [kg kg-1 s-1] + real(kind_phys) :: nliq2pr ! Rate of conversion from liquid to precipitation [kg kg-1 s-1] + real(kind_phys) :: nliq2snow ! Rate of conversion from liquid to snow [kg kg-1 s-1] + real(kind_phys) :: precab(ncol) ! Rate of precipitation entering layer [kg m-2 s-1] + real(kind_phys) :: prprov(ncol) ! Provisional precipitation at bottom of layer [kg m-2 s-1] + + ! Cloud properties + real(kind_phys) :: cldm(ncol) ! Mean cloud fraction over timestep [1] + real(kind_phys) :: cldmax(ncol) ! Maximum cloud fraction above current level [1] + real(kind_phys) :: icwc(ncol) ! In-cloud water content [kg kg-1] + real(kind_phys) :: cwm(ncol) ! Cloud water mixing ratio at midpoint of timestep [kg kg-1] + real(kind_phys) :: cwn(ncol) ! Cloud water mixing ratio at end of timestep [kg kg-1] + real(kind_phys) :: coef(ncol) ! Conversion timescale for condensate to rain [s-1] + + ! Thermodynamic variables + real(kind_phys) :: t(ncol,pver) ! Temperature before timestep [K] + real(kind_phys) :: tsp(ncol,pver) ! Saturation point temperature [K] + real(kind_phys) :: es(ncol) ! Saturation vapor pressure [Pa] + real(kind_phys) :: q(ncol,pver) ! Water vapor mixing ratio [kg kg-1] + real(kind_phys) :: qs(ncol) ! Saturation specific humidity [kg kg-1] + real(kind_phys) :: qsp(ncol,pver) ! Saturation point mixing ratio [kg kg-1] + real(kind_phys) :: relhum(ncol) ! Relative humidity [1] + real(kind_phys) :: relhum1(ncol) ! Updated relative humidity [1] + real(kind_phys) :: rhu_adj(ncol,pver) ! Adjusted relative humidity limit for strat dehydration [1] + + ! Cloud-Evaporation scheme variables + real(kind_phys) :: calpha(ncol) ! Alpha term in C-E formulation [1] + real(kind_phys) :: cbeta(ncol) ! Beta term in C-E formulation [1] + real(kind_phys) :: cbetah(ncol) ! Beta-hat at saturation portion [1] + real(kind_phys) :: cgamma(ncol) ! Gamma term in C-E formulation [1] + real(kind_phys) :: cgamah(ncol) ! Gamma-hat at saturation portion [1] + real(kind_phys) :: rcgama(ncol) ! Ratio of gamma to gamma-hat [1] + real(kind_phys) :: csigma(ncol) ! Sigma term in C-E formulation [1] + real(kind_phys) :: qmeres(ncol) ! Residual condensation after C-E and evapprec [kg kg-1 s-1] + real(kind_phys) :: qmec1(ncol) ! Cloud condensation coefficient 1 C-E formulation [1] + real(kind_phys) :: qmec2(ncol) ! Cloud condensation coefficient 2 C-E formulation [1] + real(kind_phys) :: qmec3(ncol) ! Cloud condensation coefficient 3 C-E formulation [1] + real(kind_phys) :: qmec4(ncol) ! Cloud condensation coefficient 4 C-E formulation [1] + + ! Diagnostic arrays for cloud water budget + ! Hardcoded 2 here = number of iterations (iter below) + real(kind_phys) :: rcwn(ncol,2,pver) ! Cloud water ratio evolution [kg kg-1] + real(kind_phys) :: rliq(ncol,2,pver) ! Liquid water ratio evolution [kg kg-1] + real(kind_phys) :: rice(ncol,2,pver) ! Ice ratio evolution [kg kg-1] + + ! Constants and parameters + real(kind_phys), parameter :: omsm = 0.99999_kind_phys ! unity minus small number for rounding + real(kind_phys) :: mincld ! Minimum cloud fraction [1] + real(kind_phys) :: cpohl ! Ratio of specific heat to latent heat [K-1] + real(kind_phys) :: hlocp ! Ratio of latent heat to specific heat [K] + real(kind_phys) :: clrh2o ! Ratio of latent heat to water vapor gas constant [K] + real(kind_phys) :: dto2 ! Half timestep [s] + + ! Work variables + real(kind_phys) :: denom ! Denominator work variable [1] + real(kind_phys) :: dqsdt ! Change in saturation specific humidity with temperature [kg kg-1 K-1] + real(kind_phys) :: gamma(ncol) ! Temperature derivative of saturation specific humidity [kg kg-1 K-1] + real(kind_phys) :: qtl(ncol) ! Saturation tendency [kg kg-1 s-1] + real(kind_phys) :: qtmp ! Temporary mixing ratio [kg kg-1] + real(kind_phys) :: ttmp ! Temporary temperature [K] + real(kind_phys) :: qsn ! Updated saturation specific humidity [kg kg-1] + real(kind_phys) :: esn ! Updated saturation vapor pressure [Pa] + real(kind_phys) :: prtmp ! Temporary precipitation rate [kg m-2 s-1] + real(kind_phys) :: ctmp ! Temporary condensation rate [kg kg-1 s-1] + real(kind_phys) :: cdt ! Timestep factor [1] + real(kind_phys) :: wtthick ! Layer thickness weight [1] + real(kind_phys) :: pol ! Production over loss ratio [1] + + errmsg = '' + errflg = 0 + error_found = .false. + + clrh2o = latvap/rh2o + cpohl = cpair/latvap + hlocp = latvap/cpair + dto2 = 0.5_kind_phys * deltat + +#ifdef PERGRO + mincld = 1.e-4_kind_phys + iter = 1 ! number of times to iterate the precipitation calculation +#else + mincld = 1.e-4_kind_phys + iter = 2 +#endif + + ! initialize total cloud water [kg kg-1] + cwat(:ncol,:pver) = cldice(:ncol,:pver) + cldliq(:ncol,:pver) + + ! initialize single level and multi level fields + do i = 1, ncol + precip(i) = 0.0_kind_phys + precab(i) = 0.0_kind_phys + snowab(i) = 0.0_kind_phys + cldmax(i) = 0.0_kind_phys + end do + + do k = 1, pver + do i = 1, ncol + q(i,k) = qn(i,k) + t(i,k) = tn(i,k) + ! q(i,k)=qn(i,k)-qtend(i,k)*deltat + ! t(i,k)=tn(i,k)-ttend(i,k)*deltat + end do + end do + + qme (:ncol,:) = 0._kind_phys + evapprec(:ncol,:) = 0._kind_phys + prodprec(:ncol,:) = 0._kind_phys + evapsnow(:ncol,:) = 0._kind_phys + prodsnow(:ncol,:) = 0._kind_phys + evapheat(:ncol,:) = 0._kind_phys + meltheat(:ncol,:) = 0._kind_phys + prfzheat(:ncol,:) = 0._kind_phys + ice2pr(:ncol,:) = 0._kind_phys + liq2pr(:ncol,:) = 0._kind_phys + liq2snow(:ncol,:) = 0._kind_phys + fwaut(:ncol,:) = 0._kind_phys + fsaut(:ncol,:) = 0._kind_phys + fracw(:ncol,:) = 0._kind_phys + fsacw(:ncol,:) = 0._kind_phys + fsaci(:ncol,:) = 0._kind_phys + lsflxprc(:ncol,:) = 0._kind_phys + lsflxsnw(:ncol,:) = 0._kind_phys + pracwo(:ncol,:) = 0._kind_phys + psacwo(:ncol,:) = 0._kind_phys + psacio(:ncol,:) = 0._kind_phys + + ! reset diagnostic arrays + rcwn(:,:,:) = 0._kind_phys + rliq(:,:,:) = 0._kind_phys + rice(:,:,:) = 0._kind_phys + + ! find the wet bulb temp and saturation value + ! for the provisional t and q without condensation + mp_level_loop: do k = top_lev, pver + ! "True" means that ice will be taken into account. + call findsp_vc(qn(:ncol,k), tn(:ncol,k), pmid(:ncol,k), .true., & + tsp(:ncol,k), qsp(:ncol,k)) + + call qsat(t(1:ncol,k), pmid(1:ncol,k), es(1:ncol), qs(1:ncol), ncol, gam=gamma(1:ncol)) + + do i = 1, ncol + relhum(i) = q(i,k)/qs(i) + cldm(i) = max(cldn(i,k),mincld) + + ! the max cloud fraction above this level + cldmax(i) = max(cldmax(i), cldm(i)) + + ! define the coefficients for C - E calculation + calpha(i) = 1.0_kind_phys/qs(i) + cbeta (i) = q(i,k)/qs(i)**2*gamma(i)*cpohl + cbetah(i) = 1.0_kind_phys/qs(i)*gamma(i)*cpohl + cgamma(i) = calpha(i)+latvap*cbeta(i)/cpair + cgamah(i) = calpha(i)+latvap*cbetah(i)/cpair + rcgama(i) = cgamma(i)/cgamah(i) + + if(cldm(i) > mincld) then + icwc(i) = max(0._kind_phys,cwat(i,k)/cldm(i)) + else + icwc(i) = 0.0_kind_phys + endif + ! PJR the above logic give zero icwc with nonzero cwat, dont like it! + ! PJR generates problems with csigma + ! PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds + ! icwc(i) = max(1.e-8_kind_phys,cwat(i,k)/cldm(i)) + + ! initial guess of evaporation, will be updated within iteration + evapprec(i,k) = conke*(1._kind_phys - cldm(i))*sqrt(precab(i)) & + *(1._kind_phys - min(relhum(i),1._kind_phys)) + + ! zero qmeres before iteration for each level + qmeres(i) = 0.0_kind_phys + end do + + do i = 1, ncol + ! calculate the cooling due to a phase change of the rainwater from above + if (t(i,k) >= t0) then + meltheat(i,k) = -latice * snowab(i) * gravit/pdel(i,k) + snowab(i) = 0._kind_phys + else + meltheat(i,k) = 0._kind_phys + endif + end do + + ! calculate qme and formation of precip. + ! + ! The cloud microphysics is highly nonlinear and coupled with qme + ! Both rain processes and qme are calculated iteratively. + qme_iter_loop: do l = 1, iter + qme_update_loop: do i = 1, ncol + ! calculation of qme has 4 scenarios + call relhum_min_adj(ncol, pver, tropLev, dlat, rhu00, rhu_adj) + + if(relhum(i) > rhu_adj(i,k)) then + ! 1. whole grid saturation + if(relhum(i) >= 0.999_kind_phys .or. cldm(i) >= 0.999_kind_phys) then + qme(i,k) = (calpha(i)*qtend(i,k)-cbetah(i)*ttend(i,k))/cgamah(i) + ! 2. fractional saturation + else + if (rhdfda(i,k) .eq. 0._kind_phys .and. icwc(i) .eq. 0._kind_phys) then + write(iulog,*) 'prognostic_cloud_water: empty rh cloud @ ', i, k + write(iulog,*) 'relhum, iter ', relhum(i), l, rhu_adj(i,k), cldm(i), cldn(i,k) + + errflg = 1 + errmsg = 'prognostic_cloud_water: empty RH cloud' + return + endif + + csigma(i) = 1.0_kind_phys/(rhdfda(i,k)+cgamma(i)*icwc(i)) + qmec1(i) = (1.0_kind_phys-cldm(i))*csigma(i)*rhdfda(i,k) + qmec2(i) = cldm(i)*calpha(i)/cgamah(i)+(1.0_kind_phys-rcgama(i)*cldm(i))* & + csigma(i)*calpha(i)*icwc(i) + qmec3(i) = cldm(i)*cbetah(i)/cgamah(i) + & + (cbeta(i)-rcgama(i)*cldm(i)*cbetah(i))*csigma(i)*icwc(i) + qmec4(i) = csigma(i)*cgamma(i)*icwc(i) + + ! Q = C-E = -C1*Al + C2*Aq - C3*At + C4*Er + qme(i,k) = -qmec1(i)*ltend(i,k) + qmec2(i)*qtend(i,k) & + -qmec3(i)*ttend(i,k) + qmec4(i)*evapprec(i,k) + endif + ! 3. when rh < rhu00, evaporate existing cloud water + else if(cwat(i,k) > 0.0_kind_phys) then + ! liquid water should be evaporated but not to exceed + ! saturation point. if qn > qsp, not to evaporate cwat + qme(i,k) = -min(max(0._kind_phys,qsp(i,k)-qn(i,k)),cwat(i,k))/deltat + ! 4. no condensation nor evaporation + else + qme(i,k) = 0.0_kind_phys + endif + end do qme_update_loop + + ! Because of the finite time step, place a bound here not to exceed wet bulb point + ! and not to evaporate more than available water + qme_bound_loop: do i = 1, ncol + qtmp = qn(i,k) - qme(i,k)*deltat + + ! possibilities to have qtmp > qsp + ! + ! 1. if qn > qs(tn), it condenses; + ! if after applying qme, qtmp > qsp, more condensation is applied. + ! + ! 2. if qn < qs, evaporation should not exceed qsp, + if(qtmp > qsp(i,k)) then + qme(i,k) = qme(i,k) + (qtmp-qsp(i,k))/deltat + endif + + ! if net evaporation, it should not exceed available cwat + if(qme(i,k) < -cwat(i,k)/deltat) then + qme(i,k) = -cwat(i,k)/deltat + endif + + ! addition of residual condensation from previous step of iteration + qme(i,k) = qme(i,k) + qmeres(i) + + ! limit qme for roundoff errors (multiply by slightly less than unity) + qme(i,k) = qme(i,k) * omsm + end do qme_bound_loop + + do i = 1, ncol + ! as a safe limit, condensation should not reduce grid mean rh below rhu00 + if(qme(i,k) > 0.0_kind_phys .and. relhum(i) > rhu_adj(i,k)) then + qme(i,k) = min(qme(i,k), (qn(i,k)-qs(i)*rhu_adj(i,k))/deltat) + endif + + ! initial guess for cwm (mean cloud water over time step) if 1st iteration + if(l < 2) then + cwm(i) = max(cwat(i,k)+qme(i,k)*dto2, 0._kind_phys) + endif + enddo + + ! provisional precipitation falling through model layer + prprov_update_loop: do i = 1, ncol + ! prprov(i) = precab(i) + prodprec(i,k)*pdel(i,k)/gravit + ! rain produced in this layer not too effective in collection process + wtthick = max(0._kind_phys,min(0.5_kind_phys,((zi(i,k)-zi(i,k+1))/1000._kind_phys)**2)) + prprov(i) = precab(i) + wtthick*prodprec(i,k)*pdel(i,k)/gravit + end do prprov_update_loop + + ! calculate conversion of condensate to precipitation by cloud microphysics + call findmcnew( & + ncol = ncol, & + pver = pver, & + pi = pi, & + k = k, & + precab = prprov(:ncol), & + snowab = snowab(:ncol), & + t = t(:ncol,:), & + p = pmid(:ncol,:), & + cwm = cwm(:ncol), & + cldm = cldm(:ncol), & + cldmax = cldmax(:ncol), & + fice = fice(:ncol,k), & + landm = landm(:ncol), & + seaicef = seaicef(:ncol), & + snowh = snowh(:ncol), & + ! below output + coef = coef(:ncol), & + fwaut = fwaut(:ncol,k), & + fsaut = fsaut(:ncol,k), & + fracw = fracw(:ncol,k), & + fsacw = fsacw(:ncol,k), & + fsaci = fsaci(:ncol,k), & + pracwo = pracwo(:ncol,k),& + psacwo = psacwo(:ncol,k),& + psacio = psacio(:ncol,k)) + + ! calculate the precip rate + error_found = .false. + precip_update_loop: do i = 1, ncol + if (cldm(i) > 0) then + ! first predict the cloud water + cdt = coef(i)*deltat + if(cdt > 0.01_kind_phys) then + pol = qme(i,k)/coef(i) ! production over loss + cwn(i) = max(0._kind_phys,(cwat(i,k)-pol)*exp(-cdt)+ pol) + else + cwn(i) = max(0._kind_phys,(cwat(i,k) + qme(i,k)*deltat)/(1+cdt)) + endif + + ! now back out the tendency of net rain production + prodprec(i,k) = max(0._kind_phys,qme(i,k)-(cwn(i)-cwat(i,k))/deltat) + else + prodprec(i,k) = 0.0_kind_phys + cwn(i) = 0._kind_phys + endif + + ! provisional calculation of conversion terms + ice2pr(i,k) = prodprec(i,k)*(fsaut(i,k)+fsaci(i,k)) + liq2pr(i,k) = prodprec(i,k)*(fwaut(i,k)+fsacw(i,k)+fracw(i,k)) + + ! revision suggested by Jim McCaa + ! it controls the amount of snow hitting the sfc + ! by forcing a lot of conversion of cloud liquid to snow phase + ! it might be better done later by an explicit representation of + ! rain accreting ice (and freezing), or by an explicit freezing of raindrops + ! + ! old: + ! liq2snow(i,k) = prodprec(i,k)*fsacw(i,k) + liq2snow(i,k) = max(prodprec(i,k)*fsacw(i,k), fsnow(i,k)*liq2pr(i,k)) + + ! bounds + nice2pr = min(ice2pr(i,k),(cwat(i,k)+qme(i,k)*deltat)*fice(i,k)/deltat) + nliq2pr = min(liq2pr(i,k),(cwat(i,k)+qme(i,k)*deltat)*(1._kind_phys-fice(i,k))/deltat) + + if (liq2pr(i,k) .ne. 0._kind_phys) then + nliq2snow = liq2snow(i,k)*nliq2pr/liq2pr(i,k) ! correction + else + nliq2snow = liq2snow(i,k) + endif + + ! avoid roundoff problems generating negatives + nliq2snow = nliq2snow*omsm + nliq2pr = nliq2pr*omsm + nice2pr = nice2pr*omsm + + ! final estimates of conversion to precip and snow + prodprec(i,k) = (nliq2pr + nice2pr) + prodsnow(i,k) = (nice2pr + nliq2snow) + + ! compute diagnostics and sanity checks + rcwn(i,l,k) = cwat(i,k) + (qme(i,k) - prodprec(i,k))*deltat + rliq(i,l,k) = (cwat(i,k) + qme(i,k)*deltat)*(1._kind_phys-fice(i,k)) - nliq2pr*deltat + rice(i,l,k) = (cwat(i,k) + qme(i,k)*deltat) * fice(i,k) - nice2pr*deltat + + ! Sanity checks + if(abs(rcwn(i,l,k)) < 1.e-300_kind_phys) rcwn(i,l,k) = 0._kind_phys + if(abs(rliq(i,l,k)) < 1.e-300_kind_phys) rliq(i,l,k) = 0._kind_phys + if(abs(rice(i,l,k)) < 1.e-300_kind_phys) rice(i,l,k) = 0._kind_phys + if(rcwn(i,l,k) < 0._kind_phys) error_found = .true. + if(rliq(i,l,k) < 0._kind_phys) error_found = .true. + if(rice(i,l,k) < 0._kind_phys) error_found = .true. + if(error_found) then + if(rcwn(i,l,k) < 0._kind_phys) then + write(iulog,*) 'prognostic_cloud_water: prob with neg rcwn1 ', rcwn(i,l,k), cwn(i) + write(iulog,*) ' cwat, qme*deltat, prodprec*deltat ', & + cwat(i,k), qme(i,k)*deltat, & + prodprec(i,k)*deltat, & + (qme(i,k)-prodprec(i,k))*deltat + + errflg = 1 + errmsg = 'prognostic_cloud_water: negative rcwn1' + return + endif + + if (rliq(i,l,k) < 0._kind_phys) then + write(iulog,*) 'prognostic_cloud_water: prob with neg rliq1 ', rliq(i,l,k) + errflg = 1 + errmsg = 'prognostic_cloud_water: negative rliq1' + return + endif + + if (rice(i,l,k) < 0._kind_phys) then + write(iulog,*) 'prognostic_cloud_water: prob with neg rice ', rice(i,l,k) + errflg = 1 + errmsg = 'prognostic_cloud_water: negative rice' + return + endif + endif + + ! Final version of condensate to precip terms + liq2pr(i,k) = nliq2pr + liq2snow(i,k) = nliq2snow + ice2pr(i,k) = nice2pr + cwn(i) = rcwn(i,l,k) + + ! update any remaining provisional values + cwm(i) = (cwn(i) + cwat(i,k))*0.5_kind_phys + + ! update in cloud water + if(cldm(i) > mincld) then + icwc(i) = cwm(i)/cldm(i) + else + icwc(i) = 0.0_kind_phys + endif + ! PJR the above logic give zero icwc with nonzero cwat, dont like it! + ! PJR generates problems with csigma + ! PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds + ! icwc(i) = max(1.e-8_kind_phys,cwm(i)/cldm(i)) + end do precip_update_loop + + ! calculate provisional value of cloud water for + ! evaporation of precipitate (evapprec) calculation + qtl_update_loop: do i = 1, ncol + qtmp = qn(i,k) - qme(i,k)*deltat + ttmp = tn(i,k) + deltat/cpair * ( meltheat(i,k) + (latvap + latice*fice(i,k)) * qme(i,k) ) + esn = estblf(ttmp) + qsn = svp_to_qsat(esn, pmid(i,k)) + qtl(i) = max((qsn - qtmp)/deltat,0._kind_phys) + relhum1(i) = qtmp/qsn + end do qtl_update_loop + + evap_update_loop: do i = 1, ncol +#ifdef PERGRO + evapprec(i,k) = conke*(1._kind_phys - max(cldm(i),mincld))* & + sqrt(precab(i))*(1._kind_phys - min(relhum1(i),1._kind_phys)) +#else + evapprec(i,k) = conke*(1._kind_phys - cldm(i))*sqrt(precab(i)) & + *(1._kind_phys - min(relhum1(i),1._kind_phys)) +#endif + + ! limit the evaporation to the amount which is entering the box + ! or saturates the box + prtmp = precab(i)*gravit/pdel(i,k) + evapprec(i,k) = min(evapprec(i,k), prtmp, qtl(i))*omsm +#ifdef PERGRO + ! zeroing needed for pert growth + evapprec(i,k) = 0._kind_phys +#endif + + ! Partition evaporation of precipitate between rain and snow using + ! the fraction of snow falling into the box. Determine the heating + ! due to evaporation. Note that evaporation is positive (loss of precip, + ! gain of vapor) and that heating is negative. + if (evapprec(i,k) > 0._kind_phys) then + evapsnow(i,k) = evapprec(i,k) * snowab(i) / precab(i) + evapheat(i,k) = -latvap * evapprec(i,k) - latice * evapsnow(i,k) + else + evapsnow(i,k) = 0._kind_phys + evapheat(i,k) = 0._kind_phys + end if + + ! Account for the latent heat of fusion for liquid drops collected by falling snow + prfzheat(i,k) = latice * liq2snow(i,k) + end do evap_update_loop + + ! now remove the residual of any over-saturation. Normally, + ! the oversaturated water vapor should have been removed by + ! qme formulation plus constraints by wet bulb tsp/qsp + ! as computed above. However, because of non-linearity, + ! addition of (qme-evapprec) to update t and q may still cause + ! a very small amount of over saturation. It is called a + ! residual of over-saturation because theoretically, qme + ! should have taken care of all of large scale condensation. + qmeres_update_loop: do i = 1, ncol + qtmp = qn(i,k)-(qme(i,k)-evapprec(i,k))*deltat + ttmp = tn(i,k) + deltat/cpair * ( meltheat(i,k) + evapheat(i,k) + prfzheat(i,k) & + + (latvap + latice*fice(i,k)) * qme(i,k) ) + + call qsat(ttmp, pmid(i,k), esn, qsn, dqsdt=dqsdt) + + if(qtmp > qsn) then + ! now extra condensation to bring air to just saturation + ctmp = (qtmp-qsn)/(1._kind_phys+hlocp*dqsdt)/deltat + qme(i,k) = qme(i,k)+ctmp + ! save residual on qmeres to addtion to qme on entering next iteration + ! qme exit here contain the residual but overrided if back to iteration + qmeres(i) = ctmp + else + qmeres(i) = 0.0_kind_phys + endif + end do qmeres_update_loop + end do qme_iter_loop ! loop over l (iteration) + + ! precipitation + precip_snow_loop: do i = 1, ncol + precip(i) = precip(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k)) + precab(i) = precab(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k)) + if(precab(i) < 0._kind_phys) then + precab(i) = 0._kind_phys + endif + + ! snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodprec(i,k)*fice(i,k) - evapsnow(i,k)) + snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodsnow(i,k) - evapsnow(i,k)) + + ! If temperature above freezing, all precip is rain flux. if temperature below freezing, all precip is snow flux. + lsflxprc(i,k+1) = precab(i) !! making this consistent with other precip fluxes. prc = rain + snow + ! lsflxprc(i,k+1) = precab(i) - snowab(i) + lsflxsnw(i,k+1) = snowab(i) + end do precip_snow_loop + end do mp_level_loop ! loop over k (level) + + ! Convert precip (lwe_stratiform_precipitation_rate_at_surface) + ! and snowab (lwe_snow_precipitation_rate_at_surface_due_to_microphysics) + ! from kg m-2 s-1 to m s-1 (precipitation units) for surface exchange. + ! + ! If this conversion is removed in the future, the metadata needs to + ! be updated. + precip(:ncol) = precip(:ncol)/1000._kind_phys + snowab(:ncol) = snowab(:ncol)/1000._kind_phys + end subroutine prognostic_cloud_water_run + + ! Calculate the conversion of condensate to precipitate + ! Rasch, P. J., and J. E. Kristjánsson, 1998: A Comparison of the CCM3 Model Climate Using Diagnosed and Predicted Condensate Parameterizations. J. Climate, 11, 1587–1614 + ! https://doi.org/10.1175/1520-0442(1998)011<1587:ACOTCM>2.0.CO;2 + ! + ! Original author: P. Rasch + subroutine findmcnew( & + ncol, pver, & + pi, & + k, & + precab, snowab, t, p, cwm, cldm, cldmax, fice, & + landm, seaicef, snowh, & + coef, fwaut, fsaut, & + fracw, fsacw, fsaci, & + pracwo, psacwo, psacio) + + ! input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of levels + real(kind_phys), intent(in) :: pi ! pi_constant [1] + integer, intent(in) :: k ! level index + + real(kind_phys), intent(in) :: precab(:) ! rate of precipitation from above [kg m-2 s-1] + real(kind_phys), intent(in) :: t(:,:) ! temperature [K] + real(kind_phys), intent(in) :: p(:,:) ! air_pressure [Pa] + real(kind_phys), intent(in) :: cwm(:) ! condensate mixing ratio [kg kg-1] + real(kind_phys), intent(in) :: cldm(:) ! cloud fraction + real(kind_phys), intent(in) :: cldmax(:) ! max cloud fraction above this level + real(kind_phys), intent(in) :: fice(:) ! fraction of cwat that is ice + real(kind_phys), intent(in) :: landm(:) ! Land fraction ramped over water + real(kind_phys), intent(in) :: seaicef(:) ! sea ice fraction + real(kind_phys), intent(in) :: snowab(:) ! rate of snow from above [kg m-2 s-1] + real(kind_phys), intent(in) :: snowh(:) ! Snow depth over land, water equivalent [m] + + ! output arguments + real(kind_phys), intent(out) :: coef(:) ! conversion rate [s-1] + real(kind_phys), intent(out) :: fwaut(:) ! relative importance of liquid autoconversion (a diagnostic) + real(kind_phys), intent(out) :: fsaut(:) ! relative importance of ice autoconversion (a diagnostic) + real(kind_phys), intent(out) :: fracw(:) ! relative importance of rain accreting liquid (a diagnostic) + real(kind_phys), intent(out) :: fsacw(:) ! relative importance of snow accreting liquid (a diagnostic) + real(kind_phys), intent(out) :: fsaci(:) ! relative importance of snow accreting ice (a diagnostic) + real(kind_phys), intent(out) :: pracwo(:) ! accretion of cloud water by rain [s-1] + real(kind_phys), intent(out) :: psacwo(:) ! accretion of cloud water by snow [s-1] + real(kind_phys), intent(out) :: psacio(:) ! accretion of cloud ice by snow [s-1] + + ! local variables + integer :: i, ii ! Loop index [index] + integer :: ncols ! Number of active columns for microphysics (different from ncol!!) [count] + integer :: ind(ncol) ! Active column indices [index] + real(kind_phys) :: capn ! Local cloud particle number concentration [cm-3] + real(kind_phys) :: capnoice ! Cloud particle concentration excluding sea ice [cm-3] + real(kind_phys) :: cldloc(ncol) ! Non-zero cloud fraction [1] + real(kind_phys) :: cldpr(ncol) ! Cloud fraction for precipitation [1] + real(kind_phys) :: totmr(ncol) ! In-cloud total water mixing ratio [kg kg-1] + real(kind_phys) :: icemr(ncol) ! In-cloud ice mixing ratio [kg kg-1] + real(kind_phys) :: liqmr(ncol) ! In-cloud liquid water mixing ratio [kg kg-1] + real(kind_phys) :: rainmr(ncol) ! In-cloud rain mixing ratio [kg kg-1] + real(kind_phys) :: ciaut ! Ice autoconversion coefficient [s-1] + real(kind_phys) :: pracw ! Rate of rain collecting cloud water [s-1] + real(kind_phys) :: psaci ! Rate of snow collecting cloud ice [s-1] + real(kind_phys) :: psacw ! Rate of snow collecting cloud water [s-1] + real(kind_phys) :: psaut ! Rate of ice autoconversion [s-1] + real(kind_phys) :: pwaut ! Rate of liquid water autoconversion [s-1] + real(kind_phys) :: ptot ! Total conversion rate [s-1] + real(kind_phys) :: prlloc(ncol) ! Local rain flux [mm day-1] + real(kind_phys) :: prscgs(ncol) ! Local snow amount [g cm-2] + real(kind_phys) :: snowfr ! Snow fraction of precipitation [1] + real(kind_phys) :: vfallw ! Fall speed of liquid precipitation [m s-1] + real(kind_phys) :: rho(ncol) ! Air density [kg m-3] + real(kind_phys) :: rhocgs ! Air density in CGS units [g cm-3] + real(kind_phys) :: r3l ! Cloud droplet volume radius [m] + real(kind_phys) :: icrit ! Ice autoconversion threshold [kg kg-1] + real(kind_phys) :: wsi ! Sea ice weight factor [1] + real(kind_phys) :: wt ! Ice fraction weight [1] + real(kind_phys) :: wland ! Land fraction weight [1] + real(kind_phys) :: wp ! Pressure dependence weight [1] + real(kind_phys) :: ftot ! Total fraction for conversion processes [1] + real(kind_phys) :: con1 ! Work constant for radius calculation [m] + real(kind_phys) :: con2 ! Work constant for density ratios [1] + real(kind_phys) :: csacx ! Constant used for snow accreting liquid or ice [??] + real(kind_phys) :: rat1 ! Density ratio work variable [1] + real(kind_phys) :: rat2 ! Mixing ratio ratio work variable [1] + + ! find all the points where we need to do the microphysics + ! and set the output variables to zero + ncols = 0 + do i = 1, ncol + coef(i) = 0._kind_phys + fwaut(i) = 0._kind_phys + fsaut(i) = 0._kind_phys + fracw(i) = 0._kind_phys + fsacw(i) = 0._kind_phys + fsaci(i) = 0._kind_phys + liqmr(i) = 0._kind_phys + rainmr(i) = 0._kind_phys + if (cwm(i) > 1.e-20_kind_phys) then + ncols = ncols + 1 + ind(ncols) = i + endif + end do + + do ii = 1, ncols + ! map from active microphysics columns to physics columns + i = ind(ii) + + ! the local cloudiness at this level + cldloc(i) = max(cldmin,cldm(i)) + + ! a weighted mean between max cloudiness above, and this layer + cldpr(i) = max(cldmin,(cldmax(i)+cldm(i))*0.5_kind_phys) + + ! decompose the suspended condensate into + ! an incloud liquid and ice phase component + totmr(i) = cwm(i)/cldloc(i) + icemr(i) = totmr(i)*fice(i) + liqmr(i) = totmr(i)*(1._kind_phys-fice(i)) + + ! density + rho(i) = p(i,k)/(287._kind_phys*t(i,k)) + rhocgs = rho(i)*1.e-3_kind_phys ! density in cgs units + + ! decompose the precipitate into a liquid and ice phase + if (t(i,k) > t0) then + vfallw = convfw/sqrt(rho(i)) + rainmr(i) = precab(i)/(rho(i)*vfallw*cldpr(i)) + snowfr = 0 + else + snowfr = 1 + rainmr(i) = 0._kind_phys + endif + + ! local snow amount in cgs units + prscgs(i) = precab(i)/cldpr(i)*0.1_kind_phys*snowfr + + ! local rain amount in mm/day + prlloc(i) = precab(i)*86400._kind_phys/cldpr(i) + end do + + con1 = 1._kind_phys/(1.333_kind_phys*pi)**0.333_kind_phys * 0.01_kind_phys ! meters + + ! calculate the conversion terms + do ii = 1, ncols + i = ind(ii) + rhocgs = rho(i)*1.e-3_kind_phys ! density in cgs units + + ! some temperature dependent constants + wt = fice(i) + icrit = icritc*wt + icritw*(1-wt) + + ! jrm Reworked droplet number concentration algorithm + ! Start with pressure-dependent value appropriate for continental air + ! Note: reltab has a temperature dependence here + capn = capnw + (capnc-capnw) * min(1._kind_phys,max(0._kind_phys,1.0_kind_phys-(p(i,k)-0.8_kind_phys*p(i,pver))/(0.2_kind_phys*p(i,pver)))) + ! Modify for snow depth over land + capn = capn + (capnc-capn) * min(1.0_kind_phys,max(0.0_kind_phys,snowh(i)*10._kind_phys)) + ! Ramp between polluted value over land to clean value over ocean. + capn = capn + (capnc-capn) * min(1.0_kind_phys,max(0.0_kind_phys,1.0_kind_phys-landm(i))) + ! Ramp between the resultant value and a sea ice value in the presence of ice. + capn = capn + (capnsi-capn) * min(1.0_kind_phys,max(0.0_kind_phys,seaicef(i))) + ! end jrm + + ! useful terms in following calculations + rat1 = rhocgs/rhow + rat2 = liqmr(i)/capn + con2 = (rat1*rat2)**0.333_kind_phys + + ! volume radius + r3l = con1*con2 + + ! critical threshold for autoconversion if modified for mixed phase + ! clouds to mimic a bergeron findeisen process + ! r3lc2 = r3lcrit*(1.-0.5*fice(i)*(1-fice(i))) + ! + ! autoconversion of liquid + ! + ! cwaut = 2.e-4 + ! cwaut = 1.e-3 + ! lcrit = 2.e-4 + ! lcrit = 5.e-4 + ! pwaut = max(0._kind_phys,liqmr(i)-lcrit)*cwaut + ! + ! pwaut is following tripoli and cotton (and many others) + ! we reduce the autoconversion below critpr, because these are regions where + ! the drop size distribution is likely to imply much smaller collector drops than + ! those relevant for a cloud distribution corresponding to the value of effc = 0.55 + ! suggested by cotton (see austin 1995 JAS, baker 1993) + + ! easy to follow form + ! pwaut = capc*liqmr(i)**2*rhocgs/rhow + ! *(liqmr(i)*rhocgs/(rhow*capn))**(.333) + ! *heavy(r3l,r3lcrit) + ! *max(0.10_kind_phys,min(1._kind_phys,prlloc(i)/critpr)) + ! somewhat faster form + + ! using modified Heaviside function: + pwaut = capc*liqmr(i)**2*rat1*con2*heavymp(r3l,r3lcrit) * & + max(0.10_kind_phys,min(1._kind_phys,prlloc(i)/critpr)) + + ! not using modified Heaviside function: + ! pwaut = capc*liqmr(i)**2*rat1*con2*heavym(r3l,r3lcrit)* & + ! max(0.10_kind_phys,min(1._kind_phys,prlloc(i)/critpr)) + + ! autoconversion of ice + ciaut = ciautb + + ! autoconversion of ice condensate +#ifdef PERGRO + psaut = heavyp(icemr(i),icrit)*icemr(i)*ciaut +#else + psaut = max(0._kind_phys,icemr(i)-icrit)*ciaut +#endif + + ! collection of liquid by rain + ! pracw = cracw*rho(i)*liqmr(i)*rainmr(i) !(beheng 1994) + pracw = cracw*rho(i)*sqrt(rho(i))*liqmr(i)*rainmr(i) !(tripoli and cotton) + pracwo(i) = pracw + + ! the following lines calculate the slope parameter and snow mixing ratio + ! from the precip rate using the equations found in lin et al 83 + ! in the most natural form, but it is expensive, so after some tedious + ! algebraic manipulation you can use the cheaper form found below + ! vfalls = c*gam4pd/(6*lamdas**d)*sqrt(rhonot/rhocgs) + ! $ *0.01 ! convert from cm/s to m/s + ! snowmr(i) = snowfr*precab(i)/(rho(i)*vfalls*cldpr(i)) + ! snowmr(i) = ( prscgs(i)*mcon02 * (rhocgs**mcon03) )**mcon04 + ! lamdas = (prhonos/max(rhocgs*snowmr(i),small))**0.25 + ! csacw = mcon01*sqrt(rhonot/rhocgs)/(lamdas**thrpd) + ! + ! coefficient for collection by snow independent of phase + csacx = mcon07*rhocgs**mcon08*prscgs(i)**mcon05 + + ! collection of liquid by snow (lin et al 1983) + psacw = csacx*liqmr(i)*esw + +#ifdef PERGRO + ! this is necessary for pergro + psacw = 0._kind_phys +#endif + + psacwo(i) = psacw + + ! collection of ice by snow (lin et al 1983) + psaci = csacx*icemr(i)*esi + psacio(i) = psaci + + ! total conversion of condensate to precipitate + ptot = pwaut + psaut + pracw + psacw + psaci + + ! the reciprocal of cloud water amount (or zero if no cloud water) + ! rcwm = totmr(i)/(max(totmr(i),small)**2) + + ! turn the tendency back into a loss rate [s-1] + if (totmr(i) > 0._kind_phys) then + coef(i) = ptot/totmr(i) + else + coef(i) = 0._kind_phys + endif + + if (ptot.gt.0._kind_phys) then + fwaut(i) = pwaut/ptot + fsaut(i) = psaut/ptot + fracw(i) = pracw/ptot + fsacw(i) = psacw/ptot + fsaci(i) = psaci/ptot + else + fwaut(i) = 0._kind_phys + fsaut(i) = 0._kind_phys + fracw(i) = 0._kind_phys + fsacw(i) = 0._kind_phys + fsaci(i) = 0._kind_phys + endif + + ftot = fwaut(i)+fsaut(i)+fracw(i)+fsacw(i)+fsaci(i) + end do + + end subroutine findmcnew + + ! For special WACCM/CAM-chem cases with CAM4 physics: + ! Sets rhu to a different value poleward of +/- 50 deg latitude and + ! levels above the tropopause if polstrat_rhmin is specified + subroutine relhum_min_adj(ncol, pver, tropLev, dlat, rhu, rhu_adj) + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: tropLev(:) + real(kind_phys), intent(in) :: dlat(:) ! latitudes in degrees + real(kind_phys), intent(in) :: rhu(:,:) + real(kind_phys), intent(out) :: rhu_adj(:,:) + + integer :: i, k + + rhu_adj(:,:) = rhu(:,:) + if (.not. do_psrhmin) return + + do k = 1, pver + do i = 1, ncol + ! assumption that up is lower index + if ((k .lt. tropLev(i)) .and. (abs(dlat(i)) .gt. 50._kind_phys)) then + rhu_adj(i,k) = psrhmin + endif + enddo + enddo + end subroutine relhum_min_adj + + ! Heavyside step functions used in findmcnew + ! Heavyside step function + pure function heavy(a1, a2) result(h) + ! h: 1 if a1 > a2, 0 otherwise + real(kind_phys), intent(in) :: a1, a2 + real(kind_phys) :: h + h = max(0.0_kind_phys, sign(1.0_kind_phys, a1-a2)) + end function heavy + + ! Modified Heavyside function with minimum value of 0.01 + pure function heavym(a1, a2) result(h) + ! h: 1 if a1 > a2, 0.01 otherwise + real(kind_phys), intent(in) :: a1, a2 + real(kind_phys) :: h + h = max(0.01_kind_phys, sign(1.0_kind_phys, a1-a2)) + end function heavym + + ! Smoothed Heavyside function using ratio formulation + pure function heavyp(a1, a2) result(h) + ! h: Ratio of a1/(a1+a2+eps) + real(kind_phys), intent(in) :: a1, a2 + real(kind_phys) :: h + real(kind_phys), parameter :: eps = 1.0e-36_kind_phys + h = a1/(a2 + a1 + eps) + end function heavyp + + ! Modified smooth Heavyside with offset term + pure function heavymp(a1, a2) result(h) + ! h: Ratio (a1 + 0.01*a2)/(a1+a2+eps) + real(kind_phys), intent(in) :: a1, a2 + real(kind_phys) :: h + real(kind_phys), parameter :: eps = 1.0e-36_kind_phys + h = (a1 + 0.01_kind_phys*a2)/(a2 + a1 + eps) + end function heavymp + +end module prognostic_cloud_water diff --git a/schemes/rasch_kristjansson/prognostic_cloud_water.meta b/schemes/rasch_kristjansson/prognostic_cloud_water.meta new file mode 100644 index 00000000..6f59d935 --- /dev/null +++ b/schemes/rasch_kristjansson/prognostic_cloud_water.meta @@ -0,0 +1,437 @@ +[ccpp-table-properties] + name = prognostic_cloud_water + type = scheme + +[ccpp-arg-table] + name = prognostic_cloud_water_init + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ tmelt ] + standard_name = freezing_point_of_water + units = K + type = real | kind = kind_phys + dimensions = () + intent = in +[ rhodair ] + standard_name = density_of_dry_air_at_stp + units = kg m-3 + type = real | kind = kind_phys + dimensions = () + intent = in +[ pi ] + standard_name = pi_constant + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ icritc_in ] + standard_name = tunable_parameter_for_autoconversion_of_cold_ice_for_rk_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ icritw_in ] + standard_name = tunable_parameter_for_autoconversion_of_warm_ice_for_rk_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ conke_in ] + standard_name = tunable_parameter_for_precipitation_evaporation_for_rk_microphysics + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ r3lcrit_in ] + standard_name = tunable_parameter_for_cloud_water_autoconversion_for_rk_microphysics + units = m + type = real | kind = kind_phys + dimensions = () + intent = in +[ do_psrhmin_in ] + standard_name = flag_for_relative_humidity_threshold_for_cloud_formation_in_polar_stratosphere_for_rk_microphysics + units = flag + type = logical + dimensions = () + intent = in +[ psrhmin_in ] + standard_name = relative_humidity_threshold_for_cloud_formation_in_polar_stratosphere_for_rk_microphysics + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = prognostic_cloud_water_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ top_lev ] + standard_name = vertical_layer_index_of_troposphere_cloud_physics_top + units = index + type = integer + dimensions = () + intent = in +[ deltat ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ pi ] + standard_name = pi_constant + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rh2o ] + standard_name = gas_constant_of_water_vapor + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ epsilo ] + standard_name = ratio_of_water_vapor_to_dry_air_molecular_weights + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ latvap ] + standard_name = latent_heat_of_vaporization_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ latice ] + standard_name = latent_heat_of_fusion_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ cpair ] + standard_name = specific_heat_of_dry_air_at_constant_pressure + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ dlat ] + standard_name = latitude_degrees_north + units = degrees + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ pmid ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ pdel ] + standard_name = air_pressure_thickness + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ zi ] + standard_name = geopotential_height_wrt_surface_at_interface + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ tropLev ] + standard_name = tropopause_vertical_layer_index + units = index + type = integer + dimensions = (horizontal_loop_extent) + intent = in +[ ttend ] + standard_name = tendency_of_air_temperature_not_due_to_microphysics_tbd + units = K s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ tn ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ qtend ] + standard_name = tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_not_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ qn ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ ltend ] + standard_name = tendency_of_cloud_water_mixing_ratio_wrt_moist_air_and_condensed_water_not_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldice ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldliq ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ omega ] + standard_name = lagrangian_tendency_of_air_pressure + units = Pa s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldn ] + standard_name = cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ fice ] + standard_name = mass_fraction_of_ice_content_within_stratiform_cloud + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ fsnow ] + standard_name = mass_fraction_of_snow_content_within_stratiform_cloud + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rhdfda ] + standard_name = derivative_of_relative_humidity_wrt_cloud_fraction_tbd + units = percent + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rhu00 ] + standard_name = relative_humidity_threshold_for_prognostic_cloud_water_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ landm ] + standard_name = smoothed_land_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ seaicef ] + standard_name = sea_ice_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snowh ] + standard_name = lwe_surface_snow_depth_over_land + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ qme ] + standard_name = net_condensation_rate_due_to_microphysics + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ prodprec ] + standard_name = precipitation_production_due_to_microphysics + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ prodsnow ] + standard_name = snow_production_due_to_microphysics + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ evapprec ] + standard_name = precipitation_evaporation_due_to_microphysics + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ evapsnow ] + standard_name = rate_of_evaporation_of_falling_snow_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ evapheat ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_precipitation_evaporation_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ prfzheat ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_precipitation_freezing_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ meltheat ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_precipitation_phase_change_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ precip ] + standard_name = lwe_stratiform_precipitation_rate_at_surface + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ snowab ] + standard_name = lwe_snow_precipitation_rate_at_surface_due_to_microphysics + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ ice2pr ] + standard_name = tendency_of_cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ liq2pr ] + standard_name = tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ liq2snow ] + standard_name = tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water_due_to_conversion_to_snow_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ lsflxprc ] + standard_name = stratiform_rain_and_snow_flux_at_interface + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = out +[ lsflxsnw ] + standard_name = stratiform_snow_flux_at_interface + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = out +[ pracwo ] + standard_name = accretion_of_cloud_liquid_water_by_rain_tbd + units = s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ psacwo ] + standard_name = accretion_of_cloud_liquid_water_by_snow_tbd + units = s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ psacio ] + standard_name = accretion_of_cloud_ice_by_snow_tbd + units = s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ fwaut ] + standard_name = relative_importance_of_liquid_autoconversion_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ fsaut ] + standard_name = relative_importance_of_ice_autoconversion_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ fracw ] + standard_name = relative_importance_of_rain_accreting_liquid_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ fsacw ] + standard_name = relative_importance_of_snow_accreting_liquid_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ fsaci ] + standard_name = relative_importance_of_snow_accreting_ice_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/rasch_kristjansson/prognostic_cloud_water_namelist.xml b/schemes/rasch_kristjansson/prognostic_cloud_water_namelist.xml new file mode 100644 index 00000000..7e2250af --- /dev/null +++ b/schemes/rasch_kristjansson/prognostic_cloud_water_namelist.xml @@ -0,0 +1,167 @@ + + + + + + + + + real + conv + rk_stratiform_nl + tunable_parameter_for_autoconversion_of_cold_ice_for_rk_microphysics + kg kg-1 + + Threshold for autoconversion of cold ice in RK microphysics scheme. + + + 5.0e-6 + 18.0e-6 + + + + real + conv + rk_stratiform_nl + tunable_parameter_for_autoconversion_of_warm_ice_for_rk_microphysics + kg kg-1 + + Threshold for autoconversion of cold ice in RK microphysics scheme. + + + 4.0e-4 + 2.0e-4 + + + + real + conv + rk_stratiform_nl + tunable_parameter_for_precipitation_evaporation_for_rk_microphysics + 1 + + + Tunable constant for evaporation of precip in RK microphysics scheme. + + + 10.0e-6 + 5.0e-6 + + + + real + conv + rk_stratiform_nl + tunable_parameter_for_cloud_water_autoconversion_for_rk_microphysics + m + + Critical radius at which autoconversion become efficient in RK microphysics scheme. + + + 10.0e-6 + + + + real + conv + rk_stratiform_nl + relative_humidity_threshold_for_cloud_formation_in_polar_stratosphere_for_rk_microphysics + 1 + + Relative humidity threshold for stratospheric cloud water condensation in RK microphysics + poleward of 50 degrees. Used if rk_strat_polstrat_dorhmin is set to true. + + + 1.200D0 + + + + logical + conv + rk_stratiform_nl + flag_for_relative_humidity_threshold_for_cloud_formation_in_polar_stratosphere_for_rk_microphysics + flag + + Flag to set rhu to a different value in the stratosphere poleward of 50 degrees for stratospheric + cloud water condensation. Used for special WACCM/CAM-chem cases with CAM4 physics. + + + .false. + + + diff --git a/schemes/rasch_kristjansson/rk_stratiform.F90 b/schemes/rasch_kristjansson/rk_stratiform.F90 new file mode 100644 index 00000000..64da34dd --- /dev/null +++ b/schemes/rasch_kristjansson/rk_stratiform.F90 @@ -0,0 +1,556 @@ +! Rasch and Kristjansson prognostic cloud microphysics and CAM4 macrophysics +! CCPP-ized: Haipeng Lin, January 2025 +module rk_stratiform + use ccpp_kinds, only: kind_phys + + implicit none + private + save + + ! public CCPP-compliant subroutines + ! note: cloud_fraction_perturbation_run calls the compute_cloud_fraction + ! scheme run phase for perturbation with a modified rhpert_flag = .true. + ! + ! refer to test SDF suite_rasch_kristjansson.xml for total order of operations, + ! as the full RK-stratiform requires other schemes not included in this module. + public :: rk_stratiform_check_qtlcwat_run + ! -- cloud_particle_sedimentation -- + public :: rk_stratiform_sedimentation_run + public :: rk_stratiform_detrain_convective_condensate_run + public :: rk_stratiform_cloud_fraction_perturbation_run ! see note. + public :: rk_stratiform_external_forcings_run + public :: rk_stratiform_condensate_repartioning_run + ! -- prognostic_cloud_water -- + public :: rk_stratiform_prognostic_cloud_water_tendencies_run + public :: rk_stratiform_cloud_optical_properties_run + public :: rk_stratiform_save_qtlcwat_run + + ! + +contains + +!> \section arg_table_rk_stratiform_check_qtlcwat_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_check_qtlcwat_run.html + subroutine rk_stratiform_check_qtlcwat_run( & + ncol, pver, & + t, q_wv, cldice, cldliq, & + qcwat, tcwat, lcwat, & ! from end of last microphysics/macrophysics call. + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(kind_phys), intent(in) :: t(:,:) ! air_temperature [K] + real(kind_phys), intent(in) :: q_wv(:, :) ! adv: water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldice(:,:) ! adv: cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldliq(:,:) ! adv: cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + + ! Input/output arguments + real(kind_phys), intent(inout) :: qcwat(:,:) ! [kg kg-1] + real(kind_phys), intent(inout) :: tcwat(:,:) ! [K] + real(kind_phys), intent(inout) :: lcwat(:,:) ! [kg kg-1] + + ! Output arguments + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + errmsg = '' + errflg = 0 + + ! Check that qcwat and tcwat were initialized - if not then initialize + ! this is made as a separate "run" scheme so it does not have to be used in current CAM + + ! lcwat is initialized from initial conditions of cldice, cldliq + ! TODO: check if this is always done (appears to be from physpkg.F90) or should be read from snapshot. + if(qcwat(1,1) < 0._kind_phys .or. tcwat(1,1) < 0._kind_phys) then + lcwat(:ncol,:) = cldice(:ncol,:) + cldliq(:ncol,:) + endif + + ! The registry will set default negative values if not read from snapshot. + if(qcwat(1,1) < 0._kind_phys) then + qcwat(:ncol,:) = q_wv(:ncol,:) + endif + + if(tcwat(1,1) < 0._kind_phys) then + tcwat(:ncol,:) = t(:ncol,:) + endif + + end subroutine rk_stratiform_check_qtlcwat_run + +!> \section arg_table_rk_stratiform_sedimentation_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_sedimentation_run.html + subroutine rk_stratiform_sedimentation_run( & + ncol, & + sfliq, snow_sed, & + prec_sed, & + prec_str, snow_str, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + real(kind_phys), intent(in) :: sfliq(:) ! stratiform_rain_flux_at_surface_due_to_sedimentation [kg m-2 s-1] + real(kind_phys), intent(in) :: snow_sed(:) ! sfice = lwe_cloud_ice_sedimentation_rate_at_surface_due_to_microphysics [m s-1] + + ! Output arguments + real(kind_phys), intent(out) :: prec_sed(:) ! stratiform_cloud_water_surface_flux_due_to_sedimentation [m s-1] + real(kind_phys), intent(out) :: prec_str(:) ! lwe_large_scale_precipitation_rate_at_surface [m s-1] + real(kind_phys), intent(out) :: snow_str(:) ! lwe_snow_and_cloud_ice_precipitation_rate_at_surface_due_to_microphysics [m s-1] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + errmsg = '' + errflg = 0 + + ! Convert rain flux to precip units from mass units + ! and create cloud water surface flux (rain + snow) + prec_sed(:ncol) = sfliq(:ncol)/1000._kind_phys + snow_sed(:ncol) + + ! Start accumulation of precipitation and snow flux [m s-1] + prec_str(:ncol) = 0._kind_phys + prec_sed(:ncol) + snow_str(:ncol) = 0._kind_phys + snow_sed(:ncol) + + end subroutine rk_stratiform_sedimentation_run + +!> \section arg_table_rk_stratiform_detrain_convective_condensate_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_detrain_convective_condensate_run.html + subroutine rk_stratiform_detrain_convective_condensate_run( & + ncol, & + dlf, & + rliq, & + prec_str, & + tend_cldliq, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + real(kind_phys), intent(in) :: dlf(:,:) ! detrainment_of_cloud_liquid_water_wrt_moist_air_and_condensed_water_due_to_all_convection [kg kg-1 s-1] + real(kind_phys), intent(in) :: rliq(:) ! vertically_integrated_cloud_liquid_water_tendency_due_to_all_convection_to_be_applied_later_in_time_loop [m s-1] + + ! Input/output arguments + real(kind_phys), intent(inout) :: prec_str(:) ! lwe_large_scale_precipitation_rate_at_surface [m s-1] + + ! Output arguments + real(kind_phys), intent(out) :: tend_cldliq(:,:) ! tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1 s-1] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + errmsg = '' + errflg = 0 + + ! Apply detrainment tendency to cloud liquid water + tend_cldliq(:ncol,:) = dlf(:ncol,:) + + ! Accumulate precipitation and snow after + ! reserved liquid (vertical integral) has now been used + ! (snow contribution is zero) + prec_str(:ncol) = prec_str(:ncol) - rliq(:ncol) + + end subroutine rk_stratiform_detrain_convective_condensate_run + + ! Call perturbed cloud fraction and compute perturbation threshold criteria + ! necessary for prognostic_cloud_water scheme +!> \section arg_table_rk_stratiform_cloud_fraction_perturbation_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_cloud_fraction_perturbation_run.html + subroutine rk_stratiform_cloud_fraction_perturbation_run( & + ncol, pver, & + cappa, gravit, rair, tmelt, pref, lapse_rate, & + top_lev_cloudphys, & + pmid, ps, temp, sst, & + q_wv, cldice, & + phis, & + shallowcu, deepcu, concld, & ! inputs from convective_cloud_cover + landfrac, ocnfrac, snowh, & + cloud, relhum, rhu00, & ! inputs from unperturbed compute_cloud_fraction + rhdfda, & ! output for prognostic_cloud_water + errmsg, errflg) + + ! Dependency: compute_cloud_fraction CCPPized scheme run phase. + ! this scheme is called with an altered rhpert_flag = .true. + ! then the outputs are combined with the "regular" output of the compute_cloud_fraction + ! CCPP scheme to get the perturbed quantities for prognostic_cloud_water. + use compute_cloud_fraction, only: compute_cloud_fraction_run + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(kind_phys), intent(in) :: cappa + real(kind_phys), intent(in) :: gravit + real(kind_phys), intent(in) :: rair + real(kind_phys), intent(in) :: tmelt + real(kind_phys), intent(in) :: pref + real(kind_phys), intent(in) :: lapse_rate + integer, intent(in) :: top_lev_cloudphys ! vertical_layer_index_of_cloud_fraction_top [index] + real(kind_phys), intent(in) :: pmid(:, :) ! air_pressure [Pa] + real(kind_phys), intent(in) :: ps(:) ! surface_air_pressure [Pa] + real(kind_phys), intent(in) :: temp(:, :) ! air_temperature [K] + real(kind_phys), intent(in) :: sst(:) ! sea_surface_temperature [K] + real(kind_phys), intent(in) :: q_wv(:, :) ! adv: water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldice(:, :) ! adv: cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: phis(:) ! surface_geopotential [m2 s-2] + real(kind_phys), intent(in) :: shallowcu(:, :) ! shallow convective cloud fraction + real(kind_phys), intent(in) :: deepcu(:, :) ! deep convective cloud fraction + real(kind_phys), intent(in) :: concld(:, :) ! convective_cloud_area_fraction [fraction] + real(kind_phys), intent(in) :: landfrac(:) ! land_area_fraction [fraction] + real(kind_phys), intent(in) :: ocnfrac(:) ! ocean_area_fraction [fraction] + real(kind_phys), intent(in) :: snowh(:) ! lwe_surface_snow_depth_over_land [m] + + real(kind_phys), intent(in) :: cloud(:, :) ! cloud_area_fraction [fraction] + real(kind_phys), intent(in) :: relhum(:, :) ! RH for prognostic cldwat [percent] + + ! Input/output arguments + ! Note: CAM4 intentionally mutates the top level of rhu00 so this is inout here. + real(kind_phys), intent(inout) :: rhu00(:, :) ! RH threshold for cloud [fraction] + + ! Output arguments + real(kind_phys), intent(out) :: rhdfda(:, :) ! derivative of RH w.r.t. cloud fraction for prognostic cloud water [percent] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local variables + integer :: i, k + + ! Local variables (outputs from perturbed compute_cloud_fraction) + real(kind_phys) :: cloud2(ncol, pver) + + ! Dummy outputs (unused) + real(kind_phys) :: rhcloud2(ncol, pver) + real(kind_phys) :: cldst2(ncol, pver) + real(kind_phys) :: rhu002(ncol, pver) + real(kind_phys) :: icecldf2(ncol, pver) + real(kind_phys) :: liqcldf2(ncol, pver) + real(kind_phys) :: relhum2(ncol, pver) + + errmsg = '' + errflg = 0 + + ! Call perturbed version of compute_cloud_fraction scheme + call compute_cloud_fraction_run( & + ncol = ncol, & + pver = pver, & + cappa = cappa, & + gravit = gravit, & + rair = rair, & + tmelt = tmelt, & + pref = pref, & + lapse_rate = lapse_rate, & + top_lev_cloudphys = top_lev_cloudphys, & ! CAM4 macrophysics - top lev is 1 + pmid = pmid(:ncol,:), & + ps = ps(:ncol), & + temp = temp(:ncol,:), & + sst = sst(:ncol), & + q = q_wv(:ncol,:), & + cldice = cldice(:ncol,:), & + phis = phis(:ncol), & + shallowcu = shallowcu(:ncol,:), & + deepcu = deepcu(:ncol,:), & + concld = concld(:ncol,:), & + landfrac = landfrac(:ncol), & + ocnfrac = ocnfrac(:ncol), & + snowh = snowh(:ncol), & + rhpert_flag = .true., & ! ** apply perturbation here ** + cloud = cloud2(:ncol, :), & + rhcloud = rhcloud2(:ncol, :), & + cldst = cldst2(:ncol,:), & + rhu00 = rhu002(:ncol,:), & + icecldf = icecldf2(:ncol,:), & + liqcldf = liqcldf2(:ncol,:), & + relhum = relhum2(:ncol,:), & + errmsg = errmsg, & + errflg = errflg) + + ! Compute rhdfda (derivative of RH w.r.t. cloud fraction) + ! for use in the prognostic_cloud_water scheme. + rhu00(:ncol,1) = 2.0_kind_phys ! arbitrary number larger than 1 (100%) + do k = 1, pver + do i = 1, ncol + if( relhum(i,k) < rhu00(i,k) ) then + rhdfda(i,k) = 0.0_kind_phys + elseif( relhum(i,k) >= 1.0_kind_phys ) then + rhdfda(i,k) = 0.0_kind_phys + else + ! Under certain circumstances, rh+ cause cld not to changed + ! when at an upper limit, or w/ strong subsidence + if( ( cloud2(i,k) - cloud(i,k) ) < 1.e-4_kind_phys ) then + rhdfda(i,k) = 0.01_kind_phys*relhum(i,k)*1.e+4_kind_phys + else + rhdfda(i,k) = 0.01_kind_phys*relhum(i,k)/(cloud2(i,k)-cloud(i,k)) + endif + endif + enddo + enddo + + end subroutine rk_stratiform_cloud_fraction_perturbation_run + + + ! Compute non-micro and non-macrophysical external forcings + ! for computing of net condensation rate. + ! Note: advective forcing of condensate is aggregated into liquid phase. +!> \section arg_table_rk_stratiform_external_forcings_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_external_forcings_run.html + subroutine rk_stratiform_external_forcings_run( & + ncol, pver, & + dtime, & + t, & + q_wv, cldice, cldliq, & + qcwat, tcwat, lcwat, & ! from end of last physics timestep. + qtend, ttend, ltend, & ! output for prognostic_cloud_water + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(kind_phys), intent(in) :: dtime + real(kind_phys), intent(in) :: t(:,:) ! air_temperature [K] + real(kind_phys), intent(in) :: q_wv(:,:) ! adv: water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldice(:,:) ! adv: cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldliq(:,:) ! adv: cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + + real(kind_phys), intent(in) :: qcwat(:,:) ! [kg kg-1] + real(kind_phys), intent(in) :: tcwat(:,:) ! [K] + real(kind_phys), intent(in) :: lcwat(:,:) ! [kg kg-1] + + ! Output arguments (for prognostic_cloud_water) + real(kind_phys), intent(out) :: qtend(:,:) ! not due to micro/macrophysics [kg kg-1 s-1] + real(kind_phys), intent(out) :: ttend(:,:) ! not due to micro/macrophysics [K s-1] + real(kind_phys), intent(out) :: ltend(:,:) ! not due to micro/macrophysics [kg kg-1 s-1] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local arguments + real(kind_phys) :: totcw(ncol, pver) + + errmsg = '' + errflg = 0 + + totcw(:ncol,:) = cldice(:ncol,:) + cldliq(:ncol,:) + + qtend(:ncol,:pver) = 1.0_kind_phys / dtime * (q_wv (:ncol,:pver) - qcwat(:ncol,:pver)) + ttend(:ncol,:pver) = 1.0_kind_phys / dtime * (t (:ncol,:pver) - tcwat(:ncol,:pver)) + ltend(:ncol,:pver) = 1.0_kind_phys / dtime * (totcw (:ncol,:pver) - lcwat(:ncol,:pver)) + + end subroutine rk_stratiform_external_forcings_run + + ! Repartitioning of stratiform condensate, + ! and compute repartition heating from change in cloud ice +!> \section arg_table_rk_stratiform_condensate_repartioning_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_condensate_repartioning_run.html + subroutine rk_stratiform_condensate_repartioning_run( & + ncol, pver, & + dtime, & + latice, & + cldice, cldliq, & + fice, & ! from cloud_fraction_fice + repartht, & + tend_cldice, & + tend_cldliq, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(kind_phys), intent(in) :: dtime + real(kind_phys), intent(in) :: latice + real(kind_phys), intent(in) :: cldice(:,:) ! adv: cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldliq(:,:) ! adv: cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: fice(:,:) ! mass_fraction_of_ice_content_within_stratiform_cloud [fraction] + + ! Input/output arguments + + ! Output arguments + real(kind_phys), intent(out) :: repartht(:,:) ! [J kg-1 s-1] + real(kind_phys), intent(out) :: tend_cldice(:,:) ! tendency_of_cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1 s-1] + real(kind_phys), intent(out) :: tend_cldliq(:,:) ! tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1 s-1] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local arguments + real(kind_phys) :: totcw(ncol, pver) + + errmsg = '' + errflg = 0 + + totcw(:ncol,:) = cldice(:ncol,:) + cldliq(:ncol,:) + tend_cldice(:ncol,:) = 1.0_kind_phys / dtime * ( totcw(:ncol,:)*fice(:ncol,:) - cldice(:ncol,:) ) + tend_cldliq(:ncol,:) = 1.0_kind_phys / dtime * ( totcw(:ncol,:)*(1.0_kind_phys-fice(:ncol,:)) - cldliq(:ncol,:) ) + + repartht(:ncol,:pver) = latice * tend_cldice(:ncol,:pver) + + end subroutine rk_stratiform_condensate_repartioning_run + + ! Determine tendencies from prognostic cloud water +!> \section arg_table_rk_stratiform_prognostic_cloud_water_tendencies_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_prognostic_cloud_water_tendencies_run.html + subroutine rk_stratiform_prognostic_cloud_water_tendencies_run( & + ncol, pver, & + dtime, & + latvap, latice, & + qme, fice, & + evapheat, prfzheat, meltheat, & + repartht, & + evapprec, & + ice2pr, liq2pr, & + prec_pcw, snow_pcw, & + prec_str, snow_str, & + cmeheat, cmeice, cmeliq, & + tend_s, & + tend_q, & + tend_cldice, & + tend_cldliq, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(kind_phys), intent(in) :: dtime + real(kind_phys), intent(in) :: latvap + real(kind_phys), intent(in) :: latice + real(kind_phys), intent(in) :: qme(:,:) ! net_condensation_rate_due_to_microphysics [s-1] + real(kind_phys), intent(in) :: fice(:,:) ! mass_fraction_of_ice_content_within_stratiform_cloud [fraction] + real(kind_phys), intent(in) :: evapheat(:,:) ! + real(kind_phys), intent(in) :: prfzheat(:,:) ! + real(kind_phys), intent(in) :: meltheat(:,:) ! + real(kind_phys), intent(in) :: repartht(:,:) ! (from microphysical tend) + real(kind_phys), intent(in) :: evapprec(:,:) ! + real(kind_phys), intent(in) :: ice2pr(:,:) ! + real(kind_phys), intent(in) :: liq2pr(:,:) ! + real(kind_phys), intent(in) :: prec_pcw(:) ! lwe_stratiform_precipitation_rate_at_surface [m s-1] + real(kind_phys), intent(in) :: snow_pcw(:) ! lwe_snow_precipitation_rate_at_surface_due_to_microphysics [m s-1] + + ! Input/output arguments + real(kind_phys), intent(inout) :: prec_str(:) ! lwe_large_scale_precipitation_rate_at_surface [m s-1] + real(kind_phys), intent(inout) :: snow_str(:) ! lwe_snow_and_cloud_ice_precipitation_rate_at_surface_due_to_microphysics [m s-1] + + ! Output arguments + real(kind_phys), intent(out) :: cmeheat(:,:) ! ... [J kg-1 s-1] + real(kind_phys), intent(out) :: cmeice(:,:) ! ... [kg kg-1 s-1] + real(kind_phys), intent(out) :: cmeliq(:,:) ! ... [kg kg-1 s-1] + real(kind_phys), intent(out) :: tend_s(:,:) ! tendency_of_dry_air_enthalpy_at_constant_pressure [J kg-1 s-1] + real(kind_phys), intent(out) :: tend_q(:,:) ! tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1 s-1] + real(kind_phys), intent(out) :: tend_cldice(:,:) ! tendency_of_cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1 s-1] + real(kind_phys), intent(out) :: tend_cldliq(:,:) ! tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1 s-1] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + integer :: i, k + errmsg = '' + errflg = 0 + + do k = 1, pver + do i = 1, ncol + ! Heating from cond-evap within the cloud [J kg-1 s-1] + cmeheat(i,k) = qme(i,k)*(latvap + latice*fice(i,k)) + + ! Rate of cond-evap of ice within the cloud [kg kg-1 s-1] + cmeice(i,k) = qme(i,k)*fice(i,k) + + ! Rate of cond-evap of liq within the cloud [kg kg-1 s-1] + cmeliq(i,k) = qme(i,k)*(1._kind_phys-fice(i,k)) + + ! Tendencies from after prognostic_cloud_water... + tend_s(i,k) = cmeheat(i,k) + & + evapheat(i,k) + prfzheat(i,k) + meltheat(i,k) + repartht(i,k) + tend_q(i,k) = - qme(i,k) + evapprec(i,k) + tend_cldice(i,k) = cmeice(i,k) - ice2pr(i,k) + tend_cldliq(i,k) = cmeliq(i,k) - liq2pr(i,k) + end do + end do + + prec_str(:ncol) = prec_str(:ncol) + prec_pcw(:ncol) + snow_str(:ncol) = snow_str(:ncol) + snow_pcw(:ncol) + + end subroutine rk_stratiform_prognostic_cloud_water_tendencies_run + + ! Save Q, T, cloud water at end of stratiform microphysics for use in next timestep + ! to determine non-microphysical/macrophysical tendencies +!> \section arg_table_rk_stratiform_save_qtlcwat_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_save_qtlcwat_run.html + subroutine rk_stratiform_save_qtlcwat_run( & + ncol, pver, & + t, & + q_wv, cldice, cldliq, & + qcwat, tcwat, lcwat, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + + real(kind_phys), intent(in) :: t(:,:) ! air_temperature [K] + real(kind_phys), intent(in) :: q_wv(:, :) ! adv: water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldice(:,:) ! adv: cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + real(kind_phys), intent(in) :: cldliq(:,:) ! adv: cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water [kg kg-1] + + ! Output arguments + real(kind_phys), intent(out) :: qcwat(:,:) ! [kg kg-1] + real(kind_phys), intent(out) :: tcwat(:,:) ! [K] + real(kind_phys), intent(out) :: lcwat(:,:) ! [kg kg-1] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local variables + integer :: k + + errmsg = '' + errflg = 0 + + do k = 1, pver + qcwat(:ncol,k) = q_wv(:ncol,k) + tcwat(:ncol,k) = t(:ncol,k) + lcwat(:ncol,k) = cldice(:ncol,k) + cldliq(:ncol,k) + enddo + + end subroutine rk_stratiform_save_qtlcwat_run + + ! Compute and save cloud water and ice particle sizes for radiation +!> \section arg_table_rk_stratiform_cloud_optical_properties_run Argument Table +!! \htmlinclude arg_table_rk_stratiform_cloud_optical_properties_run.html + subroutine rk_stratiform_cloud_optical_properties_run( & + ncol, pver, & + tmelt, & + landfrac, icefrac, snowh, landm, & + t, ps, pmid, & + rel, rei, & + errmsg, errflg) + + ! Dependency: to_be_ccppized + use cloud_optical_properties, only: cldefr + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + + real(kind_phys), intent(in) :: tmelt + real(kind_phys), intent(in) :: landfrac(:) ! land_area_fraction [fraction] + real(kind_phys), intent(in) :: icefrac(:) ! sea_ice_area_fraction [fraction] + real(kind_phys), intent(in) :: snowh(:) ! lwe_surface_snow_depth_over_land [m] + real(kind_phys), intent(in) :: landm(:) ! smoothed_land_area_fraction [fraction] + real(kind_phys), intent(in) :: t(:,:) ! air_temperature [K] + real(kind_phys), intent(in) :: ps(:) ! surface_air_pressure [Pa] + real(kind_phys), intent(in) :: pmid(:,:) ! air_pressure [Pa] + + ! Output arguments + real(kind_phys), intent(out) :: rel(:,:) ! effective_radius_of_stratiform_cloud_liquid_water_particle [um] + real(kind_phys), intent(out) :: rei(:,:) ! effective_radius_of_stratiform_cloud_ice_particle [um] + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + errmsg = '' + errflg = 0 + + call cldefr( & + ncol = ncol, & + pver = pver, & + tmelt = tmelt, & + landfrac = landfrac(:ncol), & + icefrac = icefrac(:ncol), & + snowh = snowh(:ncol), & + landm = landm(:ncol), & + t = t(:ncol,:), & + ps = ps(:ncol), & + pmid = pmid(:ncol,:), & ! below output: + rel = rel(:ncol,:), & + rei = rei(:ncol,:)) + + end subroutine rk_stratiform_cloud_optical_properties_run + +end module rk_stratiform diff --git a/schemes/rasch_kristjansson/rk_stratiform.meta b/schemes/rasch_kristjansson/rk_stratiform.meta new file mode 100644 index 00000000..f994c28b --- /dev/null +++ b/schemes/rasch_kristjansson/rk_stratiform.meta @@ -0,0 +1,883 @@ +[ccpp-table-properties] + name = rk_stratiform + type = scheme + +[ccpp-table-properties] + name = rk_stratiform_check_qtlcwat + type = scheme + +[ccpp-arg-table] + name = rk_stratiform_check_qtlcwat_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ q_wv ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = true +[ cldice ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = true +[ cldliq ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in + advected = true +[ qcwat ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_at_end_of_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = inout +[ tcwat ] + standard_name = air_temperature_at_end_of_microphysics + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = inout +[ lcwat ] + standard_name = cloud_water_mixing_ratio_wrt_moist_air_and_condensed_water_at_end_of_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = inout +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = rk_stratiform_sedimentation + type = scheme + +[ccpp-arg-table] + name = rk_stratiform_sedimentation_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ sfliq ] + standard_name = stratiform_rain_flux_at_surface_due_to_sedimentation + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snow_sed ] + standard_name = lwe_cloud_ice_sedimentation_rate_at_surface_due_to_microphysics + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ prec_sed ] + standard_name = stratiform_cloud_water_surface_flux_due_to_sedimentation + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ prec_str ] + standard_name = lwe_large_scale_precipitation_rate_at_surface + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ snow_str ] + standard_name = lwe_snow_and_cloud_ice_precipitation_rate_at_surface_due_to_microphysics + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = rk_stratiform_detrain_convective_condensate + type = scheme + +[ccpp-arg-table] + name = rk_stratiform_detrain_convective_condensate_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ dlf ] + standard_name = detrainment_of_cloud_liquid_water_wrt_moist_air_and_condensed_water_due_to_all_convection + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rliq ] + standard_name = vertically_integrated_cloud_liquid_water_tendency_due_to_all_convection_to_be_applied_later_in_time_loop + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ prec_str ] + standard_name = lwe_large_scale_precipitation_rate_at_surface + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ tend_cldliq ] + standard_name = tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = true +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = rk_stratiform_cloud_fraction_perturbation + type = scheme + +[ccpp-arg-table] + name = rk_stratiform_cloud_fraction_perturbation_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ cappa ] + standard_name = ratio_of_dry_air_gas_constant_to_specific_heat_of_dry_air_at_constant_pressure + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rair ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ tmelt ] + standard_name = freezing_point_of_water + units = K + type = real | kind = kind_phys + dimensions = () + intent = in +[ pref ] + standard_name = surface_reference_pressure + units = Pa + type = real | kind = kind_phys + dimensions = () + intent = in +[ lapse_rate ] + standard_name = reference_temperature_lapse_rate + units = K m-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ top_lev_cloudphys ] + standard_name = vertical_layer_index_of_cloud_fraction_top + units = index + type = integer + dimensions = () + intent = in +[ pmid ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ ps ] + standard_name = surface_air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ temp ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ sst ] + standard_name = sea_surface_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ q_wv ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldice ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ phis ] + standard_name = surface_geopotential + units = m2 s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ shallowcu ] + standard_name = shallow_convective_cloud_area_fraction_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ deepcu ] + standard_name = deep_convective_cloud_area_fraction_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ concld ] + standard_name = convective_cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ landfrac ] + standard_name = land_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ ocnfrac ] + standard_name = ocean_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snowh ] + standard_name = lwe_surface_snow_depth_over_land + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ cloud ] + standard_name = cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ relhum ] + standard_name = relative_humidity_for_prognostic_cloud_water_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rhu00 ] + standard_name = relative_humidity_threshold_for_prognostic_cloud_water_tbd + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = inout +[ rhdfda ] + standard_name = derivative_of_relative_humidity_wrt_cloud_fraction_tbd + units = percent + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = rk_stratiform_external_forcings + type = scheme + +[ccpp-arg-table] + name = rk_stratiform_external_forcings_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ dtime ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ q_wv ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldice ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldliq ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ qcwat ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_at_end_of_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ tcwat ] + standard_name = air_temperature_at_end_of_microphysics + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ lcwat ] + standard_name = cloud_water_mixing_ratio_wrt_moist_air_and_condensed_water_at_end_of_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ qtend ] + standard_name = tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_not_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ ttend ] + standard_name = tendency_of_air_temperature_not_due_to_microphysics_tbd + units = K s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ ltend ] + standard_name = tendency_of_cloud_water_mixing_ratio_wrt_moist_air_and_condensed_water_not_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = rk_stratiform_condensate_repartioning + type = scheme + +[ccpp-arg-table] + name = rk_stratiform_condensate_repartioning_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ dtime ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ latice ] + standard_name = latent_heat_of_fusion_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ cldice ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldliq ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ fice ] + standard_name = mass_fraction_of_ice_content_within_stratiform_cloud + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ repartht ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_precipitation_repartitioning_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ tend_cldice ] + standard_name = tendency_of_cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = true +[ tend_cldliq ] + standard_name = tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = true +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = rk_stratiform_prognostic_cloud_water_tendencies + type = scheme + +[ccpp-arg-table] + name = rk_stratiform_prognostic_cloud_water_tendencies_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ dtime ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ latvap ] + standard_name = latent_heat_of_vaporization_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ latice ] + standard_name = latent_heat_of_fusion_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ qme ] + standard_name = net_condensation_rate_due_to_microphysics + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ fice ] + standard_name = mass_fraction_of_ice_content_within_stratiform_cloud + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ evapheat ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_precipitation_evaporation_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ prfzheat ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_precipitation_freezing_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ meltheat ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_precipitation_phase_change_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ repartht ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_precipitation_repartitioning_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ evapprec ] + standard_name = precipitation_evaporation_due_to_microphysics + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ ice2pr ] + standard_name = tendency_of_cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ liq2pr ] + standard_name = tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water_due_to_microphysics_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ prec_pcw ] + standard_name = lwe_stratiform_precipitation_rate_at_surface + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snow_pcw ] + standard_name = lwe_snow_precipitation_rate_at_surface_due_to_microphysics + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ prec_str ] + standard_name = lwe_large_scale_precipitation_rate_at_surface + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ snow_str ] + standard_name = lwe_snow_and_cloud_ice_precipitation_rate_at_surface_due_to_microphysics + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ cmeheat ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_condensation_evaporation_within_stratiform_cloud_tbd + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ cmeice ] + standard_name = rate_of_condensation_evaporation_of_cloud_ice_within_stratiform_cloud_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ cmeliq ] + standard_name = rate_of_condensation_evaporation_of_cloud_liquid_water_within_stratiform_cloud_tbd + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ tend_s ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ tend_q ] + standard_name = tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = true +[ tend_cldice ] + standard_name = tendency_of_cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = true +[ tend_cldliq ] + standard_name = tendency_of_cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out + constituent = true +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = rk_stratiform_save_qtlcwat + type = scheme + +[ccpp-arg-table] + name = rk_stratiform_save_qtlcwat_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ q_wv ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldice ] + standard_name = cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cldliq ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ qcwat ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_at_end_of_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ tcwat ] + standard_name = air_temperature_at_end_of_microphysics + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ lcwat ] + standard_name = cloud_water_mixing_ratio_wrt_moist_air_and_condensed_water_at_end_of_microphysics + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-table-properties] + name = rk_stratiform_cloud_optical_properties + type = scheme + dependencies = ../../to_be_ccppized/cloud_optical_properties.F90 + +[ccpp-arg-table] + name = rk_stratiform_cloud_optical_properties_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ tmelt ] + standard_name = freezing_point_of_water + units = K + type = real | kind = kind_phys + dimensions = () + intent = in +[ landfrac ] + standard_name = land_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ icefrac ] + standard_name = sea_ice_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snowh ] + standard_name = lwe_surface_snow_depth_over_land + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ landm ] + standard_name = smoothed_land_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ ps ] + standard_name = surface_air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ pmid ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rel ] + standard_name = effective_radius_of_stratiform_cloud_liquid_water_particle + units = um + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ rei ] + standard_name = effective_radius_of_stratiform_cloud_ice_particle + units = um + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/sima_diagnostics/convect_shallow_diagnostics.F90 b/schemes/sima_diagnostics/convect_shallow_diagnostics.F90 index f4fc68d3..1504a1c9 100644 --- a/schemes/sima_diagnostics/convect_shallow_diagnostics.F90 +++ b/schemes/sima_diagnostics/convect_shallow_diagnostics.F90 @@ -1,6 +1,3 @@ -! Copyright (C) 2024-2025 National Science Foundation-National Center for Atmospheric Research -! SPDX-License-Identifier: Apache-2.0 -! ! Diagnostics for shallow convection and merged deep + shallow convection ! Haipeng Lin, December 2024 module convect_shallow_diagnostics diff --git a/schemes/zhang_mcfarlane/zm_conv_options_namelist.xml b/schemes/zhang_mcfarlane/zm_conv_options_namelist.xml index b744c1f1..001fde77 100644 --- a/schemes/zhang_mcfarlane/zm_conv_options_namelist.xml +++ b/schemes/zhang_mcfarlane/zm_conv_options_namelist.xml @@ -39,7 +39,6 @@ The number of negative buoyancy regions that are allowed before the convection top and CAPE calculations are completed. 5 - 1 1 diff --git a/test/docker/Dockerfile.musica b/test/docker/Dockerfile.musica index 53331e46..89b5f0eb 100644 --- a/test/docker/Dockerfile.musica +++ b/test/docker/Dockerfile.musica @@ -85,7 +85,7 @@ RUN cd atmospheric_physics/test \ ENV CCPP_SRC_PATH="lib/ccpp-framework/src" ENV CCPP_FORTRAN_TOOLS_PATH="lib/ccpp-framework/scripts/fortran_tools" -# Make the CCPPStandardNames available +# Make the ESMStandardNames available RUN cd atmospheric_physics/test/lib \ && git clone --depth 1 https://github.com/ESCOMP/ESMStandardNames.git ENV CCPP_STD_NAMES_PATH="lib/ESMStandardNames" diff --git a/test/test_suites/suite_rasch_kristjansson.xml b/test/test_suites/suite_rasch_kristjansson.xml new file mode 100644 index 00000000..1281df89 --- /dev/null +++ b/test/test_suites/suite_rasch_kristjansson.xml @@ -0,0 +1,96 @@ + + + + + to_be_ccppized_temporary + + + tropopause_find + + + rk_stratiform_diagnostics + + + rk_stratiform_check_qtlcwat + + + cloud_particle_sedimentation + + apply_constituent_tendencies + apply_heating_rate + qneg + geopotential_temp + rk_stratiform_sedimentation + + + rk_stratiform_detrain_convective_condensate + apply_constituent_tendencies + qneg + geopotential_temp + + + convective_cloud_cover + + compute_cloud_fraction + rk_stratiform_cloud_fraction_perturbation + + + + rk_stratiform_external_forcings + + + cloud_fraction_fice + + + prognostic_cloud_water + + + rk_stratiform_condensate_repartioning + + apply_constituent_tendencies + qneg + geopotential_temp + + + rk_stratiform_prognostic_cloud_water_tendencies + + apply_constituent_tendencies + apply_heating_rate + qneg + geopotential_temp + + + compute_cloud_fraction + + + + rk_stratiform_cloud_optical_properties + + + + rk_stratiform_save_qtlcwat + + diff --git a/to_be_ccppized/cloud_optical_properties.F90 b/to_be_ccppized/cloud_optical_properties.F90 new file mode 100644 index 00000000..340b73c9 --- /dev/null +++ b/to_be_ccppized/cloud_optical_properties.F90 @@ -0,0 +1,309 @@ +! Cloud optical properties (to-be-ccppized utility module) +! Computes liquid and ice particle size and emissivity +! Author: Byron Boville, Sept 2002 assembled from existing subroutines +module cloud_optical_properties + + use shr_kind_mod, only: r8 => shr_kind_r8 + use ppgrid, only: pcols, pver, pverp + + implicit none + private + save + + public :: cldefr, cldovrlap, cldclw, reitab, reltab + public :: cldems_rk, cldems + +contains + + ! Compute cloud water and ice particle size [um] + ! using empirical formulas to construct effective radii + ! Original author: J.T. Kiehl, B.A. Boville, P. Rasch + subroutine cldefr( & + ncol, pver, & + tmelt, & + landfrac, t, rel, rei, ps, pmid, landm, icefrac, snowh) + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver + real(r8), intent(in) :: tmelt + + real(r8), intent(in) :: landfrac(:) ! Land fraction + real(r8), intent(in) :: icefrac(:) ! Ice fraction + real(r8), intent(in) :: t(:, :) ! Temperature + real(r8), intent(in) :: ps(:) ! Surface pressure + real(r8), intent(in) :: pmid(:, :) ! Midpoint pressures + real(r8), intent(in) :: landm(:) + real(r8), intent(in) :: snowh(:) ! Snow depth over land, water equivalent [m] + + ! Output arguments + real(r8), intent(out) :: rel(:, :) ! Liquid effective drop size [um] + real(r8), intent(out) :: rei(:, :) ! Ice effective drop size [um] + + ! following Kiehl + call reltab(ncol, pver, tmelt, t(:ncol,:), landfrac(:ncol), landm(:ncol), icefrac(:ncol), snowh(:), rel(:ncol,:)) + + ! following Kristjansson and Mitchell + call reitab(ncol, pver, t(:ncol,:), rei(:ncol,:)) + end subroutine cldefr + + ! Compute cloud emissivity using cloud liquid water path [g m-2] + ! Original author: J.T. Kiehl + ! + ! Variant 1 used for RK microphysics + subroutine cldems_rk(ncol, pver, clwp, fice, rei, emis, cldtau) + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of vertical levels + real(r8), intent(in) :: clwp(pcols, pver) ! cloud liquid water path (g/m**2) + real(r8), intent(in) :: rei(pcols, pver) ! ice effective drop size (microns) + real(r8), intent(in) :: fice(pcols, pver) ! fractional ice content within cloud + + real(r8), intent(out) :: emis(pcols, pver) ! cloud emissivity (fraction) + real(r8), intent(out) :: cldtau(pcols, pver) ! cloud optical depth + + integer :: i, k ! longitude, level indices + real(r8) :: kabs ! longwave absorption coeff (m**2/g) + real(r8) :: kabsi ! ice absorption coefficient + real(r8) :: kabsl ! longwave liquid absorption coeff (m**2/g) + parameter(kabsl=0.090361_r8) + + do k = 1, pver + do i = 1, ncol + ! note that optical properties for ice valid only + ! in range of 13 > rei > 130 micron (Ebert and Curry 92) + kabsi = 0.005_r8 + 1._r8/rei(i, k) + kabs = kabsl*(1._r8 - fice(i, k)) + kabsi*fice(i, k) + emis(i, k) = 1._r8 - exp(-1.66_r8*kabs*clwp(i, k)) + cldtau(i, k) = kabs*clwp(i, k) + end do + end do + end subroutine cldems_rk + + ! Variant 2 used for other microphysical schemes + subroutine cldems(ncol, pver, clwp, fice, rei, emis, cldtau) + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of vertical levels + real(r8), intent(in) :: clwp(pcols, pver) ! cloud liquid water path (g/m**2) + real(r8), intent(in) :: rei(pcols, pver) ! ice effective drop size (microns) + real(r8), intent(in) :: fice(pcols, pver) ! fractional ice content within cloud + + real(r8), intent(out) :: emis(pcols, pver) ! cloud emissivity (fraction) + real(r8), intent(out) :: cldtau(pcols, pver) ! cloud optical depth + + integer :: i, k ! longitude, level indices + real(r8) :: kabs ! longwave absorption coeff (m**2/g) + real(r8) :: kabsi ! ice absorption coefficient + real(r8) :: kabsl ! longwave liquid absorption coeff (m**2/g) + parameter(kabsl=0.090361_r8) + + do k = 1, pver + do i = 1, ncol + ! note that optical properties for ice valid only + ! in range of 13 > rei > 130 micron (Ebert and Curry 92) + kabsi = 0.005_r8 + 1._r8/min(max(13._r8, rei(i, k)), 130._r8) + kabs = kabsl*(1._r8 - fice(i, k)) + kabsi*fice(i, k) + emis(i, k) = 1._r8 - exp(-1.66_r8*kabs*clwp(i, k)) + cldtau(i, k) = kabs*clwp(i, k) + end do + end do + end subroutine cldems + + ! Partitions each column into regions with clouds in neighboring layers. + ! This information is used to implement maximum overlap in these regions + ! with random overlap between them. + ! On output, + ! nmxrgn contains the number of regions in each column + ! pmxrgn contains the interface pressures for the lower boundaries of + ! each region! + subroutine cldovrlap(ncol, pver, pverp, pint, cld, nmxrgn, pmxrgn) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: pverp + + real(r8), intent(in) :: pint(:, :) ! Interface pressure + real(r8), intent(in) :: cld(:, :) ! Fractional cloud cover + + ! Output arguments + integer, intent(out) :: nmxrgn(:) ! Number of maximally overlapped regions + real(r8), intent(out) :: pmxrgn(:, :) ! Maximum values of pressure for each + ! maximally overlapped region. + ! 0->pmxrgn(i,1) is range of pressure for + ! 1st region,pmxrgn(i,1)->pmxrgn(i,2) for + ! 2nd region, etc + ! (ncol, pverp) + + integer :: i, k + integer :: n ! Max-overlap region counter + real(r8) :: pnm(ncol, pverp) ! Interface pressure + logical :: cld_found ! Flag for detection of cloud + logical :: cld_layer(pver) ! Flag for cloud in layer + + do i = 1, ncol + cld_found = .false. + cld_layer(:) = cld(i, :) > 0.0_r8 + pmxrgn(i, :) = 0.0_r8 + pnm(i, :) = pint(i, :)*10._r8 + n = 1 + do k = 1, pver + if (cld_layer(k) .and. .not. cld_found) then + cld_found = .true. + else if (.not. cld_layer(k) .and. cld_found) then + cld_found = .false. + if (count(cld_layer(k:pver)) == 0) then + exit + end if + pmxrgn(i, n) = pnm(i, k) + n = n + 1 + end if + end do + pmxrgn(i, n) = pnm(i, pverp) + nmxrgn(i) = n + end do + end subroutine cldovrlap + + ! Evaluate cloud liquid water path clwp [g m-2] + ! Original author: Author: J.T. Kiehl + subroutine cldclw(ncol, zi, clwp, tpw, hl) + + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + + real(r8), intent(in) :: zi(pcols, pverp) ! height at layer interfaces(m) + real(r8), intent(in) :: tpw(pcols) ! total precipitable water (mm) + + ! Output arguments + real(r8), intent(out) :: clwp(pcols, pver) ! cloud liquid water path (g/m**2) + real(r8), intent(out) :: hl(pcols) ! liquid water scale height + + integer :: i, k ! longitude, level indices + real(r8) :: clwc0 ! reference liquid water concentration (g/m**3) + real(r8) :: emziohl(pcols, pverp) ! exp(-zi/hl) + real(r8) :: rhl(pcols) ! 1/hl + + ! Set reference liquid water concentration + clwc0 = 0.21_r8 + + ! Diagnose liquid water scale height from precipitable water + do i = 1, ncol + hl(i) = 700.0_r8*log(max(tpw(i) + 1.0_r8, 1.0_r8)) + rhl(i) = 1.0_r8/hl(i) + end do + + ! Evaluate cloud liquid water path (vertical integral of exponential fn) + do k = 1, pverp + do i = 1, ncol + emziohl(i, k) = exp(-zi(i, k)*rhl(i)) + end do + end do + do k = 1, pver + do i = 1, ncol + clwp(i, k) = clwc0*hl(i)*(emziohl(i, k + 1) - emziohl(i, k)) + end do + end do + end subroutine cldclw + + + ! Compute cloud water size + ! analytic formula following the formulation originally developed by J. T. Kiehl + ! Author: Phil Rasch + subroutine reltab(ncol, pver, tmelt, t, landfrac, landm, icefrac, snowh, rel) + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(r8), intent(in) :: tmelt + real(r8), intent(in) :: landfrac(:) ! Land fraction + real(r8), intent(in) :: landm(:) ! Land fraction ramping to zero over ocean + real(r8), intent(in) :: icefrac(:) ! Ice fraction + real(r8), intent(in) :: t(:, :) ! Temperature [K] + real(r8), intent(in) :: snowh(:) ! Snow depth over land, water equivalent [m] + + ! Output arguments + real(r8), intent(out) :: rel(:, :) ! Liquid effective drop size (microns) + + integer i, k + real(r8) :: rliqland ! liquid drop size if over land + real(r8) :: rliqocean ! liquid drop size if over ocean + real(r8) :: rliqice ! liquid drop size if over sea ice + + rliqocean = 14.0_r8 + rliqice = 14.0_r8 + rliqland = 8.0_r8 + + do k = 1, pver + do i = 1, ncol + ! jrm Reworked effective radius algorithm + ! Start with temperature-dependent value appropriate for continental air + ! Note: findmcnew has a pressure dependence here + rel(i, k) = rliqland + (rliqocean - rliqland)*min(1.0_r8, max(0.0_r8, (tmelt - t(i, k))*0.05_r8)) + ! Modify for snow depth over land + rel(i, k) = rel(i, k) + (rliqocean - rel(i, k))*min(1.0_r8, max(0.0_r8, snowh(i)*10._r8)) + ! Ramp between polluted value over land to clean value over ocean. + rel(i, k) = rel(i, k) + (rliqocean - rel(i, k))*min(1.0_r8, max(0.0_r8, 1.0_r8 - landm(i))) + ! Ramp between the resultant value and a sea ice value in the presence of ice. + rel(i, k) = rel(i, k) + (rliqice - rel(i, k))*min(1.0_r8, max(0.0_r8, icefrac(i))) + ! end jrm + end do + end do + end subroutine reltab + + subroutine reitab(ncol, pver, t, re) + + integer, intent(in) :: ncol + integer, intent(in) :: pver + real(r8), intent(in) :: t(:, :) + real(r8), intent(out) :: re(:, :) + integer, parameter :: len_retab = 138 + real(r8), parameter :: min_retab = 136._r8 + real(r8) :: retab(len_retab) + real(r8) :: corr + integer :: i + integer :: k + integer :: index + ! + ! Tabulated values of re(T) in the temperature interval + ! 180 K -- 274 K; hexagonal columns assumed: + ! + ! Modified for pmc formation: 136K -- 274K + ! + data retab / & + 0.05_r8, 0.05_r8, 0.05_r8, 0.05_r8, 0.05_r8, 0.05_r8, & + 0.055_r8, 0.06_r8, 0.07_r8, 0.08_r8, 0.09_r8, 0.1_r8, & + 0.2_r8, 0.3_r8, 0.40_r8, 0.50_r8, 0.60_r8, 0.70_r8, & + 0.8_r8, 0.9_r8, 1.0_r8, 1.1_r8, 1.2_r8, 1.3_r8, & + 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 2.0_r8, 2.2_r8, & + 2.4_r8, 2.6_r8, 2.8_r8, 3.0_r8, 3.2_r8, 3.5_r8, & + 3.8_r8, 4.1_r8, 4.4_r8, 4.7_r8, 5.0_r8, 5.3_r8, & + 5.6_r8, & + 5.92779_r8, 6.26422_r8, 6.61973_r8, 6.99539_r8, 7.39234_r8, & + 7.81177_r8, 8.25496_r8, 8.72323_r8, 9.21800_r8, 9.74075_r8, 10.2930_r8, & + 10.8765_r8, 11.4929_r8, 12.1440_r8, 12.8317_r8, 13.5581_r8, 14.2319_r8, & + 15.0351_r8, 15.8799_r8, 16.7674_r8, 17.6986_r8, 18.6744_r8, 19.6955_r8, & + 20.7623_r8, 21.8757_r8, 23.0364_r8, 24.2452_r8, 25.5034_r8, 26.8125_r8, & + 27.7895_r8, 28.6450_r8, 29.4167_r8, 30.1088_r8, 30.7306_r8, 31.2943_r8, & + 31.8151_r8, 32.3077_r8, 32.7870_r8, 33.2657_r8, 33.7540_r8, 34.2601_r8, & + 34.7892_r8, 35.3442_r8, 35.9255_r8, 36.5316_r8, 37.1602_r8, 37.8078_r8, & + 38.4720_r8, 39.1508_r8, 39.8442_r8, 40.5552_r8, 41.2912_r8, 42.0635_r8, & + 42.8876_r8, 43.7863_r8, 44.7853_r8, 45.9170_r8, 47.2165_r8, 48.7221_r8, & + 50.4710_r8, 52.4980_r8, 54.8315_r8, 57.4898_r8, 60.4785_r8, 63.7898_r8, & + 65.5604_r8, 71.2885_r8, 75.4113_r8, 79.7368_r8, 84.2351_r8, 88.8833_r8, & + 93.6658_r8, 98.5739_r8, 103.603_r8, 108.752_r8, 114.025_r8, 119.424_r8, & + 124.954_r8, 130.630_r8, 136.457_r8, 142.446_r8, 148.608_r8, 154.956_r8, & + 161.503_r8, 168.262_r8, 175.248_r8, 182.473_r8, 189.952_r8, 197.699_r8, & + 205.728_r8, 214.055_r8, 222.694_r8, 231.661_r8, 240.971_r8, 250.639_r8/ + save retab + + do k = 1, pver + do i = 1, ncol + index = int(t(i, k) - min_retab) + index = min(max(index, 1), len_retab - 1) + corr = t(i, k) - int(t(i, k)) + re(i, k) = retab(index)*(1._r8 - corr) & + + retab(index + 1)*corr + ! re(i,k) = amax1(amin1(re(i,k),30.),10.) + end do + end do + end subroutine reitab + +end module cloud_optical_properties