diff --git a/physics/MP/GFDL/v1/fv_sat_adj.F90 b/physics/MP/GFDL/v1/fv_sat_adj.F90 new file mode 100644 index 000000000..53543485b --- /dev/null +++ b/physics/MP/GFDL/v1/fv_sat_adj.F90 @@ -0,0 +1,1431 @@ +!>\file fv_sat_adj.F90 +!! This file contains the GFDL in-core fast saturation adjustment. +!! and it is an "intermediate physics" implemented in the remapping Lagrangian to +!! Eulerian loop of FV3 solver. +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the GFDL Cloud Microphysics. +!* +!* The GFDL Cloud Microphysics is free software: you can +!8 redistribute it and/or modify it under the terms of the +!* GNU Lesser General Public License as published by the +!* Free Software Foundation, either version 3 of the License, or +!* (at your option) any later version. +!* +!* The GFDL Cloud Microphysics is distributed in the hope it will be +!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty +!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +!* See the GNU General Public License for more details. +!* +!* You should have received a copy of the GNU Lesser General Public +!* License along with the GFDL Cloud Microphysics. +!* If not, see . +!*********************************************************************** + +!> This module contains the GFDL in-core fast saturation adjustment +!! called in FV3 dynamics solver. +module fv_sat_adj +! Modules Included: +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +! +!
Module NameFunctions Included
constants_modrvgas, rdgas, grav, hlv, hlf, cp_air
fv_arrays_mod r_grid
fv_mp_modis_master
gfdl_cloud_microphys_modql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt, +! tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r, +! rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs
+ ! DH* TODO - MAKE THIS INPUT ARGUMENTS *DH + use physcons, only : rdgas => con_rd_dyn, & + rvgas => con_rv_dyn, & + grav => con_g_dyn, & + hlv => con_hvap_dyn, & + hlf => con_hfus_dyn, & + cp_air => con_cp_dyn + ! *DH + use machine, only: kind_grid, kind_dyn + use gfdl_cloud_microphys_mod, only: ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt + use gfdl_cloud_microphys_mod, only: icloud_f, sat_adj0, t_sub, cld_min + use gfdl_cloud_microphys_mod, only: tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r + use gfdl_cloud_microphys_mod, only: rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs +#ifdef MULTI_GASES + use ccpp_multi_gases_mod, only: multi_gases_init, & + multi_gases_finalize, & + virq_qpz, vicpqd_qpz, & + vicvqd_qpz, num_gas +#endif + + implicit none + + private + + public fv_sat_adj_init, fv_sat_adj_run, fv_sat_adj_finalize + + logical :: is_initialized = .false. + + real(kind=kind_dyn), parameter :: rrg = -rdgas/grav + ! real, parameter :: cp_air = cp_air ! 1004.6, heat capacity of dry air at constant pressure, come from constants_mod + real(kind=kind_dyn), parameter :: cp_vap = 4.0 * rvgas !< 1846.0, heat capacity of water vapor at constant pressure + real(kind=kind_dyn), parameter :: cv_air = cp_air - rdgas !< 717.55, heat capacity of dry air at constant volume + real(kind=kind_dyn), parameter :: cv_vap = 3.0 * rvgas !< 1384.5, heat capacity of water vapor at constant volume + ! http: / / www.engineeringtoolbox.com / ice - thermal - properties - d_576.html + ! c_ice = 2050.0 at 0 deg c + ! c_ice = 1972.0 at - 15 deg c + ! c_ice = 1818.0 at - 40 deg c + ! http: / / www.engineeringtoolbox.com / water - thermal - properties - d_162.html + ! c_liq = 4205.0 at 4 deg c + ! c_liq = 4185.5 at 15 deg c + ! c_liq = 4178.0 at 30 deg c + ! real, parameter :: c_ice = 2106.0 ! ifs: heat capacity of ice at 0 deg c + ! real, parameter :: c_liq = 4218.0 ! ifs: heat capacity of liquid at 0 deg c + real(kind=kind_dyn), parameter :: c_ice = 1972.0 !< gfdl: heat capacity of ice at - 15 deg c + real(kind=kind_dyn), parameter :: c_liq = 4185.5 !< gfdl: heat capacity of liquid at 15 deg c + real(kind=kind_dyn), parameter :: dc_vap = cp_vap - c_liq !< - 2339.5, isobaric heating / cooling + real(kind=kind_dyn), parameter :: dc_ice = c_liq - c_ice !< 2213.5, isobaric heating / colling + real(kind=kind_dyn), parameter :: tice = 273.16 !< freezing temperature + real(kind=kind_dyn), parameter :: t_wfr = tice - 40. !< homogeneous freezing temperature + real(kind=kind_dyn), parameter :: lv0 = hlv - dc_vap * tice !< 3.13905782e6, evaporation latent heat coefficient at 0 deg k + real(kind=kind_dyn), parameter :: li00 = hlf - dc_ice * tice !< - 2.7105966e5, fusion latent heat coefficient at 0 deg k + ! real (kind_grid), parameter :: e00 = 610.71 ! gfdl: saturation vapor pressure at 0 deg c + real (kind_grid), parameter :: e00 = 611.21 !< ifs: saturation vapor pressure at 0 deg c + real (kind_grid), parameter :: d2ice = dc_vap + dc_ice !< - 126, isobaric heating / cooling + real (kind_grid), parameter :: li2 = lv0 + li00 !< 2.86799816e6, sublimation latent heat coefficient at 0 deg k + real(kind=kind_dyn), parameter :: lat2 = (hlv + hlf) ** 2 !< used in bigg mechanism + real(kind=kind_dyn) :: d0_vap !< the same as dc_vap, except that cp_vap can be cp_vap or cv_vap + real(kind=kind_dyn) :: lv00 !< the same as lv0, except that cp_vap can be cp_vap or cv_vap + real(kind=kind_dyn), allocatable :: table (:), table2 (:), tablew (:), des2 (:), desw (:) + +contains + +!>\brief The subroutine 'fv_sat_adj_init' initializes lookup tables for the saturation mixing ratio. +!! \section arg_table_fv_sat_adj_init Argument Table +!! \htmlinclude fv_sat_adj_init.html +!! +subroutine fv_sat_adj_init(do_sat_adj, kmp, nwat, ngas, rilist, cpilist, & + mpirank, mpiroot, errmsg, errflg) + + implicit none + + ! Interface variables + logical, intent(in ) :: do_sat_adj + integer, intent(in ) :: kmp + integer, intent(in ) :: nwat + integer, intent(in ) :: ngas + real(kind_dyn), intent(in ) :: rilist(0:ngas) + real(kind_dyn), intent(in ) :: cpilist(0:ngas) + integer, intent(in ) :: mpirank + integer, intent(in ) :: mpiroot + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Local variables + integer, parameter :: length = 2621 + integer :: i + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + ! If saturation adjustment is not used, return immediately + if (.not.do_sat_adj) then + write(errmsg,'(a)') 'Logic error: fv_sat_adj_init is called but do_sat_adj is set to false' + errflg = 1 + return + end if + + if (.not.nwat==6) then + write(errmsg,'(a)') 'Logic error: fv_sat_adj requires six water species (nwat=6)' + errflg = 1 + return + end if + + if (is_initialized) return + + ! generate es table (dt = 0.1 deg c) + + allocate (table (length)) + allocate (table2 (length)) + allocate (tablew (length)) + allocate (des2 (length)) + allocate (desw (length)) + + call qs_table (length) + call qs_table2 (length) + call qs_tablew (length) + + do i = 1, length - 1 + des2 (i) = max (0., table2 (i + 1) - table2 (i)) + desw (i) = max (0., tablew (i + 1) - tablew (i)) + enddo + des2 (length) = des2 (length - 1) + desw (length) = desw (length - 1) + +#ifdef MULTI_GASES + call multi_gases_init(ngas,nwat,rilist,cpilist,mpirank==mpiroot) +#endif + + is_initialized = .true. + +end subroutine fv_sat_adj_init + +!\ingroup fast_sat_adj +!>\brief The subroutine 'fv_sat_adj_finalize' deallocates lookup tables for the saturation mixing ratio. +!! \section arg_table_fv_sat_adj_finalize Argument Table +!! \htmlinclude fv_sat_adj_finalize.html +!! +subroutine fv_sat_adj_finalize (errmsg, errflg) + + implicit none + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + if (.not.is_initialized) return + + if (allocated(table )) deallocate(table ) + if (allocated(table2)) deallocate(table2) + if (allocated(tablew)) deallocate(tablew) + if (allocated(des2 )) deallocate(des2 ) + if (allocated(desw )) deallocate(desw ) + +#ifdef MULTI_GASES + call multi_gases_finalize() +#endif + + is_initialized = .false. + +end subroutine fv_sat_adj_finalize + +!>\defgroup fast_sat_adj GFDL In-Core Fast Saturation Adjustment Module +!> @{ +!! The subroutine 'fv_sat_adj' implements the fast processes in the GFDL +!! Cloud MP. It is part of the GFDL Cloud MP. +!>\author Shian-Jiann Lin, Linjiong Zhou +!! +!>\brief The subroutine 'fv_sat_adj' performs the fast processes in the GFDL microphysics. +!>\details This is designed for single-moment 6-class cloud microphysics schemes. +!! It handles the heat release due to in situ phase changes. +!! +!! \section arg_table_fv_sat_adj_run Argument Table +!! \htmlinclude fv_sat_adj_run.html +!! +subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, jsd, jed, & + ng, hydrostatic, fast_mp_consv, te0_2d, te0, ngas, qvi, qv, ql, qi, qr, & + qs, qg, hs, peln, delz, delp, pt, pkz, q_con, akap, cappa, area, dtdt, & + out_dt, last_step, do_qa, qa, & + nthreads, errmsg, errflg) + + implicit none + + ! Interface variables + real(kind=kind_dyn), intent(in) :: mdt + real(kind=kind_dyn), intent(in) :: zvir + integer, intent(in) :: is + integer, intent(in) :: ie + integer, intent(in) :: isd + integer, intent(in) :: ied + integer, intent(in) :: kmp + integer, intent(in) :: km + integer, intent(in) :: kmdelz + integer, intent(in) :: js + integer, intent(in) :: je + integer, intent(in) :: jsd + integer, intent(in) :: jed + integer, intent(in) :: ng + logical, intent(in) :: hydrostatic + logical, intent(in) :: fast_mp_consv + real(kind=kind_dyn), intent(inout) :: te0_2d(is:ie, js:je) + real(kind=kind_dyn), intent( out) :: te0(isd:ied, jsd:jed, 1:km) + ! If multi-gases physics are not used, ngas is one and qvi identical to qv + integer, intent(in) :: ngas + real(kind=kind_dyn), intent(inout) :: qvi(isd:ied, jsd:jed, 1:km, 1:ngas) + real(kind=kind_dyn), intent(inout) :: qv(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: ql(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: qi(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: qr(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: qs(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: qg(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(in) :: hs(isd:ied, jsd:jed) + real(kind=kind_dyn), intent(in) :: peln(is:ie, 1:km+1, js:je) + ! For hydrostatic build, kmdelz=1, otherwise kmdelz=km (see fv_arrays.F90) + real(kind=kind_dyn), intent(in) :: delz(is:ie, js:je, 1:kmdelz) + real(kind=kind_dyn), intent(in) :: delp(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: pt(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: pkz(is:ie, js:je, 1:km) +#ifdef USE_COND + real(kind=kind_dyn), intent(inout) :: q_con(isd:ied, jsd:jed, 1:km) +#else + real(kind=kind_dyn), intent(inout) :: q_con(isd:isd, jsd:jsd, 1) +#endif + real(kind=kind_dyn), intent(in) :: akap +#ifdef MOIST_CAPPA + real(kind=kind_dyn), intent(inout) :: cappa(isd:ied, jsd:jed, 1:km) +#else + real(kind=kind_dyn), intent(inout) :: cappa(isd:ied, jsd:jed, 1) +#endif + ! DH* WARNING, allocation in fv_arrays.F90 is area(isd_2d:ied_2d, jsd_2d:jed_2d), + ! where normally isd_2d = isd etc, but for memory optimization, these can be set + ! to isd_2d=0, ied_2d=-1 etc. I don't believe this optimization is actually used, + ! as it would break a whole lot of code (including the one below)! + ! Assume thus that isd_2d = isd etc. + real(kind_grid), intent(in) :: area(isd:ied, jsd:jed) + real(kind=kind_dyn), intent(inout) :: dtdt(is:ie, js:je, 1:km) + logical, intent(in) :: out_dt + logical, intent(in) :: last_step + logical, intent(in) :: do_qa + real(kind=kind_dyn), intent( out) :: qa(isd:ied, jsd:jed, 1:km) + integer, intent(in) :: nthreads + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Local variables + real(kind=kind_dyn), dimension(is:ie,js:je) :: dpln + integer :: kdelz + integer :: k, j, i + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + +#ifndef FV3 +! Open parallel region if not already opened by host model +!$OMP parallel num_threads(nthreads) default(none) & +!$OMP shared(kmp,km,js,je,is,ie,peln,mdt, & +!$OMP isd,jsd,delz,q_con,cappa,qa, & +!$OMP do_qa,last_step,out_dt,dtdt, & +!$OMP area,delp,pt,hs,qg,qs,qr,qi, & +!$OMP ql,qv,te0,fast_mp_consv, & +!$OMP hydrostatic,ng,zvir,pkz, & +!$OMP akap,te0_2d,ngas,qvi) & +!$OMP private(k,j,i,kdelz,dpln) +#endif + +!$OMP do + do k=kmp,km + do j=js,je + do i=is,ie + dpln(i,j) = peln(i,k+1,j) - peln(i,k,j) + enddo + enddo + if (hydrostatic) then + kdelz = 1 + else + kdelz = k + end if + call fv_sat_adj_work(abs(mdt), zvir, is, ie, js, je, ng, hydrostatic, fast_mp_consv, & + te0(isd,jsd,k), & +#ifdef MULTI_GASES + qvi(isd,jsd,k,1:ngas), & +#else + qv(isd,jsd,k), & +#endif + ql(isd,jsd,k), qi(isd,jsd,k), & + qr(isd,jsd,k), qs(isd,jsd,k), qg(isd,jsd,k), & + hs, dpln, delz(is:,js:,kdelz), pt(isd,jsd,k), delp(isd,jsd,k),& + q_con(isd:,jsd:,k), cappa(isd:,jsd:,k), area, dtdt(is,js,k), & + out_dt, last_step, do_qa, qa(isd,jsd,k)) + if ( .not. hydrostatic ) then + do j=js,je + do i=is,ie +#ifdef MOIST_CAPPA + pkz(i,j,k) = exp(cappa(i,j,k)*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k))) +#else +#ifdef MULTI_GASES + pkz(i,j,k) = exp(akap*(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k))) +#else + pkz(i,j,k) = exp(akap*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k))) +#endif +#endif + enddo + enddo + endif + enddo +!$OMP end do + + if ( fast_mp_consv ) then +!$OMP do + do j=js,je + do i=is,ie + do k=kmp,km + te0_2d(i,j) = te0_2d(i,j) + te0(i,j,k) + enddo + enddo + enddo +!$OMP end do + endif + +#ifndef FV3 +!$OMP end parallel +#endif + + return + +end subroutine fv_sat_adj_run + +!>\ingroup fast_sat_adj +!> This subroutine includes the entity of the fast saturation adjustment processes. +!>\section fast_gen GFDL Cloud Fast Physics General Algorithm +!> @{ +subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, te0, & +#ifdef MULTI_GASES + qvi, & +#else + qv, & +#endif + ql, qi, qr, qs, qg, hs, dpln, delz, pt, dp, q_con, cappa, & + area, dtdt, out_dt, last_step, do_qa, qa) + + implicit none + + ! Interface variables + integer, intent (in) :: is, ie, js, je, ng + logical, intent (in) :: hydrostatic, consv_te, out_dt, last_step, do_qa + real(kind=kind_dyn), intent (in) :: zvir, mdt ! remapping time step + real(kind=kind_dyn), intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: dp, hs + real(kind=kind_dyn), intent (in), dimension (is:ie, js:je) :: dpln, delz + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: pt +#ifdef MULTI_GASES + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng, 1:1, 1:num_gas) :: qvi +#else + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: qv +#endif + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: ql, qi, qr, qs, qg + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: q_con, cappa + real(kind=kind_dyn), intent (inout), dimension (is:ie, js:je) :: dtdt + real(kind=kind_dyn), intent (out), dimension (is - ng:ie + ng, js - ng:je + ng) :: qa, te0 + real (kind_grid), intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: area + + ! Local variables +#ifdef MULTI_GASES + real, dimension (is - ng:ie + ng, js - ng:je + ng) :: qv +#endif + real(kind=kind_dyn), dimension (is:ie) :: wqsat, dq2dt, qpz, cvm, t0, pt1, qstar + real(kind=kind_dyn), dimension (is:ie) :: icp2, lcp2, tcp2, tcp3 + real(kind=kind_dyn), dimension (is:ie) :: den, q_liq, q_sol, q_cond, src, sink, hvar + real(kind=kind_dyn), dimension (is:ie) :: mc_air, lhl, lhi + real(kind=kind_dyn) :: qsw, rh + real(kind=kind_dyn) :: tc, qsi, dqsdt, dq, dq0, pidep, qi_crt, tmp, dtmp + real(kind=kind_dyn) :: tin, rqi, q_plus, q_minus + real(kind=kind_dyn) :: sdt, dt_bigg, adj_fac + real(kind=kind_dyn) :: fac_smlt, fac_r2g, fac_i2s, fac_imlt, fac_l2r, fac_v2l, fac_l2v + real(kind=kind_dyn) :: factor, qim, tice0, c_air, c_vap, dw + integer :: i, j + +#ifdef MULTI_GASES + qv(:,:) = qvi(:,:,1,1) +#endif + sdt = 0.5 * mdt ! half remapping time step + dt_bigg = mdt ! bigg mechinism time step + + tice0 = tice - 0.01 ! 273.15, standard freezing temperature + + ! ----------------------------------------------------------------------- + !> - Define conversion scalar / factor. + ! ----------------------------------------------------------------------- + + fac_i2s = 1. - exp (- mdt / tau_i2s) + fac_v2l = 1. - exp (- sdt / tau_v2l) + fac_r2g = 1. - exp (- mdt / tau_r2g) + fac_l2r = 1. - exp (- mdt / tau_l2r) + + fac_l2v = 1. - exp (- sdt / tau_l2v) + fac_l2v = min (sat_adj0, fac_l2v) + + fac_imlt = 1. - exp (- sdt / tau_imlt) + fac_smlt = 1. - exp (- mdt / tau_smlt) + + ! ----------------------------------------------------------------------- + !> - Define heat capacity of dry air and water vapor based on hydrostatical property. + ! ----------------------------------------------------------------------- + + if (hydrostatic) then + c_air = cp_air + c_vap = cp_vap + else + c_air = cv_air + c_vap = cv_vap + endif + d0_vap = c_vap - c_liq + lv00 = hlv - d0_vap * tice + ! dc_vap = cp_vap - c_liq ! - 2339.5 + ! d0_vap = cv_vap - c_liq ! - 2801.0 + + do j = js, je ! start j loop + + do i = is, ie + q_liq (i) = ql (i, j) + qr (i, j) + q_sol (i) = qi (i, j) + qs (i, j) + qg (i, j) + qpz (i) = q_liq (i) + q_sol (i) +#ifdef MULTI_GASES + pt1 (i) = pt (i, j) / virq_qpz(qvi(i,j,1,1:num_gas),qpz(i)) +#else +#ifdef USE_COND + pt1 (i) = pt (i, j) / ((1 + zvir * qv (i, j)) * (1 - qpz (i))) +#else + pt1 (i) = pt (i, j) / (1 + zvir * qv (i, j)) +#endif +#endif + t0 (i) = pt1 (i) ! true temperature + qpz (i) = qpz (i) + qv (i, j) ! total_wat conserved in this routine + enddo + + ! ----------------------------------------------------------------------- + !> - Define air density based on hydrostatical property. + ! ----------------------------------------------------------------------- + + if (hydrostatic) then + do i = is, ie + den (i) = dp (i, j) / (dpln (i, j) * rdgas * pt (i, j)) + enddo + else + do i = is, ie + den (i) = - dp (i, j) / (grav * delz (i, j)) ! moist_air density + enddo + endif + + ! ----------------------------------------------------------------------- + !> - Define heat capacity and latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie +#ifdef MULTI_GASES + if (hydrostatic) then + c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i)) + else + c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i)) + endif +#endif + mc_air (i) = (1. - qpz (i)) * c_air ! constant + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + lhi (i) = li00 + dc_ice * pt1 (i) + icp2 (i) = lhi (i) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Fix energy conservation. + ! ----------------------------------------------------------------------- + + if (consv_te) then + if (hydrostatic) then + do i = is, ie +#ifdef MULTI_GASES + c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i)) +#endif + te0 (i, j) = - c_air * t0 (i) + enddo + else + do i = is, ie +#ifdef USE_COND + te0 (i, j) = - cvm (i) * t0 (i) +#else +#ifdef MULTI_GASES + c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i)) +#endif + te0 (i, j) = - c_air * t0 (i) +#endif + enddo + endif + endif + + ! ----------------------------------------------------------------------- + !> - Fix negative cloud ice with snow. + ! ----------------------------------------------------------------------- + + do i = is, ie + if (qi (i, j) < 0.) then + qs (i, j) = qs (i, j) + qi (i, j) + qi (i, j) = 0. + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Melting of cloud ice to cloud water and rain. + ! ----------------------------------------------------------------------- + + do i = is, ie + if (qi (i, j) > 1.e-8 .and. pt1 (i) > tice) then + sink (i) = min (qi (i, j), fac_imlt * (pt1 (i) - tice) / icp2 (i)) + qi (i, j) = qi (i, j) - sink (i) + ! sjl, may 17, 2017 + ! tmp = min (sink (i), dim (ql_mlt, ql (i, j))) ! max ql amount + ! ql (i, j) = ql (i, j) + tmp + ! qr (i, j) = qr (i, j) + sink (i) - tmp + ! sjl, may 17, 2017 + ql (i, j) = ql (i, j) + sink (i) + q_liq (i) = q_liq (i) + sink (i) + q_sol (i) = q_sol (i) - sink (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) - sink (i) * lhi (i) / cvm (i) + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhi (i) = li00 + dc_ice * pt1 (i) + icp2 (i) = lhi (i) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Fix negative snow with graupel or graupel with available snow. + ! ----------------------------------------------------------------------- + + do i = is, ie + if (qs (i, j) < 0.) then + qg (i, j) = qg (i, j) + qs (i, j) + qs (i, j) = 0. + elseif (qg (i, j) < 0.) then + tmp = min (- qg (i, j), max (0., qs (i, j))) + qg (i, j) = qg (i, j) + tmp + qs (i, j) = qs (i, j) - tmp + endif + enddo + + ! after this point cloud ice & snow are positive definite + + ! ----------------------------------------------------------------------- + !> - Fix negative cloud water with rain or rain with available cloud water. + ! ----------------------------------------------------------------------- + + do i = is, ie + if (ql (i, j) < 0.) then + tmp = min (- ql (i, j), max (0., qr (i, j))) + ql (i, j) = ql (i, j) + tmp + qr (i, j) = qr (i, j) - tmp + elseif (qr (i, j) < 0.) then + tmp = min (- qr (i, j), max (0., ql (i, j))) + ql (i, j) = ql (i, j) - tmp + qr (i, j) = qr (i, j) + tmp + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Enforce complete freezing of cloud water to cloud ice below - 48 c. + ! ----------------------------------------------------------------------- + + do i = is, ie + dtmp = tice - 48. - pt1 (i) + if (ql (i, j) > 0. .and. dtmp > 0.) then + sink (i) = min (ql (i, j), dtmp / icp2 (i)) + ql (i, j) = ql (i, j) - sink (i) + qi (i, j) = qi (i, j) + sink (i) + q_liq (i) = q_liq (i) - sink (i) + q_sol (i) = q_sol (i) + sink (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i) + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhl (i) = lv00 + d0_vap * pt1 (i) + lhi (i) = li00 + dc_ice * pt1 (i) + lcp2 (i) = lhl (i) / cvm (i) + icp2 (i) = lhi (i) / cvm (i) + tcp3 (i) = lcp2 (i) + icp2 (i) * min (1., dim (tice, pt1 (i)) /48.) + enddo + + ! ----------------------------------------------------------------------- + !> - Condensation/evaporation between water vapor and cloud water. + ! ----------------------------------------------------------------------- + + call wqs2_vect (is, ie, pt1, den, wqsat, dq2dt) + + adj_fac = sat_adj0 + do i = is, ie + dq0 = (qv (i, j) - wqsat (i)) / (1. + tcp3 (i) * dq2dt (i)) + if (dq0 > 0.) then ! whole grid - box saturated + src (i) = min (adj_fac * dq0, max (ql_gen - ql (i, j), fac_v2l * dq0)) + else ! evaporation of ql + ! sjl 20170703 added ql factor to prevent the situation of high ql and rh<1 + ! factor = - min (1., fac_l2v * sqrt (max (0., ql (i, j)) / 1.e-5) * 10. * (1. - qv (i, j) / wqsat (i))) + ! factor = - fac_l2v + ! factor = - 1 + factor = - min (1., fac_l2v * 10. * (1. - qv (i, j) / wqsat (i))) ! the rh dependent factor = 1 at 90% + src (i) = - min (ql (i, j), factor * dq0) + endif + qv (i, j) = qv (i, j) - src (i) +#ifdef MULTI_GASES + qvi(i,j,1,1) = qv (i, j) +#endif + ql (i, j) = ql (i, j) + src (i) + q_liq (i) = q_liq (i) + src (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) + src (i) * lhl (i) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhl (i) = lv00 + d0_vap * pt1 (i) + lhi (i) = li00 + dc_ice * pt1 (i) + lcp2 (i) = lhl (i) / cvm (i) + icp2 (i) = lhi (i) / cvm (i) + tcp3 (i) = lcp2 (i) + icp2 (i) * min (1., dim (tice, pt1 (i)) / 48.) + enddo + + if (last_step) then + + ! ----------------------------------------------------------------------- + !> - condensation/evaporation between water vapor and cloud water, last time step + !! enforce upper (no super_sat) & lower (critical rh) bounds. + ! final iteration: + ! ----------------------------------------------------------------------- + + call wqs2_vect (is, ie, pt1, den, wqsat, dq2dt) + + do i = is, ie + dq0 = (qv (i, j) - wqsat (i)) / (1. + tcp3 (i) * dq2dt (i)) + if (dq0 > 0.) then ! remove super - saturation, prevent super saturation over water + src (i) = dq0 + else ! evaporation of ql + ! factor = - min (1., fac_l2v * sqrt (max (0., ql (i, j)) / 1.e-5) * 10. * (1. - qv (i, j) / wqsat (i))) ! the rh dependent factor = 1 at 90% + ! factor = - fac_l2v + ! factor = - 1 + factor = - min (1., fac_l2v * 10. * (1. - qv (i, j) / wqsat (i))) ! the rh dependent factor = 1 at 90% + src (i) = - min (ql (i, j), factor * dq0) + endif + adj_fac = 1. + qv (i, j) = qv (i, j) - src (i) +#ifdef MULTI_GASES + qvi(i,j,1,1) = qv(i,j) +#endif + ql (i, j) = ql (i, j) + src (i) + q_liq (i) = q_liq (i) + src (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) + src (i) * lhl (i) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhl (i) = lv00 + d0_vap * pt1 (i) + lhi (i) = li00 + dc_ice * pt1 (i) + lcp2 (i) = lhl (i) / cvm (i) + icp2 (i) = lhi (i) / cvm (i) + enddo + + endif + + ! ----------------------------------------------------------------------- + !> - Homogeneous freezing of cloud water to cloud ice. + ! ----------------------------------------------------------------------- + + do i = is, ie + dtmp = t_wfr - pt1 (i) ! [ - 40, - 48] + if (ql (i, j) > 0. .and. dtmp > 0.) then + sink (i) = min (ql (i, j), ql (i, j) * dtmp * 0.125, dtmp / icp2 (i)) + ql (i, j) = ql (i, j) - sink (i) + qi (i, j) = qi (i, j) + sink (i) + q_liq (i) = q_liq (i) - sink (i) + q_sol (i) = q_sol (i) + sink (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i) + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhi (i) = li00 + dc_ice * pt1 (i) + icp2 (i) = lhi (i) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - bigg mechanism (heterogeneous freezing of cloud water to cloud ice). + ! ----------------------------------------------------------------------- + + do i = is, ie + tc = tice0 - pt1 (i) + if (ql (i, j) > 0.0 .and. tc > 0.) then + sink (i) = 3.3333e-10 * dt_bigg * (exp (0.66 * tc) - 1.) * den (i) * ql (i, j) ** 2 + sink (i) = min (ql (i, j), tc / icp2 (i), sink (i)) + ql (i, j) = ql (i, j) - sink (i) + qi (i, j) = qi (i, j) + sink (i) + q_liq (i) = q_liq (i) - sink (i) + q_sol (i) = q_sol (i) + sink (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i) + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhi (i) = li00 + dc_ice * pt1 (i) + icp2 (i) = lhi (i) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Freezing of rain to graupel. + ! ----------------------------------------------------------------------- + + do i = is, ie + dtmp = (tice - 0.1) - pt1 (i) + if (qr (i, j) > 1.e-7 .and. dtmp > 0.) then + tmp = min (1., (dtmp * 0.025) ** 2) * qr (i, j) ! no limit on freezing below - 40 deg c + sink (i) = min (tmp, fac_r2g * dtmp / icp2 (i)) + qr (i, j) = qr (i, j) - sink (i) + qg (i, j) = qg (i, j) + sink (i) + q_liq (i) = q_liq (i) - sink (i) + q_sol (i) = q_sol (i) + sink (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i) + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhi (i) = li00 + dc_ice * pt1 (i) + icp2 (i) = lhi (i) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Melting of snow to rain or cloud water. + ! ----------------------------------------------------------------------- + + do i = is, ie + dtmp = pt1 (i) - (tice + 0.1) + if (qs (i, j) > 1.e-7 .and. dtmp > 0.) then + tmp = min (1., (dtmp * 0.1) ** 2) * qs (i, j) ! no limter on melting above 10 deg c + sink (i) = min (tmp, fac_smlt * dtmp / icp2 (i)) + tmp = min (sink (i), dim (qs_mlt, ql (i, j))) ! max ql due to snow melt + qs (i, j) = qs (i, j) - sink (i) + ql (i, j) = ql (i, j) + tmp + qr (i, j) = qr (i, j) + sink (i) - tmp + ! qr (i, j) = qr (i, j) + sink (i) + q_liq (i) = q_liq (i) + sink (i) + q_sol (i) = q_sol (i) - sink (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) - sink (i) * lhi (i) / cvm (i) + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Autoconversion from cloud water to rain. + ! ----------------------------------------------------------------------- + + do i = is, ie + if (ql (i, j) > ql0_max) then + sink (i) = fac_l2r * (ql (i, j) - ql0_max) + qr (i, j) = qr (i, j) + sink (i) + ql (i, j) = ql (i, j) - sink (i) + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhi (i) = li00 + dc_ice * pt1 (i) + lhl (i) = lv00 + d0_vap * pt1 (i) + lcp2 (i) = lhl (i) / cvm (i) + icp2 (i) = lhi (i) / cvm (i) + tcp2 (i) = lcp2 (i) + icp2 (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Sublimation/deposition between water vapor and cloud ice. + ! ----------------------------------------------------------------------- + + do i = is, ie + src (i) = 0. + if (pt1 (i) < t_sub) then ! too cold to be accurate; freeze qv as a fix + src (i) = dim (qv (i, j), 1.e-6) + elseif (pt1 (i) < tice0) then + qsi = iqs2 (pt1 (i), den (i), dqsdt) + dq = qv (i, j) - qsi + sink (i) = adj_fac * dq / (1. + tcp2 (i) * dqsdt) + if (qi (i, j) > 1.e-8) then + pidep = sdt * dq * 349138.78 * exp (0.875 * log (qi (i, j) * den (i))) & + / (qsi * den (i) * lat2 / (0.0243 * rvgas * pt1 (i) ** 2) + 4.42478e4) + else + pidep = 0. + endif + if (dq > 0.) then ! vapor - > ice + tmp = tice - pt1 (i) + qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den (i) + src (i) = min (sink (i), max (qi_crt - qi (i, j), pidep), tmp / tcp2 (i)) + else + pidep = pidep * min (1., dim (pt1 (i), t_sub) * 0.2) + src (i) = max (pidep, sink (i), - qi (i, j)) + endif + endif + qv (i, j) = qv (i, j) - src (i) +#ifdef MULTI_GASES + qvi(i,j,1,1) = qv(i,j) +#endif + qi (i, j) = qi (i, j) + src (i) + q_sol (i) = q_sol (i) + src (i) + cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + pt1 (i) = pt1 (i) + src (i) * (lhl (i) + lhi (i)) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Virtual temperature updated. + ! ----------------------------------------------------------------------- + + do i = is, ie +#ifdef USE_COND + q_con (i, j) = q_liq (i) + q_sol (i) +#ifdef MULTI_GASES + pt (i, j) = pt1 (i) * virq_qpz(qvi(i,j,1,1:num_gas),q_con(i,j)) +#else + tmp = 1. + zvir * qv (i, j) + pt (i, j) = pt1 (i) * tmp * (1. - q_con (i, j)) +#endif + tmp = rdgas * tmp + cappa (i, j) = tmp / (tmp + cvm (i)) +#else +#ifdef MULTI_GASES + q_con (i, j) = q_liq (i) + q_sol (i) + pt (i, j) = pt1 (i) * virq_qpz(qvi(i,j,1,1:num_gas),q_con(i,j)) * (1. - q_con(i,j)) +#else + pt (i, j) = pt1 (i) * (1. + zvir * qv (i, j)) +#endif +#endif + enddo + + ! ----------------------------------------------------------------------- + !> - Fix negative graupel with available cloud ice. + ! ----------------------------------------------------------------------- + + do i = is, ie + if (qg (i, j) < 0.) then + tmp = min (- qg (i, j), max (0., qi (i, j))) + qg (i, j) = qg (i, j) + tmp + qi (i, j) = qi (i, j) - tmp + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Autoconversion from cloud ice to snow. + ! ----------------------------------------------------------------------- + + do i = is, ie + qim = qi0_max / den (i) + if (qi (i, j) > qim) then + sink (i) = fac_i2s * (qi (i, j) - qim) + qi (i, j) = qi (i, j) - sink (i) + qs (i, j) = qs (i, j) + sink (i) + endif + enddo + + if (out_dt) then + do i = is, ie + dtdt (i, j) = dtdt (i, j) + pt1 (i) - t0 (i) + enddo + endif + + ! ----------------------------------------------------------------------- + !> - Fix energy conservation. + ! ----------------------------------------------------------------------- + + if (consv_te) then + do i = is, ie + if (hydrostatic) then +#ifdef MULTI_GASES + c_air = cp_air * vicpqd_qpz(qvi(i,j,1,1:num_gas),qpz(i)) +#endif + te0 (i, j) = dp (i, j) * (te0 (i, j) + c_air * pt1 (i)) + else +#ifdef USE_COND + te0 (i, j) = dp (i, j) * (te0 (i, j) + cvm (i) * pt1 (i)) +#else +#ifdef MULTI_GASES + c_air = cv_air * vicvqd_qpz(qvi(i,j,1,1:num_gas),qpz(i)) +#endif + te0 (i, j) = dp (i, j) * (te0 (i, j) + c_air * pt1 (i)) +#endif + endif + enddo + endif + + ! ----------------------------------------------------------------------- + !> - Update latend heat coefficient. + ! ----------------------------------------------------------------------- + + do i = is, ie + lhi (i) = li00 + dc_ice * pt1 (i) + lhl (i) = lv00 + d0_vap * pt1 (i) + cvm (i) = mc_air (i) + (qv (i, j) + q_liq (i) + q_sol (i)) * c_vap + lcp2 (i) = lhl (i) / cvm (i) + icp2 (i) = lhi (i) / cvm (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Compute cloud fraction. + ! ----------------------------------------------------------------------- + + if (do_qa .and. last_step) then + + ! ----------------------------------------------------------------------- + !> - If it is the last step, combine water species. + ! ----------------------------------------------------------------------- + + if (rad_snow) then + if (rad_graupel) then + do i = is, ie + q_sol (i) = qi (i, j) + qs (i, j) + qg (i, j) + enddo + else + do i = is, ie + q_sol (i) = qi (i, j) + qs (i, j) + enddo + endif + else + do i = is, ie + q_sol (i) = qi (i, j) + enddo + endif + if (rad_rain) then + do i = is, ie + q_liq (i) = ql (i, j) + qr (i, j) + enddo + else + do i = is, ie + q_liq (i) = ql (i, j) + enddo + endif + do i = is, ie + q_cond (i) = q_sol (i) + q_liq (i) + enddo + + ! ----------------------------------------------------------------------- + !> - Use the "liquid - frozen water temperature" (tin) to compute saturated specific humidity. + ! ----------------------------------------------------------------------- + + do i = is, ie + + if(tintqs) then + tin = pt1(i) + else + tin = pt1 (i) - (lcp2 (i) * q_cond (i) + icp2 (i) * q_sol (i)) ! minimum temperature + ! tin = pt1 (i) - ((lv00 + d0_vap * pt1 (i)) * q_cond (i) + & + ! (li00 + dc_ice * pt1 (i)) * q_sol (i)) / (mc_air (i) + qpz (i) * c_vap) + endif + + ! ----------------------------------------------------------------------- + ! determine saturated specific humidity + ! ----------------------------------------------------------------------- + + if (tin <= t_wfr) then + ! ice phase: + qstar (i) = iqs1 (tin, den (i)) + elseif (tin >= tice) then + ! liquid phase: + qstar (i) = wqs1 (tin, den (i)) + else + ! mixed phase: + qsi = iqs1 (tin, den (i)) + qsw = wqs1 (tin, den (i)) + if (q_cond (i) > 1.e-6) then + rqi = q_sol (i) / q_cond (i) + else + ! mostly liquid water clouds at initial cloud development stage + rqi = ((tice - tin) / (tice - t_wfr)) + endif + qstar (i) = rqi * qsi + (1. - rqi) * qsw + endif + !> - higher than 10 m is considered "land" and will have higher subgrid variability + dw = dw_ocean + (dw_land - dw_ocean) * min (1., abs (hs (i, j)) / (10. * grav)) + !> - "scale - aware" subgrid variability: 100 - km as the base + hvar (i) = min (0.2, max (0.01, dw * sqrt (sqrt (area (i, j)) / 100.e3))) + + ! ----------------------------------------------------------------------- + !> - calculate partial cloudiness by pdf; + !! assuming subgrid linear distribution in horizontal; this is effectively a smoother for the + !! binary cloud scheme; qa = 0.5 if qstar (i) == qpz + ! ----------------------------------------------------------------------- + + rh = qpz (i) / qstar (i) + + ! ----------------------------------------------------------------------- + ! icloud_f = 0: bug - fixed + ! icloud_f = 1: old fvgfs gfdl) mp implementation + ! icloud_f = 2: binary cloud scheme (0 / 1) + ! ----------------------------------------------------------------------- + + if (rh > 0.75 .and. qpz (i) > 1.e-8) then + dq = hvar (i) * qpz (i) + q_plus = qpz (i) + dq + q_minus = qpz (i) - dq + if (icloud_f == 2) then + if (qpz (i) > qstar (i)) then + qa (i, j) = 1. + elseif (qstar (i) < q_plus .and. q_cond (i) > 1.e-8) then + qa (i, j) = ((q_plus - qstar (i)) / dq) ** 2 + qa (i, j) = min (1., qa (i, j)) + else + qa (i, j) = 0. + endif + else + if (qstar (i) < q_minus) then + qa (i, j) = 1. + else + if (qstar (i) < q_plus) then + if (icloud_f == 0) then + qa (i, j) = (q_plus - qstar (i)) / (dq + dq) + else + qa (i, j) = (q_plus - qstar (i)) / (2. * dq * (1. - q_cond (i))) + endif + else + qa (i, j) = 0. + endif + ! impose minimum cloudiness if substantial q_cond (i) exist + if (q_cond (i) > 1.e-8) then + qa (i, j) = max (cld_min, qa (i, j)) + endif + qa (i, j) = min (1., qa (i, j)) + endif + endif + else + qa (i, j) = 0. + endif + + enddo + + endif + + enddo ! end j loop + +end subroutine fv_sat_adj_work +!> @} + +! ======================================================================= +!>\ingroup fast_sat_adj +!>\brief the function 'wqs1' computes the +!! saturated specific humidity for table ii. +! ======================================================================= +real(kind=kind_dyn) function wqs1 (ta, den) + + implicit none + + ! pure water phase; universal dry / moist formular using air density + ! input "den" can be either dry or moist air density + + real(kind=kind_dyn), intent (in) :: ta, den + + real(kind=kind_dyn) :: es, ap1, tmin + + integer :: it + + tmin = tice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = tablew (it) + (ap1 - it) * desw (it) + wqs1 = es / (rvgas * ta * den) + +end function wqs1 + +! ======================================================================= +!>\ingroup fast_sat_adj +!>\brief the function 'wqs1' computes the saturated specific humidity +!! for table iii +! ======================================================================= +real(kind=kind_dyn) function iqs1 (ta, den) + + implicit none + + ! water - ice phase; universal dry / moist formular using air density + ! input "den" can be either dry or moist air density + + real(kind=kind_dyn), intent (in) :: ta, den + + real(kind=kind_dyn) :: es, ap1, tmin + + integer :: it + + tmin = tice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = table2 (it) + (ap1 - it) * des2 (it) + iqs1 = es / (rvgas * ta * den) + +end function iqs1 + +! ======================================================================= +!>\ingroup fast_sat_adj +!>\brief The function 'wqs2'computes the gradient of saturated specific +!! humidity for table ii +! ======================================================================= +real(kind=kind_dyn) function wqs2 (ta, den, dqdt) + + implicit none + + ! pure water phase; universal dry / moist formular using air density + ! input "den" can be either dry or moist air density + + real(kind=kind_dyn), intent (in) :: ta, den + + real(kind=kind_dyn), intent (out) :: dqdt + + real(kind=kind_dyn) :: es, ap1, tmin + + integer :: it + + tmin = tice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = tablew (it) + (ap1 - it) * desw (it) + wqs2 = es / (rvgas * ta * den) + it = ap1 - 0.5 + ! finite diff, del_t = 0.1: + dqdt = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta * den) + +end function wqs2 + +! ======================================================================= +!>\ingroup fast_sat_adj +!>\brief The function wqs2_vect computes the gradient of saturated +!! specific humidity for table ii. +!! It is the same as "wqs2", but written as vector function. +! ======================================================================= +subroutine wqs2_vect (is, ie, ta, den, wqsat, dqdt) + + implicit none + + ! pure water phase; universal dry / moist formular using air density + ! input "den" can be either dry or moist air density + + integer, intent (in) :: is, ie + + real(kind=kind_dyn), intent (in), dimension (is:ie) :: ta, den + + real(kind=kind_dyn), intent (out), dimension (is:ie) :: wqsat, dqdt + + real(kind=kind_dyn) :: es, ap1, tmin + + integer :: i, it + + tmin = tice - 160. + + do i = is, ie + ap1 = 10. * dim (ta (i), tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = tablew (it) + (ap1 - it) * desw (it) + wqsat (i) = es / (rvgas * ta (i) * den (i)) + it = ap1 - 0.5 + ! finite diff, del_t = 0.1: + dqdt (i) = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta (i) * den (i)) + enddo + +end subroutine wqs2_vect + +! ======================================================================= +!>\ingroup fast_sat_adj +!>\brief The function 'iqs2' computes the gradient of saturated specific +!! humidity for table iii. +! ======================================================================= +real(kind=kind_dyn) function iqs2 (ta, den, dqdt) + + implicit none + + ! water - ice phase; universal dry / moist formular using air density + ! input "den" can be either dry or moist air density + + real(kind=kind_dyn), intent (in) :: ta, den + + real(kind=kind_dyn), intent (out) :: dqdt + + real(kind=kind_dyn) :: es, ap1, tmin + + integer :: it + + tmin = tice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = table2 (it) + (ap1 - it) * des2 (it) + iqs2 = es / (rvgas * ta * den) + it = ap1 - 0.5 + ! finite diff, del_t = 0.1: + dqdt = 10. * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) / (rvgas * ta * den) + +end function iqs2 + +! ======================================================================= +!>\ingroup fast_sat_adj +!! saturation water vapor pressure table i +! 3 - phase table +! ======================================================================= + +subroutine qs_table (n) + + implicit none + + integer, intent (in) :: n + real (kind_grid) :: delt = 0.1 + real (kind_grid) :: tmin, tem, esh20 + real (kind_grid) :: wice, wh2o, fac0, fac1, fac2 + real (kind_grid) :: esupc (200) + integer :: i + + tmin = tice - 160. + + ! ----------------------------------------------------------------------- + ! compute es over ice between - 160 deg c and 0 deg c. + ! ----------------------------------------------------------------------- + + do i = 1, 1600 + tem = tmin + delt * real (i - 1) + fac0 = (tem - tice) / (tem * tice) + fac1 = fac0 * li2 + fac2 = (d2ice * log (tem / tice) + fac1) / rvgas + table (i) = e00 * exp (fac2) + enddo + + ! ----------------------------------------------------------------------- + ! compute es over water between - 20 deg c and 102 deg c. + ! ----------------------------------------------------------------------- + + do i = 1, 1221 + tem = 253.16 + delt * real (i - 1) + fac0 = (tem - tice) / (tem * tice) + fac1 = fac0 * lv0 + fac2 = (dc_vap * log (tem / tice) + fac1) / rvgas + esh20 = e00 * exp (fac2) + if (i <= 200) then + esupc (i) = esh20 + else + table (i + 1400) = esh20 + endif + enddo + + ! ----------------------------------------------------------------------- + ! derive blended es over ice and supercooled water between - 20 deg c and 0 deg c + ! ----------------------------------------------------------------------- + + do i = 1, 200 + tem = 253.16 + delt * real (i - 1) + wice = 0.05 * (tice - tem) + wh2o = 0.05 * (tem - 253.16) + table (i + 1400) = wice * table (i + 1400) + wh2o * esupc (i) + enddo + +end subroutine qs_table + +! ======================================================================= +!>\ingroup fast_sat_adj +!! saturation water vapor pressure table ii. +! 1 - phase table +! ======================================================================= + +subroutine qs_tablew (n) + + implicit none + + integer, intent (in) :: n + real (kind_grid) :: delt = 0.1 + real (kind_grid) :: tmin, tem, fac0, fac1, fac2 + integer :: i + + tmin = tice - 160. + + ! ----------------------------------------------------------------------- + ! compute es over water + ! ----------------------------------------------------------------------- + + do i = 1, n + tem = tmin + delt * real (i - 1) + fac0 = (tem - tice) / (tem * tice) + fac1 = fac0 * lv0 + fac2 = (dc_vap * log (tem / tice) + fac1) / rvgas + tablew (i) = e00 * exp (fac2) + enddo + +end subroutine qs_tablew + +! ======================================================================= +!>\ingroup fast_sat_adj +!! saturation water vapor pressure table iii. +! 2 - phase table +! ======================================================================= + +subroutine qs_table2 (n) + + implicit none + + integer, intent (in) :: n + real (kind_grid) :: delt = 0.1 + real (kind_grid) :: tmin, tem0, tem1, fac0, fac1, fac2 + integer :: i, i0, i1 + + tmin = tice - 160. + + do i = 1, n + tem0 = tmin + delt * real (i - 1) + fac0 = (tem0 - tice) / (tem0 * tice) + if (i <= 1600) then + ! ----------------------------------------------------------------------- + ! compute es over ice between - 160 deg c and 0 deg c. + ! ----------------------------------------------------------------------- + fac1 = fac0 * li2 + fac2 = (d2ice * log (tem0 / tice) + fac1) / rvgas + else + ! ----------------------------------------------------------------------- + ! compute es over water between 0 deg c and 102 deg c. + ! ----------------------------------------------------------------------- + fac1 = fac0 * lv0 + fac2 = (dc_vap * log (tem0 / tice) + fac1) / rvgas + endif + table2 (i) = e00 * exp (fac2) + enddo + + ! ----------------------------------------------------------------------- + ! smoother around 0 deg c + ! ----------------------------------------------------------------------- + + i0 = 1600 + i1 = 1601 + tem0 = 0.25 * (table2 (i0 - 1) + 2. * table (i0) + table2 (i0 + 1)) + tem1 = 0.25 * (table2 (i1 - 1) + 2. * table (i1) + table2 (i1 + 1)) + table2 (i0) = tem0 + table2 (i1) = tem1 + +end subroutine qs_table2 + +end module fv_sat_adj +!> @} diff --git a/physics/MP/GFDL/v1/fv_sat_adj.meta b/physics/MP/GFDL/v1/fv_sat_adj.meta new file mode 100644 index 000000000..c79eb6a25 --- /dev/null +++ b/physics/MP/GFDL/v1/fv_sat_adj.meta @@ -0,0 +1,441 @@ +[ccpp-table-properties] + name = fv_sat_adj + type = scheme + dependencies = ../../../hooks/machine.F,../../../hooks/physcons.F90 + dependencies = module_gfdl_cloud_microphys.F90,../multi_gases.F90 + dependencies = ../../module_mp_radar.F90 + +######################################################################## +[ccpp-arg-table] + name = fv_sat_adj_init + type = scheme +[do_sat_adj] + standard_name = flag_for_saturation_adjustment_for_microphysics_in_dynamics + long_name = flag for saturation adjustment for microphysics in dynamics + units = none + dimensions = () + type = logical + intent = in +[kmp] + standard_name = top_layer_index_for_fast_physics + long_name = top_layer_inder_for_gfdl_mp + units = index + dimensions = () + type = integer + intent = in +[nwat] + standard_name = number_of_water_species + long_name = number of water species + units = count + dimensions = () + type = integer + intent = in +[ngas] + standard_name = number_of_gases_for_multi_gases_physics + long_name = number of gases for multi gases physics + units = count + dimensions = () + type = integer + intent = in +[rilist] + standard_name = gas_constants_for_multi_gases_physics + long_name = gas constants for multi gases physics + units = J kg-1 K-1 + dimensions = (0:number_of_gases_for_multi_gases_physics) + type = real + kind = kind_dyn + intent = in +[cpilist] + standard_name = specific_heat_capacities_for_multi_gases_physics + long_name = specific heat capacities for multi gases physics + units = J kg-1 K-1 + dimensions = (0:number_of_gases_for_multi_gases_physics) + type = real + kind = kind_dyn + intent = in +[mpirank] + standard_name = mpi_rank_for_fast_physics + long_name = current MPI-rank for fast physics schemes + units = index + dimensions = () + type = integer + intent = in +[mpiroot] + standard_name = mpi_root_for_fast_physics + long_name = master MPI-rank for fast physics schemes + units = index + dimensions = () + type = integer + intent = in +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out + +######################################################################## +[ccpp-arg-table] + name = fv_sat_adj_finalize + type = scheme +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out + +######################################################################## +[ccpp-arg-table] + name = fv_sat_adj_run + type = scheme +[mdt] + standard_name = time_step_for_remapping_for_fast_physics + long_name = remapping time step for fast physics + units = s + dimensions = () + type = real + kind = kind_dyn + intent = in +[zvir] + standard_name = ratio_of_vapor_to_dry_air_gas_constants_minus_one_default_kind + long_name = zvir=rv/rd-1.0 + units = none + dimensions = () + type = real + kind = kind_dyn + intent = in +[is] + standard_name = starting_x_direction_index + long_name = starting X direction index + units = count + dimensions = () + type = integer + intent = in +[ie] + standard_name = ending_x_direction_index + long_name = ending X direction index + units = count + dimensions = () + type = integer + intent = in +[isd] + standard_name = starting_x_direction_index_domain + long_name = starting X direction index for domain + units = count + dimensions = () + type = integer + intent = in +[ied] + standard_name = ending_x_direction_index_domain + long_name = ending X direction index for domain + units = count + dimensions = () + type = integer + intent = in +[kmp] + standard_name = top_layer_index_for_fast_physics + long_name = top layer index for GFDL mp + units = index + dimensions = () + type = integer + intent = in +[km] + standard_name = vertical_dimension_for_fast_physics + long_name = number of vertical levels + units = count + dimensions = () + type = integer + intent = in +[kmdelz] + standard_name = vertical_dimension_for_thickness_at_Lagrangian_surface + long_name = vertical dimension for thickness at Lagrangian surface + units = count + dimensions = () + type = integer + intent = in +[js] + standard_name = starting_y_direction_index + long_name = starting Y direction index + units = count + dimensions = () + type = integer + intent = in +[je] + standard_name = ending_y_direction_index + long_name = ending Y direction index + units = count + dimensions = () + type = integer + intent = in +[jsd] + standard_name = starting_y_direction_index_domain + long_name = starting X direction index for domain + units = count + dimensions = () + type = integer + intent = in +[jed] + standard_name = ending_y_direction_index_domain + long_name = ending X direction index for domain + units = count + dimensions = () + type = integer + intent = in +[ng] + standard_name = number_of_ghost_zones + long_name = number of ghost zones defined in fv_mp + units = count + dimensions = () + type = integer + intent = in +[hydrostatic] + standard_name = flag_for_hydrostatic_solver_for_fast_physics + long_name = flag for use the hydrostatic or nonhydrostatic solver + units = flag + dimensions = () + type = logical + intent = in +[fast_mp_consv] + standard_name = flag_for_fast_microphysics_energy_conservation + long_name = flag for fast microphysics energy conservation + units = flag + dimensions = () + type = logical + intent = in +[te0_2d] + standard_name = atmosphere_energy_content_in_column + long_name = atmosphere total energy in columns + units = J m-2 + dimensions = (starting_x_direction_index:ending_x_direction_index,starting_y_direction_index:ending_y_direction_index) + type = real + kind = kind_dyn + intent = inout +[te0] + standard_name = atmosphere_energy_content_at_Lagrangian_surface + long_name = atmosphere total energy at Lagrangian surface + units = J m-2 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = out +[ngas] + standard_name = number_of_gases_for_multi_gases_physics + long_name = number of gases for multi gases physics + units = count + dimensions = () + type = integer + intent = in +[qvi] + standard_name = gas_tracers_for_multi_gas_physics_at_Lagrangian_surface + long_name = gas tracers for multi gas physics at Lagrangian surface + units = kg kg-1 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics,1:number_of_gases_for_multi_gases_physics) + type = real + kind = kind_dyn + intent = inout +[qv] + standard_name = water_vapor_specific_humidity_at_Lagrangian_surface + long_name = water vapor specific humidity updated by fast physics at Lagrangian surface + units = kg kg-1 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[ql] + standard_name = cloud_liquid_water_specific_humidity_at_Lagrangian_surface + long_name = cloud liquid water specific humidity updated by fast physics at Lagrangian surface + units = kg kg-1 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[qi] + standard_name = cloud_ice_specific_humidity_at_Lagrangian_surface + long_name = cloud ice specific humidity updated by fast physics at Lagrangian surface + units = kg kg-1 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[qr] + standard_name = cloud_rain_specific_humidity_at_Lagrangian_surface + long_name = cloud rain specific humidity updated by fast physics at Lagrangian surface + units = kg kg-1 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[qs] + standard_name = cloud_snow_specific_humidity_at_Lagrangian_surface + long_name = cloud snow specific humidity updated by fast physics at Lagrangian surface + units = kg kg-1 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[qg] + standard_name = cloud_graupel_specific_humidity_at_Lagrangian_surface + long_name = cloud graupel specific humidity updated by fast physics at Lagrangian surface + units = kg kg-1 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[hs] + standard_name = surface_geopotential_at_Lagrangian_surface + long_name = surface geopotential at Lagrangian surface + units = m2 s-2 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain) + type = real + kind = kind_dyn + intent = in +[peln] + standard_name = log_pressure_at_Lagrangian_surface + long_name = logarithm of pressure at Lagrangian surface + units = Pa + dimensions = (starting_x_direction_index:ending_x_direction_index,1:vertical_dimension_for_fast_physics_plus_one,starting_y_direction_index:ending_y_direction_index) + type = real + kind = kind_dyn + intent = in +[delz] + standard_name = thickness_at_Lagrangian_surface + long_name = thickness at Lagrangian_surface + units = m + dimensions = (starting_x_direction_index:ending_x_direction_index,starting_y_direction_index:ending_y_direction_index,1:vertical_dimension_for_thickness_at_Lagrangian_surface) + type = real + kind = kind_dyn + intent = in +[delp] + standard_name = pressure_thickness_at_Lagrangian_surface + long_name = pressure thickness at Lagrangian surface + units = Pa + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = in +[pt] + standard_name = virtual_temperature_at_Lagrangian_surface + long_name = virtual temperature at Lagrangian surface + units = K + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[pkz] + standard_name = finite_volume_mean_edge_pressure_raised_to_the_power_of_kappa + long_name = finite-volume mean edge pressure in Pa raised to the power of kappa + units = 1 + dimensions = (starting_x_direction_index:ending_x_direction_index,starting_y_direction_index:ending_y_direction_index,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[q_con] + standard_name = cloud_condensed_water_specific_humidity_at_Lagrangian_surface + long_name = cloud condensed water specific humidity updated by fast physics at Lagrangian surface + units = kg kg-1 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_condensed_water_at_Lagrangian_surface) + type = real + kind = kind_dyn + intent = inout +[akap] + standard_name = kappa_dry_for_fast_physics + long_name = modified kappa for dry air, fast physics + units = none + dimensions = () + type = real + kind = kind_dyn + intent = in +[cappa] + standard_name = cappa_moist_gas_constant_at_Lagrangian_surface + long_name = cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/(1.+r_vir*q(i,j,k,sphum)) ) + units = none + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_cappa_at_Lagrangian_surface) + type = real + kind = kind_dyn + intent = inout +[area] + standard_name = cell_area_for_fast_physics + long_name = area of the grid cell for fast physics + units = m2 + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain) + type = real + kind = kind_grid + intent = in +[dtdt] + standard_name = tendency_of_air_temperature_at_Lagrangian_surface + long_name = air temperature tendency due to fast physics at Lagrangian surface + units = K s-1 + dimensions = (starting_x_direction_index:ending_x_direction_index,starting_y_direction_index:ending_y_direction_index,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = inout +[out_dt] + standard_name = flag_for_tendency_of_air_temperature_at_Lagrangian_surface + long_name = flag for calculating tendency of air temperature due to fast physics + units = flag + dimensions = () + type = logical + intent = in +[last_step] + standard_name = flag_for_the_last_step_of_k_split_remapping + long_name = flag for the last step of k-split remapping + units = flag + dimensions = () + type = logical + intent = in +[do_qa] + standard_name = flag_for_inline_cloud_fraction_calculation + long_name = flag for the inline cloud fraction calculation + units = flag + dimensions = () + type = logical + intent = in +[qa] + standard_name = cloud_fraction_at_Lagrangian_surface + long_name = cloud fraction at Lagrangian surface + units = none + dimensions = (starting_x_direction_index_domain:ending_x_direction_index_domain,starting_y_direction_index_domain:ending_y_direction_index_domain,1:vertical_dimension_for_fast_physics) + type = real + kind = kind_dyn + intent = out +[nthreads] + standard_name = omp_threads_for_fast_physics + long_name = number of OpenMP threads available for fast physics schemes + units = count + dimensions = () + type = integer + intent = in +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out diff --git a/physics/MP/GFDL/v1/gfdl_cloud_microphys.F90 b/physics/MP/GFDL/v1/gfdl_cloud_microphys.F90 new file mode 100644 index 000000000..0fd84c7ea --- /dev/null +++ b/physics/MP/GFDL/v1/gfdl_cloud_microphys.F90 @@ -0,0 +1,331 @@ +!> \file gfdl_cloud_microphys.F90 +!! This file contains the CCPP entry point for the column GFDL cloud microphysics ( Chen and Lin (2013) +!! \cite chen_and_lin_2013 ). +module gfdl_cloud_microphys + + use gfdl_cloud_microphys_mod, only: gfdl_cloud_microphys_mod_init, & + gfdl_cloud_microphys_mod_driver, & + gfdl_cloud_microphys_mod_end, & + cloud_diagnosis + + implicit none + + private + + public gfdl_cloud_microphys_run, gfdl_cloud_microphys_init, gfdl_cloud_microphys_finalize + + logical :: is_initialized = .false. + +contains + +! ----------------------------------------------------------------------- +! CCPP entry points for gfdl cloud microphysics +! ----------------------------------------------------------------------- + +!>\brief The subroutine initializes the GFDL +!! cloud microphysics. +!! +!> \section arg_table_gfdl_cloud_microphys_init Argument Table +!! \htmlinclude gfdl_cloud_microphys_init.html +!! + subroutine gfdl_cloud_microphys_init (me, master, nlunit, input_nml_file, logunit, fn_nml, & + imp_physics, imp_physics_gfdl, do_shoc, errmsg, errflg) + + implicit none + + integer, intent (in) :: me + integer, intent (in) :: master + integer, intent (in) :: nlunit + integer, intent (in) :: logunit + character(len=*), intent (in) :: fn_nml + character(len=*), intent (in) :: input_nml_file(:) + integer, intent( in) :: imp_physics + integer, intent( in) :: imp_physics_gfdl + logical, intent( in) :: do_shoc + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (is_initialized) return + + if (imp_physics/=imp_physics_gfdl) then + write(errmsg,'(*(a))') 'Namelist option for microphysics does not match choice in suite definition file' + errflg = 1 + return + end if + + if (do_shoc) then + write(errmsg,'(*(a))') 'SHOC is not currently compatible with GFDL MP' + errflg = 1 + return + endif + + call gfdl_cloud_microphys_mod_init(me, master, nlunit, input_nml_file, logunit, fn_nml, errmsg, errflg) + + is_initialized = .true. + + end subroutine gfdl_cloud_microphys_init + +! ======================================================================= +!>\brief The subroutine 'gfdl_cloud_microphys_finalize' terminates the GFDL +!! cloud microphysics. +!! +!! \section arg_table_gfdl_cloud_microphys_finalize Argument Table +!! \htmlinclude gfdl_cloud_microphys_finalize.html +!! + subroutine gfdl_cloud_microphys_finalize(errmsg, errflg) + + implicit none + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (.not.is_initialized) return + + call gfdl_cloud_microphys_mod_end() + + is_initialized = .false. + + end subroutine gfdl_cloud_microphys_finalize + +!>\defgroup gfdlmp GFDL Cloud Microphysics Module +!! This is cloud microphysics package for GFDL global cloud resolving model. +!! The algorithms are originally derived from Lin et al. (1983) \cite lin_et_al_1983. +!! Most of the key elements have been simplified/improved. This code at this stage +!! bears little to no similarity to the original Lin MP. +!! Therefore, it is best to be called GFDL microphysics (GFDL MP) . +!! +!>\brief The module contains the GFDL cloud +!! microphysics (Chen and Lin (2013) \cite chen_and_lin_2013 ). +!> The module is paired with \ref fast_sat_adj, which performs the "fast" +!! processes. +!! +!>\brief The subroutine executes the full GFDL cloud microphysics. +!! \section arg_table_gfdl_cloud_microphys_run Argument Table +!! \htmlinclude gfdl_cloud_microphys_run.html +!! + subroutine gfdl_cloud_microphys_run( & + levs, im, rainmin, con_g, con_fvirt, con_rd, con_eps, frland, garea, islmsk, & + gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, gq0_ntsw, gq0_ntgl, gq0_ntclamt, & + gt0, gu0, gv0, vvl, prsl, phii, del, & + rain0, ice0, snow0, graupel0, prcp0, sr, & + dtp, hydrostatic, phys_hydrostatic, lradar, refl_10cm, & + reset, effr_in, rew, rei, rer, res, reg, & + cplchm, pfi_lsan, pfl_lsan, errmsg, errflg) + + use machine, only: kind_phys + + implicit none + + ! DH* TODO: CLEANUP, all of these should be coming in through the argument list + ! parameters + real(kind=kind_phys), parameter :: one = 1.0d0 + real(kind=kind_phys), parameter :: con_p001= 0.001d0 + real(kind=kind_phys), parameter :: con_day = 86400.d0 + !real(kind=kind_phys), parameter :: rainmin = 1.0d-13 + ! *DH + + ! interface variables + integer, intent(in ) :: levs, im + real(kind=kind_phys), intent(in ) :: con_g, con_fvirt, con_rd, con_eps, rainmin + real(kind=kind_phys), intent(in ), dimension(:) :: frland, garea + integer, intent(in ), dimension(:) :: islmsk + real(kind=kind_phys), intent(inout), dimension(:,:) :: gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, & + gq0_ntsw, gq0_ntgl, gq0_ntclamt + real(kind=kind_phys), intent(inout), dimension(:,:) :: gt0, gu0, gv0 + real(kind=kind_phys), intent(in ), dimension(:,:) :: vvl, prsl, del + real(kind=kind_phys), intent(in ), dimension(:,:) :: phii + + ! rain/snow/ice/graupel/precip amounts, fraction of frozen precip + real(kind_phys), intent(out ), dimension(:) :: rain0 + real(kind_phys), intent(out ), dimension(:) :: snow0 + real(kind_phys), intent(out ), dimension(:) :: ice0 + real(kind_phys), intent(out ), dimension(:) :: graupel0 + real(kind_phys), intent(out ), dimension(:) :: prcp0 + real(kind_phys), intent(out ), dimension(:) :: sr + + real(kind_phys), intent(in) :: dtp ! physics time step + logical, intent (in) :: hydrostatic, phys_hydrostatic + + logical, intent (in) :: lradar + real(kind=kind_phys), intent(inout), dimension(:,:) :: refl_10cm + logical, intent (in) :: reset, effr_in + real(kind=kind_phys), intent(inout), dimension(:,:) :: rew, rei, rer, res, reg + logical, intent (in) :: cplchm + ! ice and liquid water 3d precipitation fluxes - only allocated if cplchm is .true. + real(kind=kind_phys), intent(inout), dimension(:,:) :: pfi_lsan, pfl_lsan + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! local variables + integer :: iis, iie, jjs, jje, kks, kke, kbot, ktop + integer :: i, k, kk + real(kind=kind_phys), dimension(1:im,1:levs) :: delp, dz, uin, vin, pt, qv1, ql1, qr1, qg1, qa1, qn1, qi1, & + qs1, pt_dt, qa_dt, u_dt, v_dt, w, qv_dt, ql_dt, qr_dt, qi_dt, & + qs_dt, qg_dt, p123, refl + real(kind=kind_phys), dimension(1:im,1,1:levs) :: pfils, pflls + real(kind=kind_phys), dimension(:,:), allocatable :: den + real(kind=kind_phys) :: onebg + real(kind=kind_phys) :: tem + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + iis = 1 + iie = im + jjs = 1 + jje = 1 + kks = 1 + kke = levs + ! flipping of vertical direction + ktop = 1 + kbot = levs + + onebg = one/con_g + + do k = 1, levs + kk = levs-k+1 + do i = 1, im + qv_dt(i,k) = 0.0 + ql_dt(i,k) = 0.0 + qr_dt(i,k) = 0.0 + qi_dt(i,k) = 0.0 + qs_dt(i,k) = 0.0 + qg_dt(i,k) = 0.0 + qa_dt(i,k) = 0.0 + pt_dt(i,k) = 0.0 + u_dt(i,k) = 0.0 + v_dt(i,k) = 0.0 + qn1(i,k) = 0.0 + pfils(i,1,k) = 0.0 + pflls(i,1,k) = 0.0 + ! flip vertical (k) coordinate + qv1(i,k) = gq0(i,kk) + ql1(i,k) = gq0_ntcw(i,kk) + qr1(i,k) = gq0_ntrw(i,kk) + qi1(i,k) = gq0_ntiw(i,kk) + qs1(i,k) = gq0_ntsw(i,kk) + qg1(i,k) = gq0_ntgl(i,kk) + qa1(i,k) = gq0_ntclamt(i,kk) + pt(i,k) = gt0(i,kk) + w(i,k) = -vvl(i,kk) * (one+con_fvirt * gq0(i,kk)) & + * gt0(i,kk) / prsl(i,kk) * (con_rd*onebg) + uin(i,k) = gu0(i,kk) + vin(i,k) = gv0(i,kk) + delp(i,k) = del(i,kk) + dz(i,k) = (phii(i,kk)-phii(i,kk+1))*onebg + p123(i,k) = prsl(i,kk) + refl(i,k) = refl_10cm(i,kk) + enddo + enddo + + ! reset precipitation amounts to zero + rain0 = 0 + ice0 = 0 + snow0 = 0 + graupel0 = 0 + + call gfdl_cloud_microphys_mod_driver(iis, iie, jjs, jje, kks, kke, ktop, kbot, & + qv1, ql1, qr1, qi1, qs1, qg1, qa1, qn1, qv_dt, ql_dt, qr_dt, qi_dt, & + qs_dt, qg_dt, qa_dt, pt_dt, pt, w, uin, vin, u_dt, v_dt, dz, delp, & + garea, dtp, frland, rain0, snow0, ice0, graupel0, hydrostatic, & + phys_hydrostatic, p123, lradar, refl, reset, pfils, pflls) + tem = dtp*con_p001/con_day + + ! fix negative values + do i = 1, im + !rain0(i) = max(con_d00, rain0(i)) + !snow0(i) = max(con_d00, snow0(i)) + !ice0(i) = max(con_d00, ice0(i)) + !graupel0(i) = max(con_d00, graupel0(i)) + if(rain0(i)*tem < rainmin) then + rain0(i) = 0.0 + endif + if(ice0(i)*tem < rainmin) then + ice0(i) = 0.0 + endif + if(snow0(i)*tem < rainmin) then + snow0(i) = 0.0 + endif + if(graupel0(i)*tem < rainmin) then + graupel0(i) = 0.0 + endif + enddo + + ! calculate fraction of frozen precipitation using unscaled + ! values of rain0, ice0, snow0, graupel0 (for bit-for-bit) + do i=1,im + prcp0(i) = (rain0(i)+snow0(i)+ice0(i)+graupel0(i)) * tem + if ( prcp0(i) > rainmin ) then + sr(i) = (snow0(i) + ice0(i) + graupel0(i)) & + / (rain0(i) + snow0(i) + ice0(i) + graupel0(i)) + else + sr(i) = 0.0 + endif + enddo + + ! convert rain0, ice0, snow0, graupel0 from mm per day to m per physics timestep + rain0 = rain0*tem + ice0 = ice0*tem + snow0 = snow0*tem + graupel0 = graupel0*tem + + ! flip vertical coordinate back + do k=1,levs + kk = levs-k+1 + do i=1,im + gq0(i,k) = qv1(i,kk) + qv_dt(i,kk) * dtp + gq0_ntcw(i,k) = ql1(i,kk) + ql_dt(i,kk) * dtp + gq0_ntrw(i,k) = qr1(i,kk) + qr_dt(i,kk) * dtp + gq0_ntiw(i,k) = qi1(i,kk) + qi_dt(i,kk) * dtp + gq0_ntsw(i,k) = qs1(i,kk) + qs_dt(i,kk) * dtp + gq0_ntgl(i,k) = qg1(i,kk) + qg_dt(i,kk) * dtp + gq0_ntclamt(i,k) = qa1(i,kk) + qa_dt(i,kk) * dtp + gt0(i,k) = gt0(i,k) + pt_dt(i,kk) * dtp + gu0(i,k) = gu0(i,k) + u_dt(i,kk) * dtp + gv0(i,k) = gv0(i,k) + v_dt(i,kk) * dtp + refl_10cm(i,k) = refl(i,kk) + enddo + enddo + + ! output ice and liquid water 3d precipitation fluxes if requested + if (cplchm) then + do k=1,levs + kk = levs-k+1 + do i=1,im + pfi_lsan(i,k) = pfils(i,1,kk) + pfl_lsan(i,k) = pflls(i,1,kk) + enddo + enddo + endif + + if(effr_in) then + allocate(den(1:im,1:levs)) + do k=1,levs + do i=1,im + den(i,k)=con_eps*prsl(i,k)/(con_rd*gt0(i,k)*(gq0(i,k)+con_eps)) + enddo + enddo + call cloud_diagnosis (1, im, 1, levs, den(1:im,1:levs), & + del(1:im,1:levs), islmsk(1:im), & + gq0_ntcw(1:im,1:levs), gq0_ntiw(1:im,1:levs), & + gq0_ntrw(1:im,1:levs), & + gq0_ntsw(1:im,1:levs) + gq0_ntgl(1:im,1:levs), & + gq0_ntgl(1:im,1:levs)*0.0, gt0(1:im,1:levs), & + rew(1:im,1:levs), rei(1:im,1:levs), rer(1:im,1:levs),& + res(1:im,1:levs), reg(1:im,1:levs)) + deallocate(den) + endif + + end subroutine gfdl_cloud_microphys_run + +end module gfdl_cloud_microphys diff --git a/physics/MP/GFDL/v1/gfdl_cloud_microphys.meta b/physics/MP/GFDL/v1/gfdl_cloud_microphys.meta new file mode 100644 index 000000000..236d81ed3 --- /dev/null +++ b/physics/MP/GFDL/v1/gfdl_cloud_microphys.meta @@ -0,0 +1,482 @@ +[ccpp-table-properties] + name = gfdl_cloud_microphys + type = scheme + dependencies = ../../../hooks/machine.F + dependencies = ../../module_mp_radar.F90 + dependencies = module_gfdl_cloud_microphys.F90 + +######################################################################## +[ccpp-arg-table] + name = gfdl_cloud_microphys_init + type = scheme +[me] + standard_name = mpi_rank + long_name = MPI rank of current process + units = index + dimensions = () + type = integer + intent = in +[master] + standard_name = mpi_root + long_name = MPI rank of master process + units = index + dimensions = () + type = integer + intent = in +[nlunit] + standard_name = iounit_of_namelist + long_name = fortran unit number for opening nameliust file + units = none + dimensions = () + type = integer + intent = in +[input_nml_file] + standard_name = filename_of_internal_namelist + long_name = character string to store full namelist contents + units = none + dimensions = (number_of_lines_in_internal_namelist) + type = character + kind = len=* + intent = in +[logunit] + standard_name = iounit_of_log + long_name = fortran unit number for writing logfile + units = none + dimensions = () + type = integer + intent = in +[fn_nml] + standard_name = filename_of_namelist + long_name = namelist filename + units = none + dimensions = () + type = character + kind = len=* + intent = in +[imp_physics] + standard_name = control_for_microphysics_scheme + long_name = choice of microphysics scheme + units = flag + dimensions = () + type = integer + intent = in +[imp_physics_gfdl] + standard_name = identifier_for_gfdl_microphysics_scheme + long_name = choice of GFDL microphysics scheme + units = flag + dimensions = () + type = integer + intent = in +[do_shoc] + standard_name = flag_for_shoc + long_name = flag to indicate use of SHOC + units = flag + dimensions = () + type = logical + intent = in +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out + +######################################################################## +[ccpp-arg-table] + name = gfdl_cloud_microphys_finalize + type = scheme +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out + +######################################################################## +[ccpp-arg-table] + name = gfdl_cloud_microphys_run + type = scheme +[levs] + standard_name = vertical_layer_dimension + long_name = number of vertical levels + units = count + dimensions = () + type = integer + intent = in +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in +[rainmin] + standard_name = lwe_thickness_of_minimum_rain_amount + long_name = minimum rain amount + units = m + dimensions = () + type = real + kind = kind_phys + intent = in +[con_g] + standard_name = gravitational_acceleration + long_name = gravitational acceleration + units = m s-2 + dimensions = () + type = real + kind = kind_phys + intent = in +[con_fvirt] + standard_name = ratio_of_vapor_to_dry_air_gas_constants_minus_one + long_name = rv/rd - 1 (rv = ideal gas constant for water vapor) + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[con_rd] + standard_name = gas_constant_of_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in +[frland] + standard_name = land_area_fraction_for_microphysics + long_name = land area fraction used in microphysics schemes + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[garea] + standard_name = cell_area + long_name = area of grid cell + units = m2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[islmsk] + standard_name = sea_land_ice_mask + long_name = sea/land/ice mask (=0/1/2) + units = flag + dimensions = (horizontal_loop_extent) + type = integer + intent = in +[gq0] + standard_name = specific_humidity_of_new_state + long_name = water vapor specific humidity updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gq0_ntcw] + standard_name = cloud_liquid_water_mixing_ratio_of_new_state + long_name = cloud condensed water mixing ratio updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gq0_ntrw] + standard_name = rain_mixing_ratio_of_new_state + long_name = moist mixing ratio of rain updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gq0_ntiw] + standard_name = cloud_ice_mixing_ratio_of_new_state + long_name = moist mixing ratio of cloud ice updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gq0_ntsw] + standard_name = snow_mixing_ratio_of_new_state + long_name = moist mixing ratio of snow updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gq0_ntgl] + standard_name = graupel_mixing_ratio_of_new_state + long_name = moist ratio of mass of graupel to mass of dry air plus vapor (without condensates) updated by physics + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gq0_ntclamt] + standard_name = cloud_area_fraction_in_atmosphere_layer_of_new_state + long_name = cloud fraction updated by physics + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gt0] + standard_name = air_temperature_of_new_state + long_name = air temperature updated by physics + units = K + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gu0] + standard_name = x_wind_of_new_state + long_name = zonal wind updated by physics + units = m s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[gv0] + standard_name = y_wind_of_new_state + long_name = meridional wind updated by physics + units = m s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[vvl] + standard_name = lagrangian_tendency_of_air_pressure + long_name = layer mean vertical velocity + units = Pa s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[prsl] + standard_name = air_pressure + long_name = mean layer pressure + units = Pa + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[phii] + standard_name = geopotential_at_interface + long_name = geopotential at model layer interfaces + units = m2 s-2 + dimensions = (horizontal_loop_extent,vertical_interface_dimension) + type = real + kind = kind_phys + intent = in +[del] + standard_name = air_pressure_difference_between_midlayers + long_name = air pressure difference between mid-layers + units = Pa + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[rain0] + standard_name = lwe_thickness_of_explicit_rain_amount + long_name = explicit rain on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[ice0] + standard_name = lwe_thickness_of_ice_amount + long_name = ice fall on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[snow0] + standard_name = lwe_thickness_of_snow_amount + long_name = snow fall on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[graupel0] + standard_name = lwe_thickness_of_graupel_amount + long_name = graupel fall on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[prcp0] + standard_name = lwe_thickness_of_explicit_precipitation_amount + long_name = explicit precipitation (rain, ice, snow, graupel) on physics timestep + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[sr] + standard_name = ratio_of_snowfall_to_rainfall + long_name = snow ratio: ratio of snow to total precipitation + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[dtp] + standard_name = timestep_for_physics + long_name = physics timestep + units = s + dimensions = () + type = real + kind = kind_phys + intent = in +[hydrostatic] + standard_name = flag_for_hydrostatic_solver + long_name = flag indicating hydrostatic solver + units = flag + dimensions = () + type = logical + intent = in +[phys_hydrostatic] + standard_name = flag_for_hydrostatic_heating_from_physics + long_name = flag indicating hydrostatic heating from physics + units = flag + dimensions = () + type = logical + intent = in +[lradar] + standard_name = flag_for_radar_reflectivity + long_name = flag for radar reflectivity + units = flag + dimensions = () + type = logical + intent = in +[refl_10cm] + standard_name = radar_reflectivity_10cm + long_name = instantaneous refl_10cm + units = dBZ + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[reset] + standard_name = flag_reset_maximum_hourly_fields + long_name = flag for resetting maximum hourly fields + units = flag + dimensions = () + type = logical + intent = in +[effr_in] + standard_name = flag_for_cloud_effective_radii + long_name = flag for cloud effective radii calculations in GFDL microphysics + units = flag + dimensions = () + type = logical + intent = in +[rew] + standard_name = effective_radius_of_stratiform_cloud_liquid_water_particle + long_name = eff. radius of cloud liquid water particle in micrometer + units = um + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[rei] + standard_name = effective_radius_of_stratiform_cloud_ice_particle + long_name = eff. radius of cloud ice water particle in micrometer + units = um + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[rer] + standard_name = effective_radius_of_stratiform_cloud_rain_particle + long_name = effective radius of cloud rain particle in micrometers + units = um + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[res] + standard_name = effective_radius_of_stratiform_cloud_snow_particle + long_name = effective radius of cloud snow particle in micrometers + units = um + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[reg] + standard_name = effective_radius_of_stratiform_cloud_graupel_particle + long_name = eff. radius of cloud graupel particle in micrometer + units = um + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[cplchm] + standard_name = flag_for_chemistry_coupling + long_name = flag controlling cplchm collection (default off) + units = flag + dimensions = () + type = logical + intent = in +[pfi_lsan] + standard_name = ice_flux_due_to_large_scale_precipitation + long_name = instantaneous 3D flux of ice from nonconvective precipitation + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[pfl_lsan] + standard_name = liquid_flux_due_to_large_scale_precipitation + long_name = instantaneous 3D flux of liquid water from nonconvective precipitation + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out diff --git a/physics/MP/GFDL/v1/module_gfdl_cloud_microphys.F90 b/physics/MP/GFDL/v1/module_gfdl_cloud_microphys.F90 new file mode 100644 index 000000000..5cab1abbc --- /dev/null +++ b/physics/MP/GFDL/v1/module_gfdl_cloud_microphys.F90 @@ -0,0 +1,5075 @@ +!> \file gfdl_cloud_microphys.F90 +!! This file contains the full GFDL cloud microphysics (Chen and Lin (2013) +!! \cite chen_and_lin_2013 and Zhou et al. 2019 \cite zhou_etal_2019). +!! The module is paired with 'gfdl_fv_sat_adj', which performs the "fast" +!! processes +!>author Shian-Jiann Lin, Linjiong Zhou +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the GFDL Cloud Microphysics. +!* +!* The GFDL Cloud Microphysics is free software: you can +!* redistribute it and/or modify it under the terms of the +!* GNU Lesser General Public License as published by the +!* Free Software Foundation, either version 3 of the License, or +!* (at your option) any later version. +!* +!* The GFDL Cloud Microphysics is distributed in the hope it will be +!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty +!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +!* See the GNU General Public License for more details. +!* +!* You should have received a copy of the GNU Lesser General Public +!* License along with the GFDL Cloud Microphysics. +!* If not, see . +!*********************************************************************** +! ======================================================================= +!>\defgroup mod_gfdl_cloud_mp GFDL Cloud MP modules +!!\ingroup gfdlmp +!! This module contains the column GFDL Cloud microphysics scheme. +module gfdl_cloud_microphys_mod + + ! use mpp_mod, only: stdlog, mpp_pe, mpp_root_pe, mpp_clock_id, & + ! mpp_clock_begin, mpp_clock_end, clock_routine, & + ! input_nml_file + ! use diag_manager_mod, only: register_diag_field, send_data + ! use time_manager_mod, only: time_type, get_time + ! use constants_mod, only: grav, rdgas, rvgas, cp_air, hlv, hlf, pi => pi_8 + ! use fms_mod, only: write_version_number, open_namelist_file, & + ! check_nml_error, file_exist, close_file + + use module_mp_radar + + implicit none + + private + + public gfdl_cloud_microphys_mod_driver, gfdl_cloud_microphys_mod_init, & + gfdl_cloud_microphys_mod_end, cloud_diagnosis +! public wqs1, wqs2, qs_blend, wqsat_moist, wqsat2_moist +! public qsmith_init, qsmith, es2_table1d, es3_table1d, esw_table1d +! public setup_con, wet_bulb + + real :: missing_value = - 1.e10 + + logical :: module_is_initialized = .false. + logical :: qsmith_tables_initialized = .false. + + character (len = 17) :: mod_name = 'gfdl_cloud_microphys' + + real, parameter :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6 + real, parameter :: rhos = 0.1e3, rhog = 0.4e3 + real, parameter :: grav = 9.80665 !< gfs: acceleration due to gravity + real, parameter :: rdgas = 287.05 !< gfs: gas constant for dry air + real, parameter :: rvgas = 461.50 !< gfs: gas constant for water vapor + real, parameter :: cp_air = 1004.6 !< gfs: heat capacity of dry air at constant pressure + real, parameter :: hlv = 2.5e6 !< gfs: latent heat of evaporation + real, parameter :: hlf = 3.3358e5 !< gfs: latent heat of fusion + real, parameter :: pi = 3.1415926535897931 !< gfs: ratio of circle circumference to diameter + + ! real, parameter :: rdgas = 287.04 !< gfdl: gas constant for dry air + + ! real, parameter :: cp_air = rdgas * 7. / 2. ! 1004.675, heat capacity of dry air at constant pressure + real, parameter :: cp_vap = 4.0 * rvgas !< 1846.0, heat capacity of water vapore at constnat pressure + ! real, parameter :: cv_air = 717.56 ! satoh value + real, parameter :: cv_air = cp_air - rdgas !< 717.55, heat capacity of dry air at constant volume + ! real, parameter :: cv_vap = 1410.0 ! emanuel value + real, parameter :: cv_vap = 3.0 * rvgas !< 1384.5, heat capacity of water vapor at constant volume + + ! the following two are from emanuel's book "atmospheric convection" + ! real, parameter :: c_ice = 2106.0 ! heat capacity of ice at 0 deg c: c = c_ice + 7.3 * (t - tice) + ! real, parameter :: c_liq = 4190.0 ! heat capacity of water at 0 deg c + + real, parameter :: c_ice = 1972.0 !< gfdl: heat capacity of ice at - 15 deg c + real, parameter :: c_liq = 4185.5 !< gfdl: heat capacity of water at 15 deg c + ! real, parameter :: c_liq = 4218.0 ! ifs: heat capacity of liquid at 0 deg c + + real, parameter :: eps = rdgas / rvgas !< 0.6219934995 + real, parameter :: zvir = rvgas / rdgas - 1. !< 0.6077338443 + + real, parameter :: t_ice = 273.16 !< freezing temperature + real, parameter :: table_ice = 273.16 !< freezing point for qs table + + ! real, parameter :: e00 = 610.71 ! gfdl: saturation vapor pressure at 0 deg c + real, parameter :: e00 = 611.21 !< ifs: saturation vapor pressure at 0 deg c + + real, parameter :: dc_vap = cp_vap - c_liq !< - 2339.5, isobaric heating / cooling + real, parameter :: dc_ice = c_liq - c_ice !< 2213.5, isobaric heating / colling + + real, parameter :: hlv0 = hlv !< gfs: evaporation latent heat coefficient at 0 deg c + ! real, parameter :: hlv0 = 2.501e6 ! emanuel appendix - 2 + real, parameter :: hlf0 = hlf !< gfs: fussion latent heat coefficient at 0 deg c + ! real, parameter :: hlf0 = 3.337e5 ! emanuel + + real, parameter :: lv0 = hlv0 - dc_vap * t_ice !< 3.13905782e6, evaporation latent heat coefficient at 0 deg k + real, parameter :: li00 = hlf0 - dc_ice * t_ice !< - 2.7105966e5, fusion latent heat coefficient at 0 deg k + + real, parameter :: d2ice = dc_vap + dc_ice !< - 126, isobaric heating/cooling + real, parameter :: li2 = lv0 + li00 !< 2.86799816e6, sublimation latent heat coefficient at 0 deg k + + real, parameter :: qrmin = 1.e-8 !< min value for rain + real, parameter :: qvmin = 1.e-20 !< min value for water vapor (treated as zero) + real, parameter :: qcmin = 1.e-12 !< min value for cloud condensates + + real, parameter :: vr_min = 1.e-3 !< min fall speed for rain + real, parameter :: vf_min = 1.e-5 !< min fall speed for cloud ice, snow, graupel + + real, parameter :: dz_min = 1.e-2 !< use for correcting flipped height + + real, parameter :: sfcrho = 1.2 !< surface air density + real, parameter :: rhor = 1.e3 !< density of rain water, lin83 + ! intercept parameters + + real, parameter :: rnzr = 8.0e6 ! lin83 + real, parameter :: rnzs = 3.0e6 ! lin83 + real, parameter :: rnzg = 4.0e6 ! rh84 + real, parameter :: rnzh = 4.0e4 ! lin83 --- lmh 29 sep 17 + + ! density parameters + + real, parameter :: rhoh = 0.917e3 ! lin83 --- lmh 29 sep 17 + + public rhor, rhos, rhog, rhoh, rnzr, rnzs, rnzg, rnzh + real :: cracs, csacr, cgacr, cgacs, csacw, craci, csaci, cgacw, cgaci, cracw !< constants for accretions + real :: acco (3, 4) !< constants for accretions + real :: cssub (5), cgsub (5), crevp (5), cgfr (2), csmlt (5), cgmlt (5) + + real :: es0, ces0 + real :: pie, rgrav, fac_rc + real :: c_air, c_vap + + real :: lati, latv, lats, lat2, lcp, icp, tcp !< used in Bigg mechanism and wet bulk + + real :: d0_vap !< the same as dc_vap, except that cp_vap can be cp_vap or cv_vap + real :: lv00 !< the same as lv0, except that cp_vap can be cp_vap or cv_vap + + ! cloud microphysics switchers + + integer :: icloud_f = 0 !< cloud scheme + integer :: irain_f = 0 !< cloud water to rain auto conversion scheme + + logical :: de_ice = .false. !< to prevent excessive build - up of cloud ice from external sources + logical :: sedi_transport = .true. !< transport of momentum in sedimentation + logical :: do_sedi_w = .false. !< transport of vertical motion in sedimentation + logical :: do_sedi_heat = .true. !< transport of heat in sedimentation + logical :: prog_ccn = .false. !< do prognostic ccn (yi ming's method) + logical :: do_qa = .true. !< do inline cloud fraction + logical :: rad_snow = .true. !< consider snow in cloud fraciton calculation + logical :: rad_graupel = .true. !< consider graupel in cloud fraction calculation + logical :: rad_rain = .true. !< consider rain in cloud fraction calculation + logical :: fix_negative = .false. !< fix negative water species + logical :: do_setup = .true. !< setup constants and parameters + logical :: p_nonhydro = .false. !< perform hydrosatic adjustment on air density + + real, allocatable :: table (:), table2 (:), table3 (:), tablew (:) + real, allocatable :: des (:), des2 (:), des3 (:), desw (:) + + logical :: tables_are_initialized = .false. + + ! logical :: master + ! integer :: id_rh, id_vtr, id_vts, id_vtg, id_vti, id_rain, id_snow, id_graupel, & + ! id_ice, id_prec, id_cond, id_var, id_droplets + real, parameter :: dt_fr = 8. !< homogeneous freezing of all cloud water at t_wfr - dt_fr + ! minimum temperature water can exist (moore & molinero nov. 2011, nature) + ! dt_fr can be considered as the error bar + + real :: p_min = 100. !< minimum pressure (pascal) for mp to operate + + ! slj, the following parameters are for cloud - resolving resolution: 1 - 5 km + + ! qi0_crt = 0.8e-4 + ! qs0_crt = 0.6e-3 + ! c_psaci = 0.1 + ! c_pgacs = 0.1 + + ! ----------------------------------------------------------------------- + ! namelist parameters + ! ----------------------------------------------------------------------- + + real :: cld_min = 0.05 !< minimum cloud fraction + real :: tice = 273.16 !< set tice = 165. to trun off ice - phase phys (kessler emulator) + + real :: t_min = 178. !< min temp to freeze - dry all water vapor + real :: t_sub = 184. !< min temp for sublimation of cloud ice + real :: mp_time = 150. !< maximum micro - physics time step (sec) + + ! relative humidity increment + + real :: rh_inc = 0.25 !< rh increment for complete evaporation of cloud water and cloud ice + real :: rh_inr = 0.25 !< rh increment for minimum evaporation of rain + real :: rh_ins = 0.25 !< rh increment for sublimation of snow + + ! conversion time scale + + real :: tau_r2g = 900. !< rain freezing during fast_sat + real :: tau_smlt = 900. !< snow melting + real :: tau_g2r = 600. !< graupel melting to rain + real :: tau_imlt = 600. !< cloud ice melting + real :: tau_i2s = 1000. !< cloud ice to snow auto-conversion + real :: tau_l2r = 900. !< cloud water to rain auto-conversion + real :: tau_v2l = 150. !< water vapor to cloud water (condensation) + real :: tau_l2v = 300. !< cloud water to water vapor (evaporation) + real :: tau_g2v = 900. !< graupel sublimation + real :: tau_v2g = 21600. !< graupel deposition -- make it a slow process + + ! horizontal subgrid variability + + real :: dw_land = 0.20 !< base value for subgrid deviation / variability over land + real :: dw_ocean = 0.10 !< base value for ocean + + ! prescribed ccn + + real :: ccn_o = 90. !< ccn over ocean (cm^ - 3) + real :: ccn_l = 270. !< ccn over land (cm^ - 3) + + real :: rthresh = 10.0e-6 !< critical cloud drop radius (micro m) + + ! ----------------------------------------------------------------------- + ! wrf / wsm6 scheme: qi_gen = 4.92e-11 * (1.e3 * exp (0.1 * tmp)) ** 1.33 + ! optimized: qi_gen = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp))) + ! qi_gen ~ 4.808e-7 at 0 c; 1.818e-6 at - 10 c, 9.82679e-5 at - 40c + ! the following value is constructed such that qc_crt = 0 at zero c and @ - 10c matches + ! wrf / wsm6 ice initiation scheme; qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den + ! ----------------------------------------------------------------------- + + real :: sat_adj0 = 0.90 !< adjustment factor (0: no, 1: full) during fast_sat_adj + + real :: qc_crt = 5.0e-8 !< mini condensate mixing ratio to allow partial cloudiness + + real :: qi_lim = 1. !< cloud ice limiter to prevent large ice build up + + real :: ql_mlt = 2.0e-3 !< max value of cloud water allowed from melted cloud ice + real :: qs_mlt = 1.0e-6 !< max cloud water due to snow melt + + real :: ql_gen = 1.0e-3 !< max cloud water generation during remapping step if fast_sat_adj = .t. + real :: qi_gen = 1.82e-6 !< max cloud ice generation during remapping step + + ! cloud condensate upper bounds: "safety valves" for ql & qi + + real :: ql0_max = 2.0e-3 !< max cloud water value (auto converted to rain) + real :: qi0_max = 1.0e-4 !< max cloud ice value (by other sources) + + real :: qi0_crt = 1.0e-4 !< cloud ice to snow autoconversion threshold (was 1.e-4); + !! qi0_crt is highly dependent on horizontal resolution + real :: qr0_crt = 1.0e-4 !< rain to snow or graupel/hail threshold + ! lfo used * mixing ratio * = 1.e-4 (hail in lfo) + real :: qs0_crt = 1.0e-3 !< snow to graupel density threshold (0.6e-3 in purdue lin scheme) + + real :: c_paut = 0.55 !< autoconversion cloud water to rain (use 0.5 to reduce autoconversion) + real :: c_psaci = 0.02 !< accretion: cloud ice to snow (was 0.1 in zetac) + real :: c_piacr = 5.0 !< accretion: rain to ice: + real :: c_cracw = 0.9 !< rain accretion efficiency + real :: c_pgacs = 2.0e-3 !< snow to graupel "accretion" eff. (was 0.1 in zetac) + + ! decreasing clin to reduce csacw (so as to reduce cloud water --- > snow) + + real :: alin = 842.0 !< "a" in lin1983 + real :: clin = 4.8 !< "c" in lin 1983, 4.8 -- > 6. (to ehance ql -- > qs) + + ! fall velocity tuning constants: + + logical :: const_vi = .false. !< if .t. the constants are specified by v * _fac + logical :: const_vs = .false. !< if .t. the constants are specified by v * _fac + logical :: const_vg = .false. !< if .t. the constants are specified by v * _fac + logical :: const_vr = .false. !< if .t. the constants are specified by v * _fac + + ! good values: + + real :: vi_fac = 1. !< if const_vi: 1 / 3 + real :: vs_fac = 1. !< if const_vs: 1. + real :: vg_fac = 1. !< if const_vg: 2. + real :: vr_fac = 1. !< if const_vr: 4. + + ! upper bounds of fall speed (with variable speed option) + + real :: vi_max = 0.5 !< max fall speed for ice + real :: vs_max = 5.0 !< max fall speed for snow + real :: vg_max = 8.0 !< max fall speed for graupel + real :: vr_max = 12. !< max fall speed for rain + + ! cloud microphysics switchers + + logical :: fast_sat_adj = .false. !< has fast saturation adjustments + logical :: z_slope_liq = .true. !< use linear mono slope for autocconversions + logical :: z_slope_ice = .false. !< use linear mono slope for autocconversions + logical :: use_ccn = .false. !< must be true when prog_ccn is false + logical :: use_ppm = .false. !< use ppm fall scheme + logical :: mono_prof = .true. !< perform terminal fall with mono ppm scheme + logical :: mp_print = .false. !< cloud microphysics debugging printout + logical :: do_hail = .false. !< use hail parameters instead of graupel + + ! real :: global_area = - 1. + + real :: log_10, tice0, t_wfr + + integer :: reiflag = 1 + ! 1: Heymsfield and Mcfarquhar, 1996 + ! 2: Wyser, 1998 + + logical :: tintqs = .false. !< use temperature in the saturation mixing in PDF + + real :: rewmin = 5.0, rewmax = 10.0 + real :: reimin = 10.0, reimax = 150.0 + real :: rermin = 10.0, rermax = 10000.0 + real :: resmin = 150.0, resmax = 10000.0 + real :: regmin = 300.0, regmax = 10000.0 + + ! ----------------------------------------------------------------------- + ! namelist + ! ----------------------------------------------------------------------- + + namelist / gfdl_cloud_microphysics_nml / & + mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, & + vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, & + vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max, & + qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, & + const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, & + tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v, & + tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, & + z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, & + rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, & + do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, & + mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, & + resmin, resmax, regmin, regmax, tintqs, do_hail + + public & + mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, & + vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, & + vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max, & + qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, & + const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, & + tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v, & + tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, & + z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, & + rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, & + do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, & + mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, & + resmin, resmax, regmin, regmax, tintqs, do_hail + +contains + +! ----------------------------------------------------------------------- +! the driver of the gfdl cloud microphysics +! ----------------------------------------------------------------------- + +!>\ingroup mod_gfdl_cloud_mp +!! This subroutine is the driver of the GFDL cloud microphysics +subroutine gfdl_cloud_microphys_mod_driver ( & + iis, iie, jjs, jje, kks, kke, ktop, kbot, & + qv, ql, qr, qi, qs, qg, qa, qn, & + qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, pt_dt, pt, w, & + uin, vin, udt, vdt, dz, delp, area, dt_in, land, & + rain, snow, ice, graupel, hydrostatic, phys_hydrostatic, & + p, lradar, refl_10cm, reset, pfils, pflls) + + implicit none + + ! Interface variables + integer, intent (in) :: iis, iie, jjs, jje ! physics window + integer, intent (in) :: kks, kke ! vertical dimension + integer, intent (in) :: ktop, kbot ! vertical compute domain + + real, intent (in) :: dt_in ! physics time step + + real, intent (in), dimension (iis:iie, jjs:jje) :: area ! cell area + real, intent (in), dimension (iis:iie, jjs:jje) :: land ! land fraction + + real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: delp, dz, uin, vin + real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: pt, qv, ql, qr, qg, qa, qn + + real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qi, qs + real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: pt_dt, qa_dt, udt, vdt, w + real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qv_dt, ql_dt, qr_dt + real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qi_dt, qs_dt, qg_dt + + real, intent (out), dimension (iis:iie, jjs:jje) :: rain, snow, ice, graupel + + logical, intent (in) :: hydrostatic, phys_hydrostatic + + !integer, intent (in) :: seconds + real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: p + logical, intent (in) :: lradar + real, intent (out), dimension (iis:iie, jjs:jje, kks:kke) :: refl_10cm + logical, intent (in) :: reset + real, intent (out), dimension (iis:iie, jjs:jje, kks:kke) :: pfils, pflls + + ! Local variables + logical :: melti = .false. + + real :: mpdt, rdt, dts, convt, tot_prec + + integer :: i, j, k + integer :: is, ie, js, je ! physics window + integer :: ks, ke ! vertical dimension + integer :: days, ntimes, kflip + + real, dimension (iie-iis+1, jje-jjs+1) :: prec_mp, prec1, cond, w_var, rh0 + + real, dimension (iie-iis+1, jje-jjs+1,kke-kks+1) :: vt_r, vt_s, vt_g, vt_i, qn2 + + real, dimension (size(pt,1), size(pt,3)) :: m2_rain, m2_sol + + real :: allmax +!+---+-----------------------------------------------------------------+ +!For 3D reflectivity calculations + real, dimension(ktop:kbot):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dBZ +!+---+-----------------------------------------------------------------+ + + is = 1 + js = 1 + ks = 1 + ie = iie - iis + 1 + je = jje - jjs + 1 + ke = kke - kks + 1 + ! call mpp_clock_begin (gfdl_mp_clock) + + ! ----------------------------------------------------------------------- + ! define heat capacity of dry air and water vapor based on hydrostatical property + ! ----------------------------------------------------------------------- + + if (phys_hydrostatic .or. hydrostatic) then + c_air = cp_air + c_vap = cp_vap + p_nonhydro = .false. + else + c_air = cv_air + c_vap = cv_vap + p_nonhydro = .true. + endif + d0_vap = c_vap - c_liq + lv00 = hlv0 - d0_vap * t_ice + + if (hydrostatic) do_sedi_w = .false. + + ! ----------------------------------------------------------------------- + ! define latent heat coefficient used in wet bulb and Bigg mechanism + ! ----------------------------------------------------------------------- + + latv = hlv + lati = hlf + lats = latv + lati + lat2 = lats * lats + + lcp = latv / cp_air + icp = lati / cp_air + tcp = (latv + lati) / cp_air + + ! tendency zero out for am moist processes should be done outside the driver + + ! ----------------------------------------------------------------------- + ! define cloud microphysics sub time step + ! ----------------------------------------------------------------------- + + mpdt = min (dt_in, mp_time) + rdt = 1. / dt_in + ntimes = nint (dt_in / mpdt) + + ! small time step: + dts = dt_in / real (ntimes) + + ! call get_time (time, seconds, days) + + ! ----------------------------------------------------------------------- + ! initialize precipitation + ! ----------------------------------------------------------------------- + + do j = js, je + do i = is, ie + graupel (i, j) = 0. + rain (i, j) = 0. + snow (i, j) = 0. + ice (i, j) = 0. + cond (i, j) = 0. + enddo + enddo + + pfils = 0. + pflls = 0. + + ! ----------------------------------------------------------------------- + ! major cloud microphysics + ! ----------------------------------------------------------------------- + + do j = js, je + call mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, qg,& + qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, & + rain (:, j), snow (:, j), graupel (:, j), ice (:, j), m2_rain, & + m2_sol, cond (:, j), area (:, j), land (:, j), udt, vdt, pt_dt, & + qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, w_var, vt_r, & + vt_s, vt_g, vt_i, qn2) + do k = ktop, kbot + do i = is, ie + pfils(i, j, k) = m2_sol (i, k) + pflls(i, j, k) = m2_rain(i, k) + enddo + enddo + enddo + + ! ----------------------------------------------------------------------- + ! no clouds allowed above ktop + ! ----------------------------------------------------------------------- + + if (ks < ktop) then + do k = ks, ktop + if (do_qa) then + do j = js, je + do i = is, ie + qa_dt (i, j, k) = 0. + enddo + enddo + else + do j = js, je + do i = is, ie + ! qa_dt (i, j, k) = - qa (i, j, k) * rdt + qa_dt (i, j, k) = 0. ! gfs + enddo + enddo + endif + enddo + endif + + ! ----------------------------------------------------------------------- + ! diagnostic output + ! ----------------------------------------------------------------------- + + ! if (id_vtr > 0) then + ! used = send_data (id_vtr, vt_r, time, is_in = iis, js_in = jjs) + ! endif + + ! if (id_vts > 0) then + ! used = send_data (id_vts, vt_s, time, is_in = iis, js_in = jjs) + ! endif + + ! if (id_vtg > 0) then + ! used = send_data (id_vtg, vt_g, time, is_in = iis, js_in = jjs) + ! endif + + ! if (id_vti > 0) then + ! used = send_data (id_vti, vt_i, time, is_in = iis, js_in = jjs) + ! endif + + ! if (id_droplets > 0) then + ! used = send_data (id_droplets, qn2, time, is_in = iis, js_in = jjs) + ! endif + + ! if (id_var > 0) then + ! used = send_data (id_var, w_var, time, is_in = iis, js_in = jjs) + ! endif + + ! convert to mm / day + + convt = 86400. * rdt * rgrav + do j = js, je + do i = is, ie + rain (i, j) = rain (i, j) * convt + snow (i, j) = snow (i, j) * convt + ice (i, j) = ice (i, j) * convt + graupel (i, j) = graupel (i, j) * convt + prec_mp (i, j) = rain (i, j) + snow (i, j) + ice (i, j) + graupel (i, j) + enddo + enddo + + ! if (id_cond > 0) then + ! do j = js, je + ! do i = is, ie + ! cond (i, j) = cond (i, j) * rgrav + ! enddo + ! enddo + ! used = send_data (id_cond, cond, time, is_in = iis, js_in = jjs) + ! endif + + ! if (id_snow > 0) then + ! used = send_data (id_snow, snow, time, iis, jjs) + ! used = send_data (id_snow, snow, time, is_in = iis, js_in = jjs) + ! if (mp_print .and. seconds == 0) then + ! tot_prec = g_sum (snow, is, ie, js, je, area, 1) + ! if (master) write (*, *) 'mean snow = ', tot_prec + ! endif + ! endif + ! + ! if (id_graupel > 0) then + ! used = send_data (id_graupel, graupel, time, iis, jjs) + ! used = send_data (id_graupel, graupel, time, is_in = iis, js_in = jjs) + ! if (mp_print .and. seconds == 0) then + ! tot_prec = g_sum (graupel, is, ie, js, je, area, 1) + ! if (master) write (*, *) 'mean graupel = ', tot_prec + ! endif + ! endif + ! + ! if (id_ice > 0) then + ! used = send_data (id_ice, ice, time, iis, jjs) + ! used = send_data (id_ice, ice, time, is_in = iis, js_in = jjs) + ! if (mp_print .and. seconds == 0) then + ! tot_prec = g_sum (ice, is, ie, js, je, area, 1) + ! if (master) write (*, *) 'mean ice_mp = ', tot_prec + ! endif + ! endif + ! + ! if (id_rain > 0) then + ! used = send_data (id_rain, rain, time, iis, jjs) + ! used = send_data (id_rain, rain, time, is_in = iis, js_in = jjs) + ! if (mp_print .and. seconds == 0) then + ! tot_prec = g_sum (rain, is, ie, js, je, area, 1) + ! if (master) write (*, *) 'mean rain = ', tot_prec + ! endif + ! endif + ! + ! if (id_rh > 0) then !not used? + ! used = send_data (id_rh, rh0, time, iis, jjs) + ! used = send_data (id_rh, rh0, time, is_in = iis, js_in = jjs) + ! endif + ! + ! + ! if (id_prec > 0) then + ! used = send_data (id_prec, prec_mp, time, iis, jjs) + ! used = send_data (id_prec, prec_mp, time, is_in = iis, js_in = jjs) + ! endif + + ! if (mp_print) then + ! prec1 (:, :) = prec1 (:, :) + prec_mp (:, :) + ! if (seconds == 0) then + ! prec1 (:, :) = prec1 (:, :) * dt_in / 86400. + ! tot_prec = g_sum (prec1, is, ie, js, je, area, 1) + ! if (master) write (*, *) 'daily prec_mp = ', tot_prec + ! prec1 (:, :) = 0. + ! endif + ! endif + + ! call mpp_clock_end (gfdl_mp_clock) + if(lradar) then + ! Only set melti to true at the output times + if (reset) then + melti=.true. + else + melti=.false. + endif + do j = js, je + do i = is, ie + do k = ktop,kbot + kflip = kbot-ktop+1-k+1 + t1d(k) = pt(i,j,kflip) + p1d(k) = p(i,j,kflip) + qv1d(k) = qv(i,j,kflip)/(1-qv(i,j,kflip)) + qr1d(k) = qr(i,j,kflip) + qs1d(k) = qs(i,j,kflip) + qg1d(k) = qg(i,j,kflip) + enddo + call refl10cm_gfdl (qv1d, qr1d, qs1d, qg1d, & + t1d, p1d, dBZ, ktop, kbot, i, j, melti) + do k = ktop,kbot + kflip = kbot-ktop+1-k+1 + refl_10cm(i,j,kflip) = MAX(-35., dBZ(k)) + enddo + enddo + enddo + endif + +end subroutine gfdl_cloud_microphys_mod_driver + +! ----------------------------------------------------------------------- +!>\ingroup mod_gfdl_cloud_mp +!>\brief GFDL cloud microphysics, major program, and is based on +!! Lin et al.(1983) \cite lin_et_al_1983 and +!! Rutledge and Hobbs (1984) \cite rutledge_and_hobbs_1984. +!! +!>\section detmpdrv GFDL Cloud mpdrv General Algorithm +subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & + qg, qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, & + rain, snow, graupel, ice, m2_rain, m2_sol, cond, area1, land, & + u_dt, v_dt, pt_dt, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, & + w_var, vt_r, vt_s, vt_g, vt_i, qn2) + + implicit none + + logical, intent (in) :: hydrostatic + + integer, intent (in) :: j, is, ie, js, je, ks, ke + integer, intent (in) :: ntimes, ktop, kbot + + real, intent (in) :: dt_in + + real, intent (in), dimension (is:) :: area1, land + + real, intent (in), dimension (is:, js:, ks:) :: uin, vin, delp, pt, dz + real, intent (in), dimension (is:, js:, ks:) :: qv, ql, qr, qg, qa, qn + + real, intent (inout), dimension (is:, js:, ks:) :: qi, qs + real, intent (inout), dimension (is:, js:, ks:) :: u_dt, v_dt, w, pt_dt, qa_dt + real, intent (inout), dimension (is:, js:, ks:) :: qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt + + real, intent (inout), dimension (is:) :: rain, snow, ice, graupel, cond + + real, intent (out), dimension (is:, js:) :: w_var + + real, intent (out), dimension (is:, js:, ks:) :: vt_r, vt_s, vt_g, vt_i, qn2 + + real, intent (out), dimension (is:, ks:) :: m2_rain, m2_sol + + real, dimension (ktop:kbot) :: qvz, qlz, qrz, qiz, qsz, qgz, qaz + real, dimension (ktop:kbot) :: vtiz, vtsz, vtgz, vtrz + real, dimension (ktop:kbot) :: dp0, dp1, dz0, dz1 + real, dimension (ktop:kbot) :: qv0, ql0, qr0, qi0, qs0, qg0, qa0 + real, dimension (ktop:kbot) :: t0, den, den0, tz, p1, denfac + real, dimension (ktop:kbot) :: ccn, c_praut, m1_rain, m1_sol, m1 + real, dimension (ktop:kbot) :: u0, v0, u1, v1, w1 + + real :: cpaut, rh_adj, rh_rain + real :: r1, s1, i1, g1, rdt, ccn0 + real :: dt_rain, dts + real :: s_leng, t_land, t_ocean, h_var + real :: cvm, tmp, omq + real :: dqi, qio, qin + + integer :: i, k, n + + dts = dt_in / real (ntimes) + dt_rain = dts * 0.5 + rdt = 1. / dt_in + + ! ----------------------------------------------------------------------- + ! use local variables + ! ----------------------------------------------------------------------- + + !GJF: assign values to intent(out) variables that are commented out + w_var = 0.0 + vt_r = 0.0 + vt_s = 0.0 + vt_g = 0.0 + vt_i = 0.0 + qn2 = 0.0 + !GJF + + do i = is, ie + + do k = ktop, kbot + qiz (k) = qi (i, j, k) + qsz (k) = qs (i, j, k) + enddo + + ! ----------------------------------------------------------------------- + !> - Prevent excessive build-up of cloud ice from external sources. + ! ----------------------------------------------------------------------- + + if (de_ice) then + do k = ktop, kbot + qio = qiz (k) - dt_in * qi_dt (i, j, k) ! original qi before phys + qin = max (qio, qi0_max) ! adjusted value + if (qiz (k) > qin) then + qsz (k) = qsz (k) + qiz (k) - qin + qiz (k) = qin + dqi = (qin - qio) * rdt ! modified qi tendency + qs_dt (i, j, k) = qs_dt (i, j, k) + qi_dt (i, j, k) - dqi + qi_dt (i, j, k) = dqi + qi (i, j, k) = qiz (k) + qs (i, j, k) = qsz (k) + endif + enddo + endif + + do k = ktop, kbot + + t0 (k) = pt (i, j, k) + tz (k) = t0 (k) + dp1 (k) = delp (i, j, k) + dp0 (k) = dp1 (k) ! moist air mass * grav + + ! ----------------------------------------------------------------------- + !> - Convert moist mixing ratios to dry mixing ratios. + ! ----------------------------------------------------------------------- + + qvz (k) = qv (i, j, k) + qlz (k) = ql (i, j, k) + qrz (k) = qr (i, j, k) + qgz (k) = qg (i, j, k) + + ! dp1: dry air_mass + ! dp1 (k) = dp1 (k) * (1. - (qvz (k) + qlz (k) + qrz (k) + qiz (k) + qsz (k) + qgz (k))) + dp1 (k) = dp1 (k) * (1. - qvz (k)) ! gfs + omq = dp0 (k) / dp1 (k) + + qvz (k) = qvz (k) * omq + qlz (k) = qlz (k) * omq + qrz (k) = qrz (k) * omq + qiz (k) = qiz (k) * omq + qsz (k) = qsz (k) * omq + qgz (k) = qgz (k) * omq + + qa0 (k) = qa (i, j, k) + qaz (k) = 0. + dz0 (k) = dz (i, j, k) + + den0 (k) = - dp1 (k) / (grav * dz0 (k)) ! density of dry air + p1 (k) = den0 (k) * rdgas * t0 (k) ! dry air pressure + + ! ----------------------------------------------------------------------- + ! save a copy of old value for computing tendencies + ! ----------------------------------------------------------------------- + + qv0 (k) = qvz (k) + ql0 (k) = qlz (k) + qr0 (k) = qrz (k) + qi0 (k) = qiz (k) + qs0 (k) = qsz (k) + qg0 (k) = qgz (k) + + ! ----------------------------------------------------------------------- + ! for sedi_momentum + ! ----------------------------------------------------------------------- + + m1 (k) = 0. + u0 (k) = uin (i, j, k) + v0 (k) = vin (i, j, k) + u1 (k) = u0 (k) + v1 (k) = v0 (k) + + enddo + + if (do_sedi_w) then + do k = ktop, kbot + w1 (k) = w (i, j, k) + enddo + ! Set to -999 if not used + else + w1(:) = -999.0 + endif + + ! ----------------------------------------------------------------------- + !> - Calculate cloud condensation nuclei (ccn), + !! following klein eq. 15 + ! ----------------------------------------------------------------------- + + cpaut = c_paut * 0.104 * grav / 1.717e-5 + + if (prog_ccn) then + do k = ktop, kbot + ! convert # / cc to # / m^3 + ccn (k) = qn (i, j, k) * 1.e6 + c_praut (k) = cpaut * (ccn (k) * rhor) ** (- 1. / 3.) + enddo + use_ccn = .false. + else + ccn0 = (ccn_l * land (i) + ccn_o * (1. - land (i))) * 1.e6 + if (use_ccn) then + ! ----------------------------------------------------------------------- + ! ccn is formulted as ccn = ccn_surface * (den / den_surface) + ! ----------------------------------------------------------------------- + ccn0 = ccn0 * rdgas * tz (kbot) / p1 (kbot) + endif + tmp = cpaut * (ccn0 * rhor) ** (- 1. / 3.) + do k = ktop, kbot + c_praut (k) = tmp + ccn (k) = ccn0 + enddo + endif + + ! ----------------------------------------------------------------------- + !> - Calculate horizontal subgrid variability, which is used in cloud + !! fraction, relative humidity calculation, evaporation and condensation + !! processes. Horizontal sub-grid variability is a function of cell area + !! and land/sea mask: + !!\n Over land: + !!\f[ + !! t_{land}=dw_{land}(\frac{A_{r}}{10^{10}})^{0.25} + !!\f] + !!\n Over ocean: + !!\f[ + !! t_{ocean}=dw_{ocean}(\frac{A_{r}}{10^{10}})^{0.25} + !!\f] + !! where \f$A_{r}\f$ is cell area. \f$dw_{land}=0.16\f$ and \f$dw_{ocean}=0.10\f$ + !! are base value for sub-grid variability over land and ocean. + !! The total horizontal sub-grid variability is: + !!\f[ + !! h_{var}=t_{land}\times fr_{land}+t_{ocean}\times (1-fr_{land}) + !!\f] + !!\f[ + !! h_{var}=min[0.2,max(0.01,h_{var})] + !!\f] + ! total water subgrid deviation in horizontal direction + ! default area dependent form: use dx ~ 100 km as the base + ! ----------------------------------------------------------------------- + + s_leng = sqrt (sqrt (area1 (i) / 1.e10)) + t_land = dw_land * s_leng + t_ocean = dw_ocean * s_leng + h_var = t_land * land (i) + t_ocean * (1. - land (i)) + h_var = min (0.20, max (0.01, h_var)) + ! if (id_var > 0) w_var (i, j) = h_var + + ! ----------------------------------------------------------------------- + !> - Calculate relative humidity increment. + ! ----------------------------------------------------------------------- + + rh_adj = 1. - h_var - rh_inc + rh_rain = max (0.35, rh_adj - rh_inr) ! rh_inr = 0.25 + + ! ----------------------------------------------------------------------- + !> - If requested, call neg_adj() and fix all negative water species. + ! ----------------------------------------------------------------------- + + if (fix_negative) & + call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz) + + m2_rain (i, :) = 0. + m2_sol (i, :) = 0. + + !> - Do loop on cloud microphysics sub time step. + do n = 1, ntimes + + ! ----------------------------------------------------------------------- + !> - Define air density based on hydrostatical property. + ! ----------------------------------------------------------------------- + + if (p_nonhydro) then + do k = ktop, kbot + dz1 (k) = dz0 (k) + den (k) = den0 (k) ! dry air density remains the same + denfac (k) = sqrt (sfcrho / den (k)) + enddo + else + do k = ktop, kbot + dz1 (k) = dz0 (k) * tz (k) / t0 (k) ! hydrostatic balance + den (k) = den0 (k) * dz0 (k) / dz1 (k) + denfac (k) = sqrt (sfcrho / den (k)) + enddo + endif + + ! ----------------------------------------------------------------------- + !> - Call warm_rain() - time-split warm rain processes: 1st pass. + ! ----------------------------------------------------------------------- + + call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, & + qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var) + + rain (i) = rain (i) + r1 + + do k = ktop, kbot + m2_rain (i, k) = m2_rain (i, k) + m1_rain (k) + m1 (k) = m1 (k) + m1_rain (k) + enddo + + ! ----------------------------------------------------------------------- + !> - Sedimentation of cloud ice, snow, and graupel. + ! ----------------------------------------------------------------------- + + !> - Call fall_speed() to calculate the fall velocity of cloud ice, snow + !! and graupel. + call fall_speed (ktop, kbot, den, qsz, qiz, qgz, qlz, tz, vtsz, vtiz, vtgz) + + !> - Call terminal_fall() to calculate the terminal fall speed. + call terminal_fall (dts, ktop, kbot, tz, qvz, qlz, qrz, qgz, qsz, qiz, & + dz1, dp1, den, vtgz, vtsz, vtiz, r1, g1, s1, i1, m1_sol, w1) + + rain (i) = rain (i) + r1 ! from melted snow & ice that reached the ground + snow (i) = snow (i) + s1 + graupel (i) = graupel (i) + g1 + ice (i) = ice (i) + i1 + + ! ----------------------------------------------------------------------- + !> - Call sedi_heat() to calculate heat transportation during sedimentation. + ! ----------------------------------------------------------------------- + + if (do_sedi_heat) & + call sedi_heat (ktop, kbot, dp1, m1_sol, dz1, tz, qvz, qlz, qrz, qiz, & + qsz, qgz, c_ice) + + ! ----------------------------------------------------------------------- + !> - Call warm_rain() to - time-split warm rain processes: 2nd pass + ! ----------------------------------------------------------------------- + + call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, & + qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var) + + rain (i) = rain (i) + r1 + + do k = ktop, kbot + m2_rain (i, k) = m2_rain (i, k) + m1_rain (k) + m2_sol (i, k) = m2_sol (i, k) + m1_sol (k) + m1 (k) = m1 (k) + m1_rain (k) + m1_sol (k) + enddo + + ! ----------------------------------------------------------------------- + !> - Call icloud(): ice-phase microphysics + ! ----------------------------------------------------------------------- + + call icloud (ktop, kbot, tz, p1, qvz, qlz, qrz, qiz, qsz, qgz, dp1, den, & + denfac, vtsz, vtgz, vtrz, qaz, rh_adj, rh_rain, dts, h_var) + + enddo + + ! convert units from Pa*kg/kg to kg/m^2/s + m2_rain (i, :) = m2_rain (i, :) * rdt * rgrav + m2_sol (i, :) = m2_sol (i, :) * rdt * rgrav + + ! ----------------------------------------------------------------------- + !> - Calculate momentum transportation during sedimentation. + ! note: dp1 is dry mass; dp0 is the old moist (total) mass + ! ----------------------------------------------------------------------- + + if (sedi_transport) then + do k = ktop + 1, kbot + u1 (k) = (dp0 (k) * u1 (k) + m1 (k - 1) * u1 (k - 1)) / (dp0 (k) + m1 (k - 1)) + v1 (k) = (dp0 (k) * v1 (k) + m1 (k - 1) * v1 (k - 1)) / (dp0 (k) + m1 (k - 1)) + u_dt (i, j, k) = u_dt (i, j, k) + (u1 (k) - u0 (k)) * rdt + v_dt (i, j, k) = v_dt (i, j, k) + (v1 (k) - v0 (k)) * rdt + enddo + endif + + if (do_sedi_w) then + do k = ktop, kbot + w (i, j, k) = w1 (k) + enddo + endif + + ! ----------------------------------------------------------------------- + !> - Update moist air mass (actually hydrostatic pressure). + ! convert to dry mixing ratios + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + omq = dp1 (k) / dp0 (k) + qv_dt (i, j, k) = qv_dt (i, j, k) + rdt * (qvz (k) - qv0 (k)) * omq + ql_dt (i, j, k) = ql_dt (i, j, k) + rdt * (qlz (k) - ql0 (k)) * omq + qr_dt (i, j, k) = qr_dt (i, j, k) + rdt * (qrz (k) - qr0 (k)) * omq + qi_dt (i, j, k) = qi_dt (i, j, k) + rdt * (qiz (k) - qi0 (k)) * omq + qs_dt (i, j, k) = qs_dt (i, j, k) + rdt * (qsz (k) - qs0 (k)) * omq + qg_dt (i, j, k) = qg_dt (i, j, k) + rdt * (qgz (k) - qg0 (k)) * omq + cvm = c_air + qvz (k) * c_vap + (qrz (k) + qlz (k)) * c_liq + (qiz (k) + qsz (k) + qgz (k)) * c_ice + pt_dt (i, j, k) = pt_dt (i, j, k) + rdt * (tz (k) - t0 (k)) * cvm / cp_air + enddo + + ! ----------------------------------------------------------------------- + !> - Update cloud fraction tendency. + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + if (do_qa) then + qa_dt (i, j, k) = 0. + else + qa_dt (i, j, k) = qa_dt (i, j, k) + rdt * (qaz (k) / real (ntimes) - qa0 (k)) + endif + enddo + + ! ----------------------------------------------------------------------- + ! fms diagnostics: + ! ----------------------------------------------------------------------- + + ! if (id_cond > 0) then + ! do k = ktop, kbot ! total condensate + ! cond (i) = cond (i) + dp1 (k) * (qlz (k) + qrz (k) + qsz (k) + qiz (k) + qgz (k)) + ! enddo + ! endif + ! + ! if (id_vtr > 0) then + ! do k = ktop, kbot + ! vt_r (i, j, k) = vtrz (k) + ! enddo + ! endif + ! + ! if (id_vts > 0) then + ! do k = ktop, kbot + ! vt_s (i, j, k) = vtsz (k) + ! enddo + ! endif + ! + ! if (id_vtg > 0) then + ! do k = ktop, kbot + ! vt_g (i, j, k) = vtgz (k) + ! enddo + ! endif + ! + ! if (id_vts > 0) then + ! do k = ktop, kbot + ! vt_i (i, j, k) = vtiz (k) + ! enddo + ! endif + ! + ! if (id_droplets > 0) then + ! do k = ktop, kbot + ! qn2 (i, j, k) = ccn (k) + ! enddo + ! endif + + enddo + +end subroutine mpdrv + +! ----------------------------------------------------------------------- +!>\ingroup mod_gfdl_cloud_mp +!>\brief This subroutine calculates sedimentation of heat. +subroutine sedi_heat (ktop, kbot, dm, m1, dz, tz, qv, ql, qr, qi, qs, qg, cw) + + implicit none + + ! input q fields are dry mixing ratios, and dm is dry air mass + + integer, intent (in) :: ktop, kbot + + real, intent (in), dimension (ktop:kbot) :: dm, m1, dz, qv, ql, qr, qi, qs, qg + + real, intent (inout), dimension (ktop:kbot) :: tz + + real, intent (in) :: cw ! heat capacity + + real, dimension (ktop:kbot) :: dgz, cvn + + real :: tmp + + integer :: k + + do k = ktop, kbot + dgz (k) = - 0.5 * grav * dz (k) ! > 0 + cvn (k) = dm (k) * (cv_air + qv (k) * cv_vap + (qr (k) + ql (k)) * & + c_liq + (qi (k) + qs (k) + qg (k)) * c_ice) + enddo + + ! ----------------------------------------------------------------------- + ! sjl, july 2014 + ! assumption: the ke in the falling condensates is negligible compared to the potential energy + ! that was unaccounted for. local thermal equilibrium is assumed, and the loss in pe is transformed + ! into internal energy (to heat the whole grid box) + ! backward time - implicit upwind transport scheme: + ! dm here is dry air mass + ! ----------------------------------------------------------------------- + + k = ktop + tmp = cvn (k) + m1 (k) * cw + tz (k) = (tmp * tz (k) + m1 (k) * dgz (k)) / tmp + + ! ----------------------------------------------------------------------- + ! implicit algorithm: can't be vectorized + ! needs an inner i - loop for vectorization + ! ----------------------------------------------------------------------- + + do k = ktop + 1, kbot + tz (k) = ((cvn (k) + cw * (m1 (k) - m1 (k - 1))) * tz (k) + m1 (k - 1) * & + cw * tz (k - 1) + dgz (k) * (m1 (k - 1) + m1 (k))) / (cvn (k) + cw * m1 (k)) + enddo + +end subroutine sedi_heat + +! ----------------------------------------------------------------------- +!>\ingroup mod_gfdl_cloud_mp +!> This subroutine includes warm rain cloud microphysics. +!>\section warm_gen GFDL Cloud warm_rain General Algorithm +subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, & + den, denfac, ccn, c_praut, rh_rain, vtr, r1, m1_rain, w1, h_var) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in) :: dt ! time step (s) + real, intent (in) :: rh_rain, h_var + + real, intent (in), dimension (ktop:kbot) :: dp, dz, den + real, intent (in), dimension (ktop:kbot) :: denfac, ccn, c_praut + + real, intent (inout), dimension (ktop:kbot) :: tz, vtr + real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qi, qs, qg + real, intent (inout), dimension (ktop:kbot) :: m1_rain, w1 + + real, intent (out) :: r1 + + real, parameter :: so3 = 7. / 3. + + real, dimension (ktop:kbot) :: dl, dm + real, dimension (ktop:kbot + 1) :: ze, zt + + real :: sink, dq, qc0, qc + real :: qden + real :: zs = 0. + real :: dt5 + + integer :: k + + ! fall velocity constants: + + real, parameter :: vconr = 2503.23638966667 + real, parameter :: normr = 25132741228.7183 + real, parameter :: thr = 1.e-8 + + logical :: no_fall + + dt5 = 0.5 * dt + + ! ----------------------------------------------------------------------- + !> - Call check_column() to check if the water species is large enough to fall. + ! ----------------------------------------------------------------------- + + m1_rain (:) = 0. + + call check_column (ktop, kbot, qr, no_fall) + + if (no_fall) then + vtr (:) = vf_min + r1 = 0. + else + + ! ----------------------------------------------------------------------- + !> - Calculate fall speed of rain. + ! ----------------------------------------------------------------------- + + if (const_vr) then + vtr (:) = vr_fac ! ifs_2016: 4.0 + else + do k = ktop, kbot + qden = qr (k) * den (k) + if (qr (k) < thr) then + vtr (k) = vr_min + else + vtr (k) = vr_fac * vconr * sqrt (min (10., sfcrho / den (k))) * & + exp (0.2 * log (qden / normr)) + vtr (k) = min (vr_max, max (vr_min, vtr (k))) + endif + enddo + endif + + ze (kbot + 1) = zs + do k = kbot, ktop, - 1 + ze (k) = ze (k + 1) - dz (k) ! dz < 0 + enddo + + ! ----------------------------------------------------------------------- + !> - Call revap_racc(), to calculate evaporation and accretion + !! of rain for the first 1/2 time step. + ! ----------------------------------------------------------------------- + + ! if (.not. fast_sat_adj) & + call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var) + + if (do_sedi_w) then + do k = ktop, kbot + dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k)) + enddo + endif + + ! ----------------------------------------------------------------------- + !> - Calculate mass flux induced by falling rain + !! ( use_ppm =.false, call implicit_fall(): time-implicit monotonic fall scheme.) + ! ----------------------------------------------------------------------- + + if (use_ppm) then + zt (ktop) = ze (ktop) + do k = ktop + 1, kbot + zt (k) = ze (k) - dt5 * (vtr (k - 1) + vtr (k)) + enddo + zt (kbot + 1) = zs - dt * vtr (kbot) + + do k = ktop, kbot + if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min + enddo + call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qr, r1, m1_rain, mono_prof) + else + call implicit_fall (dt, ktop, kbot, ze, vtr, dp, qr, r1, m1_rain) + endif + + ! ----------------------------------------------------------------------- + ! - Calculate vertical velocity transportation during sedimentation. + ! (do_sedi_w =.true. to turn on vertical motion tranport during sedimentation + ! .false. by default) + ! ----------------------------------------------------------------------- + + if (do_sedi_w) then + w1 (ktop) = (dm (ktop) * w1 (ktop) + m1_rain (ktop) * vtr (ktop)) / (dm (ktop) - m1_rain (ktop)) + do k = ktop + 1, kbot + w1 (k) = (dm (k) * w1 (k) - m1_rain (k - 1) * vtr (k - 1) + m1_rain (k) * vtr (k)) & + / (dm (k) + m1_rain (k - 1) - m1_rain (k)) + enddo + endif + + ! ----------------------------------------------------------------------- + !> - Call sedi_heat() to calculate heat transportation during sedimentation. + ! ----------------------------------------------------------------------- + + if (do_sedi_heat) & + call sedi_heat (ktop, kbot, dp, m1_rain, dz, tz, qv, ql, qr, qi, qs, qg, c_liq) + + ! ----------------------------------------------------------------------- + !> - Call revap_racc() to calculate evaporation and accretion + !! of rain for the remaing 1/2 time step. + ! ----------------------------------------------------------------------- + + call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var) + + endif + + ! ----------------------------------------------------------------------- + !> - Auto-conversion + !! (assuming linear subgrid vertical distribution of cloud water + !! following Lin et al. (1994) \cite lin_et_al_1994.) + ! ----------------------------------------------------------------------- + + if (irain_f /= 0) then + + ! ----------------------------------------------------------------------- + ! no subgrid varaibility + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + qc0 = fac_rc * ccn (k) + if (tz (k) > t_wfr) then + if (use_ccn) then + ! ----------------------------------------------------------------------- + ! ccn is formulted as ccn = ccn_surface * (den / den_surface) + ! ----------------------------------------------------------------------- + qc = qc0 + else + qc = qc0 / den (k) + endif + dq = ql (k) - qc + if (dq > 0.) then + sink = min (dq, dt * c_praut (k) * den (k) * exp (so3 * log (ql (k)))) + ql (k) = ql (k) - sink + qr (k) = qr (k) + sink + endif + endif + enddo + + else + + ! ----------------------------------------------------------------------- + !> - Call linear_prof() to calculate vertical subgrid variability of cloud water. + ! ----------------------------------------------------------------------- + + call linear_prof (kbot - ktop + 1, ql (ktop), dl (ktop), z_slope_liq, h_var) + + do k = ktop, kbot + qc0 = fac_rc * ccn (k) + if (tz (k) > t_wfr + dt_fr) then + dl (k) = min (max (1.e-6, dl (k)), 0.5 * ql (k)) + ! -------------------------------------------------------------------- + ! as in klein's gfdl am2 stratiform scheme (with subgrid variations) + ! -------------------------------------------------------------------- + if (use_ccn) then + ! -------------------------------------------------------------------- + ! ccn is formulted as ccn = ccn_surface * (den / den_surface) + ! -------------------------------------------------------------------- + qc = qc0 + else + qc = qc0 / den (k) + endif + dq = 0.5 * (ql (k) + dl (k) - qc) + ! -------------------------------------------------------------------- + ! dq = dl if qc == q_minus = ql - dl + ! dq = 0 if qc == q_plus = ql + dl + ! -------------------------------------------------------------------- + if (dq > 0.) then ! q_plus > qc + ! -------------------------------------------------------------------- + !> - Revised continuous form: linearly decays (with subgrid dl) to zero at qc == ql + dl + ! -------------------------------------------------------------------- + sink = min (1., dq / dl (k)) * dt * c_praut (k) * den (k) * exp (so3 * log (ql (k))) + ql (k) = ql (k) - sink + qr (k) = qr (k) + sink + endif + endif + enddo + endif + +end subroutine warm_rain + +! ----------------------------------------------------------------------- +!>\ingroup mod_gfdl_cloud_mp +!> This subroutine calculates evaporation of rain and accretion of rain. +!!\section gen_ravap GFDL Cloud revap_racc General Algorithm +subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in) :: dt ! time step (s) + real, intent (in) :: rh_rain, h_var + + real, intent (in), dimension (ktop:kbot) :: den, denfac + + real, intent (inout), dimension (ktop:kbot) :: tz, qv, qr, ql, qi, qs, qg + + real, dimension (ktop:kbot) :: lhl, cvm, q_liq, q_sol, lcpk + + real :: dqv, qsat, dqsdt, evap, t2, qden, q_plus, q_minus, sink + real :: qpz, dq, dqh, tin + + integer :: k + + do k = ktop, kbot + + if (tz (k) > t_wfr .and. qr (k) > qrmin) then + + ! ----------------------------------------------------------------------- + !> - Define heat capacity and latent heat coefficient. + ! ----------------------------------------------------------------------- + + lhl (k) = lv00 + d0_vap * tz (k) + q_liq (k) = ql (k) + qr (k) + q_sol (k) = qi (k) + qs (k) + qg (k) + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + lcpk (k) = lhl (k) / cvm (k) + + tin = tz (k) - lcpk (k) * ql (k) ! presence of clouds suppresses the rain evap + qpz = qv (k) + ql (k) + qsat = wqs2 (tin, den (k), dqsdt) + dqh = max (ql (k), h_var * max (qpz, qcmin)) + dqh = min (dqh, 0.2 * qpz) ! new limiter + dqv = qsat - qv (k) ! use this to prevent super - sat the gird box + q_minus = qpz - dqh + q_plus = qpz + dqh + + ! ----------------------------------------------------------------------- + !> - qsat must be > q_minus to activate evaporation; + !! qsat must be < q_plus to activate accretion + ! ----------------------------------------------------------------------- + + ! ----------------------------------------------------------------------- + !> - rain evaporation (\f$evap\f$) + ! ----------------------------------------------------------------------- + + if (dqv > qvmin .and. qsat > q_minus) then + if (qsat > q_plus) then + dq = qsat - qpz + else + ! ----------------------------------------------------------------------- + ! q_minus < qsat < q_plus + ! dq == dqh if qsat == q_minus + ! ----------------------------------------------------------------------- + dq = 0.25 * (q_minus - qsat) ** 2 / dqh + endif + qden = qr (k) * den (k) + t2 = tin * tin + evap = crevp (1) * t2 * dq * (crevp (2) * sqrt (qden) + crevp (3) * & + exp (0.725 * log (qden))) / (crevp (4) * t2 + crevp (5) * qsat * den (k)) + evap = min (qr (k), dt * evap, dqv / (1. + lcpk (k) * dqsdt)) + ! ----------------------------------------------------------------------- + ! alternative minimum evap in dry environmental air + ! sink = min (qr (k), dim (rh_rain * qsat, qv (k)) / (1. + lcpk (k) * dqsdt)) + ! evap = max (evap, sink) + ! ----------------------------------------------------------------------- + qr (k) = qr (k) - evap + qv (k) = qv (k) + evap + q_liq (k) = q_liq (k) - evap + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) - evap * lhl (k) / cvm (k) + endif + + ! ----------------------------------------------------------------------- + !> - accretion: \f$P_{racc}\f$ + ! ----------------------------------------------------------------------- + + ! if (qr (k) > qrmin .and. ql (k) > 1.e-7 .and. qsat < q_plus) then + if (qr (k) > qrmin .and. ql (k) > 1.e-6 .and. qsat < q_minus) then + sink = dt * denfac (k) * cracw * exp (0.95 * log (qr (k) * den (k))) + sink = sink / (1. + sink) * ql (k) + ql (k) = ql (k) - sink + qr (k) = qr (k) + sink + endif + + endif ! warm - rain + enddo + +end subroutine revap_racc + +! ----------------------------------------------------------------------- +!>\ingroup mod_gfdl_cloud_mp +!> Definition of vertical subgrid variability +!! used for cloud ice and cloud water autoconversion. +! qi -- > ql & ql -- > qr +! edges: qe == qbar + / - dm +!>\section gen_linear GFDL cloud linear_prof General Algorithm +subroutine linear_prof (km, q, dm, z_var, h_var) + + implicit none + + integer, intent (in) :: km + + real, intent (in) :: q (km), h_var + + real, intent (out) :: dm (km) + + logical, intent (in) :: z_var + + real :: dq (km) + + integer :: k + + if (z_var) then + do k = 2, km + dq (k) = 0.5 * (q (k) - q (k - 1)) + enddo + dm (1) = 0. + + ! ----------------------------------------------------------------------- + !> - Use twice the strength of the positive definiteness limiter (Lin et al.(1994) + !! \cite lin_et_al_1994). + ! ----------------------------------------------------------------------- + + do k = 2, km - 1 + dm (k) = 0.5 * min (abs (dq (k) + dq (k + 1)), 0.5 * q (k)) + if (dq (k) * dq (k + 1) <= 0.) then + if (dq (k) > 0.) then ! local max + dm (k) = min (dm (k), dq (k), - dq (k + 1)) + else + dm (k) = 0. + endif + endif + enddo + dm (km) = 0. + + ! ----------------------------------------------------------------------- + !> - Impose the background horizontal variability (\f$h_{var}\f$) that + !! is proportional to the value itself. + ! ----------------------------------------------------------------------- + + do k = 1, km + dm (k) = max (dm (k), qvmin, h_var * q (k)) + enddo + else + do k = 1, km + dm (k) = max (qvmin, h_var * q (k)) + enddo + endif + +end subroutine linear_prof + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!> This subroutine includes cloud ice microphysics processes. +!>\author Shian-Jiann Lin, GFDL +!! +!! This scheme is featured with: +!! - bulk cloud microphysics +!! - processes splitting with some un-split sub-grouping +!! - time implicit (when possible) accretion and autoconversion +!>\section det_icloud GFDL icloud Detailed Algorithm +subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, & + den, denfac, vts, vtg, vtr, qak, rh_adj, rh_rain, dts, h_var) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in), dimension (ktop:kbot) :: p1, dp1, den, denfac, vts, vtg, vtr + + real, intent (inout), dimension (ktop:kbot) :: tzk, qvk, qlk, qrk, qik, qsk, qgk, qak + + real, intent (in) :: rh_adj, rh_rain, dts, h_var + + real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, di, lhl, lhi + real, dimension (ktop:kbot) :: cvm, q_liq, q_sol + + real :: rdts, fac_g2v, fac_v2g, fac_i2s, fac_imlt + real :: tz, qv, ql, qr, qi, qs, qg, melt + real :: pracs, psacw, pgacw, psacr, pgacr, pgaci, praci, psaci + real :: pgmlt, psmlt, pgfr, pgaut, psaut, pgsub + real :: tc, tsq, dqs0, qden, qim, qsm + real :: dt5, factor, sink, qi_crt + real :: tmp, qsw, qsi, dqsdt, dq + real :: dtmp, qc, q_plus, q_minus + + integer :: k + + dt5 = 0.5 * dts + + rdts = 1. / dts + + ! ----------------------------------------------------------------------- + !> - Define conversion scalar/factor. + ! ----------------------------------------------------------------------- + + fac_i2s = 1. - exp (- dts / tau_i2s) + fac_g2v = 1. - exp (- dts / tau_g2v) + fac_v2g = 1. - exp (- dts / tau_v2g) + + fac_imlt = 1. - exp (- dt5 / tau_imlt) + + ! ----------------------------------------------------------------------- + !> - Define heat capacity and latend heat coefficient. + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + lhi (k) = li00 + dc_ice * tzk (k) + q_liq (k) = qlk (k) + qrk (k) + q_sol (k) = qik (k) + qsk (k) + qgk (k) + cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + icpk (k) = lhi (k) / cvm (k) + enddo + + ! ----------------------------------------------------------------------- + ! sources of cloud ice: pihom, cold rain, and the sat_adj + ! (initiation plus deposition) + ! sources of snow: cold rain, auto conversion + accretion (from cloud ice) + ! sat_adj (deposition; requires pre - existing snow) ; initial snow comes from auto conversion + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + if (tzk (k) > tice .and. qik (k) > qcmin) then + + ! ----------------------------------------------------------------------- + !> - Calculate \f$P_{imlt}\f$: instant melting of cloud ice. + ! ----------------------------------------------------------------------- + + melt = min (qik (k), fac_imlt * (tzk (k) - tice) / icpk (k)) + tmp = min (melt, dim (ql_mlt, qlk (k))) ! max ql amount + qlk (k) = qlk (k) + tmp + qrk (k) = qrk (k) + melt - tmp + qik (k) = qik (k) - melt + q_liq (k) = q_liq (k) + melt + q_sol (k) = q_sol (k) - melt + cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tzk (k) = tzk (k) - melt * lhi (k) / cvm (k) + + elseif (tzk (k) < t_wfr .and. qlk (k) > qcmin) then + + ! ----------------------------------------------------------------------- + !> - Calculate \f$P_{ihom}\f$: homogeneous freezing of cloud water into cloud ice. + !! This is the 1st occurance of liquid water freezing in the split MP process. + ! ----------------------------------------------------------------------- + + dtmp = t_wfr - tzk (k) + factor = min (1., dtmp / dt_fr) + sink = min (qlk (k) * factor, dtmp / icpk (k)) + qi_crt = qi_gen * min (qi_lim, 0.1 * (tice - tzk (k))) / den (k) + tmp = min (sink, dim (qi_crt, qik (k))) + qlk (k) = qlk (k) - sink + qsk (k) = qsk (k) + sink - tmp + qik (k) = qik (k) + tmp + q_liq (k) = q_liq (k) - sink + q_sol (k) = q_sol (k) + sink + cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tzk (k) = tzk (k) + sink * lhi (k) / cvm (k) + + endif + enddo + + ! ----------------------------------------------------------------------- + !> - Call linear_prof() to calculate vertical subgrid variability of cloud ice. + ! ----------------------------------------------------------------------- + + call linear_prof (kbot - ktop + 1, qik (ktop), di (ktop), z_slope_ice, h_var) + + ! ----------------------------------------------------------------------- + !> - Update capacity heat and latend heat coefficient. + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + lhl (k) = lv00 + d0_vap * tzk (k) + lhi (k) = li00 + dc_ice * tzk (k) + lcpk (k) = lhl (k) / cvm (k) + icpk (k) = lhi (k) / cvm (k) + tcpk (k) = lcpk (k) + icpk (k) + enddo + + do k = ktop, kbot + + ! ----------------------------------------------------------------------- + ! do nothing above p_min + ! ----------------------------------------------------------------------- + + if (p1 (k) < p_min) cycle + + tz = tzk (k) + qv = qvk (k) + ql = qlk (k) + qi = qik (k) + qr = qrk (k) + qs = qsk (k) + qg = qgk (k) + + pgacr = 0. + pgacw = 0. + tc = tz - tice + + if (tc .ge. 0.) then + + ! ----------------------------------------------------------------------- + !> - Melting of snow: + ! ----------------------------------------------------------------------- + + dqs0 = ces0 / p1 (k) - qv + + if (qs > qcmin) then + + ! ----------------------------------------------------------------------- + !> - \f$P_{sacw}\f$: accretion of cloud water by snow + ! only rate is used (for snow melt) since tc > 0. + ! ----------------------------------------------------------------------- + + if (ql > qrmin) then + factor = denfac (k) * csacw * exp (0.8125 * log (qs * den (k))) + psacw = factor / (1. + dts * factor) * ql ! rate + else + psacw = 0. + endif + + ! ----------------------------------------------------------------------- + !> - \f$P_{sacr}\f$: accretion of rain by melted snow + !> - \f$P_{racs}\f$: accretion of snow by rain + ! ----------------------------------------------------------------------- + + if (qr > qrmin) then + psacr = min (acr3d (vts (k), vtr (k), qr, qs, csacr, acco (1, 2), & + den (k)), qr * rdts) + pracs = acr3d (vtr (k), vts (k), qs, qr, cracs, acco (1, 1), den (k)) + else + psacr = 0. + pracs = 0. + endif + + ! ----------------------------------------------------------------------- + !> - Total snow sink: + !! \f$P_{smlt}\f$: snow melt (smlt(); due to rain accretion) + ! ----------------------------------------------------------------------- + + psmlt = max (0., smlt (tc, dqs0, qs * den (k), psacw, psacr, csmlt, & + den (k), denfac (k))) + sink = min (qs, dts * (psmlt + pracs), tc / icpk (k)) + qs = qs - sink + ! sjl, 20170321: + tmp = min (sink, dim (qs_mlt, ql)) ! max ql due to snow melt + ql = ql + tmp + qr = qr + sink - tmp + ! qr = qr + sink + ! sjl, 20170321: + q_liq (k) = q_liq (k) + sink + q_sol (k) = q_sol (k) - sink + cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz = tz - sink * lhi (k) / cvm (k) + tc = tz - tice + + endif + + ! ----------------------------------------------------------------------- + !> - Update capacity heat and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhi (k) = li00 + dc_ice * tz + icpk (k) = lhi (k) / cvm (k) + + ! ----------------------------------------------------------------------- + !> - Melting of graupel: + ! ----------------------------------------------------------------------- + + if (qg > qcmin .and. tc > 0.) then + + ! ----------------------------------------------------------------------- + !> - \f$P_{gacr}\f$: accretion of rain by graupel + ! ----------------------------------------------------------------------- + + if (qr > qrmin) & + pgacr = min (acr3d (vtg (k), vtr (k), qr, qg, cgacr, acco (1, 3), & + den (k)), rdts * qr) + + ! ----------------------------------------------------------------------- + !> - \f$P_{gacw}\f$: accretion of cloud water by graupel + ! ----------------------------------------------------------------------- + + qden = qg * den (k) + if (ql > qrmin) then + factor = cgacw * qden / sqrt (den (k) * sqrt (sqrt (qden))) + pgacw = factor / (1. + dts * factor) * ql ! rate + endif + + ! ----------------------------------------------------------------------- + !> - \f$P_{gmlt}\f$: graupel melt (gmlt()) + ! ----------------------------------------------------------------------- + + pgmlt = dts * gmlt (tc, dqs0, qden, pgacw, pgacr, cgmlt, den (k)) + pgmlt = min (max (0., pgmlt), qg, tc / icpk (k)) + qg = qg - pgmlt + qr = qr + pgmlt + q_liq (k) = q_liq (k) + pgmlt + q_sol (k) = q_sol (k) - pgmlt + cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz = tz - pgmlt * lhi (k) / cvm (k) + + endif + + else + + ! ----------------------------------------------------------------------- + !> - Cloud ice processes: + ! ----------------------------------------------------------------------- + + ! ----------------------------------------------------------------------- + !> - \f$P_{saci}\f$: accretion of cloud ice by snow + ! ----------------------------------------------------------------------- + + if (qi > 3.e-7) then ! cloud ice sink terms + + if (qs > 1.e-7) then + ! ----------------------------------------------------------------------- + ! sjl added (following lin eq. 23) the temperature dependency + ! to reduce accretion, use esi = exp (0.05 * tc) as in hong et al 2004 + ! ----------------------------------------------------------------------- + factor = dts * denfac (k) * csaci * exp (0.05 * tc + 0.8125 * log (qs * den (k))) + psaci = factor / (1. + factor) * qi + else + psaci = 0. + endif + + ! ----------------------------------------------------------------------- + !> - \f$P_{saut}\f$: autoconversion: cloud ice \f$\rightarrow\f$ snow. + !!\n similar to Lin et al.(1983) \cite lin_et_al_1983 eq. 21 solved implicitly; + !! threshold from wsm6 scheme, Hong et al. (2004) \cite hong_et_al_2004, + !! eq (13) : qi0_crt ~0.8e-4. + ! ----------------------------------------------------------------------- + if (qi0_crt < 0.) then + qim = - qi0_crt + else + qim = qi0_crt / den (k) + endif + ! ----------------------------------------------------------------------- + ! assuming linear subgrid vertical distribution of cloud ice + ! the mismatch computation following lin et al. 1994, mwr + ! ----------------------------------------------------------------------- + + if (const_vi) then + tmp = fac_i2s + else + tmp = fac_i2s * exp (0.025 * tc) + endif + + di (k) = max (di (k), qrmin) + q_plus = qi + di (k) + if (q_plus > (qim + qrmin)) then + if (qim > (qi - di (k))) then + dq = (0.25 * (q_plus - qim) ** 2) / di (k) + else + dq = qi - qim + endif + psaut = tmp * dq + else + psaut = 0. + endif + ! ----------------------------------------------------------------------- + ! sink is no greater than 75% of qi + ! ----------------------------------------------------------------------- + sink = min (0.75 * qi, psaci + psaut) + qi = qi - sink + qs = qs + sink + + ! ----------------------------------------------------------------------- + !> - \f$P_{gaci}\f$: accretion of cloud ice by graupel + ! ----------------------------------------------------------------------- + + if (qg > 1.e-6) then + ! ----------------------------------------------------------------------- + ! factor = dts * cgaci / sqrt (den (k)) * exp (0.05 * tc + 0.875 * log (qg * den (k))) + ! simplified form: remove temp dependency & set the exponent "0.875" -- > 1 + ! ----------------------------------------------------------------------- + factor = dts * cgaci * sqrt (den (k)) * qg + pgaci = factor / (1. + factor) * qi + qi = qi - pgaci + qg = qg + pgaci + endif + + endif + + ! ----------------------------------------------------------------------- + !> - Cold-rain processes: + ! ----------------------------------------------------------------------- + + ! ----------------------------------------------------------------------- + ! rain to ice, snow, graupel processes: + ! ----------------------------------------------------------------------- + + tc = tz - tice + + if (qr > 1.e-7 .and. tc < 0.) then + + ! ----------------------------------------------------------------------- + ! * sink * terms to qr: psacr + pgfr + ! source terms to qs: psacr + ! source terms to qg: pgfr + ! ----------------------------------------------------------------------- + + ! ----------------------------------------------------------------------- + !> - \f$P_{sacr}\f$: accretion of rain by snow (acr3d()) + ! ----------------------------------------------------------------------- + + if (qs > 1.e-7) then ! if snow exists + psacr = dts * acr3d (vts (k), vtr (k), qr, qs, csacr, acco (1, 2), den (k)) + else + psacr = 0. + endif + + ! ----------------------------------------------------------------------- + !> - \f$P_{gfr}\f$: rain freezing \f$\rightarrow\f$ graupel + ! ----------------------------------------------------------------------- + + pgfr = dts * cgfr (1) / den (k) * (exp (- cgfr (2) * tc) - 1.) * & + exp (1.75 * log (qr * den (k))) + + ! ----------------------------------------------------------------------- + !> - Calculate total sink to \f$q_r\f$ : + !!\n Sink terms to \f$q_r\f$: \f$P_{sacr}+P_{gfr}\f$ + !!\n source term to \f$q_s\f$: \f$P_{sacr}\f$ + !!\n source term to \f$q_g\f$: \f$P_{gfr}\f$ + ! ----------------------------------------------------------------------- + + sink = psacr + pgfr + factor = min (sink, qr, - tc / icpk (k)) / max (sink, qrmin) + + psacr = factor * psacr + pgfr = factor * pgfr + + sink = psacr + pgfr + qr = qr - sink + qs = qs + psacr + qg = qg + pgfr + q_liq (k) = q_liq (k) - sink + q_sol (k) = q_sol (k) + sink + cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz = tz + sink * lhi (k) / cvm (k) + + endif + + ! ----------------------------------------------------------------------- + !> - Update capacity heat and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhi (k) = li00 + dc_ice * tz + icpk (k) = lhi (k) / cvm (k) + + ! ----------------------------------------------------------------------- + !> - Graupel production terms: + ! ----------------------------------------------------------------------- + + if (qs > 1.e-7) then + + ! ----------------------------------------------------------------------- + !> - accretion: snow \f$\rightarrow\f$ graupel (acr3d()) + ! ----------------------------------------------------------------------- + + if (qg > qrmin) then + sink = dts * acr3d (vtg (k), vts (k), qs, qg, cgacs, acco (1, 4), den (k)) + else + sink = 0. + endif + + ! ----------------------------------------------------------------------- + !> - autoconversion: snow \f$\rightarrow\f$ graupel + ! ----------------------------------------------------------------------- + + qsm = qs0_crt / den (k) + if (qs > qsm) then + factor = dts * 1.e-3 * exp (0.09 * (tz - tice)) + sink = sink + factor / (1. + factor) * (qs - qsm) + endif + sink = min (qs, sink) + qs = qs - sink + qg = qg + sink + + endif ! snow existed + + if (qg > 1.e-7 .and. tz < tice0) then + + ! ----------------------------------------------------------------------- + !> - \f$P_{gacw}\f$: accretion of cloud water by graupel + ! ----------------------------------------------------------------------- + + if (ql > 1.e-6) then + qden = qg * den (k) + factor = dts * cgacw * qden / sqrt (den (k) * sqrt (sqrt (qden))) + pgacw = factor / (1. + factor) * ql + else + pgacw = 0. + endif + + ! ----------------------------------------------------------------------- + !> - \f$P_{gacr}\f$: accretion of rain by graupel (acr3d()) + ! ----------------------------------------------------------------------- + + if (qr > 1.e-6) then + pgacr = min (dts * acr3d (vtg (k), vtr (k), qr, qg, cgacr, acco (1, 3), & + den (k)), qr) + else + pgacr = 0. + endif + + sink = pgacr + pgacw + factor = min (sink, dim (tice, tz) / icpk (k)) / max (sink, qrmin) + pgacr = factor * pgacr + pgacw = factor * pgacw + + sink = pgacr + pgacw + qg = qg + sink + qr = qr - pgacr + ql = ql - pgacw + q_liq (k) = q_liq (k) - sink + q_sol (k) = q_sol (k) + sink + cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz = tz + sink * lhi (k) / cvm (k) + + endif + + endif + + tzk (k) = tz + qvk (k) = qv + qlk (k) = ql + qik (k) = qi + qrk (k) = qr + qsk (k) = qs + qgk (k) = qg + + enddo + + ! ----------------------------------------------------------------------- + !> - Call subgrid_z_proc() for subgrid cloud microphysics. + ! ----------------------------------------------------------------------- + + call subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tzk, qvk, & + qlk, qrk, qik, qsk, qgk, qak, h_var, rh_rain) + +end subroutine icloud + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!> This subroutine calculates temperature sentive high vertical resolution processes. +!>\section gen_subz GFDL Cloud subgrid_z_proc General Algorithm +subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tz, qv, & + ql, qr, qi, qs, qg, qa, h_var, rh_rain) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in), dimension (ktop:kbot) :: p1, den, denfac + + real, intent (in) :: dts, rh_adj, h_var, rh_rain + + real, intent (inout), dimension (ktop:kbot) :: tz, qv, ql, qr, qi, qs, qg, qa + + real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, tcp3, lhl, lhi + real, dimension (ktop:kbot) :: cvm, q_liq, q_sol, q_cond + + real :: fac_v2l, fac_l2v + + real :: pidep, qi_crt + + ! ----------------------------------------------------------------------- + ! qstar over water may be accurate only down to - 80 deg c with ~10% uncertainty + ! must not be too large to allow psc + ! ----------------------------------------------------------------------- + + real :: rh, rqi, tin, qsw, qsi, qpz, qstar + real :: dqsdt, dwsdt, dq, dq0, factor, tmp + real :: q_plus, q_minus, dt_evap, dt_pisub + real :: evap, sink, tc, pisub, q_adj, dtmp + real :: pssub, pgsub, tsq, qden, fac_g2v, fac_v2g + + integer :: k + + if (fast_sat_adj) then + dt_evap = 0.5 * dts + else + dt_evap = dts + endif + + ! ----------------------------------------------------------------------- + !> - Define conversion scalar/factor. + ! ----------------------------------------------------------------------- + + fac_v2l = 1. - exp (- dt_evap / tau_v2l) + fac_l2v = 1. - exp (- dt_evap / tau_l2v) + + fac_g2v = 1. - exp (- dts / tau_g2v) + fac_v2g = 1. - exp (- dts / tau_v2g) + + ! ----------------------------------------------------------------------- + !> - Define heat capacity and latend heat coefficient. + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + lhl (k) = lv00 + d0_vap * tz (k) + lhi (k) = li00 + dc_ice * tz (k) + q_liq (k) = ql (k) + qr (k) + q_sol (k) = qi (k) + qs (k) + qg (k) + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + lcpk (k) = lhl (k) / cvm (k) + icpk (k) = lhi (k) / cvm (k) + tcpk (k) = lcpk (k) + icpk (k) + tcp3 (k) = lcpk (k) + icpk (k) * min (1., dim (tice, tz (k)) / (tice - t_wfr)) + enddo + + do k = ktop, kbot + + if (p1 (k) < p_min) cycle + + ! ----------------------------------------------------------------------- + !> - Instant deposit all water vapor to cloud ice when temperature is super low. + ! ----------------------------------------------------------------------- + + if (tz (k) < t_min) then + sink = dim (qv (k), 1.e-7) + qv (k) = qv (k) - sink + qi (k) = qi (k) + sink + q_sol (k) = q_sol (k) + sink + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) + sink * (lhl (k) + lhi (k)) / cvm (k) + if (.not. do_qa) qa (k) = qa (k) + 1. ! air fully saturated; 100 % cloud cover + cycle + endif + + ! ----------------------------------------------------------------------- + !> - Update heat capacity and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhl (k) = lv00 + d0_vap * tz (k) + lhi (k) = li00 + dc_ice * tz (k) + lcpk (k) = lhl (k) / cvm (k) + icpk (k) = lhi (k) / cvm (k) + tcpk (k) = lcpk (k) + icpk (k) + tcp3 (k) = lcpk (k) + icpk (k) * min (1., dim (tice, tz (k)) / (tice - t_wfr)) + + ! ----------------------------------------------------------------------- + !> - Instant evaporation/sublimation of all clouds if rh < rh_adj \f$\rightarrow\f$ cloud free. + ! ----------------------------------------------------------------------- + + qpz = qv (k) + ql (k) + qi (k) + tin = tz (k) - (lhl (k) * (ql (k) + qi (k)) + lhi (k) * qi (k)) / (c_air + & + qpz * c_vap + qr (k) * c_liq + (qs (k) + qg (k)) * c_ice) + if (tin > t_sub + 6.) then + rh = qpz / iqs1 (tin, den (k)) + if (rh < rh_adj) then ! qpz / rh_adj < qs + tz (k) = tin + qv (k) = qpz + ql (k) = 0. + qi (k) = 0. + cycle ! cloud free + endif + endif + + ! ----------------------------------------------------------------------- + !> - cloud water \f$\Leftrightarrow\f$ vapor adjustment: + ! ----------------------------------------------------------------------- + + qsw = wqs2 (tz (k), den (k), dwsdt) + dq0 = qsw - qv (k) + if (dq0 > 0.) then + ! SJL 20170703 added ql factor to prevent the situation of high ql and low RH + ! factor = min (1., fac_l2v * sqrt (max (0., ql (k)) / 1.e-5) * 10. * dq0 / qsw) + ! factor = fac_l2v + ! factor = 1 + factor = min (1., fac_l2v * (10. * dq0 / qsw)) ! the rh dependent factor = 1 at 90% + evap = min (ql (k), factor * dq0 / (1. + tcp3 (k) * dwsdt)) + else ! condensate all excess vapor into cloud water + ! ----------------------------------------------------------------------- + ! evap = fac_v2l * dq0 / (1. + tcp3 (k) * dwsdt) + ! sjl, 20161108 + ! ----------------------------------------------------------------------- + evap = dq0 / (1. + tcp3 (k) * dwsdt) + endif + qv (k) = qv (k) + evap + ql (k) = ql (k) - evap + q_liq (k) = q_liq (k) - evap + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) - evap * lhl (k) / cvm (k) + + ! ----------------------------------------------------------------------- + !> - Update heat capacity and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhi (k) = li00 + dc_ice * tz (k) + icpk (k) = lhi (k) / cvm (k) + + ! ----------------------------------------------------------------------- + !> - Enforce complete freezing below \f$-48^oC\f$. + ! ----------------------------------------------------------------------- + + dtmp = t_wfr - tz (k) ! [ - 40, - 48] + if (dtmp > 0. .and. ql (k) > qcmin) then + sink = min (ql (k), ql (k) * dtmp * 0.125, dtmp / icpk (k)) + ql (k) = ql (k) - sink + qi (k) = qi (k) + sink + q_liq (k) = q_liq (k) - sink + q_sol (k) = q_sol (k) + sink + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) + sink * lhi (k) / cvm (k) + endif + + ! ----------------------------------------------------------------------- + !> - Update heat capacity and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhi (k) = li00 + dc_ice * tz (k) + icpk (k) = lhi (k) / cvm (k) + + ! ----------------------------------------------------------------------- + !> - Apply Bigg mechanism. + ! ----------------------------------------------------------------------- + + if (fast_sat_adj) then + dt_pisub = 0.5 * dts + else + dt_pisub = dts + tc = tice - tz (k) + if (ql (k) > qrmin .and. tc > 0.) then + sink = 3.3333e-10 * dts * (exp (0.66 * tc) - 1.) * den (k) * ql (k) * ql (k) + sink = min (ql (k), tc / icpk (k), sink) + ql (k) = ql (k) - sink + qi (k) = qi (k) + sink + q_liq (k) = q_liq (k) - sink + q_sol (k) = q_sol (k) + sink + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) + sink * lhi (k) / cvm (k) + endif ! significant ql existed + endif + + ! ----------------------------------------------------------------------- + !> - Update capacity heat and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhl (k) = lv00 + d0_vap * tz (k) + lhi (k) = li00 + dc_ice * tz (k) + lcpk (k) = lhl (k) / cvm (k) + icpk (k) = lhi (k) / cvm (k) + tcpk (k) = lcpk (k) + icpk (k) + + ! ----------------------------------------------------------------------- + !> - Sublimation/deposition of ice. + ! ----------------------------------------------------------------------- + + if (tz (k) < tice) then + qsi = iqs2 (tz (k), den (k), dqsdt) + dq = qv (k) - qsi + sink = dq / (1. + tcpk (k) * dqsdt) + if (qi (k) > qrmin) then + ! eq 9, hong et al. 2004, mwr + ! for a and b, see dudhia 1989: page 3103 eq (b7) and (b8) + pidep = dt_pisub * dq * 349138.78 * exp (0.875 * log (qi (k) * den (k))) & + / (qsi * den (k) * lat2 / (0.0243 * rvgas * tz (k) ** 2) + 4.42478e4) + else + pidep = 0. + endif + if (dq > 0.) then ! vapor - > ice + tmp = tice - tz (k) + ! 20160912: the following should produce more ice at higher altitude + ! qi_crt = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp))) / den (k) + qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den (k) + sink = min (sink, max (qi_crt - qi (k), pidep), tmp / tcpk (k)) + else ! ice -- > vapor + pidep = pidep * min (1., dim (tz (k), t_sub) * 0.2) + sink = max (pidep, sink, - qi (k)) + endif + qv (k) = qv (k) - sink + qi (k) = qi (k) + sink + q_sol (k) = q_sol (k) + sink + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) + sink * (lhl (k) + lhi (k)) / cvm (k) + endif + + ! ----------------------------------------------------------------------- + !> - Update capacity heat and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhl (k) = lv00 + d0_vap * tz (k) + lhi (k) = li00 + dc_ice * tz (k) + lcpk (k) = lhl (k) / cvm (k) + icpk (k) = lhi (k) / cvm (k) + tcpk (k) = lcpk (k) + icpk (k) + + ! ----------------------------------------------------------------------- + !> - Sublimation/deposition of snow. + ! this process happens for all temp rage + ! ----------------------------------------------------------------------- + + if (qs (k) > qrmin) then + qsi = iqs2 (tz (k), den (k), dqsdt) + qden = qs (k) * den (k) + tmp = exp (0.65625 * log (qden)) + tsq = tz (k) * tz (k) + dq = (qsi - qv (k)) / (1. + tcpk (k) * dqsdt) + pssub = cssub (1) * tsq * (cssub (2) * sqrt (qden) + cssub (3) * tmp * & + sqrt (denfac (k))) / (cssub (4) * tsq + cssub (5) * qsi * den (k)) + pssub = (qsi - qv (k)) * dts * pssub + if (pssub > 0.) then ! qs -- > qv, sublimation + pssub = min (pssub * min (1., dim (tz (k), t_sub) * 0.2), qs (k)) + else + if (tz (k) > tice) then + pssub = 0. ! no deposition + else + pssub = max (pssub, dq, (tz (k) - tice) / tcpk (k)) + endif + endif + qs (k) = qs (k) - pssub + qv (k) = qv (k) + pssub + q_sol (k) = q_sol (k) - pssub + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) - pssub * (lhl (k) + lhi (k)) / cvm (k) + endif + + ! ----------------------------------------------------------------------- + !> - Update capacity heat and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhl (k) = lv00 + d0_vap * tz (k) + lhi (k) = li00 + dc_ice * tz (k) + lcpk (k) = lhl (k) / cvm (k) + icpk (k) = lhi (k) / cvm (k) + tcpk (k) = lcpk (k) + icpk (k) + + ! ----------------------------------------------------------------------- + !> - Simplified 2-way grapuel sublimation-deposition mechanism. + ! ----------------------------------------------------------------------- + + if (qg (k) > qrmin) then + qsi = iqs2 (tz (k), den (k), dqsdt) + dq = (qv (k) - qsi) / (1. + tcpk (k) * dqsdt) + pgsub = (qv (k) / qsi - 1.) * qg (k) + if (pgsub > 0.) then ! deposition + if (tz (k) > tice) then + pgsub = 0. ! no deposition + else + pgsub = min (fac_v2g * pgsub, 0.2 * dq, ql (k) + qr (k), & + (tice - tz (k)) / tcpk (k)) + endif + else ! submilation + pgsub = max (fac_g2v * pgsub, dq) * min (1., dim (tz (k), t_sub) * 0.1) + endif + qg (k) = qg (k) + pgsub + qv (k) = qv (k) - pgsub + q_sol (k) = q_sol (k) + pgsub + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) + pgsub * (lhl (k) + lhi (k)) / cvm (k) + endif + +#ifdef USE_MIN_EVAP + ! ----------------------------------------------------------------------- + !> - Update capacity heat and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhl (k) = lv00 + d0_vap * tz (k) + lcpk (k) = lhl (k) / cvm (k) + + ! ----------------------------------------------------------------------- + !> - Minimum evap of rain in dry environmental air. + ! ----------------------------------------------------------------------- + + if (qr (k) > qcmin) then + qsw = wqs2 (tz (k), den (k), dqsdt) + sink = min (qr (k), dim (rh_rain * qsw, qv (k)) / (1. + lcpk (k) * dqsdt)) + qv (k) = qv (k) + sink + qr (k) = qr (k) - sink + q_liq (k) = q_liq (k) - sink + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) - sink * lhl (k) / cvm (k) + endif +#endif + + ! ----------------------------------------------------------------------- + !> - Update capacity heat and latend heat coefficient. + ! ----------------------------------------------------------------------- + + lhl (k) = lv00 + d0_vap * tz (k) + cvm (k) = c_air + (qv (k) + q_liq (k) + q_sol (k)) * c_vap + lcpk (k) = lhl (k) / cvm (k) + + ! ----------------------------------------------------------------------- + ! compute cloud fraction + ! ----------------------------------------------------------------------- + + ! ----------------------------------------------------------------------- + !> - Combine water species. + ! ----------------------------------------------------------------------- + + if (do_qa) cycle + + if (rad_snow) then + q_sol (k) = qi (k) + qs (k) + else + q_sol (k) = qi (k) + endif + if (rad_rain) then + q_liq (k) = ql (k) + qr (k) + else + q_liq (k) = ql (k) + endif + q_cond (k) = q_liq (k) + q_sol (k) + + qpz = qv (k) + q_cond (k) ! qpz is conserved + + ! ----------------------------------------------------------------------- + !> - Use the "liquid-frozen water temperature" (tin) to compute saturated specific humidity. + ! ----------------------------------------------------------------------- + + tin = tz (k) - (lcpk (k) * q_cond (k) + icpk (k) * q_sol (k)) ! minimum temperature + ! tin = tz (k) - ((lv00 + d0_vap * tz (k)) * q_cond (k) + & + ! (li00 + dc_ice * tz (k)) * q_sol (k)) / (c_air + qpz * c_vap) + + ! ----------------------------------------------------------------------- + ! determine saturated specific humidity + ! ----------------------------------------------------------------------- + + if (tin <= t_wfr) then + ! ice phase: + qstar = iqs1 (tin, den (k)) + elseif (tin >= tice) then + ! liquid phase: + qstar = wqs1 (tin, den (k)) + else + ! mixed phase: + qsi = iqs1 (tin, den (k)) + qsw = wqs1 (tin, den (k)) + if (q_cond (k) > 3.e-6) then + rqi = q_sol (k) / q_cond (k) + else + ! ----------------------------------------------------------------------- + ! mostly liquid water q_cond (k) at initial cloud development stage + ! ----------------------------------------------------------------------- + rqi = (tice - tin) / (tice - t_wfr) + endif + qstar = rqi * qsi + (1. - rqi) * qsw + endif + + ! ----------------------------------------------------------------------- + !> - Compute cloud fraction, assuming subgrid linear distribution in horizontal; + !! this is effectively a smoother for the binary cloud scheme. + ! ----------------------------------------------------------------------- + + if (qpz > qrmin) then + ! partial cloudiness by pdf: + dq = max (qcmin, h_var * qpz) + q_plus = qpz + dq ! cloud free if qstar > q_plus + q_minus = qpz - dq + if (qstar < q_minus) then + qa (k) = qa (k) + 1. ! air fully saturated; 100 % cloud cover + elseif (qstar < q_plus .and. q_cond (k) > qc_crt) then + qa (k) = qa (k) + (q_plus - qstar) / (dq + dq) ! partial cloud cover + ! qa (k) = sqrt (qa (k) + (q_plus - qstar) / (dq + dq)) + endif + endif + + enddo + +end subroutine subgrid_z_proc + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!> This subroutine calculates rain evaporation. +subroutine revap_rac1 (hydrostatic, is, ie, dt, tz, qv, ql, qr, qi, qs, qg, den, hvar) + + implicit none + + logical, intent (in) :: hydrostatic + + integer, intent (in) :: is, ie + + real, intent (in) :: dt ! time step (s) + + real, intent (in), dimension (is:ie) :: den, hvar, qi, qs, qg + + real, intent (inout), dimension (is:ie) :: tz, qv, qr, ql + + real, dimension (is:ie) :: lcp2, denfac, q_liq, q_sol, cvm, lhl + + real :: dqv, qsat, dqsdt, evap, qden, q_plus, q_minus, sink + real :: tin, t2, qpz, dq, dqh + + integer :: i + + ! ----------------------------------------------------------------------- + ! define latend heat coefficient + ! ----------------------------------------------------------------------- + + do i = is, ie + lhl (i) = lv00 + d0_vap * tz (i) + q_liq (i) = ql (i) + qr (i) + q_sol (i) = qi (i) + qs (i) + qg (i) + cvm (i) = c_air + qv (i) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + lcp2 (i) = lhl (i) / cvm (i) + ! denfac (i) = sqrt (sfcrho / den (i)) + enddo + + do i = is, ie + if (qr (i) > qrmin .and. tz (i) > t_wfr) then + qpz = qv (i) + ql (i) + tin = tz (i) - lcp2 (i) * ql (i) ! presence of clouds suppresses the rain evap + qsat = wqs2 (tin, den (i), dqsdt) + dqh = max (ql (i), hvar (i) * max (qpz, qcmin)) + dqv = qsat - qv (i) + q_minus = qpz - dqh + q_plus = qpz + dqh + + ! ----------------------------------------------------------------------- + ! qsat must be > q_minus to activate evaporation + ! qsat must be < q_plus to activate accretion + ! ----------------------------------------------------------------------- + + ! ----------------------------------------------------------------------- + ! rain evaporation + ! ----------------------------------------------------------------------- + + if (dqv > qvmin .and. qsat > q_minus) then + if (qsat > q_plus) then + dq = qsat - qpz + else + ! q_minus < qsat < q_plus + ! dq == dqh if qsat == q_minus + dq = 0.25 * (q_minus - qsat) ** 2 / dqh + endif + qden = qr (i) * den (i) + t2 = tin * tin + evap = crevp (1) * t2 * dq * (crevp (2) * sqrt (qden) + crevp (3) * exp (0.725 * log (qden))) & + / (crevp (4) * t2 + crevp (5) * qsat * den (i)) + evap = min (qr (i), dt * evap, dqv / (1. + lcp2 (i) * dqsdt)) + qr (i) = qr (i) - evap + qv (i) = qv (i) + evap + q_liq (i) = q_liq (i) - evap + cvm (i) = c_air + qv (i) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice + tz (i) = tz (i) - evap * lhl (i) / cvm (i) + endif + + ! ----------------------------------------------------------------------- + ! accretion: pracc + ! ----------------------------------------------------------------------- + + if (qr (i) > qrmin .and. ql (i) > 1.e-8 .and. qsat < q_plus) then + denfac (i) = sqrt (sfcrho / den (i)) + sink = dt * denfac (i) * cracw * exp (0.95 * log (qr (i) * den (i))) + sink = sink / (1. + sink) * ql (i) + ql (i) = ql (i) - sink + qr (i) = qr (i) + sink + endif + endif + enddo + +end subroutine revap_rac1 + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!>@brief The subroutine 'terminal_fall' computes terminal fall speed. +!>@details It considers cloud ice, snow, and graupel's melting during fall. +subroutine terminal_fall (dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, & + den, vtg, vts, vti, r1, g1, s1, i1, m1_sol, w1) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in) :: dtm ! time step (s) + + real, intent (in), dimension (ktop:kbot) :: vtg, vts, vti, den, dp, dz + + real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qg, qs, qi, tz, m1_sol, w1 + + real, intent (out) :: r1, g1, s1, i1 + + real, dimension (ktop:kbot + 1) :: ze, zt + + real :: qsat, dqsdt, dt5, evap, dtime + real :: factor, frac + real :: tmp, precip, tc, sink + + real, dimension (ktop:kbot) :: lcpk, icpk, cvm, q_liq, q_sol, lhl, lhi + real, dimension (ktop:kbot) :: m1, dm + + real :: zs = 0. + real :: fac_imlt + + integer :: k, k0, m + + logical :: no_fall + + dt5 = 0.5 * dtm + fac_imlt = 1. - exp (- dt5 / tau_imlt) + + ! ----------------------------------------------------------------------- + ! define heat capacity and latend heat coefficient + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + m1_sol (k) = 0. + lhl (k) = lv00 + d0_vap * tz (k) + lhi (k) = li00 + dc_ice * tz (k) + q_liq (k) = ql (k) + qr (k) + q_sol (k) = qi (k) + qs (k) + qg (k) + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + lcpk (k) = lhl (k) / cvm (k) + icpk (k) = lhi (k) / cvm (k) + enddo + + ! ----------------------------------------------------------------------- + ! find significant melting level + ! ----------------------------------------------------------------------- + + k0 = kbot + do k = ktop, kbot - 1 + if (tz (k) > tice) then + k0 = k + exit + endif + enddo + + ! ----------------------------------------------------------------------- + ! melting of cloud_ice (before fall) : + ! ----------------------------------------------------------------------- + + do k = k0, kbot + tc = tz (k) - tice + if (qi (k) > qcmin .and. tc > 0.) then + sink = min (qi (k), fac_imlt * tc / icpk (k)) + tmp = min (sink, dim (ql_mlt, ql (k))) + ql (k) = ql (k) + tmp + qr (k) = qr (k) + sink - tmp + qi (k) = qi (k) - sink + q_liq (k) = q_liq (k) + sink + q_sol (k) = q_sol (k) - sink + cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice + tz (k) = tz (k) - sink * lhi (k) / cvm (k) + tc = tz (k) - tice + endif + enddo + + ! ----------------------------------------------------------------------- + ! turn off melting when cloud microphysics time step is small + ! ----------------------------------------------------------------------- + + if (dtm < 60.) k0 = kbot + + ! sjl, turn off melting of falling cloud ice, snow and graupel + k0 = kbot + ! sjl, turn off melting of falling cloud ice, snow and graupel + + ze (kbot + 1) = zs + do k = kbot, ktop, - 1 + ze (k) = ze (k + 1) - dz (k) ! dz < 0 + enddo + + zt (ktop) = ze (ktop) + + ! ----------------------------------------------------------------------- + ! update capacity heat and latend heat coefficient + ! ----------------------------------------------------------------------- + + do k = k0, kbot + lhi (k) = li00 + dc_ice * tz (k) + icpk (k) = lhi (k) / cvm (k) + enddo + + ! ----------------------------------------------------------------------- + ! melting of falling cloud ice into rain + ! ----------------------------------------------------------------------- + + call check_column (ktop, kbot, qi, no_fall) + + if (vi_fac < 1.e-5 .or. no_fall) then + i1 = 0. + else + + do k = ktop + 1, kbot + zt (k) = ze (k) - dt5 * (vti (k - 1) + vti (k)) + enddo + zt (kbot + 1) = zs - dtm * vti (kbot) + + do k = ktop, kbot + if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min + enddo + + if (k0 < kbot) then + do k = kbot - 1, k0, - 1 + if (qi (k) > qrmin) then + do m = k + 1, kbot + if (zt (k + 1) >= ze (m)) exit + if (zt (k) < ze (m + 1) .and. tz (m) > tice) then + dtime = min (1.0, (ze (m) - ze (m + 1)) / (max (vr_min, vti (k)) * tau_imlt)) + sink = min (qi (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m)) + tmp = min (sink, dim (ql_mlt, ql (m))) + ql (m) = ql (m) + tmp + qr (m) = qr (m) - tmp + sink + tz (m) = tz (m) - sink * icpk (m) + qi (k) = qi (k) - sink * dp (m) / dp (k) + endif + enddo + endif + enddo + endif + + if (do_sedi_w) then + do k = ktop, kbot + dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k)) + enddo + endif + + if (use_ppm) then + call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qi, i1, m1_sol, mono_prof) + else + call implicit_fall (dtm, ktop, kbot, ze, vti, dp, qi, i1, m1_sol) + endif + + if (do_sedi_w) then + w1 (ktop) = (dm (ktop) * w1 (ktop) + m1_sol (ktop) * vti (ktop)) / (dm (ktop) - m1_sol (ktop)) + do k = ktop + 1, kbot + w1 (k) = (dm (k) * w1 (k) - m1_sol (k - 1) * vti (k - 1) + m1_sol (k) * vti (k)) & + / (dm (k) + m1_sol (k - 1) - m1_sol (k)) + enddo + endif + + endif + + ! ----------------------------------------------------------------------- + ! melting of falling snow into rain + ! ----------------------------------------------------------------------- + + r1 = 0. + + call check_column (ktop, kbot, qs, no_fall) + + if (no_fall) then + s1 = 0. + else + + do k = ktop + 1, kbot + zt (k) = ze (k) - dt5 * (vts (k - 1) + vts (k)) + enddo + zt (kbot + 1) = zs - dtm * vts (kbot) + + do k = ktop, kbot + if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min + enddo + + if (k0 < kbot) then + do k = kbot - 1, k0, - 1 + if (qs (k) > qrmin) then + do m = k + 1, kbot + if (zt (k + 1) >= ze (m)) exit + dtime = min (dtm, (ze (m) - ze (m + 1)) / (vr_min + vts (k))) + if (zt (k) < ze (m + 1) .and. tz (m) > tice) then + dtime = min (1.0, dtime / tau_smlt) + sink = min (qs (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m)) + tz (m) = tz (m) - sink * icpk (m) + qs (k) = qs (k) - sink * dp (m) / dp (k) + if (zt (k) < zs) then + r1 = r1 + sink * dp (m) ! precip as rain + else + ! qr source here will fall next time step (therefore, can evap) + qr (m) = qr (m) + sink + endif + endif + if (qs (k) < qrmin) exit + enddo + endif + enddo + endif + + if (do_sedi_w) then + do k = ktop, kbot + dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k)) + enddo + endif + + if (use_ppm) then + call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qs, s1, m1, mono_prof) + else + call implicit_fall (dtm, ktop, kbot, ze, vts, dp, qs, s1, m1) + endif + + do k = ktop, kbot + m1_sol (k) = m1_sol (k) + m1 (k) + enddo + + if (do_sedi_w) then + w1 (ktop) = (dm (ktop) * w1 (ktop) + m1 (ktop) * vts (ktop)) / (dm (ktop) - m1 (ktop)) + do k = ktop + 1, kbot + w1 (k) = (dm (k) * w1 (k) - m1 (k - 1) * vts (k - 1) + m1 (k) * vts (k)) & + / (dm (k) + m1 (k - 1) - m1 (k)) + enddo + endif + + endif + + ! ---------------------------------------------- + ! melting of falling graupel into rain + ! ---------------------------------------------- + + call check_column (ktop, kbot, qg, no_fall) + + if (no_fall) then + g1 = 0. + else + + do k = ktop + 1, kbot + zt (k) = ze (k) - dt5 * (vtg (k - 1) + vtg (k)) + enddo + zt (kbot + 1) = zs - dtm * vtg (kbot) + + do k = ktop, kbot + if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min + enddo + + if (k0 < kbot) then + do k = kbot - 1, k0, - 1 + if (qg (k) > qrmin) then + do m = k + 1, kbot + if (zt (k + 1) >= ze (m)) exit + dtime = min (dtm, (ze (m) - ze (m + 1)) / vtg (k)) + if (zt (k) < ze (m + 1) .and. tz (m) > tice) then + dtime = min (1., dtime / tau_g2r) + sink = min (qg (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m)) + tz (m) = tz (m) - sink * icpk (m) + qg (k) = qg (k) - sink * dp (m) / dp (k) + if (zt (k) < zs) then + r1 = r1 + sink * dp (m) + else + qr (m) = qr (m) + sink + endif + endif + if (qg (k) < qrmin) exit + enddo + endif + enddo + endif + + if (do_sedi_w) then + do k = ktop, kbot + dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k)) + enddo + endif + + if (use_ppm) then + call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qg, g1, m1, mono_prof) + else + call implicit_fall (dtm, ktop, kbot, ze, vtg, dp, qg, g1, m1) + endif + + do k = ktop, kbot + m1_sol (k) = m1_sol (k) + m1 (k) + enddo + + if (do_sedi_w) then + w1 (ktop) = (dm (ktop) * w1 (ktop) + m1 (ktop) * vtg (ktop)) / (dm (ktop) - m1 (ktop)) + do k = ktop + 1, kbot + w1 (k) = (dm (k) * w1 (k) - m1 (k - 1) * vtg (k - 1) + m1 (k) * vtg (k)) & + / (dm (k) + m1 (k - 1) - m1 (k)) + enddo + endif + + endif + +end subroutine terminal_fall + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!>@brief The subroutine 'check_column' checks +!! if the water species is large enough to fall. +subroutine check_column (ktop, kbot, q, no_fall) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in) :: q (ktop:kbot) + + logical, intent (out) :: no_fall + + integer :: k + + no_fall = .true. + + do k = ktop, kbot + if (q (k) > qrmin) then + no_fall = .false. + exit + endif + enddo + +end subroutine check_column + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!>@brief The subroutine computes the time-implicit monotonic +!! fall scheme. +!>@author Shian-Jiann Lin, 2016 +subroutine implicit_fall (dt, ktop, kbot, ze, vt, dp, q, precip, m1) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in) :: dt + + real, intent (in), dimension (ktop:kbot + 1) :: ze + + real, intent (in), dimension (ktop:kbot) :: vt, dp + + real, intent (inout), dimension (ktop:kbot) :: q + + real, intent (out), dimension (ktop:kbot) :: m1 + + real, intent (out) :: precip + + real, dimension (ktop:kbot) :: dz, qm, dd + + integer :: k + + do k = ktop, kbot + dz (k) = ze (k) - ze (k + 1) + dd (k) = dt * vt (k) + q (k) = q (k) * dp (k) + enddo + + ! ----------------------------------------------------------------------- + ! sedimentation: non - vectorizable loop + ! ----------------------------------------------------------------------- + + qm (ktop) = q (ktop) / (dz (ktop) + dd (ktop)) + do k = ktop + 1, kbot + qm (k) = (q (k) + dd (k - 1) * qm (k - 1)) / (dz (k) + dd (k)) + enddo + + ! ----------------------------------------------------------------------- + ! qm is density at this stage + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + qm (k) = qm (k) * dz (k) + enddo + + ! ----------------------------------------------------------------------- + ! output mass fluxes: non - vectorizable loop + ! ----------------------------------------------------------------------- + + m1 (ktop) = q (ktop) - qm (ktop) + do k = ktop + 1, kbot + m1 (k) = m1 (k - 1) + q (k) - qm (k) + enddo + precip = m1 (kbot) + + ! ----------------------------------------------------------------------- + ! update: + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + q (k) = qm (k) / dp (k) + enddo + +end subroutine implicit_fall + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! Lagrangian scheme +!> \author S.J. Lin +subroutine lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in) :: zs + + logical, intent (in) :: mono + + real, intent (in), dimension (ktop:kbot + 1) :: ze, zt + + real, intent (in), dimension (ktop:kbot) :: dp + + ! m1: flux + real, intent (inout), dimension (ktop:kbot) :: q, m1 + + real, intent (out) :: precip + + real, dimension (ktop:kbot) :: qm, dz + + real :: a4 (4, ktop:kbot) + + real :: pl, pr, delz, esl + + integer :: k, k0, n, m + + real, parameter :: r3 = 1. / 3., r23 = 2. / 3. + + ! ----------------------------------------------------------------------- + ! density: + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + dz (k) = zt (k) - zt (k + 1) ! note: dz is positive + q (k) = q (k) * dp (k) + a4 (1, k) = q (k) / dz (k) + qm (k) = 0. + enddo + + ! ----------------------------------------------------------------------- + ! construct vertical profile with zt as coordinate + ! ----------------------------------------------------------------------- + + call cs_profile (a4 (1, ktop), dz (ktop), kbot - ktop + 1, mono) + + k0 = ktop + do k = ktop, kbot + do n = k0, kbot + if (ze (k) <= zt (n) .and. ze (k) >= zt (n + 1)) then + pl = (zt (n) - ze (k)) / dz (n) + if (zt (n + 1) <= ze (k + 1)) then + ! entire new grid is within the original grid + pr = (zt (n) - ze (k + 1)) / dz (n) + qm (k) = a4 (2, n) + 0.5 * (a4 (4, n) + a4 (3, n) - a4 (2, n)) * (pr + pl) - & + a4 (4, n) * r3 * (pr * (pr + pl) + pl ** 2) + qm (k) = qm (k) * (ze (k) - ze (k + 1)) + k0 = n + goto 555 + else + qm (k) = (ze (k) - zt (n + 1)) * (a4 (2, n) + 0.5 * (a4 (4, n) + & + a4 (3, n) - a4 (2, n)) * (1. + pl) - a4 (4, n) * (r3 * (1. + pl * (1. + pl)))) + if (n < kbot) then + do m = n + 1, kbot + ! locate the bottom edge: ze (k + 1) + if (ze (k + 1) < zt (m + 1)) then + qm (k) = qm (k) + q (m) + else + delz = zt (m) - ze (k + 1) + esl = delz / dz (m) + qm (k) = qm (k) + delz * (a4 (2, m) + 0.5 * esl * & + (a4 (3, m) - a4 (2, m) + a4 (4, m) * (1. - r23 * esl))) + k0 = m + goto 555 + endif + enddo + endif + goto 555 + endif + endif + enddo + 555 continue + enddo + + m1 (ktop) = q (ktop) - qm (ktop) + do k = ktop + 1, kbot + m1 (k) = m1 (k - 1) + q (k) - qm (k) + enddo + precip = m1 (kbot) + + ! convert back to * dry * mixing ratio: + ! dp must be dry air_mass (because moist air mass will be changed due to terminal fall) . + + do k = ktop, kbot + q (k) = qm (k) / dp (k) + enddo + +end subroutine lagrangian_fall_ppm + +!>\ingroup mod_gfdl_cloud_mp +subroutine cs_profile (a4, del, km, do_mono) + + implicit none + + integer, intent (in) :: km ! vertical dimension + + real, intent (in) :: del (km) + + logical, intent (in) :: do_mono + + real, intent (inout) :: a4 (4, km) + + real, parameter :: qp_min = 1.e-6 + + real :: gam (km) + real :: q (km + 1) + real :: d4, bet, a_bot, grat, pmp, lac + real :: pmp_1, lac_1, pmp_2, lac_2 + real :: da1, da2, a6da + + integer :: k + + logical extm (km) + + grat = del (2) / del (1) ! grid ratio + bet = grat * (grat + 0.5) + q (1) = (2. * grat * (grat + 1.) * a4 (1, 1) + a4 (1, 2)) / bet + gam (1) = (1. + grat * (grat + 1.5)) / bet + + do k = 2, km + d4 = del (k - 1) / del (k) + bet = 2. + 2. * d4 - gam (k - 1) + q (k) = (3. * (a4 (1, k - 1) + d4 * a4 (1, k)) - q (k - 1)) / bet + gam (k) = d4 / bet + enddo + + a_bot = 1. + d4 * (d4 + 1.5) + q (km + 1) = (2. * d4 * (d4 + 1.) * a4 (1, km) + a4 (1, km - 1) - a_bot * q (km)) & + / (d4 * (d4 + 0.5) - a_bot * gam (km)) + + do k = km, 1, - 1 + q (k) = q (k) - gam (k) * q (k + 1) + enddo + + ! ----------------------------------------------------------------------- + ! apply constraints + ! ----------------------------------------------------------------------- + + do k = 2, km + gam (k) = a4 (1, k) - a4 (1, k - 1) + enddo + + ! ----------------------------------------------------------------------- + ! apply large - scale constraints to all fields if not local max / min + ! ----------------------------------------------------------------------- + + ! ----------------------------------------------------------------------- + ! top: + ! ----------------------------------------------------------------------- + + q (1) = max (q (1), 0.) + q (2) = min (q (2), max (a4 (1, 1), a4 (1, 2))) + q (2) = max (q (2), min (a4 (1, 1), a4 (1, 2)), 0.) + + ! ----------------------------------------------------------------------- + ! interior: + ! ----------------------------------------------------------------------- + + do k = 3, km - 1 + if (gam (k - 1) * gam (k + 1) > 0.) then + q (k) = min (q (k), max (a4 (1, k - 1), a4 (1, k))) + q (k) = max (q (k), min (a4 (1, k - 1), a4 (1, k))) + else + if (gam (k - 1) > 0.) then + ! there exists a local max + q (k) = max (q (k), min (a4 (1, k - 1), a4 (1, k))) + else + ! there exists a local min + q (k) = min (q (k), max (a4 (1, k - 1), a4 (1, k))) + q (k) = max (q (k), 0.0) + endif + endif + enddo + + ! ----------------------------------------------------------------------- + ! bottom : + ! ----------------------------------------------------------------------- + + q (km) = min (q (km), max (a4 (1, km - 1), a4 (1, km))) + q (km) = max (q (km), min (a4 (1, km - 1), a4 (1, km)), 0.) + ! q (km + 1) = max (q (km + 1), 0.) + + ! ----------------------------------------------------------------------- + ! f (s) = al + s * [ (ar - al) + a6 * (1 - s) ] (0 <= s <= 1) + ! ----------------------------------------------------------------------- + + do k = 1, km - 1 + a4 (2, k) = q (k) + a4 (3, k) = q (k + 1) + enddo + + do k = 2, km - 1 + if (gam (k) * gam (k + 1) > 0.0) then + extm (k) = .false. + else + extm (k) = .true. + endif + enddo + + if (do_mono) then + do k = 3, km - 2 + if (extm (k)) then + ! positive definite constraint only if true local extrema + if (a4 (1, k) < qp_min .or. extm (k - 1) .or. extm (k + 1)) then + a4 (2, k) = a4 (1, k) + a4 (3, k) = a4 (1, k) + endif + else + a4 (4, k) = 6. * a4 (1, k) - 3. * (a4 (2, k) + a4 (3, k)) + if (abs (a4 (4, k)) > abs (a4 (2, k) - a4 (3, k))) then + ! check within the smooth region if subgrid profile is non - monotonic + pmp_1 = a4 (1, k) - 2.0 * gam (k + 1) + lac_1 = pmp_1 + 1.5 * gam (k + 2) + a4 (2, k) = min (max (a4 (2, k), min (a4 (1, k), pmp_1, lac_1)), & + max (a4 (1, k), pmp_1, lac_1)) + pmp_2 = a4 (1, k) + 2.0 * gam (k) + lac_2 = pmp_2 - 1.5 * gam (k - 1) + a4 (3, k) = min (max (a4 (3, k), min (a4 (1, k), pmp_2, lac_2)), & + max (a4 (1, k), pmp_2, lac_2)) + endif + endif + enddo + else + do k = 3, km - 2 + if (extm (k)) then + if (a4 (1, k) < qp_min .or. extm (k - 1) .or. extm (k + 1)) then + a4 (2, k) = a4 (1, k) + a4 (3, k) = a4 (1, k) + endif + endif + enddo + endif + + do k = 1, km - 1 + a4 (4, k) = 6. * a4 (1, k) - 3. * (a4 (2, k) + a4 (3, k)) + enddo + + k = km - 1 + if (extm (k)) then + a4 (2, k) = a4 (1, k) + a4 (3, k) = a4 (1, k) + a4 (4, k) = 0. + else + da1 = a4 (3, k) - a4 (2, k) + da2 = da1 ** 2 + a6da = a4 (4, k) * da1 + if (a6da < - da2) then + a4 (4, k) = 3. * (a4 (2, k) - a4 (1, k)) + a4 (3, k) = a4 (2, k) - a4 (4, k) + elseif (a6da > da2) then + a4 (4, k) = 3. * (a4 (3, k) - a4 (1, k)) + a4 (2, k) = a4 (3, k) - a4 (4, k) + endif + endif + + call cs_limiters (km - 1, a4) + + ! ----------------------------------------------------------------------- + ! bottom layer: + ! ----------------------------------------------------------------------- + + a4 (2, km) = a4 (1, km) + a4 (3, km) = a4 (1, km) + a4 (4, km) = 0. + +end subroutine cs_profile + +!>\ingroup mod_gfdl_cloud_mp +!! This subroutine perform positive definite constraint. +subroutine cs_limiters (km, a4) + + implicit none + + integer, intent (in) :: km + + real, intent (inout) :: a4 (4, km) ! ppm array + + real, parameter :: r12 = 1. / 12. + + integer :: k + + ! ----------------------------------------------------------------------- + ! positive definite constraint + ! ----------------------------------------------------------------------- + + do k = 1, km + if (abs (a4 (3, k) - a4 (2, k)) < - a4 (4, k)) then + if ((a4 (1, k) + 0.25 * (a4 (3, k) - a4 (2, k)) ** 2 / a4 (4, k) + a4 (4, k) * r12) < 0.) then + if (a4 (1, k) < a4 (3, k) .and. a4 (1, k) < a4 (2, k)) then + a4 (3, k) = a4 (1, k) + a4 (2, k) = a4 (1, k) + a4 (4, k) = 0. + elseif (a4 (3, k) > a4 (2, k)) then + a4 (4, k) = 3. * (a4 (2, k) - a4 (1, k)) + a4 (3, k) = a4 (2, k) - a4 (4, k) + else + a4 (4, k) = 3. * (a4 (3, k) - a4 (1, k)) + a4 (2, k) = a4 (3, k) - a4 (4, k) + endif + endif + endif + enddo + +end subroutine cs_limiters + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine calculates vertical fall speed of snow/ice/graupel. +subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in), dimension (ktop:kbot) :: den, qs, qi, qg, ql, tk + real, intent (out), dimension (ktop:kbot) :: vts, vti, vtg + + ! fall velocity constants: + + real, parameter :: thi = 1.0e-8 !< cloud ice threshold for terminal fall + real, parameter :: thg = 1.0e-8 + real, parameter :: ths = 1.0e-8 + + ! coefficient for the parameterization of mass weighted fall velocity + ! as a function of temperature and IWC. + ! Table 1 in Deng and Mace (2008) \cite deng_and_mace_2008. + real, parameter :: aa = - 4.14122e-5 + real, parameter :: bb = - 0.00538922 + real, parameter :: cc = - 0.0516344 + real, parameter :: dd = 0.00216078 + real, parameter :: ee = 1.9714 + + ! marshall - palmer constants + + real, parameter :: vcons = 6.6280504 + real, parameter :: vcong = 87.2382675 + real, parameter :: vconh = vcong * sqrt (rhoh / rhog) + real, parameter :: norms = 942477796.076938 + real, parameter :: normg = 5026548245.74367 + real, parameter :: normh = pi * rhoh * rnzh + + real, dimension (ktop:kbot) :: qden, tc, rhof + + real :: vi0 + + integer :: k + + ! ----------------------------------------------------------------------- + ! marshall - palmer formula + ! ----------------------------------------------------------------------- + + ! ----------------------------------------------------------------------- + ! try the local air density -- for global model; the true value could be + ! much smaller than sfcrho over high mountains + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + rhof (k) = sqrt (min (10., sfcrho / den (k))) + enddo + + ! ----------------------------------------------------------------------- + !> - ice: use Deng and Mace (2008) \cite deng_and_mace_2008, which gives smaler + !! fall speed than Heymsfield and Donner (1990) \cite heymsfield_and_donner_1990. + ! ----------------------------------------------------------------------- + + if (const_vi) then + vti (:) = vi_fac + else + ! ----------------------------------------------------------------------- + ! use deng and mace (2008, grl), which gives smaller fall speed than hd90 formula + ! ----------------------------------------------------------------------- + vi0 = 0.01 * vi_fac + do k = ktop, kbot + if (qi (k) < thi) then ! this is needed as the fall - speed maybe problematic for small qi + vti (k) = vf_min + else + tc (k) = tk (k) - tice + vti (k) = (3. + log10 (qi (k) * den (k))) * (tc (k) * (aa * tc (k) + bb) + cc) + dd * tc (k) + ee + vti (k) = vi0 * exp (log_10 * vti (k)) * 0.9 + vti (k) = min (vi_max, max (vf_min, vti (k))) + endif + enddo + endif + + ! ----------------------------------------------------------------------- + !> - snow: + ! ----------------------------------------------------------------------- + + if (const_vs) then + vts (:) = vs_fac ! 1. ifs_2016 + else + do k = ktop, kbot + if (qs (k) < ths) then + vts (k) = vf_min + else + vts (k) = vs_fac * vcons * rhof (k) * exp (0.0625 * log (qs (k) * den (k) / norms)) + vts (k) = min (vs_max, max (vf_min, vts (k))) + endif + enddo + endif + + ! ----------------------------------------------------------------------- + !> - graupel: + ! ----------------------------------------------------------------------- + if (const_vg) then + vtg (:) = vg_fac ! 2. + else + if (do_hail) then + do k = ktop, kbot + if (qg (k) < thg) then + vtg (k) = vf_min + else + vtg (k) = vg_fac * vconh * rhof (k) * sqrt (sqrt (sqrt (qg (k) * den (k) / normh))) + vtg (k) = min (vg_max, max (vf_min, vtg (k))) + endif + enddo + else + do k = ktop, kbot + if (qg (k) < thg) then + vtg (k) = vf_min + else + vtg (k) = vg_fac * vcong * rhof (k) * sqrt (sqrt (sqrt (qg (k) * den (k) / normg))) + vtg (k) = min (vg_max, max (vf_min, vtg (k))) + endif + enddo + endif + endif + +end subroutine fall_speed + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine sets up gfdl cloud microphysics parameters. +subroutine setupm + + implicit none + + real :: gcon, cd, scm3, pisq, act (8) + real :: vdifu, tcond + real :: visk + real :: ch2o, hltf + real :: hlts, hltc, ri50 + + real, parameter :: gam263 = 1.456943, gam275 = 1.608355, gam290 = 1.827363, & + gam325 = 2.54925, gam350 = 3.323363, gam380 = 4.694155, & + gam425 = 8.285063, gam450 = 11.631769, gam480 = 17.837789, & + gam625 = 184.860962, gam680 = 496.604067 + + ! intercept parameters + +! real, parameter :: rnzr = 8.0e6 ! lin83 +! real, parameter :: rnzs = 3.0e6 ! lin83 +! real, parameter :: rnzg = 4.0e6 ! rh84 + + ! density parameters + +! real, parameter :: rhos = 0.1e3 !< lin83 (snow density; 1 / 10 of water) +! real, parameter :: rhog = 0.4e3 !< rh84 (graupel density) + real, parameter :: acc (3) = (/ 5.0, 2.0, 0.5 /) + + real den_rc + + integer :: i, k + + pie = 4. * atan (1.0) + + ! s. klein's formular (eq 16) from am2 + + fac_rc = (4. / 3.) * pie * rhor * rthresh ** 3 + + if (prog_ccn) then + ! if (master) write (*, *) 'prog_ccn option is .t.' + else + den_rc = fac_rc * ccn_o * 1.e6 + ! if (master) write (*, *) 'mp: for ccn_o = ', ccn_o, 'ql_rc = ', den_rc + den_rc = fac_rc * ccn_l * 1.e6 + ! if (master) write (*, *) 'mp: for ccn_l = ', ccn_l, 'ql_rc = ', den_rc + endif + + vdifu = 2.11e-5 + tcond = 2.36e-2 + + visk = 1.259e-5 + hlts = 2.8336e6 + hltc = 2.5e6 + hltf = 3.336e5 + + ch2o = 4.1855e3 + ri50 = 1.e-4 + + pisq = pie * pie + scm3 = (visk / vdifu) ** (1. / 3.) + + cracs = pisq * rnzr * rnzs * rhos + csacr = pisq * rnzr * rnzs * rhor + if (do_hail) then + cgacr = pisq * rnzr * rnzh * rhor + cgacs = pisq * rnzh * rnzs * rhos + else + cgacr = pisq * rnzr * rnzg * rhor + cgacs = pisq * rnzg * rnzs * rhos + endif + cgacs = cgacs * c_pgacs + + ! act: 1 - 2:racs (s - r) ; 3 - 4:sacr (r - s) ; + ! 5 - 6:gacr (r - g) ; 7 - 8:gacs (s - g) + + act (1) = pie * rnzs * rhos + act (2) = pie * rnzr * rhor + if (do_hail) then + act (6) = pie * rnzh * rhoh + else + act (6) = pie * rnzg * rhog + endif + act (3) = act (2) + act (4) = act (1) + act (5) = act (2) + act (7) = act (1) + act (8) = act (6) + + do i = 1, 3 + do k = 1, 4 + acco (i, k) = acc (i) / (act (2 * k - 1) ** ((7 - i) * 0.25) * act (2 * k) ** (i * 0.25)) + enddo + enddo + + gcon = 40.74 * sqrt (sfcrho) ! 44.628 + + csacw = pie * rnzs * clin * gam325 / (4. * act (1) ** 0.8125) + ! decreasing csacw to reduce cloud water --- > snow + + craci = pie * rnzr * alin * gam380 / (4. * act (2) ** 0.95) + csaci = csacw * c_psaci + + if (do_hail) then + cgacw = pie * rnzh * gam350 * gcon / (4. * act (6) ** 0.875) + else + cgacw = pie * rnzg * gam350 * gcon / (4. * act (6) ** 0.875) + endif + ! cgaci = cgacw * 0.1 + + ! sjl, may 28, 2012 + cgaci = cgacw * 0.05 + ! sjl, may 28, 2012 + + cracw = craci ! cracw = 3.27206196043822 + cracw = c_cracw * cracw + + ! subl and revp: five constants for three separate processes + + cssub (1) = 2. * pie * vdifu * tcond * rvgas * rnzs + if (do_hail) then + cgsub (1) = 2. * pie * vdifu * tcond * rvgas * rnzh + else + cgsub (1) = 2. * pie * vdifu * tcond * rvgas * rnzg + endif + crevp (1) = 2. * pie * vdifu * tcond * rvgas * rnzr + cssub (2) = 0.78 / sqrt (act (1)) + cgsub (2) = 0.78 / sqrt (act (6)) + crevp (2) = 0.78 / sqrt (act (2)) + cssub (3) = 0.31 * scm3 * gam263 * sqrt (clin / visk) / act (1) ** 0.65625 + cgsub (3) = 0.31 * scm3 * gam275 * sqrt (gcon / visk) / act (6) ** 0.6875 + crevp (3) = 0.31 * scm3 * gam290 * sqrt (alin / visk) / act (2) ** 0.725 + cssub (4) = tcond * rvgas + cssub (5) = hlts ** 2 * vdifu + cgsub (4) = cssub (4) + crevp (4) = cssub (4) + cgsub (5) = cssub (5) + crevp (5) = hltc ** 2 * vdifu + + cgfr (1) = 20.e2 * pisq * rnzr * rhor / act (2) ** 1.75 + cgfr (2) = 0.66 + + ! smlt: five constants (lin et al. 1983) + + csmlt (1) = 2. * pie * tcond * rnzs / hltf + csmlt (2) = 2. * pie * vdifu * rnzs * hltc / hltf + csmlt (3) = cssub (2) + csmlt (4) = cssub (3) + csmlt (5) = ch2o / hltf + + ! gmlt: five constants + + if (do_hail) then + cgmlt (1) = 2. * pie * tcond * rnzh / hltf + cgmlt (2) = 2. * pie * vdifu * rnzh * hltc / hltf + else + cgmlt (1) = 2. * pie * tcond * rnzg / hltf + cgmlt (2) = 2. * pie * vdifu * rnzg * hltc / hltf + endif + cgmlt (3) = cgsub (2) + cgmlt (4) = cgsub (3) + cgmlt (5) = ch2o / hltf + + es0 = 6.107799961e2 ! ~6.1 mb + ces0 = eps * es0 + +end subroutine setupm + +! ======================================================================= +! initialization of gfdl cloud microphysics +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine 'gfdl_cloud_microphys_init' initializes the GFDL +!! cloud microphysics. +subroutine gfdl_cloud_microphys_mod_init (me, master, nlunit, input_nml_file, logunit, & + fn_nml, errmsg, errflg) + + implicit none + + integer, intent (in) :: me + integer, intent (in) :: master + integer, intent (in) :: nlunit + integer, intent (in) :: logunit + + character (len = 64), intent (in) :: fn_nml + character (len = *), intent (in) :: input_nml_file(:) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: ios + logical :: exists + + ! integer, intent (in) :: id, jd, kd + ! integer, intent (in) :: axes (4) + ! type (time_type), intent (in) :: time + + ! integer :: unit, io, ierr, k, logunit + ! logical :: flag + ! real :: tmp, q1, q2 + + ! master = (mpp_pe () .eq.mpp_root_pe ()) + + ! Initialize CCPP error-handling + errflg = 0 + errmsg = '' + +#ifdef INTERNAL_FILE_NML + read (input_nml_file, nml = gfdl_cloud_microphysics_nml) +#else + inquire (file = trim (fn_nml), exist = exists) + if (.not. exists) then + write (6, *) 'gfdl - mp :: namelist file: ', trim (fn_nml), ' does not exist' + errflg = 1 + errmsg = 'ERROR(gfdl_cloud_microphys_mod_init): namelist file '//trim (fn_nml)//' does not exist' + return + else + open (unit = nlunit, file = fn_nml, action = 'read' , status = 'old', iostat = ios) + endif + rewind (nlunit) + read (nlunit, nml = gfdl_cloud_microphysics_nml) + close (nlunit) +#endif + + ! write version number and namelist to log file + if (me == master) then + write (logunit, *) " ================================================================== " + write (logunit, *) "gfdl_cloud_microphys_mod" + write (logunit, nml = gfdl_cloud_microphysics_nml) + endif + + if (do_setup) then + call setup_con + call setupm + do_setup = .false. + endif + + log_10 = log (10.) + + tice0 = tice - 0.01 + t_wfr = tice - 40.0 ! supercooled water can exist down to - 48 c, which is the "absolute" + + ! if (master) write (logunit, nml = gfdl_cloud_microphys_nml) + ! + ! id_vtr = register_diag_field (mod_name, 'vt_r', axes (1:3), time, & + ! 'rain fall speed', 'm / s', missing_value = missing_value) + ! id_vts = register_diag_field (mod_name, 'vt_s', axes (1:3), time, & + ! 'snow fall speed', 'm / s', missing_value = missing_value) + ! id_vtg = register_diag_field (mod_name, 'vt_g', axes (1:3), time, & + ! 'graupel fall speed', 'm / s', missing_value = missing_value) + ! id_vti = register_diag_field (mod_name, 'vt_i', axes (1:3), time, & + ! 'ice fall speed', 'm / s', missing_value = missing_value) + + ! id_droplets = register_diag_field (mod_name, 'droplets', axes (1:3), time, & + ! 'droplet number concentration', '# / m3', missing_value = missing_value) + ! id_rh = register_diag_field (mod_name, 'rh_lin', axes (1:2), time, & + ! 'relative humidity', 'n / a', missing_value = missing_value) + + ! id_rain = register_diag_field (mod_name, 'rain_lin', axes (1:2), time, & + ! 'rain_lin', 'mm / day', missing_value = missing_value) + ! id_snow = register_diag_field (mod_name, 'snow_lin', axes (1:2), time, & + ! 'snow_lin', 'mm / day', missing_value = missing_value) + ! id_graupel = register_diag_field (mod_name, 'graupel_lin', axes (1:2), time, & + ! 'graupel_lin', 'mm / day', missing_value = missing_value) + ! id_ice = register_diag_field (mod_name, 'ice_lin', axes (1:2), time, & + ! 'ice_lin', 'mm / day', missing_value = missing_value) + ! id_prec = register_diag_field (mod_name, 'prec_lin', axes (1:2), time, & + ! 'prec_lin', 'mm / day', missing_value = missing_value) + + ! if (master) write (*, *) 'prec_lin diagnostics initialized.', id_prec + + ! id_cond = register_diag_field (mod_name, 'cond_lin', axes (1:2), time, & + ! 'total condensate', 'kg / m ** 2', missing_value = missing_value) + ! id_var = register_diag_field (mod_name, 'var_lin', axes (1:2), time, & + ! 'subgrid variance', 'n / a', missing_value = missing_value) + + ! call qsmith_init + + ! testing the water vapor tables + + ! if (mp_debug .and. master) then + ! write (*, *) 'testing water vapor tables in gfdl_cloud_microphys' + ! tmp = tice - 90. + ! do k = 1, 25 + ! q1 = wqsat_moist (tmp, 0., 1.e5) + ! q2 = qs1d_m (tmp, 0., 1.e5) + ! write (*, *) nint (tmp - tice), q1, q2, 'dq = ', q1 - q2 + ! tmp = tmp + 5. + ! enddo + ! endif + + ! if (master) write (*, *) 'gfdl_cloud_micrphys diagnostics initialized.' + + module_is_initialized = .true. + +!+---+-----------------------------------------------------------------+ +!..Set these variables needed for computing radar reflectivity. These +!.. get used within radar_init to create other variables used in the +!.. radar module. + + xam_r = pi*rhor/6. + xbm_r = 3. + xmu_r = 0. + xam_s = pi*rhos/6. + xbm_s = 3. + xmu_s = 0. + xam_g = pi*rhog/6. + xbm_g = 3. + xmu_g = 0. + + call radar_init + +end subroutine gfdl_cloud_microphys_mod_init + +! ======================================================================= +! end of gfdl cloud microphysics +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine 'gfdl_cloud_microphys_init' terminates the GFDL +!! cloud microphysics. +subroutine gfdl_cloud_microphys_mod_end() + + implicit none + + deallocate (table) + deallocate (table2) + deallocate (table3) + deallocate (tablew) + deallocate (des) + deallocate (des2) + deallocate (des3) + deallocate (desw) + + tables_are_initialized = .false. + +end subroutine gfdl_cloud_microphys_mod_end + +! ======================================================================= +! qsmith table initialization +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine 'setup_con' sets up constants and calls 'qsmith_init'. +subroutine setup_con + + implicit none + + ! master = (mpp_pe () .eq.mpp_root_pe ()) + + rgrav = 1. / grav + + if (.not. qsmith_tables_initialized) call qsmith_init + + qsmith_tables_initialized = .true. + +end subroutine setup_con + +! ======================================================================= +!>\ingroup gfdlmp +!>@brief The function is an accretion function (Lin et al.(1983) +!! \cite lin_et_al_1983 ) +! ======================================================================= + +real function acr3d (v1, v2, q1, q2, c, cac, rho) + + implicit none + + real, intent (in) :: v1, v2, c, rho + real, intent (in) :: q1, q2 ! mixing ratio!!! + real, intent (in) :: cac (3) + + real :: t1, s1, s2 + + ! integer :: k + ! + ! real :: a + ! + ! a = 0.0 + ! do k = 1, 3 + ! a = a + cac (k) * ((q1 * rho) ** ((7 - k) * 0.25) * (q2 * rho) ** (k * 0.25)) + ! enddo + ! acr3d = c * abs (v1 - v2) * a / rho + + ! optimized + + t1 = sqrt (q1 * rho) + s1 = sqrt (q2 * rho) + s2 = sqrt (s1) ! s1 = s2 ** 2 + acr3d = c * abs (v1 - v2) * q1 * s2 * (cac (1) * t1 + cac (2) * sqrt (t1) * s2 + cac (3) * s1) + +end function acr3d + +! ======================================================================= +!>\ingroup gfdlmp +!>@brief Melting of snow function (Lin et al.(1983) \cite lin_et_al_1983) +!! note: psacw and psacr must be calc before smlt is called +! ======================================================================= + +real function smlt (tc, dqs, qsrho, psacw, psacr, c, rho, rhofac) + + implicit none + + real, intent (in) :: tc, dqs, qsrho, psacw, psacr, c (5), rho, rhofac + + smlt = (c (1) * tc / rho - c (2) * dqs) * (c (3) * sqrt (qsrho) + & + c (4) * qsrho ** 0.65625 * sqrt (rhofac)) + c (5) * tc * (psacw + psacr) + +end function smlt + +! ======================================================================= +!>@brief Melting of graupel function (Eq.(47) in Lin et al. 1983 \cite lin_et_al_1983) +!!\n note: \f$P_{gacw}\f$ and \f$P_{gacr}\f$ must be calculated before gmlt is called. +! ======================================================================= + +real function gmlt (tc, dqs, qgrho, pgacw, pgacr, c, rho) + + implicit none + + real, intent (in) :: tc, dqs, qgrho, pgacw, pgacr, c (5), rho + + gmlt = (c (1) * tc / rho - c (2) * dqs) * (c (3) * sqrt (qgrho) + & + c (4) * qgrho ** 0.6875 / rho ** 0.25) + c (5) * tc * (pgacw + pgacr) + +end function gmlt + +! ======================================================================= +! initialization +! prepare saturation water vapor pressure tables +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!>@brief The subroutine 'qsmith_init' initializes lookup tables for saturation +!! water vapor pressure for the following utility routines that are designed +!! to return qs consistent with the assumptions in FV3. +!>@details The calculations are highly accurate values based on the Clausius-Clapeyron +!! equation. +! ======================================================================= +subroutine qsmith_init + + implicit none + + integer, parameter :: length = 2621 + + integer :: i + + if (.not. tables_are_initialized) then + + ! master = (mpp_pe () .eq. mpp_root_pe ()) + ! if (master) print *, ' gfdl mp: initializing qs tables' + + ! debug code + ! print *, mpp_pe (), allocated (table), allocated (table2), & + ! allocated (table3), allocated (tablew), allocated (des), & + ! allocated (des2), allocated (des3), allocated (desw) + ! end debug code + + ! generate es table (dt = 0.1 deg. c) + + allocate (table (length)) + allocate (table2 (length)) + allocate (table3 (length)) + allocate (tablew (length)) + allocate (des (length)) + allocate (des2 (length)) + allocate (des3 (length)) + allocate (desw (length)) + + call qs_table (length) + call qs_table2 (length) + call qs_table3 (length) + call qs_tablew (length) + + do i = 1, length - 1 + des (i) = max (0., table (i + 1) - table (i)) + des2 (i) = max (0., table2 (i + 1) - table2 (i)) + des3 (i) = max (0., table3 (i + 1) - table3 (i)) + desw (i) = max (0., tablew (i + 1) - tablew (i)) + enddo + des (length) = des (length - 1) + des2 (length) = des2 (length - 1) + des3 (length) = des3 (length - 1) + desw (length) = desw (length - 1) + + tables_are_initialized = .true. + + endif + +end subroutine qsmith_init + +! ======================================================================= +! compute the saturated specific humidity for table ii +!>\ingroup mod_gfdl_cloud_mp +!>@brief The function 'wqs1' returns the saturation vapor pressure over pure +!! liquid water for a given temperature and air density. +real function wqs1 (ta, den) + + implicit none + + !> pure water phase; universal dry / moist formular using air density + !> input "den" can be either dry or moist air density + + real, intent (in) :: ta, den + + real :: es, ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = tablew (it) + (ap1 - it) * desw (it) + wqs1 = es / (rvgas * ta * den) + +end function wqs1 + +! ======================================================================= +! compute the gradient of saturated specific humidity for table ii +!>@brief The function 'wqs2' returns the saturation vapor pressure over pure +!! liquid water for a given temperature and air density, as well as the +!! analytic dqs/dT: rate of change of saturation vapor pressure WRT temperature. +! ======================================================================= + +real function wqs2 (ta, den, dqdt) + + implicit none + + ! pure water phase; universal dry / moist formular using air density + ! input "den" can be either dry or moist air density + + real, intent (in) :: ta, den + + real, intent (out) :: dqdt + + real :: es, ap1, tmin + + integer :: it + + tmin = table_ice - 160. + + if (.not. tables_are_initialized) call qsmith_init + + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = tablew (it) + (ap1 - it) * desw (it) + wqs2 = es / (rvgas * ta * den) + it = ap1 - 0.5 + ! finite diff, del_t = 0.1: + dqdt = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta * den) + +end function wqs2 + +! ======================================================================= +! compute wet buld temperature +!>@brief The function 'wet_bulb' uses 'wqs2' to compute the wet-bulb temperature +!! from the mixing ratio and the temperature. +! ======================================================================= + +real function wet_bulb (q, t, den) + + implicit none + + real, intent (in) :: t, q, den + + real :: qs, tp, dqdt + + wet_bulb = t + qs = wqs2 (wet_bulb, den, dqdt) + tp = 0.5 * (qs - q) / (1. + lcp * dqdt) * lcp + wet_bulb = wet_bulb - tp + + ! tp is negative if super - saturated + if (tp > 0.01) then + qs = wqs2 (wet_bulb, den, dqdt) + tp = (qs - q) / (1. + lcp * dqdt) * lcp + wet_bulb = wet_bulb - tp + endif + +end function wet_bulb + +! ======================================================================= +!>@brief The function 'iqs1' computes the saturated specific humidity +!! for table iii +! ======================================================================= + +real function iqs1 (ta, den) + + implicit none + + !> water - ice phase; universal dry / moist formular using air density + !> input "den" can be either dry or moist air density + + real, intent (in) :: ta, den + + real :: es, ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = table2 (it) + (ap1 - it) * des2 (it) + iqs1 = es / (rvgas * ta * den) + +end function iqs1 + +! ======================================================================= +!>@brief The function 'iqs2' computes the gradient of saturated specific +!! humidity for table iii +! ======================================================================= + +real function iqs2 (ta, den, dqdt) + + implicit none + + ! water - ice phase; universal dry / moist formular using air density + ! input "den" can be either dry or moist air density + + real, intent (in) :: ta, den + + real, intent (out) :: dqdt + + real :: es, ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = table2 (it) + (ap1 - it) * des2 (it) + iqs2 = es / (rvgas * ta * den) + it = ap1 - 0.5 + dqdt = 10. * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) / (rvgas * ta * den) + +end function iqs2 + +! ======================================================================= +!>@brief The function 'qs1d_moist' computes the gradient of saturated +!! specific humidity for table iii. +! ======================================================================= + +real function qs1d_moist (ta, qv, pa, dqdt) + + implicit none + + real, intent (in) :: ta, pa, qv + + real, intent (out) :: dqdt + + real :: es, ap1, tmin, eps10 + + integer :: it + + tmin = table_ice - 160. + eps10 = 10. * eps + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = table2 (it) + (ap1 - it) * des2 (it) + qs1d_moist = eps * es * (1. + zvir * qv) / pa + it = ap1 - 0.5 + dqdt = eps10 * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) * (1. + zvir * qv) / pa + +end function qs1d_moist + +! ======================================================================= +! compute the gradient of saturated specific humidity for table ii +!>@brief The function 'wqsat2_moist' computes the saturated specific humidity +!! for pure liquid water , as well as des/dT. +! ======================================================================= + +real function wqsat2_moist (ta, qv, pa, dqdt) + + implicit none + + real, intent (in) :: ta, pa, qv + + real, intent (out) :: dqdt + + real :: es, ap1, tmin, eps10 + + integer :: it + + tmin = table_ice - 160. + eps10 = 10. * eps + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = tablew (it) + (ap1 - it) * desw (it) + wqsat2_moist = eps * es * (1. + zvir * qv) / pa + it = ap1 - 0.5 + dqdt = eps10 * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) * (1. + zvir * qv) / pa + +end function wqsat2_moist + +! ======================================================================= +! compute the saturated specific humidity for table ii +!>@brief The function 'wqsat_moist' computes the saturated specific humidity +!! for pure liquid water. +! ======================================================================= + +real function wqsat_moist (ta, qv, pa) + + implicit none + + real, intent (in) :: ta, pa, qv + + real :: es, ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = tablew (it) + (ap1 - it) * desw (it) + wqsat_moist = eps * es * (1. + zvir * qv) / pa + +end function wqsat_moist + +! ======================================================================= +!>@brief The function 'qs1d_m' computes the saturated specific humidity +!! for table iii +! ======================================================================= + +real function qs1d_m (ta, qv, pa) + + implicit none + + real, intent (in) :: ta, pa, qv + + real :: es, ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = table2 (it) + (ap1 - it) * des2 (it) + qs1d_m = eps * es * (1. + zvir * qv) / pa + +end function qs1d_m + +! ======================================================================= +!>@brief The function 'd_sat' computes the difference in saturation +!! vapor * density * between water and ice +! ======================================================================= + +real function d_sat (ta, den) + + implicit none + + real, intent (in) :: ta, den + + real :: es_w, es_i, ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es_w = tablew (it) + (ap1 - it) * desw (it) + es_i = table2 (it) + (ap1 - it) * des2 (it) + d_sat = dim (es_w, es_i) / (rvgas * ta * den) ! take positive difference + +end function d_sat + +! ======================================================================= +!>@brief The function 'esw_table' computes the saturated water vapor +!! pressure for table ii +! ======================================================================= + +real function esw_table (ta) + + implicit none + + real, intent (in) :: ta + + real :: ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + esw_table = tablew (it) + (ap1 - it) * desw (it) + +end function esw_table + +! ======================================================================= +!>@brief The function 'es2_table' computes the saturated water +!! vapor pressure for table iii +! ======================================================================= + +real function es2_table (ta) + + implicit none + + real, intent (in) :: ta + + real :: ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (ta, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es2_table = table2 (it) + (ap1 - it) * des2 (it) + +end function es2_table + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine 'esw_table1d' computes the saturated water vapor +!! pressure for table ii. +subroutine esw_table1d (ta, es, n) + + implicit none + + integer, intent (in) :: n + + real, intent (in) :: ta (n) + + real, intent (out) :: es (n) + + real :: ap1, tmin + + integer :: i, it + + tmin = table_ice - 160. + + do i = 1, n + ap1 = 10. * dim (ta (i), tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es (i) = tablew (it) + (ap1 - it) * desw (it) + enddo + +end subroutine esw_table1d + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine 'es3_table1d' computes the saturated water vapor +!! pressure for table iii. +subroutine es2_table1d (ta, es, n) + + implicit none + + integer, intent (in) :: n + + real, intent (in) :: ta (n) + + real, intent (out) :: es (n) + + real :: ap1, tmin + + integer :: i, it + + tmin = table_ice - 160. + + do i = 1, n + ap1 = 10. * dim (ta (i), tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es (i) = table2 (it) + (ap1 - it) * des2 (it) + enddo + +end subroutine es2_table1d + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine 'es3_table1d' computes the saturated water vapor +!! pressure for table iv. +subroutine es3_table1d (ta, es, n) + + implicit none + + integer, intent (in) :: n + + real, intent (in) :: ta (n) + + real, intent (out) :: es (n) + + real :: ap1, tmin + + integer :: i, it + + tmin = table_ice - 160. + + do i = 1, n + ap1 = 10. * dim (ta (i), tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es (i) = table3 (it) + (ap1 - it) * des3 (it) + enddo + +end subroutine es3_table1d + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! saturation water vapor pressure table ii +! 1 - phase table +subroutine qs_tablew (n) + + implicit none + + integer, intent (in) :: n + + real :: delt = 0.1 + real :: tmin, tem, fac0, fac1, fac2 + + integer :: i + + tmin = table_ice - 160. + + ! ----------------------------------------------------------------------- + ! compute es over water + ! ----------------------------------------------------------------------- + + do i = 1, n + tem = tmin + delt * real (i - 1) + fac0 = (tem - t_ice) / (tem * t_ice) + fac1 = fac0 * lv0 + fac2 = (dc_vap * log (tem / t_ice) + fac1) / rvgas + tablew (i) = e00 * exp (fac2) + enddo + +end subroutine qs_tablew + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!>@brief saturation water vapor pressure table iii +! 2 - phase table +subroutine qs_table2 (n) + + implicit none + + integer, intent (in) :: n + + real :: delt = 0.1 + real :: tmin, tem0, tem1, fac0, fac1, fac2 + + integer :: i, i0, i1 + + tmin = table_ice - 160. + + do i = 1, n + tem0 = tmin + delt * real (i - 1) + fac0 = (tem0 - t_ice) / (tem0 * t_ice) + if (i <= 1600) then + ! ----------------------------------------------------------------------- + ! compute es over ice between - 160 deg c and 0 deg c. + ! ----------------------------------------------------------------------- + fac1 = fac0 * li2 + fac2 = (d2ice * log (tem0 / t_ice) + fac1) / rvgas + else + ! ----------------------------------------------------------------------- + ! compute es over water between 0 deg c and 102 deg c. + ! ----------------------------------------------------------------------- + fac1 = fac0 * lv0 + fac2 = (dc_vap * log (tem0 / t_ice) + fac1) / rvgas + endif + table2 (i) = e00 * exp (fac2) + enddo + + ! ----------------------------------------------------------------------- + ! smoother around 0 deg c + ! ----------------------------------------------------------------------- + + i0 = 1600 + i1 = 1601 + tem0 = 0.25 * (table2 (i0 - 1) + 2. * table (i0) + table2 (i0 + 1)) + tem1 = 0.25 * (table2 (i1 - 1) + 2. * table (i1) + table2 (i1 + 1)) + table2 (i0) = tem0 + table2 (i1) = tem1 + +end subroutine qs_table2 + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! saturation water vapor pressure table iv +! 2 - phase table with " - 2 c" as the transition point +subroutine qs_table3 (n) + + implicit none + + integer, intent (in) :: n + + real :: delt = 0.1 + real :: esbasw, tbasw, esbasi, tmin, tem, aa, b, c, d, e + real :: tem0, tem1 + + integer :: i, i0, i1 + + esbasw = 1013246.0 + tbasw = table_ice + 100. + esbasi = 6107.1 + tmin = table_ice - 160. + + do i = 1, n + tem = tmin + delt * real (i - 1) + ! if (i <= 1600) then + if (i <= 1580) then ! change to - 2 c + ! ----------------------------------------------------------------------- + ! compute es over ice between - 160 deg c and 0 deg c. + ! see smithsonian meteorological tables page 350. + ! ----------------------------------------------------------------------- + aa = - 9.09718 * (table_ice / tem - 1.) + b = - 3.56654 * alog10 (table_ice / tem) + c = 0.876793 * (1. - tem / table_ice) + e = alog10 (esbasi) + table3 (i) = 0.1 * 10 ** (aa + b + c + e) + else + ! ----------------------------------------------------------------------- + ! compute es over water between - 2 deg c and 102 deg c. + ! see smithsonian meteorological tables page 350. + ! ----------------------------------------------------------------------- + aa = - 7.90298 * (tbasw / tem - 1.) + b = 5.02808 * alog10 (tbasw / tem) + c = - 1.3816e-7 * (10 ** ((1. - tem / tbasw) * 11.344) - 1.) + d = 8.1328e-3 * (10 ** ((tbasw / tem - 1.) * (- 3.49149)) - 1.) + e = alog10 (esbasw) + table3 (i) = 0.1 * 10 ** (aa + b + c + d + e) + endif + enddo + + ! ----------------------------------------------------------------------- + ! smoother around - 2 deg c + ! ----------------------------------------------------------------------- + + i0 = 1580 + i1 = 1581 + tem0 = 0.25 * (table3 (i0 - 1) + 2. * table (i0) + table3 (i0 + 1)) + tem1 = 0.25 * (table3 (i1 - 1) + 2. * table (i1) + table3 (i1 + 1)) + table3 (i0) = tem0 + table3 (i1) = tem1 + +end subroutine qs_table3 + +! ======================================================================= +! compute the saturated specific humidity for table +! note: this routine is based on "moist" mixing ratio +!>\ingroup mod_gfdl_cloud_mp +!! The function 'qs_blend' computes the saturated specific humidity +!! with a blend of water and ice depending on the temperature. +real function qs_blend (t, p, q) + + implicit none + + real, intent (in) :: t, p, q + + real :: es, ap1, tmin + + integer :: it + + tmin = table_ice - 160. + ap1 = 10. * dim (t, tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es = table (it) + (ap1 - it) * des (it) + qs_blend = eps * es * (1. + zvir * q) / p + +end function qs_blend + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!! saturation water vapor pressure table i +! 3 - phase table +subroutine qs_table (n) + + implicit none + + integer, intent (in) :: n + + real :: delt = 0.1 + real :: tmin, tem, esh20 + real :: wice, wh2o, fac0, fac1, fac2 + real :: esupc (200) + + integer :: i + + tmin = table_ice - 160. + + ! ----------------------------------------------------------------------- + ! compute es over ice between - 160 deg c and 0 deg c. + ! ----------------------------------------------------------------------- + + do i = 1, 1600 + tem = tmin + delt * real (i - 1) + fac0 = (tem - t_ice) / (tem * t_ice) + fac1 = fac0 * li2 + fac2 = (d2ice * log (tem / t_ice) + fac1) / rvgas + table (i) = e00 * exp (fac2) + enddo + + ! ----------------------------------------------------------------------- + ! compute es over water between - 20 deg c and 102 deg c. + ! ----------------------------------------------------------------------- + + do i = 1, 1221 + tem = 253.16 + delt * real (i - 1) + fac0 = (tem - t_ice) / (tem * t_ice) + fac1 = fac0 * lv0 + fac2 = (dc_vap * log (tem / t_ice) + fac1) / rvgas + esh20 = e00 * exp (fac2) + if (i <= 200) then + esupc (i) = esh20 + else + table (i + 1400) = esh20 + endif + enddo + + ! ----------------------------------------------------------------------- + ! derive blended es over ice and supercooled water between - 20 deg c and 0 deg c + ! ----------------------------------------------------------------------- + + do i = 1, 200 + tem = 253.16 + delt * real (i - 1) + wice = 0.05 * (table_ice - tem) + wh2o = 0.05 * (tem - 253.16) + table (i + 1400) = wice * table (i + 1400) + wh2o * esupc (i) + enddo + +end subroutine qs_table + +! ======================================================================= +! compute the saturated specific humidity and the gradient of saturated specific humidity +! input t in deg k, p in pa; p = rho rdry tv, moist pressure +!>\ingroup mod_gfdl_cloud_mp +!! The function 'qsmith' computes the saturated specific humidity +!! with a blend of water and ice depending on the temperature in 3D. +!@details It als oincludes the option for computing des/dT. +subroutine qsmith (im, km, ks, t, p, q, qs, dqdt) + + implicit none + + integer, intent (in) :: im, km, ks + + real, intent (in), dimension (im, km) :: t, p, q + + real, intent (out), dimension (im, km) :: qs + + real, intent (out), dimension (im, km), optional :: dqdt + + real :: eps10, ap1, tmin + + real, dimension (im, km) :: es + + integer :: i, k, it + + tmin = table_ice - 160. + eps10 = 10. * eps + + if (.not. tables_are_initialized) then + call qsmith_init + endif + + do k = ks, km + do i = 1, im + ap1 = 10. * dim (t (i, k), tmin) + 1. + ap1 = min (2621., ap1) + it = ap1 + es (i, k) = table (it) + (ap1 - it) * des (it) + qs (i, k) = eps * es (i, k) * (1. + zvir * q (i, k)) / p (i, k) + enddo + enddo + + if (present (dqdt)) then + do k = ks, km + do i = 1, im + ap1 = 10. * dim (t (i, k), tmin) + 1. + ap1 = min (2621., ap1) - 0.5 + it = ap1 + dqdt (i, k) = eps10 * (des (it) + (ap1 - it) * (des (it + 1) - des (it))) * (1. + zvir * q (i, k)) / p (i, k) + enddo + enddo + endif + +end subroutine qsmith + +! ======================================================================= +!>\ingroup mod_gfdl_cloud_mp +!>@brief The subroutine 'neg_adj' fixes negative water species. +!>@details This is designed for 6-class micro-physics schemes. +subroutine neg_adj (ktop, kbot, pt, dp, qv, ql, qr, qi, qs, qg) + + implicit none + + integer, intent (in) :: ktop, kbot + + real, intent (in), dimension (ktop:kbot) :: dp + + real, intent (inout), dimension (ktop:kbot) :: pt, qv, ql, qr, qi, qs, qg + + real, dimension (ktop:kbot) :: lcpk, icpk + + real :: dq, cvm + + integer :: k + + ! ----------------------------------------------------------------------- + ! define heat capacity and latent heat coefficient + ! ----------------------------------------------------------------------- + + do k = ktop, kbot + cvm = c_air + qv (k) * c_vap + (qr (k) + ql (k)) * c_liq + (qi (k) + qs (k) + qg (k)) * c_ice + lcpk (k) = (lv00 + d0_vap * pt (k)) / cvm + icpk (k) = (li00 + dc_ice * pt (k)) / cvm + enddo + + do k = ktop, kbot + + ! ----------------------------------------------------------------------- + ! ice phase: + ! ----------------------------------------------------------------------- + + ! if cloud ice < 0, borrow from snow + if (qi (k) < 0.) then + qs (k) = qs (k) + qi (k) + qi (k) = 0. + endif + ! if snow < 0, borrow from graupel + if (qs (k) < 0.) then + qg (k) = qg (k) + qs (k) + qs (k) = 0. + endif + ! if graupel < 0, borrow from rain + if (qg (k) < 0.) then + qr (k) = qr (k) + qg (k) + pt (k) = pt (k) - qg (k) * icpk (k) ! heating + qg (k) = 0. + endif + + ! ----------------------------------------------------------------------- + ! liquid phase: + ! ----------------------------------------------------------------------- + + ! if rain < 0, borrow from cloud water + if (qr (k) < 0.) then + ql (k) = ql (k) + qr (k) + qr (k) = 0. + endif + ! if cloud water < 0, borrow from water vapor + if (ql (k) < 0.) then + qv (k) = qv (k) + ql (k) + pt (k) = pt (k) - ql (k) * lcpk (k) ! heating + ql (k) = 0. + endif + + enddo + + ! ----------------------------------------------------------------------- + ! fix water vapor; borrow from below + ! ----------------------------------------------------------------------- + + do k = ktop, kbot - 1 + if (qv (k) < 0.) then + qv (k + 1) = qv (k + 1) + qv (k) * dp (k) / dp (k + 1) + qv (k) = 0. + endif + enddo + + ! ----------------------------------------------------------------------- + ! bottom layer; borrow from above + ! ----------------------------------------------------------------------- + + if (qv (kbot) < 0. .and. qv (kbot - 1) > 0.) then + dq = min (- qv (kbot) * dp (kbot), qv (kbot - 1) * dp (kbot - 1)) + qv (kbot - 1) = qv (kbot - 1) - dq / dp (kbot - 1) + qv (kbot) = qv (kbot) + dq / dp (kbot) + endif + +end subroutine neg_adj + +! ======================================================================= +! compute global sum +!>@brief quick local sum algorithm +! ======================================================================= + +!real function g_sum (p, ifirst, ilast, jfirst, jlast, area, mode) +! +! use mpp_mod, only: mpp_sum +! +! implicit none +! +! integer, intent (in) :: ifirst, ilast, jfirst, jlast +! integer, intent (in) :: mode ! if == 1 divided by area +! +! real, intent (in), dimension (ifirst:ilast, jfirst:jlast) :: p, area +! +! integer :: i, j +! +! real :: gsum +! +! if (global_area < 0.) then +! global_area = 0. +! do j = jfirst, jlast +! do i = ifirst, ilast +! global_area = global_area + area (i, j) +! enddo +! enddo +! call mpp_sum (global_area) +! endif +! +! gsum = 0. +! do j = jfirst, jlast +! do i = ifirst, ilast +! gsum = gsum + p (i, j) * area (i, j) +! enddo +! enddo +! call mpp_sum (gsum) +! +! if (mode == 1) then +! g_sum = gsum / global_area +! else +! g_sum = gsum +! endif +! +!end function g_sum + +! ========================================================================== +!>\ingroup mod_gfdl_cloud_mp +!! The subroutine 'interpolate_z' interpolates to a prescribed height. +subroutine interpolate_z (is, ie, js, je, km, zl, hgt, a3, a2) + + implicit none + + integer, intent (in) :: is, ie, js, je, km + + real, intent (in), dimension (is:ie, js:je, km) :: a3 + + real, intent (in), dimension (is:ie, js:je, km + 1) :: hgt ! hgt (k) > hgt (k + 1) + + real, intent (in) :: zl + + real, intent (out), dimension (is:ie, js:je) :: a2 + + real, dimension (km) :: zm !< middle layer height + + integer :: i, j, k + + !$omp parallel do default (none) shared (is, ie, js, je, km, hgt, zl, a2, a3) private (zm) + + do j = js, je + do i = is, ie + do k = 1, km + zm (k) = 0.5 * (hgt (i, j, k) + hgt (i, j, k + 1)) + enddo + if (zl >= zm (1)) then + a2 (i, j) = a3 (i, j, 1) + elseif (zl <= zm (km)) then + a2 (i, j) = a3 (i, j, km) + else + do k = 1, km - 1 + if (zl <= zm (k) .and. zl >= zm (k + 1)) then + a2 (i, j) = a3 (i, j, k) + (a3 (i, j, k + 1) - a3 (i, j, k)) * (zm (k) - zl) / (zm (k) - zm (k + 1)) + exit + endif + enddo + endif + enddo + enddo + +end subroutine interpolate_z + +! ======================================================================= +!> \ingroup mod_gfdl_cloud_mp +!! The subroutine 'cloud_diagnosis' diagnoses the radius of cloud +!! species. +!>\author Linjiong Zhoum, Shian-Jiann Lin +! ======================================================================= +subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, qmg, t, & + rew, rei, rer, res, reg) + + implicit none + + integer, intent (in) :: is, ie, ks, ke + integer, intent (in), dimension (is:ie) :: lsm ! land sea mask, 0: ocean, 1: land, 2: sea ice + + real, intent (in), dimension (is:ie, ks:ke) :: den, delp, t + real, intent (in), dimension (is:ie, ks:ke) :: qmw, qmi, qmr, qms, qmg !< units: kg / kg + + real, intent (out), dimension (is:ie, ks:ke) :: rew, rei, rer, res, reg !< units: micron + + real, dimension (is:ie, ks:ke) :: qcw, qci, qcr, qcs, qcg !< units: g / m^2 + + integer :: i, k + + real :: lambdar, lambdas, lambdag + real :: dpg, rei_fac, mask, ccn, bw + real, parameter :: rho_0 = 50.e-3 + + real :: rhow = 1.0e3, rhor = 1.0e3, rhos = 1.0e2, rhog = 4.0e2 + real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6 + real :: alphar = 0.8, alphas = 0.25, alphag = 0.5 + real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769 + real :: qmin = 1.0e-12, beta = 1.22, qmin1 = 9.e-6 + + do k = ks, ke + do i = is, ie + + dpg = abs (delp (i, k)) / grav + mask = min (max (real(lsm (i)), 0.0), 2.0) + + ! ----------------------------------------------------------------------- + ! cloud water (Martin et al., 1994) + ! ----------------------------------------------------------------------- + + ccn = 0.80 * (- 1.15e-3 * (ccn_o ** 2) + 0.963 * ccn_o + 5.30) * abs (mask - 1.0) + & + 0.67 * (- 2.10e-4 * (ccn_l ** 2) + 0.568 * ccn_l - 27.9) * (1.0 - abs (mask - 1.0)) + + if (qmw (i, k) .gt. qmin) then + qcw (i, k) = dpg * qmw (i, k) * 1.0e3 + rew (i, k) = exp (1.0 / 3.0 * log ((3.0 * den (i, k) * qmw (i, k)) / (4.0 * pi * rhow * ccn))) * 1.0e4 + rew (i, k) = max (rewmin, min (rewmax, rew (i, k))) + else + qcw (i, k) = 0.0 + rew (i, k) = rewmin + endif + + if (reiflag .eq. 1) then + + ! ----------------------------------------------------------------------- + ! cloud ice (Heymsfield and Mcfarquhar, 1996) + ! ----------------------------------------------------------------------- + + if (qmi (i, k) .gt. qmin1) then + qci (i, k) = dpg * qmi (i, k) * 1.0e3 + rei_fac = log (1.0e3 * qmi (i, k) * den (i, k)) + if (t (i, k) - tice .lt. - 50) then + rei (i, k) = beta / 9.917 * exp (0.109 * rei_fac) * 1.0e3 + elseif (t (i, k) - tice .lt. - 40) then + rei (i, k) = beta / 9.337 * exp (0.080 * rei_fac) * 1.0e3 + elseif (t (i, k) - tice .lt. - 30) then + rei (i, k) = beta / 9.208 * exp (0.055 * rei_fac) * 1.0e3 + else + rei (i, k) = beta / 9.387 * exp (0.031 * rei_fac) * 1.0e3 + endif + rei (i, k) = max (reimin, min (reimax, rei (i, k))) + else + qci (i, k) = 0.0 + rei (i, k) = reimin + endif + + endif + + if (reiflag .eq. 2) then + + ! ----------------------------------------------------------------------- + ! cloud ice (Wyser, 1998) + ! ----------------------------------------------------------------------- + + if (qmi (i, k) .gt. qmin1) then + qci (i, k) = dpg * qmi (i, k) * 1.0e3 + bw = - 2. + 1.e-3 * log10 (den (i, k) * qmi (i, k) / rho_0) * max (0.0, tice - t (i, k)) ** 1.5 + rei (i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw)) + rei (i, k) = max (reimin, min (reimax, rei (i, k))) + else + qci (i, k) = 0.0 + rei (i, k) = reimin + endif + + endif + + ! ----------------------------------------------------------------------- + ! rain (Lin et al., 1983) + ! ----------------------------------------------------------------------- + + if (qmr (i, k) .gt. qmin) then + qcr (i, k) = dpg * qmr (i, k) * 1.0e3 + lambdar = exp (0.25 * log (pi * rhor * n0r / qmr (i, k) / den (i, k))) + rer (i, k) = 0.5 * exp (log (gammar / 6) / alphar) / lambdar * 1.0e6 + rer (i, k) = max (rermin, min (rermax, rer (i, k))) + else + qcr (i, k) = 0.0 + rer (i, k) = rermin + endif + + ! ----------------------------------------------------------------------- + ! snow (Lin et al., 1983) + ! ----------------------------------------------------------------------- + + if (qms (i, k) .gt. qmin1) then + qcs (i, k) = dpg * qms (i, k) * 1.0e3 + lambdas = exp (0.25 * log (pi * rhos * n0s / qms (i, k) / den (i, k))) + res (i, k) = 0.5 * exp (log (gammas / 6) / alphas) / lambdas * 1.0e6 + res (i, k) = max (resmin, min (resmax, res (i, k))) + else + qcs (i, k) = 0.0 + res (i, k) = resmin + endif + + ! ----------------------------------------------------------------------- + ! graupel (Lin et al., 1983) + ! ----------------------------------------------------------------------- + + if (qmg (i, k) .gt. qmin) then + qcg (i, k) = dpg * qmg (i, k) * 1.0e3 + lambdag = exp (0.25 * log (pi * rhog * n0g / qmg (i, k) / den (i, k))) + reg (i, k) = 0.5 * exp (log (gammag / 6) / alphag) / lambdag * 1.0e6 + reg (i, k) = max (regmin, min (regmax, reg (i, k))) + else + qcg (i, k) = 0.0 + reg (i, k) = regmin + endif + + enddo + enddo + +end subroutine cloud_diagnosis + +!+---+-----------------------------------------------------------------+ +!>\ingroup mod_gfdl_cloud_mp +!! This subroutine calculates radar reflectivity. + subroutine refl10cm_gfdl (qv1d, qr1d, qs1d, qg1d, & + t1d, p1d, dBZ, kts, kte, ii, jj, melti) + + IMPLICIT NONE + +!..Sub arguments + INTEGER, INTENT(IN):: kts, kte, ii,jj + REAL, DIMENSION(kts:kte), INTENT(IN):: & + qv1d, qr1d, qs1d, qg1d, t1d, p1d + REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ + +!..Local variables + REAL, DIMENSION(kts:kte):: temp, pres, qv, rho + REAL, DIMENSION(kts:kte):: rr, rs, rg +! REAL:: temp_C + + DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg + DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g + DOUBLE PRECISION:: lamr, lams, lamg + LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg + + REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel + DOUBLE PRECISION:: fmelt_s, fmelt_g + + INTEGER:: i, k, k_0, kbot, n + LOGICAL, INTENT(IN):: melti + DOUBLE PRECISION:: cback, x, eta, f_d +!+---+ + + do k = kts, kte + dBZ(k) = -35.0 + enddo + +!+---+-----------------------------------------------------------------+ +!..Put column of data into local arrays. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + temp(k) = t1d(k) +! temp_C = min(-0.001, temp(K)-273.15) + qv(k) = MAX(1.E-10, qv1d(k)) + pres(k) = p1d(k) + rho(k) = 0.622*pres(k)/(rdgas*temp(k)*(qv(k)+0.622)) + + if (qr1d(k) .gt. 1.E-9) then + rr(k) = qr1d(k)*rho(k) + N0_r(k) = n0r + lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1)) + ilamr(k) = 1./lamr + L_qr(k) = .true. + else + rr(k) = 1.E-12 + L_qr(k) = .false. + endif + + if (qs1d(k) .gt. 1.E-9) then + rs(k) = qs1d(k)*rho(k) + N0_s(k) = n0s + lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1)) + ilams(k) = 1./lams + L_qs(k) = .true. + else + rs(k) = 1.E-12 + L_qs(k) = .false. + endif + + if (qg1d(k) .gt. 1.E-9) then + rg(k) = qg1d(k)*rho(k) + N0_g(k) = n0g + lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1)) + ilamg(k) = 1./lamg + L_qg(k) = .true. + else + rg(k) = 1.E-12 + L_qg(k) = .false. + endif + enddo + +!+---+-----------------------------------------------------------------+ +!..Locate K-level of start of melting (k_0 is level above). +!+---+-----------------------------------------------------------------+ + k_0 = kts + K_LOOP:do k = kte-1, kts, -1 + if ( melti .and. (temp(k).gt.273.15) .and. L_qr(k) & + .and. (L_qs(k+1).or.L_qg(k+1)) ) then + k_0 = MAX(k+1, k_0) + EXIT K_LOOP + endif + enddo K_LOOP +!+---+-----------------------------------------------------------------+ +!..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps) +!.. and non-water-coated snow and graupel when below freezing are +!.. simple. Integrations of m(D)*m(D)*N(D)*dD. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + ze_rain(k) = 1.e-22 + ze_snow(k) = 1.e-22 + ze_graupel(k) = 1.e-22 + if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4) + if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & + * (xam_s/900.0)*(xam_s/900.0) & + * N0_s(k)*xcsg(4)*ilams(k)**xcse(4) + if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & + * (xam_g/900.0)*(xam_g/900.0) & + * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4) + enddo + + +!+---+-----------------------------------------------------------------+ +!..Special case of melting ice (snow/graupel) particles. Assume the +!.. ice is surrounded by the liquid water. Fraction of meltwater is +!.. extremely simple based on amount found above the melting level. +!.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting +!.. routines). +!+---+-----------------------------------------------------------------+ + + if (melti .and. k_0.ge.kts+1) then + do k = k_0-1, kts, -1 + +!..Reflectivity contributed by melting snow + if (L_qs(k) .and. L_qs(k_0) ) then + fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0)) + eta = 0.d0 + lams = 1./ilams(k) + do n = 1, nrbins + x = xam_s * xxDs(n)**xbm_s + call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), & + fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, & + CBACK, mixingrulestring_s, matrixstring_s, & + inclusionstring_s, hoststring_s, & + hostmatrixstring_s, hostinclusionstring_s) + f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n)) + eta = eta + f_d * CBACK * simpson(n) * xdts(n) + enddo + ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta) + endif + + +!..Reflectivity contributed by melting graupel + + if (L_qg(k) .and. L_qg(k_0) ) then + fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0)) + eta = 0.d0 + lamg = 1./ilamg(k) + do n = 1, nrbins + x = xam_g * xxDg(n)**xbm_g + call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), & + fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, & + CBACK, mixingrulestring_g, matrixstring_g, & + inclusionstring_g, hoststring_g, & + hostmatrixstring_g, hostinclusionstring_g) + f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n)) + eta = eta + f_d * CBACK * simpson(n) * xdtg(n) + enddo + ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta) + endif + + enddo + endif + + do k = kte, kts, -1 + dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18) + enddo + + + end subroutine refl10cm_gfdl +!+---+-----------------------------------------------------------------+ + +end module gfdl_cloud_microphys_mod