diff --git a/physics/gfdl_fv_sat_adj.F90 b/physics/gfdl_fv_sat_adj.F90 index 90040fab0..60fc7a742 100644 --- a/physics/gfdl_fv_sat_adj.F90 +++ b/physics/gfdl_fv_sat_adj.F90 @@ -54,16 +54,17 @@ module fv_sat_adj ! ! DH* TODO - MAKE THIS INPUT ARGUMENTS *DH !use constants_mod, only: rvgas, rdgas, grav, hlv, hlf, cp_air - use physcons, only : rdgas => con_rd, & - rvgas => con_rv, & - grav => con_g, & - hlv => con_hvap, & - hlf => con_hfus, & - cp_air => con_cp + 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 fv_mp_mod, only: is_master !use fv_arrays_mod, only: r_grid use machine, only: kind_grid + use CCPP_typedefs, only: kind_dyn use gfdl_cloud_microphys, only: ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt use gfdl_cloud_microphys, only: icloud_f, sat_adj0, t_sub, cld_min use gfdl_cloud_microphys, only: tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r @@ -73,11 +74,11 @@ module fv_sat_adj public fv_sat_adj_init, fv_sat_adj_run, fv_sat_adj_finalize - real, parameter :: rrg = -rdgas/grav + 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, 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 + real(kind=kind_dyn), parameter :: cp_vap = 4.0_kind_dyn * 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_kind_dyn * 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 @@ -88,22 +89,22 @@ module fv_sat_adj ! 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=kind_dyn), parameter :: c_ice = 1972.0_kind_dyn !< gfdl: heat capacity of ice at - 15 deg c + real(kind=kind_dyn), parameter :: c_liq = 4185.5_kind_dyn !< 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_kind_dyn !< freezing temperature + real(kind=kind_dyn), parameter :: t_wfr = tice - 40._kind_dyn !< 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 :: e00 = 611.21_kind_grid !< 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, 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 (:) + 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 @@ -146,6 +147,7 @@ subroutine fv_sat_adj_init(kmp, errmsg, errflg) 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)) @@ -195,8 +197,8 @@ end subroutine fv_sat_adj_finalize !! \section arg_table_fv_sat_adj_run Argument Table !! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | !! |----------------|---------------------------------------------------------------|----------------------------------------------------------------------------------------|-----------|------|-----------|-----------|--------|----------| -!! | mdt | time_step_for_remapping_for_fast_physics | remapping time step for fast physics | s | 0 | real | | in | F | -!! | zvir | ratio_of_vapor_to_dry_air_gas_constants_minus_one_default_kind| zvir=rv/rd-1.0 | none | 0 | real | | in | F | +!! | mdt | time_step_for_remapping_for_fast_physics | remapping time step for fast physics | s | 0 | real | kind_dyn | in | F | +!! | zvir | ratio_of_vapor_to_dry_air_gas_constants_minus_one_default_kind| zvir=rv/rd-1.0 | none | 0 | real | kind_dyn | in | F | !! | is | starting_x_direction_index | starting X direction index | count | 0 | integer | | in | F | !! | ie | ending_x_direction_index | ending X direction index | count | 0 | integer | | in | F | !! | isd | starting_x_direction_index_domain | starting X direction index for domain | count | 0 | integer | | in | F | @@ -211,29 +213,29 @@ end subroutine fv_sat_adj_finalize !! | ng | number_of_ghost_zones | number of ghost zones defined in fv_mp | count | 0 | integer | | in | F | !! | hydrostatic | flag_for_hydrostatic_solver | flag for use the hydrostatic or nonhydrostatic solver | flag | 0 | logical | | in | F | !! | fast_mp_consv | flag_for_fast_microphysics_energy_conservation | flag for fast microphysics energy conservation | flag | 0 | logical | | in | F | -!! | te0_2d | atmosphere_energy_content_in_column | atmosphere total energy in columns | J m-2 | 2 | real | | inout | F | -!! | te0 | atmosphere_energy_content_at_Lagrangian_surface | atmosphere total energy at Lagrangian surface | J m-2 | 3 | real | | out | F | -!! | qv | water_vapor_specific_humidity_at_Lagrangian_surface | water vapor specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | | inout | F | -!! | ql | cloud_liquid_water_specific_humidity_at_Lagrangian_surface | cloud liquid water specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | | inout | F | -!! | qi | cloud_ice_specific_humidity_at_Lagrangian_surface | cloud ice specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | | inout | F | -!! | qr | cloud_rain_specific_humidity_at_Lagrangian_surface | cloud rain specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | | inout | F | -!! | qs | cloud_snow_specific_humidity_at_Lagrangian_surface | cloud snow specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | | inout | F | -!! | qg | cloud_graupel_specific_humidity_at_Lagrangian_surface | cloud graupel specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | | inout | F | -!! | hs | surface_geopotential_at_Lagrangian_surface | surface geopotential at Lagrangian surface | m2 s-2 | 2 | real | | in | F | -!! | peln | log_pressure_at_Lagrangian_surface | logarithm of pressure at Lagrangian surface | Pa | 3 | real | | in | F | -!! | delz | thickness_at_Lagrangian_surface | thickness at Lagrangian_surface | m | 3 | real | | in | F | -!! | delp | pressure_thickness_at_Lagrangian_surface | pressure thickness at Lagrangian surface | Pa | 3 | real | | in | F | -!! | pt | virtual_temperature_at_Lagrangian_surface | virtual temperature at Lagrangian surface | K | 3 | real | | inout | F | -!! | pkz | finite-volume_mean_edge_pressure_raised_to_the_power_of_kappa | finite-volume mean edge pressure raised to the power of kappa | Pa**kappa | 3 | real | | inout | F | -!! | q_con | cloud_condensed_water_specific_humidity_at_Lagrangian_surface | cloud condensed water specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | | inout | F | -!! | akap | kappa_dry_for_fast_physics | modified kappa for dry air, fast physics | none | 0 | real | | in | F | -!! | cappa | cappa_moist_gas_constant_at_Lagrangian_surface | cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/(1.+r_vir*q(i,j,k,sphum)) ) | none | 3 | real | | inout | F | +!! | te0_2d | atmosphere_energy_content_in_column | atmosphere total energy in columns | J m-2 | 2 | real | kind_dyn | inout | F | +!! | te0 | atmosphere_energy_content_at_Lagrangian_surface | atmosphere total energy at Lagrangian surface | J m-2 | 3 | real | kind_dyn | out | F | +!! | qv | water_vapor_specific_humidity_at_Lagrangian_surface | water vapor specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | kind_dyn | inout | F | +!! | ql | cloud_liquid_water_specific_humidity_at_Lagrangian_surface | cloud liquid water specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | kind_dyn | inout | F | +!! | qi | cloud_ice_specific_humidity_at_Lagrangian_surface | cloud ice specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | kind_dyn | inout | F | +!! | qr | cloud_rain_specific_humidity_at_Lagrangian_surface | cloud rain specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | kind_dyn | inout | F | +!! | qs | cloud_snow_specific_humidity_at_Lagrangian_surface | cloud snow specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | kind_dyn | inout | F | +!! | qg | cloud_graupel_specific_humidity_at_Lagrangian_surface | cloud graupel specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | kind_dyn | inout | F | +!! | hs | surface_geopotential_at_Lagrangian_surface | surface geopotential at Lagrangian surface | m2 s-2 | 2 | real | kind_dyn | in | F | +!! | peln | log_pressure_at_Lagrangian_surface | logarithm of pressure at Lagrangian surface | Pa | 3 | real | kind_dyn | in | F | +!! | delz | thickness_at_Lagrangian_surface | thickness at Lagrangian_surface | m | 3 | real | kind_dyn | in | F | +!! | delp | pressure_thickness_at_Lagrangian_surface | pressure thickness at Lagrangian surface | Pa | 3 | real | kind_dyn | in | F | +!! | pt | virtual_temperature_at_Lagrangian_surface | virtual temperature at Lagrangian surface | K | 3 | real | kind_dyn | inout | F | +!! | pkz | finite-volume_mean_edge_pressure_raised_to_the_power_of_kappa | finite-volume mean edge pressure raised to the power of kappa | Pa**kappa | 3 | real | kind_dyn | inout | F | +!! | q_con | cloud_condensed_water_specific_humidity_at_Lagrangian_surface | cloud condensed water specific humidity updated by fast physics at Lagrangian surface | kg kg-1 | 3 | real | kind_dyn | inout | F | +!! | akap | kappa_dry_for_fast_physics | modified kappa for dry air, fast physics | none | 0 | real | kind_dyn | in | F | +!! | cappa | cappa_moist_gas_constant_at_Lagrangian_surface | cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/(1.+r_vir*q(i,j,k,sphum)) ) | none | 3 | real | kind_dyn | inout | F | !! | area | cell_area_for_fast_physics | area of the grid cell for fast physics | m2 | 2 | real | kind_grid | in | F | -!! | dtdt | tendency_of_air_temperature_at_Lagrangian_surface | air temperature tendency due to fast physics at Lagrangian surface | K s-1 | 3 | real | | inout | F | +!! | dtdt | tendency_of_air_temperature_at_Lagrangian_surface | air temperature tendency due to fast physics at Lagrangian surface | K s-1 | 3 | real | kind_dyn | inout | F | !! | out_dt | flag_for_tendency_of_air_temperature_at_Lagrangian_surface | flag for calculating tendency of air temperature due to fast physics | flag | 0 | logical | | in | F | !! | last_step | flag_for_the_last_step_of_k_split_remapping | flag for the last step of k-split remapping | flag | 0 | logical | | in | F | !! | do_qa | flag_for_inline_cloud_fraction_calculation | flag for the inline cloud fraction calculation | flag | 0 | logical | | in | F | -!! | qa | cloud_fraction_at_Lagrangian_surface | cloud fraction at Lagrangian surface | none | 3 | real | | out | F | +!! | qa | cloud_fraction_at_Lagrangian_surface | cloud fraction at Lagrangian surface | none | 3 | real | kind_dyn | out | F | !! | nthreads | omp_threads | number of OpenMP threads available for fast physics schemes | count | 0 | integer | | in | F | !! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | @@ -246,8 +248,8 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, implicit none ! Interface variables - real, intent(in) :: mdt - real, intent(in) :: zvir + 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 @@ -262,31 +264,31 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, integer, intent(in) :: ng logical, intent(in) :: hydrostatic logical, intent(in) :: fast_mp_consv - real, intent(inout) :: te0_2d(is:ie, js:je) - real, intent( out) :: te0(isd:ied, jsd:jed, 1:km) - real, intent(inout) :: qv(isd:ied, jsd:jed, 1:km) - real, intent(inout) :: ql(isd:ied, jsd:jed, 1:km) - real, intent(inout) :: qi(isd:ied, jsd:jed, 1:km) - real, intent(inout) :: qr(isd:ied, jsd:jed, 1:km) - real, intent(inout) :: qs(isd:ied, jsd:jed, 1:km) - real, intent(inout) :: qg(isd:ied, jsd:jed, 1:km) - real, intent(in) :: hs(isd:ied, jsd:jed) - real, intent(in) :: peln(is:ie, 1:km+1, js:je) + 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) + 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, intent(in) :: delz(isd:ied, jsd:jed, 1:kmdelz) - real, intent(in) :: delp(isd:ied, jsd:jed, 1:km) - real, intent(inout) :: pt(isd:ied, jsd:jed, 1:km) - real, intent(inout) :: pkz(is:ie, js:je, 1:km) + real(kind=kind_dyn), intent(in) :: delz(isd:ied, jsd:jed, 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, intent(inout) :: q_con(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: q_con(isd:ied, jsd:jed, 1:km) #else - real, intent(inout) :: q_con(isd:isd, jsd:jsd, 1) + real(kind=kind_dyn), intent(inout) :: q_con(isd:isd, jsd:jsd, 1) #endif - real, intent(in) :: akap + real(kind=kind_dyn), intent(in) :: akap #ifdef MOIST_CAPPA - real, intent(inout) :: cappa(isd:ied, jsd:jed, 1:km) + real(kind=kind_dyn), intent(inout) :: cappa(isd:ied, jsd:jed, 1:km) #else - real, intent(inout) :: cappa(isd:ied, jsd:jed, 1) + 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 @@ -294,17 +296,17 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, ! 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, intent(inout) :: dtdt(is:ie, js:je, 1:km) + 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, intent( out) :: qa(isd:ied, jsd:jed, 1:km) + 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, dimension(is:ie,js:je) :: dpln + real(kind=kind_dyn), dimension(is:ie,js:je) :: dpln integer :: kdelz integer :: k, j, i @@ -390,42 +392,44 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, ! Interface variables 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, delz, hs - real, intent (in), dimension (is:ie, js:je) :: dpln - real, intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: pt, qv, ql, qi, qr, qs, qg - real, intent (inout), dimension (is - ng:ie + ng, js - ng:je + 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=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, delz, hs + real(kind=kind_dyn), intent (in), dimension (is:ie, js:je) :: dpln + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: pt, qv, 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 - 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 + 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 - sdt = 0.5 * mdt ! half remapping time step + + + sdt = 0.5_kind_dyn * mdt ! half remapping time step dt_bigg = mdt ! bigg mechinism time step - tice0 = tice - 0.01 ! 273.15, standard freezing temperature + tice0 = tice - 0.01_kind_dyn ! 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_i2s = 1._kind_dyn - exp (- mdt / tau_i2s) + fac_v2l = 1._kind_dyn - exp (- sdt / tau_v2l) + fac_r2g = 1._kind_dyn - exp (- mdt / tau_r2g) + fac_l2r = 1._kind_dyn - exp (- mdt / tau_l2r) + fac_l2v = 1._kind_dyn - 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) + fac_imlt = 1._kind_dyn - exp (- sdt / tau_imlt) + fac_smlt = 1._kind_dyn - exp (- mdt / tau_smlt) ! ----------------------------------------------------------------------- !> - Define heat capacity of dry air and water vapor based on hydrostatical property. ! ----------------------------------------------------------------------- @@ -446,9 +450,9 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, q_sol (i) = qi (i, j) + qs (i, j) + qg (i, j) qpz (i) = q_liq (i) + q_sol (i) #ifdef USE_COND - pt1 (i) = pt (i, j) / ((1 + zvir * qv (i, j)) * (1 - qpz (i))) + pt1 (i) = pt (i, j) / ((1.0_kind_dyn + zvir * qv (i, j)) * (1.0_kind_dyn - qpz (i))) #else - pt1 (i) = pt (i, j) / (1 + zvir * qv (i, j)) + pt1 (i) = pt (i, j) / (1.0_kind_dyn + zvir * qv (i, j)) #endif t0 (i) = pt1 (i) ! true temperature qpz (i) = qpz (i) + qv (i, j) ! total_wat conserved in this routine @@ -469,7 +473,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, !> - Define heat capacity and latend heat coefficient. ! ----------------------------------------------------------------------- do i = is, ie - mc_air (i) = (1. - qpz (i)) * c_air ! constant + mc_air (i) = (1._kind_dyn - 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) @@ -496,7 +500,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, !> - Fix negative cloud ice with snow. ! ----------------------------------------------------------------------- do i = is, ie - if (qi (i, j) < 0.) then + if (qi (i, j) < 0._kind_dyn) then qs (i, j) = qs (i, j) + qi (i, j) qi (i, j) = 0. endif @@ -505,7 +509,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, !> - 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 + if (qi (i, j) > 1.e-8_kind_dyn .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 @@ -531,11 +535,11 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, !> - Fix negative snow with graupel or graupel with available snow. ! ----------------------------------------------------------------------- do i = is, ie - if (qs (i, j) < 0.) then + if (qs (i, j) < 0._kind_dyn) 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))) + qs (i, j) = 0._kind_dyn + elseif (qg (i, j) < 0._kind_dyn) then + tmp = min (- qg (i, j), max (0._kind_dyn, qs (i, j))) qg (i, j) = qg (i, j) + tmp qs (i, j) = qs (i, j) - tmp endif @@ -545,23 +549,25 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, !> - 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))) + if (ql (i, j) < 0._kind_dyn) then + tmp = min (- ql (i, j), max (0._kind_dyn, 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))) + elseif (qr (i, j) < 0._kind_dyn) then + tmp = min (- qr (i, j), max (0._kind_dyn, 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 + dtmp = tice - 48._kind_dyn - pt1 (i) + if (ql (i, j) > 0._kind_dyn .and. dtmp > 0._kind_dyn) 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) @@ -579,7 +585,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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.) + tcp3 (i) = lcp2 (i) + icp2 (i) * min (1._kind_dyn, dim (tice, pt1 (i)) /48._kind_dyn) enddo ! ----------------------------------------------------------------------- !> - Condensation/evaporation between water vapor and cloud water. @@ -587,15 +593,15 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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 + dq0 = (qv (i, j) - wqsat (i)) / (1._kind_dyn + tcp3 (i) * dq2dt (i)) + if (dq0 > 0._kind_dyn) 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% + factor = - min (1._kind_dyn, fac_l2v * 10._kind_dyn * (1._kind_dyn - 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) @@ -612,7 +618,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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.) + tcp3 (i) = lcp2 (i) + icp2 (i) * min (1., dim (tice, pt1 (i)) / 48._kind_dyn) enddo if (last_step) then ! ----------------------------------------------------------------------- @@ -622,17 +628,17 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, ! ----------------------------------------------------------------------- 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 + dq0 = (qv (i, j) - wqsat (i)) / (1._kind_dyn + tcp3 (i) * dq2dt (i)) + if (dq0 > 0._kind_dyn) 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% + factor = - min (1._kind_dyn, fac_l2v * 10._kind_dyn * (1._kind_dyn - qv (i, j) / wqsat (i))) ! the rh dependent factor = 1 at 90% src (i) = - min (ql (i, j), factor * dq0) endif - adj_fac = 1. + adj_fac = 1._kind_dyn qv (i, j) = qv (i, j) - src (i) ql (i, j) = ql (i, j) + src (i) q_liq (i) = q_liq (i) + src (i) @@ -654,8 +660,8 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, ! ----------------------------------------------------------------------- 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)) + if (ql (i, j) > 0._kind_dyn .and. dtmp > 0._kind_dyn) then + sink (i) = min (ql (i, j), ql (i, j) * dtmp * 0.125_kind_dyn, 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) @@ -676,8 +682,8 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, ! ----------------------------------------------------------------------- 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 + if (ql (i, j) > 0.0_kind_dyn .and. tc > 0._kind_dyn) then + sink (i) = 3.3333e-10_kind_dyn * dt_bigg * (exp (0.66_kind_dyn * tc) - 1._kind_dyn) * 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) @@ -698,9 +704,9 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, !> - 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 + dtmp = (tice - 0.1_kind_dyn) - pt1 (i) + if (qr (i, j) > 1.e-7_kind_dyn .and. dtmp > 0._kind_dyn) then + tmp = min (1._kind_dyn, (dtmp * 0.025_kind_dyn) ** 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) @@ -710,6 +716,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, pt1 (i) = pt1 (i) + sink (i) * lhi (i) / cvm (i) endif enddo + ! ----------------------------------------------------------------------- !> - Update latend heat coefficient. ! ----------------------------------------------------------------------- @@ -721,9 +728,9 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, !> - 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 + dtmp = pt1 (i) - (tice + 0.1_kind_dyn) + if (qs (i, j) > 1.e-7_kind_dyn .and. dtmp > 0._kind_dyn) then + tmp = min (1._kind_dyn, (dtmp * 0.1_kind_dyn) ** 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) @@ -736,6 +743,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, pt1 (i) = pt1 (i) - sink (i) * lhi (i) / cvm (i) endif enddo + ! ----------------------------------------------------------------------- !> - Autoconversion from cloud water to rain. ! ----------------------------------------------------------------------- @@ -746,6 +754,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, ql (i, j) = ql (i, j) - sink (i) endif enddo + ! ----------------------------------------------------------------------- !> - Update latend heat coefficient. ! ----------------------------------------------------------------------- @@ -760,25 +769,25 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, !> - Sublimation/deposition between water vapor and cloud ice. ! ----------------------------------------------------------------------- do i = is, ie - src (i) = 0. + src (i) = 0._kind_dyn 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) + src (i) = dim (qv (i, j), 1.e-6_kind_dyn) 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) + sink (i) = adj_fac * dq / (1._kind_dyn + tcp2 (i) * dqsdt) + if (qi (i, j) > 1.e-8_kind_dyn) then + pidep = sdt * dq * 349138.78_kind_dyn * exp (0.875_kind_dyn * log (qi (i, j) * den (i))) & + / (qsi * den (i) * lat2 / (0.0243_kind_dyn * rvgas * pt1 (i) ** 2) + 4.42478e4_kind_dyn) else - pidep = 0. + pidep = 0._kind_dyn endif - if (dq > 0.) then ! vapor - > ice + if (dq > 0._kind_dyn) then ! vapor - > ice tmp = tice - pt1 (i) - qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den (i) + qi_crt = qi_gen * min (qi_lim, 0.1_kind_dyn * 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) + pidep = pidep * min (1., dim (pt1 (i), t_sub) * 0.2_kind_dyn) src (i) = max (pidep, sink (i), - qi (i, j)) endif endif @@ -794,20 +803,20 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, do i = is, ie #ifdef USE_COND q_con (i, j) = q_liq (i) + q_sol (i) - tmp = 1. + zvir * qv (i, j) - pt (i, j) = pt1 (i) * tmp * (1. - q_con (i, j)) + tmp = 1._kind_dyn + zvir * qv (i, j) + pt (i, j) = pt1 (i) * tmp * (1._kind_dyn - q_con (i, j)) tmp = rdgas * tmp cappa (i, j) = tmp / (tmp + cvm (i)) #else - pt (i, j) = pt1 (i) * (1. + zvir * qv (i, j)) + pt (i, j) = pt1 (i) * (1._kind_dyn + zvir * qv (i, j)) #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))) + if (qg (i, j) < 0._kind_dyn) then + tmp = min (- qg (i, j), max (0._kind_dyn, qi (i, j))) qg (i, j) = qg (i, j) + tmp qi (i, j) = qi (i, j) - tmp endif @@ -908,18 +917,18 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, ! mixed phase: qsi = iqs1 (tin, den (i)) qsw = wqs1 (tin, den (i)) - if (q_cond (i) > 1.e-6) then + if (q_cond (i) > 1.e-6_kind_dyn) 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 + qstar (i) = rqi * qsi + (1._kind_dyn - 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)) + dw = dw_ocean + (dw_land - dw_ocean) * min (1._kind_dyn, abs (hs (i, j)) / (10._kind_dyn * 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))) + hvar (i) = min (0.2_kind_dyn, max (0.01_kind_dyn, dw * sqrt (sqrt (area (i, j)) / 100.e3_kind_dyn))) ! ----------------------------------------------------------------------- !> - calculate partial cloudiness by pdf; !! assuming subgrid linear distribution in horizontal; this is effectively a smoother for the @@ -931,37 +940,37 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, ! 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-6) then + if (rh > 0.75_kind_dyn .and. qpz (i) > 1.e-6_kind_dyn) 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-6) then + qa (i, j) = 1._kind_dyn + elseif (qstar (i) < q_plus .and. q_cond (i) > 1.e-6_kind_dyn) then qa (i, j) = ((q_plus - qstar (i)) / dq) ** 2 - qa (i, j) = min (1., qa (i, j)) + qa (i, j) = min (1._kind_dyn, qa (i, j)) else qa (i, j) = 0. endif else ! icloud_f if (qstar (i) < q_minus) then - qa (i, j) = 1. + qa (i, j) = 1._kind_dyn 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))) + qa (i, j) = (q_plus - qstar (i)) / (2._kind_dyn * dq * (1._kind_dyn - q_cond (i))) endif else - qa (i, j) = 0. + qa (i, j) = 0._kind_dyn endif ! impose minimum cloudiness if substantial q_cond (i) exist - if (q_cond (i) > 1.e-6) then + if (q_cond (i) > 1.e-6_kind_dyn) then qa (i, j) = max (cld_min, qa (i, j)) endif - qa (i, j) = min (1., qa (i, j)) + qa (i, j) = min (1._kind_dyn, qa (i, j)) endif endif else !rh @@ -970,6 +979,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, enddo endif enddo ! end j loop + end subroutine fv_sat_adj_work !! @} @@ -979,22 +989,22 @@ end subroutine fv_sat_adj_work !>@brief the function 'wqs1' computes the !! saturated specific humidity for table ii. ! ======================================================================= -real function wqs1 (ta, den) +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, intent (in) :: ta, den + real(kind=kind_dyn), intent (in) :: ta, den - real :: es, ap1, tmin + real(kind=kind_dyn) :: es, ap1, tmin integer :: it - tmin = tice - 160. - ap1 = 10. * dim (ta, tmin) + 1. - ap1 = min (2621., ap1) + tmin = tice - 160._kind_dyn + ap1 = 10._kind_dyn * dim (ta, tmin) + 1._kind_dyn + ap1 = min (2621._kind_dyn, ap1) it = ap1 es = tablew (it) + (ap1 - it) * desw (it) wqs1 = es / (rvgas * ta * den) @@ -1006,22 +1016,22 @@ end function wqs1 !>@brief the function 'wqs1' computes the saturated specific humidity !! for table iii ! ======================================================================= -real function iqs1 (ta, den) +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, intent (in) :: ta, den + real(kind=kind_dyn), intent (in) :: ta, den - real :: es, ap1, tmin + real(kind=kind_dyn) :: es, ap1, tmin integer :: it - tmin = tice - 160. - ap1 = 10. * dim (ta, tmin) + 1. - ap1 = min (2621., ap1) + tmin = tice - 160._kind_dyn + ap1 = 10._kind_dyn * dim (ta, tmin) + 1._kind_dyn + ap1 = min (2621._kind_dyn, ap1) it = ap1 es = table2 (it) + (ap1 - it) * des2 (it) iqs1 = es / (rvgas * ta * den) @@ -1033,30 +1043,30 @@ end function iqs1 !>@brief The function 'wqs2'computes the gradient of saturated specific !! humidity for table ii ! ======================================================================= -real function wqs2 (ta, den, dqdt) +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, intent (in) :: ta, den + real(kind=kind_dyn), intent (in) :: ta, den - real, intent (out) :: dqdt + real(kind=kind_dyn), intent (out) :: dqdt - real :: es, ap1, tmin + real(kind=kind_dyn) :: es, ap1, tmin integer :: it - tmin = tice - 160. - ap1 = 10. * dim (ta, tmin) + 1. - ap1 = min (2621., ap1) + tmin = tice - 160._kind_dyn + ap1 = 10._kind_dyn * dim (ta, tmin) + 1._kind_dyn + ap1 = min (2621._kind_dyn, ap1) it = ap1 es = tablew (it) + (ap1 - it) * desw (it) wqs2 = es / (rvgas * ta * den) - it = ap1 - 0.5 + it = ap1 - 0.5_kind_dyn ! finite diff, del_t = 0.1: - dqdt = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta * den) + dqdt = 10._kind_dyn * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta * den) end function wqs2 @@ -1075,25 +1085,25 @@ subroutine wqs2_vect (is, ie, ta, den, wqsat, dqdt) integer, intent (in) :: is, ie - real, intent (in), dimension (is:ie) :: ta, den + real(kind=kind_dyn), intent (in), dimension (is:ie) :: ta, den - real, intent (out), dimension (is:ie) :: wqsat, dqdt + real(kind=kind_dyn), intent (out), dimension (is:ie) :: wqsat, dqdt - real :: es, ap1, tmin + real(kind=kind_dyn) :: es, ap1, tmin integer :: i, it - tmin = tice - 160. + tmin = tice - 160._kind_dyn do i = is, ie - ap1 = 10. * dim (ta (i), tmin) + 1. - ap1 = min (2621., ap1) + ap1 = 10._kind_dyn * dim (ta (i), tmin) + 1._kind_dyn + ap1 = min (2621._kind_dyn, ap1) it = ap1 es = tablew (it) + (ap1 - it) * desw (it) wqsat (i) = es / (rvgas * ta (i) * den (i)) - it = ap1 - 0.5 + it = ap1 - 0.5_kind_dyn ! finite diff, del_t = 0.1: - dqdt (i) = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta (i) * den (i)) + dqdt (i) = 10._kind_dyn * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta (i) * den (i)) enddo end subroutine wqs2_vect @@ -1103,30 +1113,30 @@ end subroutine wqs2_vect !>@brief The function 'iqs2' computes the gradient of saturated specific !! humidity for table iii. ! ======================================================================= -real function iqs2 (ta, den, dqdt) +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, intent (in) :: ta, den + real(kind=kind_dyn), intent (in) :: ta, den - real, intent (out) :: dqdt + real(kind=kind_dyn), intent (out) :: dqdt - real :: es, ap1, tmin + real(kind=kind_dyn) :: es, ap1, tmin integer :: it - tmin = tice - 160. - ap1 = 10. * dim (ta, tmin) + 1. - ap1 = min (2621., ap1) + tmin = tice - 160._kind_dyn + ap1 = 10._kind_dyn * dim (ta, tmin) + 1._kind_dyn + ap1 = min (2621._kind_dyn, ap1) it = ap1 es = table2 (it) + (ap1 - it) * des2 (it) iqs2 = es / (rvgas * ta * den) - it = ap1 - 0.5 + it = ap1 - 0.5_kind_dyn ! finite diff, del_t = 0.1: - dqdt = 10. * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) / (rvgas * ta * den) + dqdt = 10._kind_dyn * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) / (rvgas * ta * den) end function iqs2 @@ -1136,17 +1146,17 @@ end function iqs2 subroutine qs_table (n) implicit none integer, intent (in) :: n - real (kind_grid) :: delt = 0.1 + real (kind_grid) :: delt = 0.1_kind_grid real (kind_grid) :: tmin, tem, esh20 real (kind_grid) :: wice, wh2o, fac0, fac1, fac2 real (kind_grid) :: esupc (200) integer :: i - tmin = tice - 160. + tmin = tice - 160._kind_grid ! ----------------------------------------------------------------------- ! compute es over ice between - 160 deg c and 0 deg c. ! ----------------------------------------------------------------------- do i = 1, 1600 - tem = tmin + delt * real (i - 1) + tem = tmin + delt * real (i - 1, kind=kind_grid) fac0 = (tem - tice) / (tem * tice) fac1 = fac0 * li2 fac2 = (d2ice * log (tem / tice) + fac1) / rvgas @@ -1156,7 +1166,7 @@ subroutine qs_table (n) ! compute es over water between - 20 deg c and 102 deg c. ! ----------------------------------------------------------------------- do i = 1, 1221 - tem = 253.16 + delt * real (i - 1) + tem = 253.16_kind_grid + delt * real (i - 1, kind=kind_grid) fac0 = (tem - tice) / (tem * tice) fac1 = fac0 * lv0 fac2 = (dc_vap * log (tem / tice) + fac1) / rvgas @@ -1171,9 +1181,9 @@ subroutine qs_table (n) ! 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) + tem = 253.16_kind_grid + delt * real (i - 1, kind=kind_grid) + wice = 0.05_kind_grid * (tice - tem) + wh2o = 0.05_kind_grid * (tem - 253.16_kind_grid) table (i + 1400) = wice * table (i + 1400) + wh2o * esupc (i) enddo end subroutine qs_table @@ -1184,15 +1194,16 @@ end subroutine qs_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 + real (kind_grid) :: delt = 0.1_kind_grid + real (kind_grid) :: tmin, tem, fac0, fac1, fac2,tem2 integer :: i - tmin = tice - 160. + tmin = tice - 160._kind_grid ! ----------------------------------------------------------------------- ! compute es over water ! ----------------------------------------------------------------------- do i = 1, n - tem = tmin + delt * real (i - 1) + tem2 = real (i - 1, kind=kind_grid) + tem = tmin + delt * real (i - 1, kind=kind_grid) fac0 = (tem - tice) / (tem * tice) fac1 = fac0 * lv0 fac2 = (dc_vap * log (tem / tice) + fac1) / rvgas @@ -1206,12 +1217,12 @@ end subroutine qs_tablew subroutine qs_table2 (n) implicit none integer, intent (in) :: n - real (kind_grid) :: delt = 0.1 + real (kind_grid) :: delt = 0.1_kind_grid real (kind_grid) :: tmin, tem0, tem1, fac0, fac1, fac2 integer :: i, i0, i1 - tmin = tice - 160. + tmin = tice - 160._kind_grid do i = 1, n - tem0 = tmin + delt * real (i - 1) + tem0 = tmin + delt * real (i - 1, kind=kind_grid) fac0 = (tem0 - tice) / (tem0 * tice) if (i <= 1600) then ! ----------------------------------------------------------------------- @@ -1233,8 +1244,8 @@ subroutine qs_table2 (n) ! ----------------------------------------------------------------------- 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)) + tem0 = 0.25_kind_grid * (table2 (i0 - 1) + 2._kind_grid * table (i0) + table2 (i0 + 1)) + tem1 = 0.25_kind_grid * (table2 (i1 - 1) + 2._kind_grid * table (i1) + table2 (i1 + 1)) table2 (i0) = tem0 table2 (i1) = tem1 end subroutine qs_table2 diff --git a/physics/gfdl_fv_sat_adj_pre.F90 b/physics/gfdl_fv_sat_adj_pre.F90 deleted file mode 100644 index 3fe1fa5d1..000000000 --- a/physics/gfdl_fv_sat_adj_pre.F90 +++ /dev/null @@ -1,39 +0,0 @@ -!> \file gfdl_fv_sat_adj_pre.F90 -!! Contains code related to preparing the gfdl_fv_sat_adj runs. - - module fv_sat_adj_pre - - contains - - subroutine fv_sat_adj_pre_init () - end subroutine fv_sat_adj_pre_init - - subroutine fv_sat_adj_pre_finalize() - end subroutine fv_sat_adj_pre_finalize - -!> \section arg_table_fv_sat_adj_pre_run Argument Table -!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | -!! |----------------|--------------------------------------------------------|---------------------------------------------------------|---------------|------|------------------------|-----------|--------|----------| -!! | Interstitial | CCPP_Interstitial_type | derived type CCPP_interstitial_type | DDT | 0 | CCPP_interstitial_type | | inout | F | -!! | errmsg | ccpp_error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | -!! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | -!! - subroutine fv_sat_adj_pre_run (Interstitial, errmsg, errflg) - - use CCPP_typedefs, only: CCPP_interstitial_type - - implicit none - - ! interface variables - type(CCPP_interstitial_type), intent(inout) :: Interstitial - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - - errmsg = '' - errflg = 0 - - call Interstitial%reset() - - end subroutine fv_sat_adj_pre_run - - end module fv_sat_adj_pre diff --git a/physics/physcons.F90 b/physics/physcons.F90 index 06f9309c9..937aab9f6 100644 --- a/physics/physcons.F90 +++ b/physics/physcons.F90 @@ -36,23 +36,13 @@ !! constants for GCM models. module physcons ! - use machine, only : kind_phys + use machine, only : kind_phys + use CCPP_typedefs, only : kind_dyn ! implicit none ! public -! \section arg_table_physcons Argument Table -! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | -! |-----------------------|--------------------------------------------------------|---------------------------------------------------------|---------------|------|-------------------|-----------|--------|----------| -! | con_cp | specific_heat_of_dry_air_at_constant_pressure | specific heat of dry air at constant pressure | J kg-1 K-1 | 0 | real | kind_phys | none | F | -! | con_fvirt | ratio_of_vapor_to_dry_air_gas_constants_minus_one | rv/rd - 1 (rv = ideal gas constant for water vapor) | none | 0 | real | kind_phys | none | F | -! | con_g | gravitational_acceleration | gravitational acceleration | m s-2 | 0 | real | kind_phys | none | F | -! | con_pi | pi | ratio of a circle's circumference to its diameter | radians | 0 | real | kind_phys | none | F | -! | con_rd | gas_constant_dry_air | ideal gas constant for dry air | J kg-1 K-1 | 0 | real | kind_phys | none | F | -! | con_rv | gas_constant_water_vapor | ideal gas constant for water vapor | J kg-1 K-1 | 0 | real | kind_phys | none | F | -! - !> \name Math constants real(kind=kind_phys),parameter:: con_pi =3.1415926535897931 !< pi real(kind=kind_phys),parameter:: con_sqrt2 =1.414214e+0 !< square root of 2 @@ -67,6 +57,8 @@ module physcons real(kind=kind_phys),parameter:: con_solr_old =1.3660e+3 !< solar constant (\f$W/m^{2}\f$)-Liu(2002) real(kind=kind_phys),parameter:: con_solr =1.3608e+3 !< solar constant (\f$W/m^{2}\f$)-nasa-sorce Tim(2008) ! real(kind=kind_phys),parameter:: con_solr =1.36742732e+3 ! solar constant (W/m2)-gfdl(1989) - OPR as of Jan 2006 + ! Selected geophysics/astronomy constants with kind=kind_dyn + real(kind=kind_dyn), parameter:: con_g_dyn =9.80665e+0 !< gravity (\f$m/s^{2}\f$) !> \name Thermodynamics constants real(kind=kind_phys),parameter:: con_rgas =8.314472 !< molar gas constant (\f$J/mol/K\f$) @@ -86,6 +78,12 @@ module physcons real(kind=kind_phys),parameter:: con_jcal =4.1855E+0 !< joules per calorie real(kind=kind_phys),parameter:: con_rhw0 =1022.0 !< sea water reference density (\f$kg/m^{3}\f$) real(kind=kind_phys),parameter:: con_epsq =1.0E-12 !< min q for computing precip type + ! Selected thermodynamics constants with kind=kind_dyn + real(kind=kind_dyn), parameter:: con_rd_dyn =2.8705e+2 !< gas constant air (\f$J/kg/K\f$) + real(kind=kind_dyn), parameter:: con_rv_dyn =4.6150e+2 !< gas constant H2O (\f$J/kg/K\f$) + real(kind=kind_dyn), parameter:: con_cp_dyn =1.0046e+3 !< spec heat air at p (\f$J/kg/K\f$) + real(kind=kind_dyn), parameter:: con_hvap_dyn =2.5000e+6 !< lat heat H2O cond (\f$J/kg\f$) + real(kind=kind_dyn), parameter:: con_hfus_dyn =3.3358e+5 !< lat heat H2O fusion (\f$J/kg\f$) !> \name Secondary constants real(kind=kind_phys),parameter:: con_rocp =con_rd/con_cp @@ -105,11 +103,11 @@ module physcons real(kind=kind_phys),parameter:: con_sbc =5.670400e-8 !< stefan-boltzmann (\f$W/m^{2}/K^{4}\f$) real(kind=kind_phys),parameter:: con_avgd =6.0221415e23 !< avogadro constant (\f$mol^{-1}\f$) real(kind=kind_phys),parameter:: con_gasv =22413.996e-6 !< vol of ideal gas at 273.15K, 101.325kPa (\f$m^{3}/mol\f$) -! real(kind=kind_phys),parameter:: con_amd =28.970 ! molecular wght of dry air (g/mol) +! real(kind=kind_phys),parameter:: con_amd =28.970 !< molecular wght of dry air (g/mol) real(kind=kind_phys),parameter:: con_amd =28.9644 !< molecular wght of dry air (\f$g/mol\f$) real(kind=kind_phys),parameter:: con_amw =18.0154 !< molecular wght of water vapor (\f$g/mol\f$) real(kind=kind_phys),parameter:: con_amo3 =47.9982 !< molecular wght of o3 (\f$g/mol\f$) -! real(kind=kind_phys),parameter:: con_amo3 =48.0 ! molecular wght of o3 (g/mol) +! real(kind=kind_phys),parameter:: con_amo3 =48.0 !< molecular wght of o3 (g/mol) real(kind=kind_phys),parameter:: con_amco2 =44.011 !< molecular wght of co2 (\f$g/mol\f$) real(kind=kind_phys),parameter:: con_amo2 =31.9999 !< molecular wght of o2 (\f$g/mol\f$) real(kind=kind_phys),parameter:: con_amch4 =16.043 !< molecular wght of ch4 (\f$g/mol\f$)