diff --git a/physics/MP/GFDL/fv_sat_adj.F90 b/physics/MP/GFDL/fv_sat_adj.F90 deleted file mode 100644 index 53543485b..000000000 --- a/physics/MP/GFDL/fv_sat_adj.F90 +++ /dev/null @@ -1,1431 +0,0 @@ -!>\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/fv_sat_adj.meta b/physics/MP/GFDL/fv_sat_adj.meta deleted file mode 100644 index c91e438b7..000000000 --- a/physics/MP/GFDL/fv_sat_adj.meta +++ /dev/null @@ -1,441 +0,0 @@ -[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/gfdl_cloud_microphys.F90 b/physics/MP/GFDL/gfdl_cloud_microphys.F90 deleted file mode 100644 index 0fd84c7ea..000000000 --- a/physics/MP/GFDL/gfdl_cloud_microphys.F90 +++ /dev/null @@ -1,331 +0,0 @@ -!> \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/gfdl_cloud_microphys.meta b/physics/MP/GFDL/gfdl_cloud_microphys.meta deleted file mode 100644 index 719a340e5..000000000 --- a/physics/MP/GFDL/gfdl_cloud_microphys.meta +++ /dev/null @@ -1,482 +0,0 @@ -[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/module_gfdl_cloud_microphys.F90 b/physics/MP/GFDL/module_gfdl_cloud_microphys.F90 deleted file mode 100644 index 5cab1abbc..000000000 --- a/physics/MP/GFDL/module_gfdl_cloud_microphys.F90 +++ /dev/null @@ -1,5075 +0,0 @@ -!> \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