diff --git a/physics/gfdl_cloud_microphys.F90 b/physics/gfdl_cloud_microphys.F90 index 31736d6b9..e6e781fcf 100644 --- a/physics/gfdl_cloud_microphys.F90 +++ b/physics/gfdl_cloud_microphys.F90 @@ -230,6 +230,10 @@ subroutine gfdl_cloud_microphys_run( & real(kind=kind_phys), dimension(:,:), allocatable :: den real(kind=kind_phys) :: onebg real(kind=kind_phys) :: tem +#ifdef TRANSITION + real(kind=kind_phys), volatile :: volatile_var +#endif + ! Initialize CCPP error handling variables errmsg = '' @@ -317,10 +321,17 @@ subroutine gfdl_cloud_microphys_run( & ! calculate fraction of frozen precipitation using unscaled ! values of rain0, ice0, snow0, graupel0 (for bit-for-bit) do i=1,im +#ifdef TRANSITION + volatile_var = rain0(i)+snow0(i)+ice0(i)+graupel0(i) + prcp0(i) = volatile_var * tem + if ( volatile_var * tem > rainmin ) then + sr(i) = (snow0(i) + ice0(i) + graupel0(i)) / volatile_var +#else prcp0(i) = (rain0(i)+snow0(i)+ice0(i)+graupel0(i)) * tem if ( prcp0(i) > rainmin ) then sr(i) = (snow0(i) + ice0(i) + graupel0(i)) & / (rain0(i) + snow0(i) + ice0(i) + graupel0(i)) +#endif else sr(i) = 0.0 endif diff --git a/physics/gfdl_fv_sat_adj.F90 b/physics/gfdl_fv_sat_adj.F90 index 4bc8ba717..3be08107e 100644 --- a/physics/gfdl_fv_sat_adj.F90 +++ b/physics/gfdl_fv_sat_adj.F90 @@ -61,18 +61,27 @@ module fv_sat_adj 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, kind_dyn + use machine, only: kind_grid, kind_dyn use gfdl_cloud_microphys_mod, only: ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt use gfdl_cloud_microphys_mod, only: icloud_f, sat_adj0, t_sub, cld_min use gfdl_cloud_microphys_mod, only: tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r use gfdl_cloud_microphys_mod, only: rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land + +#ifdef MULTI_GASES + use ccpp_multi_gases_mod, only: multi_gases_init, & + multi_gases_finalize, & + virq_qpz, vicpqd_qpz, & + vicvqd_qpz, num_gas +#endif + implicit none + private public fv_sat_adj_init, fv_sat_adj_run, fv_sat_adj_finalize + logical :: is_initialized = .false. + real(kind=kind_dyn), parameter :: rrg = -rdgas/grav ! real, parameter :: cp_air = cp_air ! 1004.6, heat capacity of dry air at constant pressure, come from constants_mod real(kind=kind_dyn), parameter :: cp_vap = 4.0 * rvgas !< 1846.0, heat capacity of water vapor at constant pressure @@ -109,22 +118,35 @@ module fv_sat_adj !>\brief The subroutine 'fv_sat_adj_init' initializes lookup tables for the saturation mixing ratio. !! \section arg_table_fv_sat_adj_init Argument Table -!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | -!! |----------------|---------------------------------------------------------------|----------------------------------------------------------------------------------------|---------|------|-----------|-----------|--------|----------| -!! | do_sat_adj | flag_for_saturation_adjustment_for_microphysics_in_dynamics | flag for saturation adjustment for microphysics in dynamics | none | 0 | logical | | in | F | -!! | kmp | top_layer_index_for_fast_physics | top_layer_inder_for_gfdl_mp | index | 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 | +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |----------------|---------------------------------------------------------------|----------------------------------------------------------------------------------------|------------|------|-----------|-----------|--------|----------| +!! | do_sat_adj | flag_for_saturation_adjustment_for_microphysics_in_dynamics | flag for saturation adjustment for microphysics in dynamics | none | 0 | logical | | in | F | +!! | kmp | top_layer_index_for_fast_physics | top_layer_inder_for_gfdl_mp | index | 0 | integer | | in | F | +!! | nwat | number_of_water_species | number of water species | count | 0 | integer | | in | F | +!! | ngas | number_of_gases_for_multi_gases_physics | number of gases for multi gases physics | count | 0 | integer | | in | F | +!! | rilist | gas_constants_for_multi_gases_physics | gas constants for multi gases physics | J kg-1 K-1 | 1 | real | kind_dyn | in | F | +!! | cpilist | specific_heat_capacities_for_multi_gases_physics | specific heat capacities for multi gases physics | J kg-1 K-1 | 1 | real | kind_dyn | in | F | +!! | mpirank | mpi_rank_for_fast_physics | current MPI-rank for fast physics schemes | index | 0 | integer | | in | F | +!! | mpiroot | mpi_root_for_fast_physics | master MPI-rank for fast physics schemes | index | 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 | !! -subroutine fv_sat_adj_init(do_sat_adj, kmp, errmsg, errflg) +subroutine fv_sat_adj_init(do_sat_adj, kmp, nwat, ngas, rilist, cpilist, & + mpirank, mpiroot, errmsg, errflg) implicit none ! Interface variables - logical, intent (in) :: do_sat_adj - integer, intent (in) :: kmp - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg + logical, intent(in ) :: do_sat_adj + integer, intent(in ) :: kmp + integer, intent(in ) :: nwat + integer, intent(in ) :: ngas + real(kind_dyn), intent(in ) :: rilist(0:ngas) + real(kind_dyn), intent(in ) :: cpilist(0:ngas) + integer, intent(in ) :: mpirank + integer, intent(in ) :: mpiroot + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg ! Local variables integer, parameter :: length = 2621 @@ -141,7 +163,7 @@ subroutine fv_sat_adj_init(do_sat_adj, kmp, errmsg, errflg) return end if - if (allocated(table)) return + if (is_initialized) return ! generate es table (dt = 0.1 deg c) @@ -162,6 +184,12 @@ subroutine fv_sat_adj_init(do_sat_adj, kmp, errmsg, errflg) des2 (length) = des2 (length - 1) desw (length) = desw (length - 1) +#ifdef MULTI_GASES + call multi_gases_init(ngas,nwat,rilist,cpilist,mpirank==mpiroot) +#endif + + is_initialized = .true. + end subroutine fv_sat_adj_init !\ingroup fast_sat_adj @@ -183,12 +211,20 @@ subroutine fv_sat_adj_finalize (errmsg, errflg) errmsg = '' errflg = 0 + if (.not.is_initialized) return + if (allocated(table )) deallocate(table ) if (allocated(table2)) deallocate(table2) if (allocated(tablew)) deallocate(tablew) if (allocated(des2 )) deallocate(des2 ) if (allocated(desw )) deallocate(desw ) +#ifdef MULTI_GASES + call multi_gases_finalize() +#endif + + is_initialized = .false. + end subroutine fv_sat_adj_finalize !>\defgroup fast_sat_adj GFDL MP Fast Physics @@ -222,6 +258,8 @@ end subroutine fv_sat_adj_finalize !! | 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 | 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 | +!! | ngas | number_of_gases_for_multi_gases_physics | number of gases for multi gases physics | count | 0 | integer | | in | F | +!! | qvi | gas_tracers_for_multi_gas_physics_at_Lagrangian_surface | gas tracers for multi gas physics at Lagrangian surface | kg kg-1 | 4 | real | kind_dyn | inout | 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 | @@ -248,9 +286,10 @@ end subroutine fv_sat_adj_finalize !! | errflg | ccpp_error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | !! subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, jsd, jed, & - ng, hydrostatic, fast_mp_consv, te0_2d, te0, qv, ql, qi, qr, qs, qg, & - hs, peln, delz, delp, pt, pkz, q_con, akap, cappa, area, dtdt, & - out_dt, last_step, do_qa, qa, nthreads, errmsg, errflg) + ng, hydrostatic, fast_mp_consv, te0_2d, te0, ngas, qvi, qv, ql, qi, qr, & + qs, qg, hs, peln, delz, delp, pt, pkz, q_con, akap, cappa, area, dtdt, & + out_dt, last_step, do_qa, qa, & + nthreads, errmsg, errflg) implicit none @@ -273,6 +312,9 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, logical, intent(in) :: fast_mp_consv real(kind=kind_dyn), intent(inout) :: te0_2d(is:ie, js:je) real(kind=kind_dyn), intent( out) :: te0(isd:ied, jsd:jed, 1:km) + ! If multi-gases physics are not used, ngas is one and qvi identical to qv + integer, intent(in) :: ngas + real(kind=kind_dyn), intent(inout) :: qvi(isd:ied, jsd:jed, 1:km, 1:ngas) real(kind=kind_dyn), intent(inout) :: qv(isd:ied, jsd:jed, 1:km) real(kind=kind_dyn), intent(inout) :: ql(isd:ied, jsd:jed, 1:km) real(kind=kind_dyn), intent(inout) :: qi(isd:ied, jsd:jed, 1:km) @@ -334,7 +376,7 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, !$OMP area,delp,pt,hs,qg,qs,qr,qi, & !$OMP ql,qv,te0,fast_mp_consv, & !$OMP hydrostatic,ng,zvir,pkz, & -!$OMP akap,te0_2d) & +!$OMP akap,te0_2d,ngas,qvi) & #ifdef TRANSITION !$OMP private(volatile_var) & #endif @@ -354,7 +396,13 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, kdelz = k end if call fv_sat_adj_work(abs(mdt), zvir, is, ie, js, je, ng, hydrostatic, fast_mp_consv, & - te0(isd,jsd,k), qv(isd,jsd,k), ql(isd,jsd,k), qi(isd,jsd,k), & + te0(isd,jsd,k), & +#ifdef MULTI_GASES + qvi(isd,jsd,k,1:ngas), & +#else + qv(isd,jsd,k), & +#endif + ql(isd,jsd,k), qi(isd,jsd,k), & qr(isd,jsd,k), qs(isd,jsd,k), qg(isd,jsd,k), & hs, dpln, delz(isd:,jsd:,kdelz), pt(isd,jsd,k), delp(isd,jsd,k),& q_con(isd:,jsd:,k), cappa(isd:,jsd:,k), area, dtdt(is,js,k), & @@ -371,11 +419,20 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je, #endif #else #ifdef TRANSITION +#ifdef MULTI_GASES + volatile_var = log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)) + pkz(i,j,k) = exp(akap*(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*volatile_var) +#else volatile_var = log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)) pkz(i,j,k) = exp(akap*volatile_var) +#endif +#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 #endif enddo enddo @@ -407,8 +464,13 @@ end subroutine fv_sat_adj_run !> This subroutine includes the entity of the fast saturation adjustment processes. !>\section fast_gen GFDL Cloud Fast Physics General Algorithm !! @{ -subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, & - te0, qv, ql, qi, qr, qs, qg, hs, dpln, delz, pt, dp, q_con, cappa, & +subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, te0, & +#ifdef MULTI_GASES + qvi, & +#else + qv, & +#endif + ql, qi, qr, qs, qg, hs, dpln, delz, pt, dp, q_con, cappa, & area, dtdt, out_dt, last_step, do_qa, qa) implicit none @@ -419,13 +481,22 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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) :: pt +#ifdef MULTI_GASES + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng, 1:1, 1:num_gas) :: qvi +#else + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: qv +#endif + real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: ql, qi, qr, qs, qg real(kind=kind_dyn), intent (inout), dimension (is - ng:ie + ng, js - ng:je + ng) :: q_con, cappa real(kind=kind_dyn), intent (inout), dimension (is:ie, js:je) :: dtdt real(kind=kind_dyn), intent (out), dimension (is - ng:ie + ng, js - ng:je + ng) :: qa, te0 real (kind_grid), intent (in), dimension (is - ng:ie + ng, js - ng:je + ng) :: area ! Local variables +#ifdef MULTI_GASES + real, dimension (is - ng:ie + ng, js - ng:je + ng) :: qv +#endif real(kind=kind_dyn), dimension (is:ie) :: wqsat, dq2dt, qpz, cvm, t0, pt1, qstar real(kind=kind_dyn), dimension (is:ie) :: icp2, lcp2, tcp2, tcp3 real(kind=kind_dyn), dimension (is:ie) :: den, q_liq, q_sol, q_cond, src, sink, hvar @@ -438,23 +509,33 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, real(kind=kind_dyn) :: factor, qim, tice0, c_air, c_vap, dw integer :: i, j +#ifdef MULTI_GASES + qv(:,:) = qvi(:,:,1,1) +#endif sdt = 0.5 * mdt ! half remapping time step dt_bigg = mdt ! bigg mechinism time step + tice0 = tice - 0.01 ! 273.15, standard freezing temperature + ! ----------------------------------------------------------------------- !> - Define conversion scalar / factor. ! ----------------------------------------------------------------------- + fac_i2s = 1. - exp (- mdt / tau_i2s) fac_v2l = 1. - exp (- sdt / tau_v2l) fac_r2g = 1. - exp (- mdt / tau_r2g) fac_l2r = 1. - exp (- mdt / tau_l2r) + fac_l2v = 1. - exp (- sdt / tau_l2v) fac_l2v = min (sat_adj0, fac_l2v) + fac_imlt = 1. - exp (- sdt / tau_imlt) fac_smlt = 1. - exp (- mdt / tau_smlt) + ! ----------------------------------------------------------------------- !> - Define heat capacity of dry air and water vapor based on hydrostatical property. ! ----------------------------------------------------------------------- + if (hydrostatic) then c_air = cp_air c_vap = cp_vap @@ -466,22 +547,30 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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)) @@ -491,21 +580,35 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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 @@ -513,23 +616,30 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, #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)) @@ -546,16 +656,20 @@ 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. ! ----------------------------------------------------------------------- + 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) @@ -566,10 +680,13 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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))) @@ -581,6 +698,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, qr (i, j) = qr (i, j) + tmp endif enddo + ! ----------------------------------------------------------------------- !> - Enforce complete freezing of cloud water to cloud ice below - 48 c. ! ----------------------------------------------------------------------- @@ -597,9 +715,11 @@ 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. ! ----------------------------------------------------------------------- + do i = is, ie lhl (i) = lv00 + d0_vap * pt1 (i) lhi (i) = li00 + dc_ice * pt1 (i) @@ -607,10 +727,13 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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)) @@ -625,14 +748,19 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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) @@ -640,13 +768,17 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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 @@ -660,24 +792,32 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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 @@ -690,16 +830,20 @@ 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. ! ----------------------------------------------------------------------- + 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 @@ -713,16 +857,20 @@ 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. ! ----------------------------------------------------------------------- + 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 @@ -736,16 +884,20 @@ 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. ! ----------------------------------------------------------------------- + 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 @@ -762,9 +914,11 @@ 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. ! ----------------------------------------------------------------------- + do i = is, ie if (ql (i, j) > ql0_max) then sink (i) = fac_l2r * (ql (i, j) - ql0_max) @@ -772,9 +926,11 @@ 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. ! ----------------------------------------------------------------------- + do i = is, ie lhi (i) = li00 + dc_ice * pt1 (i) lhl (i) = lv00 + d0_vap * pt1 (i) @@ -782,9 +938,11 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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 @@ -809,28 +967,44 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, endif endif qv (i, j) = qv (i, j) - src (i) +#ifdef MULTI_GASES + qvi(i,j,1,1) = qv(i,j) +#endif qi (i, j) = qi (i, j) + src (i) q_sol (i) = q_sol (i) + src (i) cvm (i) = mc_air (i) + qv (i, j) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice pt1 (i) = pt1 (i) + src (i) * (lhl (i) + lhi (i)) / cvm (i) enddo + ! ----------------------------------------------------------------------- !> - Virtual temperature updated. ! ----------------------------------------------------------------------- + do i = is, ie #ifdef USE_COND q_con (i, j) = q_liq (i) + q_sol (i) +#ifdef MULTI_GASES + pt (i, j) = pt1 (i) * virq_qpz(qvi(i,j,1,1:num_gas),q_con(i,j)) +#else tmp = 1. + zvir * qv (i, j) pt (i, j) = pt1 (i) * tmp * (1. - q_con (i, j)) +#endif tmp = rdgas * tmp cappa (i, j) = tmp / (tmp + cvm (i)) +#else +#ifdef MULTI_GASES + q_con (i, j) = q_liq (i) + q_sol (i) + pt (i, j) = pt1 (i) * virq_qpz(qvi(i,j,1,1:num_gas),q_con(i,j)) * (1. - q_con(i,j)) #else pt (i, j) = pt1 (i) * (1. + zvir * qv (i, j)) +#endif #endif enddo + ! ----------------------------------------------------------------------- !> - Fix negative graupel with available cloud ice. ! ----------------------------------------------------------------------- + do i = is, ie if (qg (i, j) < 0.) then tmp = min (- qg (i, j), max (0., qi (i, j))) @@ -838,9 +1012,11 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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 @@ -849,30 +1025,41 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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) @@ -880,13 +1067,17 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, lcp2 (i) = lhl (i) / cvm (i) icp2 (i) = lhi (i) / cvm (i) enddo + ! ----------------------------------------------------------------------- !> - Compute cloud fraction. ! ----------------------------------------------------------------------- + if (do_qa .and. last_step) then + ! ----------------------------------------------------------------------- !> - If it is the last step, combine water species. ! ----------------------------------------------------------------------- + if (rad_snow) then if (rad_graupel) then do i = is, ie @@ -914,16 +1105,21 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, 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 + 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) + ! ----------------------------------------------------------------------- ! determine saturated specific humidity ! ----------------------------------------------------------------------- + if (tin <= t_wfr) then ! ice phase: qstar (i) = iqs1 (tin, den (i)) @@ -946,17 +1142,21 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, dw = dw_ocean + (dw_land - dw_ocean) * min (1., abs (hs (i, j)) / (10. * grav)) !> - "scale - aware" subgrid variability: 100 - km as the base hvar (i) = min (0.2, max (0.01, dw * sqrt (sqrt (area (i, j)) / 100.e3))) + ! ----------------------------------------------------------------------- !> - calculate partial cloudiness by pdf; !! assuming subgrid linear distribution in horizontal; this is effectively a smoother for the !! binary cloud scheme; qa = 0.5 if qstar (i) == qpz ! ----------------------------------------------------------------------- + rh = qpz (i) / qstar (i) + ! ----------------------------------------------------------------------- ! icloud_f = 0: bug - fixed ! icloud_f = 1: old fvgfs gfdl) mp implementation ! icloud_f = 2: binary cloud scheme (0 / 1) ! ----------------------------------------------------------------------- + if (rh > 0.75 .and. qpz (i) > 1.e-6) then dq = hvar (i) * qpz (i) q_plus = qpz (i) + dq @@ -970,7 +1170,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, else qa (i, j) = 0. endif - else ! icloud_f + else if (qstar (i) < q_minus) then qa (i, j) = 1. else @@ -990,11 +1190,14 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, qa (i, j) = min (1., qa (i, j)) endif endif - else !rh + else qa (i, j) = 0. endif + enddo + endif + enddo ! end j loop end subroutine fv_sat_adj_work @@ -1156,21 +1359,29 @@ real(kind=kind_dyn) function iqs2 (ta, den, dqdt) end function iqs2 +! ======================================================================= !>\ingroup fast_sat_adj !! saturation water vapor pressure table i ! 3 - phase table +! ======================================================================= + subroutine qs_table (n) + implicit none + integer, intent (in) :: n real (kind_grid) :: delt = 0.1 real (kind_grid) :: tmin, tem, esh20 real (kind_grid) :: wice, wh2o, fac0, fac1, fac2 real (kind_grid) :: esupc (200) integer :: i + tmin = tice - 160. + ! ----------------------------------------------------------------------- ! compute es over ice between - 160 deg c and 0 deg c. ! ----------------------------------------------------------------------- + do i = 1, 1600 tem = tmin + delt * real (i - 1) fac0 = (tem - tice) / (tem * tice) @@ -1178,9 +1389,11 @@ subroutine qs_table (n) 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) @@ -1193,30 +1406,41 @@ subroutine qs_table (n) table (i + 1400) = esh20 endif enddo + ! ----------------------------------------------------------------------- ! derive blended es over ice and supercooled water between - 20 deg c and 0 deg c ! ----------------------------------------------------------------------- + do i = 1, 200 tem = 253.16 + delt * real (i - 1) wice = 0.05 * (tice - tem) wh2o = 0.05 * (tem - 253.16) table (i + 1400) = wice * table (i + 1400) + wh2o * esupc (i) enddo + end subroutine qs_table +! ======================================================================= !>\ingroup fast_sat_adj !! saturation water vapor pressure table ii. ! 1 - phase table +! ======================================================================= + subroutine qs_tablew (n) + implicit none + integer, intent (in) :: n real (kind_grid) :: delt = 0.1 real (kind_grid) :: tmin, tem, fac0, fac1, fac2 integer :: i + tmin = tice - 160. + ! ----------------------------------------------------------------------- ! compute es over water ! ----------------------------------------------------------------------- + do i = 1, n tem = tmin + delt * real (i - 1) fac0 = (tem - tice) / (tem * tice) @@ -1224,18 +1448,26 @@ subroutine qs_tablew (n) fac2 = (dc_vap * log (tem / tice) + fac1) / rvgas tablew (i) = e00 * exp (fac2) enddo + end subroutine qs_tablew +! ======================================================================= !>\ingroup fast_sat_adj !! saturation water vapor pressure table iii. ! 2 - phase table +! ======================================================================= + subroutine qs_table2 (n) + implicit none + integer, intent (in) :: n real (kind_grid) :: delt = 0.1 real (kind_grid) :: tmin, tem0, tem1, fac0, fac1, fac2 integer :: i, i0, i1 + tmin = tice - 160. + do i = 1, n tem0 = tmin + delt * real (i - 1) fac0 = (tem0 - tice) / (tem0 * tice) @@ -1254,15 +1486,18 @@ subroutine qs_table2 (n) endif table2 (i) = e00 * exp (fac2) enddo + ! ----------------------------------------------------------------------- ! smoother around 0 deg c ! ----------------------------------------------------------------------- + i0 = 1600 i1 = 1601 tem0 = 0.25 * (table2 (i0 - 1) + 2. * table (i0) + table2 (i0 + 1)) tem1 = 0.25 * (table2 (i1 - 1) + 2. * table (i1) + table2 (i1 + 1)) table2 (i0) = tem0 table2 (i1) = tem1 + end subroutine qs_table2 end module fv_sat_adj diff --git a/physics/gfs_phy_tracer_config.f b/physics/gfs_phy_tracer_config.f index 8795a921d..8ed7443d3 100644 --- a/physics/gfs_phy_tracer_config.f +++ b/physics/gfs_phy_tracer_config.f @@ -1,3 +1,5 @@ +#undef MULTI_GASES + ! !! ! Module: gfs_phy_tracer_config ! @@ -143,8 +145,14 @@ subroutine tracer_config_init (ntrac,ntoz,ntcw,ncld, if(ntrnc > 0) gfs_phy_tracer%vname(ntrnc) = 'rnc' if(ntsnc > 0) gfs_phy_tracer%vname(ntsnc) = 'snc' if(ntke > 0) gfs_phy_tracer%vname(ntke) = 'tke' +#ifdef MULTI_GASES + print *,' ++++ nto nto2 ',nto,nto2 + if(nto > 0) gfs_phy_tracer%vname(nto) = 'spfo' + if(nto2 > 0) gfs_phy_tracer%vname(nto2) = 'spfo2' +#else if(nto > 0) gfs_phy_tracer%vname(nto) = 'o' if(nto2 > 0) gfs_phy_tracer%vname(nto2) = 'o2' +#endif gfs_phy_tracer%fscav(1:gfs_phy_tracer%ntrac_met) = 0. diff --git a/physics/module_gfdl_cloud_microphys.F90 b/physics/module_gfdl_cloud_microphys.F90 index 009a50f39..6ead015df 100644 --- a/physics/module_gfdl_cloud_microphys.F90 +++ b/physics/module_gfdl_cloud_microphys.F90 @@ -855,8 +855,8 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & if (fix_negative) & call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz) - m2_rain (:, :) = 0. - m2_sol (:, :) = 0. + m2_rain (i, :) = 0. + m2_sol (i, :) = 0. !> - Do loop on cloud microphysics sub time step. do n = 1, ntimes @@ -942,6 +942,10 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & enddo + ! convert units from Pa*kg/kg to kg/m^2/s + m2_rain (i, :) = m2_rain (i, :) * rdt * rgrav + m2_sol (i, :) = m2_sol (i, :) * rdt * rgrav + ! ----------------------------------------------------------------------- !> - Calculate momentum transportation during sedimentation. ! note: dp1 is dry mass; dp0 is the old moist (total) mass diff --git a/physics/multi_gases.F90 b/physics/multi_gases.F90 new file mode 100644 index 000000000..1e62c89a1 --- /dev/null +++ b/physics/multi_gases.F90 @@ -0,0 +1,331 @@ +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the FV3 dynamical core. +!* +!* The FV3 dynamical core is free software: you can redistribute it +!* and/or modify it under the terms of the +!* GNU Lesser General Public License as published by the +!* Free Software Foundation, either version 3 of the License, or +!* (at your option) any later version. +!* +!* The FV3 dynamical core is distributed in the hope that 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 FV3 dynamical core. +!* If not, see . +!*********************************************************************** + +!>@brief The module 'multi_gases' peforms multi constitutents computations. +!>@author H.-M. H. Juang, NOAA/NWS/NCEP/EMC + +module ccpp_multi_gases_mod + +! Modules Included: +! +! +! +! +! +! +! +! +! +!
Module NameFunctions Included
constants_modrdgas, cp_air
+ use machine, only: kind_dyn + ! DH* TODO - MAKE THIS INPUT ARGUMENTS *DH + use physcons, only : rdgas => con_rd_dyn, & + cp_air => con_cp_dyn + ! *DH + + implicit none + integer num_gas + integer ind_gas + integer num_wat + integer sphum, sphump1 + real(kind_dyn), allocatable :: vir(:) + real(kind_dyn), allocatable :: vicp(:) + real(kind_dyn), allocatable :: vicv(:) + + private num_wat, sphum, sphump1 + public vir, vicp, vicv, ind_gas, num_gas + public multi_gases_init + public virq + public virq_max + public virqd + public vicpqd + public virq_nodq + public virq_qpz + public vicpqd_qpz + public vicvqd + public vicvqd_qpz + + CONTAINS +! -------------------------------------------------------- + subroutine multi_gases_init(ngas, nwat, ri, cpi, is_master) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: vir(i): ri/rdgas - r0/rdgas +! vir(0): r0/rdgas +! vicp(i): cpi/cp_air - cp0/cp_air +! vicp(0): cp0/cp_air +! cv_air = cp_air - rdgas +! vicv(i): cvi/cv_air - cv0/cv_air +! vicv(0): cv0/cv_air +!-------------------------------------------- + integer, intent(in):: ngas, nwat + real, intent(in):: ri(0:ngas) + real, intent(in):: cpi(0:ngas) + logical, intent(in):: is_master +! Local: + integer n + real cvi(0:ngas) + real cv_air + logical :: default_gas=.false. + + sphum = 1 + sphump1 = sphum+1 + num_wat = nwat + ind_gas = num_wat+1 + do n=0,ngas + if( ri(n).ne.0.0 .or. cpi(n).ne.0.0 ) num_gas=n + enddo + if ( num_gas.eq.1 ) default_gas=.true. + allocate( vir (0:num_gas) ) + allocate( vicp(0:num_gas) ) + allocate( vicv(0:num_gas) ) + + cv_air = cp_air - rdgas + do n=0,num_gas + cvi(n) = cpi(n) - ri(n) + enddo + + vir (0) = ri(0)/rdgas + vicp(0) = cpi(0)/cp_air + vicv(0) = cvi(0)/cv_air + if( default_gas ) then + vir (0) = 1.0 + vicp(0) = 1.0 + vicv(0) = 1.0 + endif + do n=1,num_gas + vir(n) = 0.0 + if( ri(n).gt.0.0 ) vir (n) = ri(n)/rdgas - vir (0) + vicp(n) = 0.0 + if( cpi(n).gt.0.0 ) vicp(n) = cpi(n)/cp_air - vicp(0) + vicv(n) = 0.0 + if( cvi(n).gt.0.0 ) vicv(n) = cvi(n)/cv_air - vicv(0) + enddo + + if( is_master ) then + write(*,*) ' multi_gases_init with ind_gas=',ind_gas + write(*,*) ' multi_gases_init with num_gas=',num_gas + write(*,*) ' multi_gases_init with vir =',vir + write(*,*) ' multi_gases_init with vicp=',vicp + write(*,*) ' multi_gases_init with vicv=',vicv + endif + + return + end subroutine multi_gases_init +! ---------------------------------------------------------------- + +! ---------------------------------------------------------------- + subroutine multi_gases_finalize() + + if(allocated(vir )) deallocate(vir ) + if(allocated(vicv)) deallocate(vicp) + if(allocated(vicp)) deallocate(vicv) + + return + end subroutine multi_gases_finalize +! ---------------------------------------------------------------- + +! -------------------------------------------------------- + pure real function virq(q) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas 1+zvir/(1-qc) +!-------------------------------------------- + real, intent(in) :: q(num_gas) +! Local: + integer :: n + + virq = vir(sphum)*q(sphum) + do n=ind_gas,num_gas + virq = virq+vir(n)*q(sphum+n-1) + end do + virq = vir(0)+virq/(1.0-sum(q(sphump1:sphum+num_wat-1))) + + return + end function +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function virq_nodq(q) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas 1+zvir without dividing by 1-qv or 1-qv-qc +!-------------------------------------------- + real, intent(in) :: q(num_gas) +! Local: + integer :: n + + virq_nodq = vir(0)+vir(sphum)*q(sphum) + do n=ind_gas,num_gas + virq_nodq = virq_nodq+vir(n)*q(sphum+n-1) + end do + + return + end function +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function virq_max(q, qmin) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas 1+zvir using max(qmin,q(sphum)) +!-------------------------------------------- + real, intent(in) :: q(num_gas) + real, intent(in) :: qmin +! Local: + integer :: n + + virq_max = vir(sphum)*max(qmin,q(sphum)) + do n=ind_gas,num_gas + virq_max = virq_max+vir(n)*q(sphum+n-1) + end do + virq_max = vir(0)+virq_max/(1.0-sum(q(sphump1:sphum+num_wat-1))) + + + return + end function +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function virq_qpz(q, qpz) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas 1+zvir/(1.-qpz): qpz in place of qv+qc from q +!-------------------------------------------- + real, intent(in) :: q(num_gas) + real, intent(in) :: qpz +! Local: + integer :: n + + virq_qpz = vir(sphum)*q(sphum) + do n=ind_gas,num_gas + virq_qpz = virq_qpz+vir(n)*q(sphum+n-1) + end do + virq_qpz = vir(0)+virq_qpz/(1.0-qpz) + + + return + end function +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function virqd(q) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas 1+zvir/(1-(qv+qc)) (dry) +!-------------------------------------------- + real, intent(in) :: q(num_gas) +! Local: + integer :: n + + virqd = 0.0 + do n=ind_gas,num_gas + virqd = virqd+vir(n)*q(sphum+n-1) + end do + virqd = vir(0)+virqd/(1.0-sum(q(sphum:sphum+num_wat-1))) + + return + end function +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function vicpqd(q) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas cp (dry) +!-------------------------------------------- + real, intent(in) :: q(num_gas) +! Local: + integer :: n + + vicpqd = 0.0 + do n=ind_gas,num_gas + vicpqd = vicpqd+vicp(n)*q(sphum+n-1) + end do + vicpqd = vicp(0)+vicpqd/(1.0-sum(q(sphum:sphum+num_wat-1))) + + return + end function +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function vicpqd_qpz(q, qpz) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas cp (dry) with qpz in place of qv+qc from q +!-------------------------------------------- + real, intent(in) :: q(num_gas) + real, intent(in) :: qpz +! Local: + integer :: n + + vicpqd_qpz = 0.0 + do n=ind_gas,num_gas + vicpqd_qpz = vicpqd_qpz+vicp(n)*q(sphum+n-1) + end do + vicpqd_qpz = vicp(0)+vicpqd_qpz/(1.0-qpz) + + return + end function +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function vicvqd(q) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas cv (dry) +!-------------------------------------------- + real, intent(in) :: q(num_gas) +! Local: + integer :: n + + vicvqd = 0.0 + do n=ind_gas,num_gas + vicvqd = vicvqd+vicv(n)*q(sphum+n-1) + end do + vicvqd = vicv(0)+vicvqd/(1.0-sum(q(sphum:sphum+num_wat-1))) + + return + end function +!-------------------------------------------- + +! -------------------------------------------------------- + pure real function vicvqd_qpz(q,qpz) +!-------------------------------------------- +! !OUTPUT PARAMETERS +! Ouput: variable gas cv (dry) with qpz in place of qv+qc from q +!-------------------------------------------- + real, intent(in) :: q(num_gas) + real, intent(in) :: qpz +! Local: + integer :: n + + vicvqd_qpz = 0.0 + do n=ind_gas,num_gas + vicvqd_qpz = vicvqd_qpz+vicv(n)*q(sphum+n-1) + end do + vicvqd_qpz = vicv(0)+vicvqd_qpz/(1.0-qpz) + + return + end function +!-------------------------------------------- + +end module ccpp_multi_gases_mod