diff --git a/docs/Doxyfile b/docs/Doxyfile index 1e2293dd9..18c928cd2 100644 --- a/docs/Doxyfile +++ b/docs/Doxyfile @@ -90,7 +90,6 @@ INPUT = ./documentation.dox \ ../model/boundary.F90 \ ../model/dyn_core.F90 \ ../model/fv_arrays.F90 \ - ../model/fv_cmp.F90 \ ../model/fv_control.F90 \ ../model/fv_dynamics.F90 \ ../model/fv_fill.F90 \ diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 37f8bb64d..3145998c9 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -167,7 +167,11 @@ module atmosphere_mod use tracer_manager_mod, only: get_tracer_index, get_number_tracers, & NO_TRACER, get_tracer_names use DYCORE_typedefs, only: DYCORE_data_type +#ifdef GFS_TYPES +use GFS_typedefs, only: IPD_data_type => GFS_data_type, kind_phys +#else use IPD_typedefs, only: IPD_data_type, kind_phys => IPD_kind_phys +#endif use fv_iau_mod, only: IAU_external_data_type !----------------- @@ -243,10 +247,7 @@ module atmosphere_mod logical :: cold_start = .false. ! used in initial condition integer, dimension(:), allocatable :: id_tracerdt_dyn - integer :: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel ! condensate species tracer indices -#ifdef CCPP - integer :: cld_amt -#endif + integer :: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, cld_amt ! condensate species tracer indices integer :: mygrid = 1 integer :: p_split = 1 @@ -271,7 +272,7 @@ module atmosphere_mod !! including the grid structures, memory, initial state (self-initialization or restart), !! and diagnostics. subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) -#ifdef CCPP + use ccpp_static_api, only: ccpp_physics_init use CCPP_data, only: ccpp_suite, & cdata => cdata_tile, & @@ -279,7 +280,7 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) #ifdef OPENMP use omp_lib #endif -#endif + type (time_type), intent(in) :: Time_init, Time, Time_step type(grid_box_type), intent(inout) :: Grid_box real(kind=kind_phys), pointer, dimension(:,:), intent(inout) :: area @@ -289,11 +290,7 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) logical :: do_atmos_nudge character(len=32) :: tracer_name, tracer_units real :: ps1, ps2 -#ifdef CCPP - integer :: nthreads - integer :: ierr -#endif - + integer :: nthreads, ierr integer :: nlunit = 9999 character (len = 64) :: fn_nml = 'input.nml' @@ -354,9 +351,7 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) rainwat = get_tracer_index (MODEL_ATMOS, 'rainwat' ) snowwat = get_tracer_index (MODEL_ATMOS, 'snowwat' ) graupel = get_tracer_index (MODEL_ATMOS, 'graupel' ) -#ifdef CCPP cld_amt = get_tracer_index (MODEL_ATMOS, 'cld_amt') -#endif if (max(sphum,liq_wat,ice_wat,rainwat,snowwat,graupel) > Atm(mygrid)%flagstruct%nwat) then call mpp_error (FATAL,' atmosphere_init: condensate species are not first in the list of & @@ -425,7 +420,6 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) call timing_off('ATMOS_INIT') -#ifdef CCPP ! Do CCPP fast physics initialization before call to adiabatic_init (since this calls fv_dynamics) ! For fast physics running over the entire domain, block @@ -474,7 +468,6 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) call mpp_error (FATAL, cdata%errmsg) end if end if -#endif ! --- initiate the start for a restarted regional forecast if ( Atm(mygrid)%gridstruct%regional .and. Atm(mygrid)%flagstruct%warm_start ) then @@ -693,17 +686,16 @@ end subroutine atmosphere_dynamics !>@brief The subroutine 'atmosphere_end' is an API for the termination of the !! FV3 dynamical core responsible for writing out a restart and final diagnostic state. subroutine atmosphere_end (Time, Grid_box, restart_endfcst) -#ifdef CCPP + use ccpp_static_api, only: ccpp_physics_finalize use CCPP_data, only: ccpp_suite use CCPP_data, only: cdata => cdata_tile -#endif + type (time_type), intent(in) :: Time type(grid_box_type), intent(inout) :: Grid_box logical, intent(in) :: restart_endfcst - -#ifdef CCPP integer :: ierr + if (Atm(mygrid)%flagstruct%do_sat_adj) then ! Finalize fast physics call ccpp_physics_finalize(cdata, suite_name=trim(ccpp_suite), group_name="fast_physics", ierr=ierr) @@ -712,7 +704,6 @@ subroutine atmosphere_end (Time, Grid_box, restart_endfcst) call mpp_error (FATAL, cdata%errmsg) end if end if -#endif ! initialize domains for writing global physics data call set_domain ( Atm(mygrid)%domain ) diff --git a/makefile b/makefile index e35a9baba..b27bf4367 100644 --- a/makefile +++ b/makefile @@ -15,12 +15,6 @@ else $(info ) endif -ifneq (,$(findstring CCPP,$(CPPDEFS))) - FAST_PHYSICS_SRCS_F90 = -else - FAST_PHYSICS_SRCS_F90 = ./model/fv_cmp.F90 -endif - LIBRARY = libfv3core.a FFLAGS += -I$(FMS_DIR) -I../gfsphysics -I../ipd -I../io -I. @@ -37,7 +31,6 @@ SRCS_F90 = \ ./model/boundary.F90 \ ./model/dyn_core.F90 \ ./model/fv_arrays.F90 \ - $(FAST_PHYSICS_SRCS_F90) \ ./model/fv_control.F90 \ ./model/fv_dynamics.F90 \ ./model/fv_fill.F90 \ diff --git a/model/fv_cmp.F90 b/model/fv_cmp.F90 deleted file mode 100644 index 34e4b3479..000000000 --- a/model/fv_cmp.F90 +++ /dev/null @@ -1,1209 +0,0 @@ - -!*********************************************************************** -!* 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 . -!*********************************************************************** - -!>@brief The module 'fv_cmp' implements the fast procesesses in the GFDL -!! microphysics -!!>@author Shian-Jiann Lin, Linjiong Zhou -! Fast saturation adjustment is part of the gfdl cloud microphysics -! ======================================================================= - -module fv_cmp_mod -! 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
- - use constants_mod, only: rvgas, rdgas, grav, hlv, hlf, cp_air - use fv_mp_mod, only: is_master - use fv_arrays_mod, only: r_grid - 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 multi_gases_mod, only: virq_qpz, vicpqd_qpz, vicvqd_qpz, num_gas -#endif - - implicit none - - private - - public fv_sat_adj, qs_init - - ! real, parameter :: cp_air = cp_air ! 1004.6, heat capacity of dry air at constant pressure, come from constants_mod - real, parameter :: cp_vap = 4.0 * rvgas !< 1846.0, heat capacity of water vapor at constant pressure - real, parameter :: cv_air = cp_air - rdgas !< 717.55, heat capacity of dry air at constant volume - real, 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, parameter :: c_ice = 1972.0 !< gfdl: heat capacity of ice at - 15 deg c - real, parameter :: c_liq = 4185.5 !< gfdl: heat capacity of liquid at 15 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 :: tice = 273.16 !< freezing temperature - real, parameter :: t_wfr = tice - 40. !< homogeneous freezing temperature - - real, parameter :: lv0 = hlv - dc_vap * tice !< 3.13905782e6, evaporation latent heat coefficient at 0 deg k - real, parameter :: li00 = hlf - dc_ice * tice !< - 2.7105966e5, fusion latent heat coefficient at 0 deg k - - ! real (kind = r_grid), parameter :: e00 = 610.71 ! gfdl: saturation vapor pressure at 0 deg c - real (kind = r_grid), parameter :: e00 = 611.21 !< ifs: saturation vapor pressure at 0 deg c - - real (kind = r_grid), parameter :: d2ice = dc_vap + dc_ice !< - 126, isobaric heating / cooling - real (kind = r_grid), parameter :: li2 = lv0 + li00 !< 2.86799816e6, sublimation latent heat coefficient at 0 deg k - - real, parameter :: lat2 = (hlv + hlf) ** 2 !< used in bigg mechanism - - 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 - - real, allocatable :: table (:), table2 (:), tablew (:), des2 (:), desw (:) - - logical :: mp_initialized = .false. - -contains - -!>@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. -subroutine fv_sat_adj (mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, te0,& -#ifdef MULTI_GASES - km, 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 - - integer, intent (in) :: is, ie, js, je, ng - - logical, intent (in) :: hydrostatic, consv_te, out_dt, last_step, do_qa - - real, intent (in) :: zvir, mdt ! remapping time step - - real, intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: dp, hs - real, intent (in), dimension (is:ie, js:je) :: dpln, delz - - -#ifdef MULTI_GASES - integer, intent(in) :: km - real, intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng,km,*) :: qvi -#else - real, intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: qv -#endif - real, intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: pt, ql, qi, qr, qs, qg - real, intent (inout), dimension (is - ng:, js - ng:) :: q_con, cappa - real, intent (inout), dimension (is:ie, js:je) :: dtdt - - real, intent (out), dimension (is - ng:ie + ng, js - ng:je + ng) :: qa, te0 - - real (kind = r_grid), intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: area - -#ifdef MULTI_GASES - real, dimension (is - ng:ie + ng, js - ng:je + ng) :: qv -#endif - real, dimension (is:ie) :: wqsat, dq2dt, qpz, cvm, t0, pt1, qstar - real, dimension (is:ie) :: icp2, lcp2, tcp2, tcp3 - real, dimension (is:ie) :: den, q_liq, q_sol, q_cond, src, sink, hvar - real, dimension (is:ie) :: mc_air, lhl, lhi - - real :: qsw, rh - real :: tc, qsi, dqsdt, dq, dq0, pidep, qi_crt, tmp, dtmp - real :: tin, rqi, q_plus, q_minus - real :: sdt, dt_bigg, adj_fac - real :: fac_smlt, fac_r2g, fac_i2s, fac_imlt, fac_l2r, fac_v2l, fac_l2v - real :: 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 temp 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 - - ! ----------------------------------------------------------------------- - ! 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))) - - ! ----------------------------------------------------------------------- - ! 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 - -! ======================================================================= -!>@brief the function 'wqs1' computes the -!! saturated specific humidity for table ii -! ======================================================================= -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 = 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 - -! ======================================================================= -!>@brief the function 'wqs1' 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 = 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 - -! ======================================================================= -!>@brief The function 'wqs2'computes the gradient of saturated specific -!! humidity for table ii -! ======================================================================= -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 = 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 - -! ======================================================================= -!>@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, intent (in), dimension (is:ie) :: ta, den - - real, intent (out), dimension (is:ie) :: wqsat, dqdt - - real :: 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 - -! ======================================================================= -!>@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 = 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 - -! ======================================================================= -! initialization -! prepare saturation water vapor pressure tables -! ======================================================================= -!>@brief The subroutine 'qs_init' initializes lookup tables for the saturation mixing ratio. -subroutine qs_init (kmp) - - implicit none - - integer, intent (in) :: kmp - - integer, parameter :: length = 2621 - - integer :: i - - if (mp_initialized) return - - if (is_master ()) write (*, *) 'top layer for gfdl_mp = ', kmp - - ! 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) - - mp_initialized = .true. - -end subroutine qs_init - -! ======================================================================= -! saturation water vapor pressure table i -! 3 - phase table -! ======================================================================= - -subroutine qs_table (n) - - implicit none - - integer, intent (in) :: n - - real (kind = r_grid) :: delt = 0.1 - real (kind = r_grid) :: tmin, tem, esh20 - real (kind = r_grid) :: wice, wh2o, fac0, fac1, fac2 - real (kind = r_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 - -! ======================================================================= -! saturation water vapor pressure table ii -! 1 - phase table -! ======================================================================= - -subroutine qs_tablew (n) - - implicit none - - integer, intent (in) :: n - - real (kind = r_grid) :: delt = 0.1 - real (kind = r_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 - -! ======================================================================= -! saturation water vapor pressure table iii -! 2 - phase table -! ======================================================================= - -subroutine qs_table2 (n) - - implicit none - - integer, intent (in) :: n - - real (kind = r_grid) :: delt = 0.1 - real (kind = r_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_cmp_mod diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index cd9a34f73..5561559ff 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -179,14 +179,12 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, gridstruct, flagstruct, neststruct, idiag, bd, & parent_grid, domain, diss_est, time_total) -#ifdef CCPP use mpp_mod, only: FATAL, mpp_error use ccpp_static_api, only: ccpp_physics_timestep_init, & ccpp_physics_timestep_finalize use CCPP_data, only: ccpp_suite use CCPP_data, only: cdata => cdata_tile use CCPP_data, only: CCPP_interstitial -#endif real, intent(IN) :: bdt !< Large time-step real, intent(IN) :: consv_te @@ -259,17 +257,11 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, ! Local Arrays real :: ws(bd%is:bd%ie,bd%js:bd%je) -#ifndef CCPP - real :: te_2d(bd%is:bd%ie,bd%js:bd%je) -#endif - real :: teq(bd%is:bd%ie,bd%js:bd%je) + real :: teq(bd%is:bd%ie,bd%js:bd%je) real :: ps2(bd%isd:bd%ied,bd%jsd:bd%jed) real :: m_fac(bd%is:bd%ie,bd%js:bd%je) real :: pfull(npz) real, dimension(bd%is:bd%ie):: cvm -#ifndef CCPP - real, allocatable :: dp1(:,:,:), dtdt_m(:,:,:), cappa(:,:,:) -#endif #ifdef MULTI_GASES real, allocatable :: kapad(:,:,:) #endif @@ -280,11 +272,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, integer :: sphum, liq_wat = -999, ice_wat = -999 ! GFDL physics integer :: rainwat = -999, snowwat = -999, graupel = -999, cld_amt = -999 integer :: theta_d = -999 -#ifdef CCPP logical used, do_omega -#else - logical used, last_step, do_omega -#endif #ifdef MULTI_GASES integer, parameter :: max_packs=13 #else @@ -294,17 +282,13 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, integer :: is, ie, js, je integer :: isd, ied, jsd, jed real :: dt2 -#ifdef CCPP integer :: ierr -#endif -#ifdef CCPP ccpp_associate: associate( cappa => CCPP_interstitial%cappa, & dp1 => CCPP_interstitial%te0, & dtdt_m => CCPP_interstitial%dtdt, & last_step => CCPP_interstitial%last_step, & te_2d => CCPP_interstitial%te0_2d ) -#endif is = bd%is ie = bd%ie @@ -326,7 +310,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, nr = nq_tot - flagstruct%dnrts rdg = -rdgas * agrav -#ifdef CCPP ! Call CCPP timestep init call ccpp_physics_timestep_init(cdata, suite_name=trim(ccpp_suite), group_name="fast_physics", ierr=ierr) ! Reset all interstitial variables for CCPP version @@ -336,21 +319,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, CCPP_interstitial%out_dt = (idiag%id_mdt > 0) end if -#else - te_2d = 0. - - allocate ( dp1(isd:ied, jsd:jed, 1:npz) ) - call init_ijk_mem(isd,ied, jsd,jed, npz, dp1, 0.) - -#ifdef MOIST_CAPPA - allocate ( cappa(isd:ied,jsd:jed,npz) ) - call init_ijk_mem(isd,ied, jsd,jed, npz, cappa, 0.) -#else - allocate ( cappa(isd:isd,jsd:jsd,1) ) - cappa = 0. -#endif -#endif - #ifdef MULTI_GASES allocate ( kapad(isd:ied, jsd:jed, npz) ) call init_ijk_mem(isd,ied, jsd,jed, npz, kapad, kappa) @@ -434,7 +402,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, enddo if ( hydrostatic ) then -#if defined(CCPP) && defined(__GFORTRAN__) +#ifdef __GFORTRAN__ !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,zvir,nwat,q,q_con,sphum,liq_wat, & #else !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,nwat,q,q_con,sphum,liq_wat, & @@ -452,7 +420,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, enddo enddo else -#if defined(CCPP) && defined(__GFORTRAN__) +#ifdef __GFORTRAN__ !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,zvir,q,q_con,sphum,liq_wat, & #else !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,q,q_con,sphum,liq_wat, & @@ -461,7 +429,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, #ifdef MULTI_GASES !$OMP kapad, & #endif -#if defined(CCPP) && defined(__GFORTRAN__) +#ifdef __GFORTRAN__ !$OMP kappa,rdg,delp,pt,delz,nwat) & #else !$OMP cappa,kappa,rdg,delp,pt,delz,nwat) & @@ -590,7 +558,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, pt_initialized = .true. endif else -#if defined(CCPP) && defined(__GFORTRAN__) +#ifdef __GFORTRAN__ !$OMP parallel do default(none) shared(is,ie,js,je,npz,pt,pkz,q_con) #else !$OMP parallel do default(none) shared(is,ie,js,je,npz,pt,dp1,pkz,q_con) @@ -617,20 +585,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, last_step = .false. mdt = bdt / real(k_split) -#ifndef CCPP - if ( idiag%id_mdt > 0 .and. (.not. do_adiabatic_init) ) then - allocate ( dtdt_m(is:ie,js:je,npz) ) -!$OMP parallel do default(none) shared(is,ie,js,je,npz,dtdt_m) - do k=1,npz - do j=js,je - do i=is,ie - dtdt_m(i,j,k) = 0. - enddo - enddo - enddo - endif -#endif - call timing_on('FV_DYN_LOOP') do n_map=1, k_split ! first level of time-split k_step = n_map @@ -650,7 +604,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE) #endif call timing_off('COMM_TOTAL') -#if defined(CCPP) && defined(__GFORTRAN__) +#ifdef __GFORTRAN__ !$OMP parallel do default(none) shared(isd,ied,jsd,jed,npz,delp) #else !$OMP parallel do default(none) shared(isd,ied,jsd,jed,npz,dp1,delp) @@ -829,7 +783,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, call timing_off('FV_DYN_LOOP') if ( idiag%id_mdt > 0 .and. (.not.do_adiabatic_init) ) then ! Output temperature tendency due to inline moist physics: -#if defined(CCPP) && defined(__GFORTRAN__) +#ifdef __GFORTRAN__ !$OMP parallel do default(none) shared(is,ie,js,je,npz,bdt) #else !$OMP parallel do default(none) shared(is,ie,js,je,npz,dtdt_m,bdt) @@ -843,9 +797,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, enddo ! call prt_mxm('Fast DTDT (deg/Day)', dtdt_m, is, ie, js, je, 0, npz, 1., gridstruct%area_64, domain) used = send_data(idiag%id_mdt, dtdt_m, fv_time) -#ifndef CCPP - deallocate ( dtdt_m ) -#endif endif if( nwat == 6 ) then @@ -927,7 +878,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, endif if( (flagstruct%consv_am.or.idiag%id_amdt>0) .and. (.not.do_adiabatic_init) ) then -#if defined(CCPP) && defined(__GFORTRAN__) +#ifdef __GFORTRAN__ !$OMP parallel do default(none) shared(is,ie,js,je,teq,dt2,ps2,ps,idiag) #else !$OMP parallel do default(none) shared(is,ie,js,je,te_2d,teq,dt2,ps2,ps,idiag) @@ -978,10 +929,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, #ifdef MULTI_GASES deallocate(kapad) #endif -#ifndef CCPP - deallocate(dp1) - deallocate(cappa) -#endif if ( flagstruct%fv_debug ) then call prt_mxm('UA', ua, is, ie, js, je, ng, npz, 1., gridstruct%area_64, domain) @@ -1002,12 +949,10 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, -50., 100., bad_range, fv_time) endif -#ifdef CCPP ! Call CCPP timestep finalize call ccpp_physics_timestep_finalize(cdata, suite_name=trim(ccpp_suite), group_name="fast_physics", ierr=ierr) end associate ccpp_associate -#endif end subroutine fv_dynamics diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90 index f4c507fa8..1479ec27c 100644 --- a/model/fv_mapz.F90 +++ b/model/fv_mapz.F90 @@ -45,10 +45,6 @@ module fv_mapz_mod ! fv_grid_type ! ! -! fv_cmp_mod -! qs_init, fv_sat_adj -! -! ! fv_fill_mod ! fillz ! @@ -61,6 +57,14 @@ module fv_mapz_mod ! is_master ! ! +! ccpp_static_api +! ccpp_physics_run +! +! +! CCPP_data +! ccpp_suite, cdata_tile, CCPP_interstitial +! +! ! fv_timing_mod ! timing_on, timing_off ! @@ -92,14 +96,11 @@ module fv_mapz_mod use fv_arrays_mod, only: fv_grid_type, fv_grid_bounds_type, R_GRID use fv_timing_mod, only: timing_on, timing_off use fv_mp_mod, only: is_master, mp_reduce_min, mp_reduce_max -#ifndef CCPP - use fv_cmp_mod, only: qs_init, fv_sat_adj -#else + ! CCPP fast physics use ccpp_static_api, only: ccpp_physics_run use CCPP_data, only: ccpp_suite use CCPP_data, only: cdata => cdata_tile use CCPP_data, only: CCPP_interstitial -#endif #ifdef MULTI_GASES use multi_gases_mod, only: virq, virqd, vicpqd, vicvqd, num_gas #endif @@ -214,11 +215,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & !----------------------------------------------------------------------- real, allocatable, dimension(:,:,:) :: dp0, u0, v0 real, allocatable, dimension(:,:,:) :: u_dt, v_dt -#ifdef CCPP real, dimension(is:ie,js:je):: te_2d, zsum0, zsum1 -#else - real, dimension(is:ie,js:je):: te_2d, zsum0, zsum1, dpln -#endif real, dimension(is:ie,km) :: q2, dp2, t0, w2 real, dimension(is:ie,km+1):: pe1, pe2, pk1, pk2, pn2, phis real, dimension(isd:ied,jsd:jed,km):: pe4 @@ -226,22 +223,13 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & real, dimension(is:ie):: gsize, gz, cvm, qv real rcp, rg, rrg, bkh, dtmp, k1k -#ifndef CCPP - logical:: fast_mp_consv -#endif integer:: i,j,k integer:: kdelz -#ifdef CCPP integer:: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, ccn_cm3, iq, n, kp, k_next integer :: ierr -#else - integer:: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, ccn_cm3, iq, n, kmp, kp, k_next -#endif -#ifdef CCPP ccpp_associate: associate( fast_mp_consv => CCPP_interstitial%fast_mp_consv, & kmp => CCPP_interstitial%kmp ) -#endif k1k = rdgas/cv_air ! akap / (1.-akap) = rg/Cv=0.4 rg = rdgas @@ -258,13 +246,6 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & if ( do_adiabatic_init .or. do_sat_adj ) then fast_mp_consv = (.not.do_adiabatic_init) .and. consv>consv_min -#ifndef CCPP - do k=1,km - kmp = k - if ( pfull(k) > 10.E2 ) exit - enddo - call qs_init(kmp) -#endif endif !$OMP parallel do default(none) shared(is,ie,js,je,km,pe,ptop,kord_tm,hydrostatic, & @@ -627,7 +608,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & 1000 continue -#if defined(CCPP) && defined(__GFORTRAN__) +#ifdef __GFORTRAN__ !$OMP parallel default(none) shared(is,ie,js,je,km,ptop,u,v,pe,ua,va,isd,ied,jsd,jed,kord_mt, & !$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, & !$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, & @@ -642,7 +623,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & !$OMP shared(num_gas) & #endif !$OMP private(q2,pe0,pe1,pe2,pe3,qv,cvm,gz,gsize,phis,kdelz,dp2,t0, ierr) -#elif defined(CCPP) +#else !$OMP parallel default(none) shared(is,ie,js,je,km,kmp,ptop,u,v,pe,ua,va,isd,ied,jsd,jed,kord_mt, & !$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, & !$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, & @@ -656,23 +637,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & #ifdef MULTI_GASES !$OMP shared(num_gas) & #endif - !$OMP private(q2,pe0,pe1,pe2,pe3,qv,cvm,gz,gsize,phis,kdelz,dp2,t0, ierr) -#else -!$OMP parallel default(none) shared(is,ie,js,je,km,kmp,ptop,u,v,pe,ua,va,isd,ied,jsd,jed,kord_mt, & -!$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln,adiabatic, & -!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, & -!$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, & -!$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, & -!$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, & -!$OMP mdt,cld_amt,cappa,dtdt,out_dt,rrg,akap,do_sat_adj, & -!$OMP fast_mp_consv,kord_tm,pe4,npx,npy,ccn_cm3, & -!$OMP u_dt,v_dt,c2l_ord,bd,dp0,ps) & - -#ifdef MULTI_GASES -!$OMP shared(num_gas) & -#endif -!$OMP private(q2,pe0,pe1,pe2,pe3,qv,cvm,gz,gsize,phis,kdelz,dpln,dp2,t0) #endif !$OMP do @@ -817,65 +782,13 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & if ( do_sat_adj ) then call timing_on('sat_adj2') -#ifdef CCPP + ! Call to CCPP fast_physics group if (cdata%initialized()) then call ccpp_physics_run(cdata, suite_name=trim(ccpp_suite), group_name='fast_physics', ierr=ierr) if (ierr/=0) call mpp_error(FATAL, "Call to ccpp_physics_run for group 'fast_physics' failed") else call mpp_error (FATAL, 'Lagrangian_to_Eulerian: can not call CCPP fast physics because CCPP not initialized') endif -#else -!$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(abs(mdt), r_vir, is, ie, js, je, ng, hydrostatic, fast_mp_consv, & - te(isd,jsd,k), & -#ifdef MULTI_GASES - km, & -#endif - q(isd,jsd,k,sphum), q(isd,jsd,k,liq_wat), & - q(isd,jsd,k,ice_wat), q(isd,jsd,k,rainwat), & - q(isd,jsd,k,snowwat), q(isd,jsd,k,graupel), & - hs, dpln, delz(is:ie,js:je,kdelz), pt(isd,jsd,k), delp(isd,jsd,k), q_con(isd:,jsd:,k), & - cappa(isd:,jsd:,k), gridstruct%area_64, dtdt(is,js,k), out_dt, last_step, cld_amt>0, q(isd,jsd,k,cld_amt)) - 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 ! OpenMP k-loop - - - 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) + te(i,j,k) - enddo - enddo - enddo - endif -#endif call timing_off('sat_adj2') endif ! do_sat_adj @@ -943,9 +856,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & endif !$OMP end parallel -#ifdef CCPP end associate ccpp_associate -#endif end subroutine Lagrangian_to_Eulerian @@ -3473,15 +3384,6 @@ subroutine moist_cv(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai enddo case(4) ! K_warm_rain with fake ice do i=is,ie -#ifndef CCPP - qv(i) = q(i,j,k,sphum) - qd(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat) -#ifdef MULTI_GASES - cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + qd(i)*c_liq -#else - cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + qd(i)*c_liq -#endif -#else qv(i) = q(i,j,k,sphum) ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat) qs(i) = q(i,j,k,ice_wat) @@ -3490,8 +3392,6 @@ subroutine moist_cv(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice #else cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice -#endif - #endif enddo case(5) @@ -3595,15 +3495,6 @@ subroutine moist_cp(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai enddo case(4) ! K_warm_rain scheme with fake ice do i=is,ie -#ifndef CCPP - qv(i) = q(i,j,k,sphum) - qd(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat) -#ifdef MULTI_GASES - cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + qd(i)*c_liq -#else - cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + qd(i)*c_liq -#endif -#else qv(i) = q(i,j,k,sphum) ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat) qs(i) = q(i,j,k,ice_wat) @@ -3612,9 +3503,6 @@ subroutine moist_cp(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice #else cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice -#endif - - #endif enddo case(5) diff --git a/model/fv_sg.F90 b/model/fv_sg.F90 index 992ad22f9..bfac101bf 100644 --- a/model/fv_sg.F90 +++ b/model/fv_sg.F90 @@ -299,17 +299,6 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & enddo elseif ( nwat==4 ) then do i=is,ie -#ifndef CCPP - q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat) -#ifdef MULTI_GASES - cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq - cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq -#else - cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq - cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq -#endif - -#else q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat) q_sol = q0(i,k,ice_wat) #ifdef MULTI_GASES @@ -318,9 +307,6 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & #else cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice -#endif - - #endif enddo elseif ( nwat==5 ) then @@ -397,11 +383,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & elseif ( nwat==4 ) then do k=1,kbot do i=is,ie -#ifndef CCPP - qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat) -#else qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat) + q0(i,k,ice_wat) -#endif enddo enddo elseif ( nwat==5 ) then @@ -470,12 +452,8 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & elseif ( nwat==3 ) then ! AM3/AM4 qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) elseif ( nwat==4 ) then ! K_warm_rain scheme with fake ice -#ifndef CCPP - qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat) -#else qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + & q0(i,km1,rainwat) -#endif elseif ( nwat==5 ) then ! K_warm_rain scheme with fake ice qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + & q0(i,km1,snowwat) + q0(i,km1,rainwat) @@ -596,16 +574,6 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & enddo elseif ( nwat == 4 ) then do i=is,ie -#ifndef CCPP - q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat) -#ifdef MULTI_GASES - cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq - cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq -#else - cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq - cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq -#endif -#else q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat) q_sol = q0(i,kk,ice_wat) #ifdef MULTI_GASES @@ -614,9 +582,6 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & #else cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice -#endif - - #endif enddo elseif ( nwat == 5 ) then @@ -888,16 +853,6 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & enddo elseif ( nwat==4 ) then do i=is,ie -#ifndef CCPP - q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat) -#ifdef MULTI_GASES - cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq - cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq -#else - cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq - cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq -#endif -#else q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat) q_sol = q0(i,k,ice_wat) #ifdef MULTI_GASES @@ -906,8 +861,6 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & #else cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice -#endif - #endif enddo elseif ( nwat==5 ) then @@ -984,11 +937,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & elseif ( nwat==4 ) then do k=1,kbot do i=is,ie -#ifndef CCPP - qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat) -#else qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat) + q0(i,k,ice_wat) -#endif enddo enddo elseif ( nwat==5 ) then @@ -1053,11 +1002,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & elseif ( nwat==3 ) then ! AM3/AM4 qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) elseif ( nwat==4 ) then ! K_warm_rain scheme with fake ice -#ifndef CCPP - qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat) -#else qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat) + q0(i,km1,ice_wat) -#endif elseif ( nwat==5 ) then qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + & q0(i,km1,snowwat) + q0(i,km1,rainwat) @@ -1177,16 +1122,6 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & enddo elseif ( nwat == 4 ) then do i=is,ie -#ifndef CCPP - q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat) -#ifdef MULTI_GASES - cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq - cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq -#else - cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq - cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq -#endif -#else q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat) q_sol = q0(i,kk,ice_wat) #ifdef MULTI_GASES @@ -1195,8 +1130,6 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & #else cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice -#endif - #endif enddo elseif ( nwat == 5 ) then diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 69ef847bc..cb328bd3a 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -225,9 +225,7 @@ subroutine get_external_ic( Atm, fv_domain, cold_start, dt_atmos ) integer :: is, ie, js, je integer :: isd, ied, jsd, jed, ng integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel -#ifdef CCPP integer :: liq_aero, ice_aero -#endif #ifdef MULTI_GASES integer :: spfo, spfo2, spfo3 #else @@ -327,10 +325,8 @@ subroutine get_external_ic( Atm, fv_domain, cold_start, dt_atmos ) #else o3mr = get_tracer_index(MODEL_ATMOS, 'o3mr') #endif -#ifdef CCPP liq_aero = get_tracer_index(MODEL_ATMOS, 'liq_aero') ice_aero = get_tracer_index(MODEL_ATMOS, 'ice_aero') -#endif if ( liq_wat > 0 ) & call prt_maxmin('liq_wat', Atm%q(:,:,:,liq_wat), is, ie, js, je, ng, Atm%npz, 1.) @@ -353,12 +349,10 @@ subroutine get_external_ic( Atm, fv_domain, cold_start, dt_atmos ) if ( o3mr > 0 ) & call prt_maxmin('O3MR', Atm%q(:,:,:,o3mr), is, ie, js, je, ng, Atm%npz, 1.) #endif -#ifdef CCPP if ( liq_aero > 0) & call prt_maxmin('liq_aero',Atm%q(:,:,:,liq_aero),is, ie, js, je, ng, Atm%npz, 1.) if ( ice_aero > 0) & call prt_maxmin('ice_aero',Atm%q(:,:,:,ice_aero),is, ie, js, je, ng, Atm%npz, 1.) -#endif endif !Now in fv_restart diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index 9bb20c293..d86a490c5 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -53,8 +53,15 @@ module fv_iau_mod get_var1_double, & get_var3_r4, & get_var1_real, check_var_exists +#ifdef GFS_TYPES + use GFS_typedefs, only: IPD_init_type => GFS_init_type, & + IPD_control_type => GFS_control_type, & + kind_phys +#else use IPD_typedefs, only: IPD_init_type, IPD_control_type, & kind_phys => IPD_kind_phys +#endif + use block_control_mod, only: block_control_type use fv_treat_da_inc_mod, only: remap_coef use tracer_manager_mod, only: get_tracer_names,get_tracer_index, get_number_tracers diff --git a/tools/fv_mp_mod.F90 b/tools/fv_mp_mod.F90 index 87c7508b3..ce3a39423 100644 --- a/tools/fv_mp_mod.F90 +++ b/tools/fv_mp_mod.F90 @@ -136,9 +136,7 @@ module fv_mp_mod integer, allocatable :: grids_master_procs(:) integer, dimension(MAX_NNEST) :: tile_fine = 0 !Global index of LAST tile in a mosaic type(nest_domain_type) :: global_nest_domain !ONE structure for ALL levels of nesting -#ifdef CCPP public commglobal -#endif public mp_start, mp_assign_gid, mp_barrier, mp_stop!, npes public domain_decomp, mp_bcst, mp_reduce_max, mp_reduce_sum, mp_gather public mp_reduce_min