From 00cb3f6991df8816517b8e968e2fa461d17bafe3 Mon Sep 17 00:00:00 2001 From: Aaron Donahue Date: Wed, 29 May 2024 16:15:36 -0700 Subject: [PATCH] Remove unused interface for ZM deep convection The ZM interface we currently have is out-of-date and no longer used. If we decide to start using ZM again it is likely that we need to rewrite the interface to reflect the current paradigm. As a result, this commit removes the old interface to avoid confusion. --- components/eamxx/src/physics/CMakeLists.txt | 3 - .../eamxx/src/physics/zm/CMakeLists.txt | 28 - .../physics/zm/eamxx_zm_process_interface.cpp | 249 - .../physics/zm/eamxx_zm_process_interface.hpp | 71 - .../src/physics/zm/scream_zm_interface.F90 | 198 - .../src/physics/zm/scream_zm_interface.hpp | 39 - components/eamxx/src/physics/zm/zm_conv.F90 | 4348 ----------------- 7 files changed, 4936 deletions(-) delete mode 100644 components/eamxx/src/physics/zm/CMakeLists.txt delete mode 100644 components/eamxx/src/physics/zm/eamxx_zm_process_interface.cpp delete mode 100644 components/eamxx/src/physics/zm/eamxx_zm_process_interface.hpp delete mode 100644 components/eamxx/src/physics/zm/scream_zm_interface.F90 delete mode 100644 components/eamxx/src/physics/zm/scream_zm_interface.hpp delete mode 100644 components/eamxx/src/physics/zm/zm_conv.F90 diff --git a/components/eamxx/src/physics/CMakeLists.txt b/components/eamxx/src/physics/CMakeLists.txt index 475fc588cd4..006426d0144 100644 --- a/components/eamxx/src/physics/CMakeLists.txt +++ b/components/eamxx/src/physics/CMakeLists.txt @@ -10,9 +10,6 @@ else() message(STATUS "WARNING: RRTMGP and COSP only supported for double precision builds; skipping") endif() add_subdirectory(shoc) -if (SCREAM_TEST_LEVEL GREATER_EQUAL SCREAM_TEST_LEVEL_EXPERIMENTAL) - add_subdirectory(zm) -endif() add_subdirectory(cld_fraction) if (SCREAM_ENABLE_ML_CORRECTION) add_subdirectory(ml_correction) diff --git a/components/eamxx/src/physics/zm/CMakeLists.txt b/components/eamxx/src/physics/zm/CMakeLists.txt deleted file mode 100644 index 9867681599f..00000000000 --- a/components/eamxx/src/physics/zm/CMakeLists.txt +++ /dev/null @@ -1,28 +0,0 @@ -set(ZM_SRCS - ${SCREAM_BASE_DIR}/../eam/src/physics/cam/physics_utils.F90 - ${SCREAM_BASE_DIR}/../eam/src/physics/cam/scream_abortutils.F90 - zm_conv.F90 - eamxx_zm_process_interface.cpp - scream_zm_interface.F90 -) - -set(ZM_HEADERS - zm.hpp - eamxx_zm_process_interface.hpp - scream_zm_interface.hpp -) - -# Add ETI source files if not on CUDA/HIP -if (NOT EAMXX_ENABLE_GPU) - list(APPEND ZM_SRCS - ) -endif() - -add_library(zm ${ZM_SRCS}) -set_target_properties(zm PROPERTIES - Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/modules -) -target_include_directories(zm PUBLIC - ${CMAKE_CURRENT_BINARY_DIR}/modules -) -target_link_libraries(zm physics_share scream_share) diff --git a/components/eamxx/src/physics/zm/eamxx_zm_process_interface.cpp b/components/eamxx/src/physics/zm/eamxx_zm_process_interface.cpp deleted file mode 100644 index f36aac6133d..00000000000 --- a/components/eamxx/src/physics/zm/eamxx_zm_process_interface.cpp +++ /dev/null @@ -1,249 +0,0 @@ -#include "physics/zm/scream_zm_interface.hpp" -#include "physics/zm/eamxx_zm_process_interface.hpp" - -#include "ekat/ekat_assert.hpp" - -namespace { -// A helper struct and fcn; -struct GridOpts { - GridOpts(const std::string& n, bool out, const ekat::units::Units& u, int idx) - : name(n), isOut(out), unit(u), field_idx(idx) {} - std::string name; - bool isOut; - const ekat::units::Units& unit; - int field_idx; -}; - -void set_grid_opts(std::map& opt_map); -} - -namespace scream -{ - -ZMDeepConvection::ZMDeepConvection (const ekat::Comm& comm,const ekat::ParameterList& params ) - : AtmosphereProcess(comm, params) -{ - // Nothing to do here -} - -void ZMDeepConvection::set_grids(const std::shared_ptr grids_manager) -{ - using namespace std; - using namespace ekat::units; - - constexpr int NVL = 72; /* TODO THIS NEEDS TO BE CHANGED TO A CONFIGURABLE */ - constexpr int QSZ = 35; /* TODO THIS NEEDS TO BE CHANGED TO A CONFIGURABLE */ - auto grid = grids_manager->get_grid("Physics"); - const int num_dofs = grid->get_num_local_dofs(); - const int nc = num_dofs; - - using namespace ShortFieldTagsNames; - - FieldLayout scalar3d_layout_mid { {COL,LEV}, {nc,NVL} }; // Note that C++ and Fortran read array dimensions in reverse - FieldLayout scalar3d_layout_int { {COL,ILEV}, {nc,NVL+1} }; // Note that C++ and Fortran read array dimensions in reverse - FieldLayout vector3d_layout_mid{ {COL,CMP,LEV}, {nc,QSZ,NVL} }; - FieldLayout tracers_layout { {COL,CMP,LEV}, {nc,QSZ,NVL} }; - FieldLayout scalar2d_layout{ {COL}, {nc} }; - - std::vector layout_opts = {scalar3d_layout_mid, scalar3d_layout_int, - vector3d_layout_mid, tracers_layout, scalar2d_layout}; - - std::map opt_map; - set_grid_opts(opt_map); - - for ( auto i = opt_map.begin(); i != opt_map.end(); ++i) { - add_required_field((i->second).name, layout_opts[((i->second).field_idx)], kg/kg, grid->name()); - if ( (i->second).isOut == true ) { - add_computed_field((i->second).name, layout_opts[((i->second).field_idx)], kg/kg, grid->name()); - } - } - -} - -// ========================================================================================= -void ZMDeepConvection::initialize_impl (const RunType /* run_type */) -{ - zm_init_f90 (*m_raw_ptrs_in["limcnv_in"], m_raw_ptrs_in["no_deep_pbl_in"]); -} -// ========================================================================================= -void ZMDeepConvection::run_impl (const double dt) -{ - std::vector in; - std::vector out; - - // Copy inputs to host. Copy also outputs, cause we might "update" them, rather than overwrite them. - for (auto& it : m_zm_fields_in) { - it.second.sync_to_host(); - } - for (auto& it : m_zm_fields_out) { - it.second.sync_to_host(); - } - - Real** temp = &m_raw_ptrs_out["fracis"]; - Real*** fracis = &temp; - - zm_main_f90(*m_raw_ptrs_out["lchnk"], *m_raw_ptrs_out["ncol"], m_raw_ptrs_out["t"], - m_raw_ptrs_out["qh"], m_raw_ptrs_out["prec"], m_raw_ptrs_out["jctop"], - m_raw_ptrs_out["jcbot"], m_raw_ptrs_out["pblh"], m_raw_ptrs_out["zm"], - m_raw_ptrs_out["geos"], m_raw_ptrs_out["zi"], m_raw_ptrs_out["qtnd"], - m_raw_ptrs_out["heat"], m_raw_ptrs_out["pap"], m_raw_ptrs_out["paph"], - m_raw_ptrs_out["dpp"], *m_raw_ptrs_out["delt"], m_raw_ptrs_out["mcon"], - m_raw_ptrs_out["cme"], m_raw_ptrs_out["cape"], m_raw_ptrs_out["tpert"], - m_raw_ptrs_out["dlf"], m_raw_ptrs_out["plfx"], m_raw_ptrs_out["zdu"], - m_raw_ptrs_out["rprd"], m_raw_ptrs_out["mu"], m_raw_ptrs_out["md"], - m_raw_ptrs_out["du"], m_raw_ptrs_out["eu"], m_raw_ptrs_out["ed"], - m_raw_ptrs_out["dp"], m_raw_ptrs_out["dsubcld"], m_raw_ptrs_out["jt"], - m_raw_ptrs_out["maxg"], m_raw_ptrs_out["ideep"], *m_raw_ptrs_out["lengath"], - m_raw_ptrs_out["ql"], m_raw_ptrs_out["rliq"], m_raw_ptrs_out["landfrac"], - m_raw_ptrs_out["hu_nm1"], m_raw_ptrs_out["cnv_nm1"], m_raw_ptrs_out["tm1"], - m_raw_ptrs_out["qm1"], &m_raw_ptrs_out["t_star"], &m_raw_ptrs_out["q_star"], - m_raw_ptrs_out["dcape"], m_raw_ptrs_out["qv"], &m_raw_ptrs_out["tend_s"], - &m_raw_ptrs_out["tend_q"], &m_raw_ptrs_out["cld"], m_raw_ptrs_out["snow"], - m_raw_ptrs_out["ntprprd"], m_raw_ptrs_out["ntsnprd"], - &m_raw_ptrs_out["flxprec"], &m_raw_ptrs_out["flxsnow"], - *m_raw_ptrs_out["ztodt"], m_raw_ptrs_out["pguall"], m_raw_ptrs_out["pgdall"], - m_raw_ptrs_out["icwu"], *m_raw_ptrs_out["ncnst"], fracis); - auto ts = timestamp(); - ts += dt; - for (auto& it : m_zm_fields_out) { - it.second.get_header().get_tracking().update_time_stamp(ts); - } -} -// ========================================================================================= -void ZMDeepConvection::finalize_impl() -{ - zm_finalize_f90 (); -} -// ========================================================================================= -void ZMDeepConvection:: -register_fields (const std::map>>& field_mgrs) const { - auto& field_mgr = *field_mgrs.at("Physics"); - for (const auto& fid : get_required_fields()) { - field_mgr.register_field(fid); - } - for (const auto& fid : get_computed_fields()) { - field_mgr.register_field(fid); - } - } - -void ZMDeepConvection::set_required_field_impl (const Field& f) { - // @Meredith: Diff between add_me_as_a_customer and get_tracking().add_customer? - - // Store a copy of the field. We need this in order to do some tracking checks - // at the beginning of the run call. Other than that, there would be really - // no need to store a scream field here; we could simply set the view ptr - // in the Homme's view, and be done with it. - const auto& name = f.get_header().get_identifier().name(); - m_zm_fields_in.emplace(name,f); - m_zm_host_views_in[name] = f.get_view(); - m_raw_ptrs_in[name] = m_zm_host_views_in[name].data(); - - // Add myself as customer to the field - add_me_as_customer(f); - -} - -void ZMDeepConvection::set_computed_field_impl (const Field< Real>& f) { - // Store a copy of the field. We need this in order to do some tracking updates - // at the end of the run call. Other than that, there would be really - // no need to store a scream field here; we could simply set the view ptr - // in the Homme's view, and be done with it. - const auto& name = f.get_header().get_identifier().name(); - m_zm_fields_out.emplace(name,f); - m_zm_host_views_out[name] = f.get_view(); - m_raw_ptrs_out[name] = m_zm_host_views_out[name].data(); - - // Add myself as provider for the field - add_me_as_provider(f); -} - -} // namespace scream - -namespace { -//Initializes struct GridOpts fields -void set_grid_opts_helper(std::map& opt_map, - const std::string& n, bool out, int field_idx, - const ekat::units::Units& unit = ekat::units::Units::nondimensional()) -{ - opt_map.emplace(std::piecewise_construct, - std::forward_as_tuple(n), - std::forward_as_tuple(n,out,unit,field_idx)); -} - -void set_grid_opts(std::map& opt_map){ - //Layout options are set as an int - //to be passed into the GridOpts struct - const int SCALAR_3D_MID = 0; - const int SCALAR_3D_INT = 1; - const int VECTOR_3D_MID = 2; - const int LINEAR = 4; - - const auto Pa = ekat::units::Pa; - - // Sets the value of the grid opt struct categories - set_grid_opts_helper(opt_map, "limcnv_in", true, 0); - set_grid_opts_helper(opt_map, "no_deep_pbl_in", true, VECTOR_3D_MID); - set_grid_opts_helper(opt_map, "lchnk", true, SCALAR_3D_MID, Pa); //temperature(K) - set_grid_opts_helper(opt_map, "ncol", true, SCALAR_3D_MID, Pa); //temperature(K) - set_grid_opts_helper(opt_map, "t", true, SCALAR_3D_MID, Pa); //temperature(K) - set_grid_opts_helper(opt_map, "qh", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "prec", true, LINEAR); - set_grid_opts_helper(opt_map, "jctop", true, LINEAR); - set_grid_opts_helper(opt_map, "jcbot", true, LINEAR); - set_grid_opts_helper(opt_map, "pblh", true, LINEAR); - set_grid_opts_helper(opt_map, "zm", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "geos", true, LINEAR); - set_grid_opts_helper(opt_map, "zi", true, SCALAR_3D_INT); - set_grid_opts_helper(opt_map, "qtnd", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "heat", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "pap", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "paph", true, SCALAR_3D_INT); - set_grid_opts_helper(opt_map, "dpp", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "delt", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "mcon", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "cme", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "cape", true, LINEAR); - set_grid_opts_helper(opt_map, "tpert", true, LINEAR); - set_grid_opts_helper(opt_map, "dlf", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "pflx", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "zdu", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "rprd", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "mu", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "md", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "du", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "eu", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "ed", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "dp", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "dsubcld", true, LINEAR); - set_grid_opts_helper(opt_map, "jt", true, LINEAR); - set_grid_opts_helper(opt_map, "maxg", true, LINEAR); - set_grid_opts_helper(opt_map, "ideep", true, LINEAR); - set_grid_opts_helper(opt_map, "lengath", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "ql", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "rliq", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "landfrac", true, LINEAR); - set_grid_opts_helper(opt_map, "hu_nm1", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "cnv_nm1", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "tm1", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "qm1", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "t_star", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "q_star", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "dcape", true, LINEAR); - set_grid_opts_helper(opt_map, "qv", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "tend_s", true, LINEAR); - set_grid_opts_helper(opt_map, "tend_q", true, LINEAR); - set_grid_opts_helper(opt_map, "cld", true, LINEAR); - set_grid_opts_helper(opt_map, "snow", true, LINEAR); - set_grid_opts_helper(opt_map, "ntprprd", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "ntsnprd", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "flxprec", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "flxsnow", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "ztodt", true, SCALAR_3D_MID); - set_grid_opts_helper(opt_map, "pguall", true, VECTOR_3D_MID); - set_grid_opts_helper(opt_map, "pgdall", true, VECTOR_3D_MID); - set_grid_opts_helper(opt_map, "icwu", true, VECTOR_3D_MID); - set_grid_opts_helper(opt_map, "ncnst", true, VECTOR_3D_MID); - set_grid_opts_helper(opt_map, "fracis", true, VECTOR_3D_MID); -} -} // anonymous namespace - diff --git a/components/eamxx/src/physics/zm/eamxx_zm_process_interface.hpp b/components/eamxx/src/physics/zm/eamxx_zm_process_interface.hpp deleted file mode 100644 index 8b25a4b09ec..00000000000 --- a/components/eamxx/src/physics/zm/eamxx_zm_process_interface.hpp +++ /dev/null @@ -1,71 +0,0 @@ -#ifndef SCREAM_ZM_DEEPCONVECTION_HPP -#define SCREAM_ZM_DEEPCONVECTION_HPP - -#include "share/atm_process/atmosphere_process.hpp" -#include "ekat/ekat_parameter_list.hpp" - -#include - -namespace scream -{ - -/* - * The class responsible to handle the atmosphere deep convection - * - * The AD should store exactly ONE instance of this class stored - * in its list of subcomponents (the AD should make sure of this). - * -*/ - -class ZMDeepConvection : public AtmosphereProcess -{ -public: - using field_type = Field< Real>; - using const_field_type = Field; - - // Constructors - ZMDeepConvection (const ekat::Comm& comm, const ekat::ParameterList& params); - - // The type of subcomponent - AtmosphereProcessType type () const { return AtmosphereProcessType::Physics; } - - // The name of the subcomponent - std::string name () const { return "DeepConvection"; } - - // Set the grid - void set_grids (const std::shared_ptr grids_manager); - - // The three main interfaces for the subcomponent - void initialize_impl (const RunType run_type); - void run_impl (const double dt); - void finalize_impl (); - - // Register all fields in the proper field manager(s). - // Note: field_mgrs[grid_name] is the FM on grid $grid_name - void register_fields (const std::map>>& field_mgrs) const; - -protected: - - std::map m_zm_fields_in; - std::map m_zm_fields_out; - - template - using view_type = field_type::view_type; - - template - using host_view_type = field_type::get_view_type,Host>; - - using host_view_in_type = host_view_type; - using host_view_out_type = host_view_type< field_type::RT>; - - std::map m_zm_host_views_in; - std::map m_zm_host_views_out; - - std::map m_raw_ptrs_in; - std::map m_raw_ptrs_out; - -}; // class ZMDeepConvection - -} // namespace scream - -#endif // SCREAM_ZM_DEEPCONVECTION_HPP diff --git a/components/eamxx/src/physics/zm/scream_zm_interface.F90 b/components/eamxx/src/physics/zm/scream_zm_interface.F90 deleted file mode 100644 index 61840a8e7ec..00000000000 --- a/components/eamxx/src/physics/zm/scream_zm_interface.F90 +++ /dev/null @@ -1,198 +0,0 @@ -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif - -module scream_zm_interface_mod - - use iso_c_binding, only: c_ptr, c_f_pointer, c_int, c_double, c_bool,C_NULL_CHAR, c_float - use physics_utils, only: r8 => rtype, rtype8, itype, btype - use zm_conv, only: zm_convr, zm_conv_evap, zm_convi, convtran, momtran,zmconv_readnl - - implicit none -#include "scream_config.f" -#ifdef SCREAM_DOUBLE_PRECISION -# define c_real c_double - integer,parameter,public :: rtype = c_double ! 8 byte real, compatible with c type double -#else -# define c_real c_float - integer,parameter,public :: rtype = c_float ! 4 byte real, compatible with c type #endif -#endif - - - - public :: zm_init_f90 - public :: zm_main_f90 - public :: zm_finalize_f90 - - real :: test - - integer(kind=c_int) :: pcols = 32 - integer(kind=c_int) :: pver = 73 - integer(kind=c_int) :: pverp = 75 - - ! real(kind=c_real) :: pcols = 32 - ! real(kind=c_real) :: pver = 75 - real(kind=c_real) :: cpair != 1004.64000000000 - real(kind=c_real) :: gravit != 9.80616000000000 - real(kind=c_real) :: latvap != 2501000.00000000 - real(kind=c_real) :: plevp = 0.0 - -! real(kind=c_real) :: pverp = 0 - -contains - - !====================================================================! - - subroutine zm_init_f90 (limcnv_in, no_deep_pbl_in) bind(c) - integer, intent(in) :: limcnv_in - logical, intent(in), optional :: no_deep_pbl_in - call zm_convi(limcnv_in, no_deep_pbl_in) - call zmconv_readnl() - - end subroutine zm_init_f90 - !====================================================================! -subroutine zm_main_f90(lchnk ,ncol , & - t ,qh ,prec ,jctop ,jcbot , & - pblh ,zm ,geos ,zi ,qtnd , & - heat ,pap ,paph ,dpp , & - delt ,mcon ,cme ,cape , & - tpert ,dlf ,pflx ,zdu ,rprd, & - mu ,md ,du ,eu ,ed , & - dp ,dsubcld ,jt ,maxg ,ideep , & - lengath ,ql ,rliq ,landfrac,hu_nm1 , & - cnv_nm1 ,tm1 ,qm1 ,t_star ,q_star , & - dcape ,q ,tend_s ,tend_q ,cld , & - snow ,ntprprd ,ntsnprd ,flxprec ,flxsnow, & - ztodt , pguall , pgdall , icwu , ncnst, fracis) bind(c) - integer :: lchnk - integer lengath - integer, intent(in) :: ncol - real(kind=c_real), intent(inout) :: t(pcols,pver) ! State array kg/kg - real(kind=c_real) :: qh(pcols,pver) - real(kind=c_real), intent(out) :: prec(pcols) - real(kind=c_real), intent(out) :: jctop(pcols) ! o row of top-of-deep-convection indices passed out. - real(kind=c_real), intent(out) :: jcbot(pcols) ! o row of base of cloud indices passed out. - integer jt(pcols) ! wg top level index of deep cumulus convection. - integer maxg(pcols) ! wg gathered values of maxi. - integer ideep(pcols) ! w holds position of gathered points vs longitude index. - integer mx(pcols) -! -! real(r8) landfracg(pcols) ! wg grid slice of landfrac - real(kind=c_real) delt ! length of model time-step in seconds. -! real(r8) :: pcont(pcols), pconb(pcols), freqzm(pcols) -! - real(kind=c_real), intent(in) :: pblh(pcols) - real(kind=c_real), intent(in) :: geos(pcols) - real(kind=c_real), intent(in) :: zi(pcols,pver+1) - real(kind=c_real), intent(in) :: zm(pcols,pver) - real(kind=c_real), intent(in) :: pap(pcols,pver) - real(kind=c_real), intent(in) :: paph(pcols,pver+1) - real(kind=c_real), intent(in) :: dpp(pcols,pver) ! local sigma half-level thickness (i.e. dshj). - real(kind=c_real), intent(in) :: tpert(pcols) - real(kind=c_real), intent(in) :: tm1(pcols,pver) ! grid slice of temperature at mid-layer. - real(kind=c_real), intent(in) :: qm1(pcols,pver) ! grid slice of specific humidity. - real(kind=c_real), intent(in) :: landfrac(pcols) ! RBN Landfrac - - - real(kind=c_real), intent(in) :: t_star(pcols,pver) ! intermediate T between n and n-1 time step - real(kind=c_real), intent(in) :: q_star(pcols,pver) ! intermediate q between n and n-1 time step - - real(kind=c_real), intent(out) :: qtnd(pcols,pver) ! specific humidity tendency (kg/kg/s) - real(kind=c_real), intent(out) :: heat(pcols,pver) ! heating rate (dry static energy tendency, W/kg) - real(kind=c_real), intent(out) :: mcon(pcols,pverp) - real(kind=c_real), intent(out) :: dlf(pcols,pver) ! scattrd version of the detraining cld h2o tend - real(kind=c_real), intent(out) :: pflx(pcols,pverp) ! scattered precip flux at each level - real(kind=c_real), intent(out) :: cme(pcols,pver) - real(kind=c_real), intent(out) :: cape(pcols) ! w convective available potential energy. - real(kind=c_real), intent(out) :: zdu(pcols,pver) - real(kind=c_real), intent(out) :: rprd(pcols,pver) ! rain production rate - real(kind=c_real), intent(out) :: mu(pcols,pver) - real(kind=c_real), intent(out) :: md(pcols,pver) - real(kind=c_real), intent(out) :: du(pcols,pver) ! detrainement rate of updraft - real(kind=c_real), intent(out) :: ed(pcols,pver) ! entrainment rate of downdraft - real(kind=c_real), intent(out) :: eu(pcols,pver) ! entrainment rate of updraft - real(kind=c_real), intent(out) :: dp(pcols,pver) ! wg layer thickness in mbs (between upper/lower interface). - real(kind=c_real), intent(out) :: dsubcld(pcols) ! wg layer thickness in mbs between lcl and maxi. - real(kind=c_real), intent(out) :: ql(pcols,pver) - real(kind=c_real), intent(out) :: rliq(pcols) ! reserved liquid (not yet in cldliq) for energy integrals - real(kind=c_real), intent(out) :: dcape(pcols) ! output dynamical CAPE -! - real(kind=c_real), intent(inout) :: hu_nm1 (pcols,pver) - real(kind=c_real), intent(inout) :: cnv_nm1 (pcols,pver) -! -! - real(kind=c_real) q(pcols,pver) -! -!! real(r8), pointer, dimension(:,:) :: rprd ! rain production rate -!!Used for convtran exclusively - real(kind=c_real), intent(in) :: fracis(pcols,pver,ncnst) - real(kind=c_real) :: fake_dpdry(pcols,pver) - real(kind=c_real) :: fake_dqdt(pcols,pver,ncnst) ! Tracer tendency array -! - integer :: il1g - integer :: nstep -!! -! - real(kind=c_real) :: tend_s_snwprd (pcols,pver) ! Heating rate of snow production - real(kind=c_real) :: tend_s_snwevmlt(pcols,pver) ! Heating rate of evap/melting of snow - real(kind=c_real),intent(inout), dimension(pcols,pver) :: tend_s ! heating rate (J/kg/s) - real(kind=c_real),intent(inout), dimension(pcols,pver) :: tend_q ! heating rate (J/kg/s) - real(kind=c_real), intent(inout), dimension(pcols,pver) :: cld - real(kind=c_real), intent(inout), dimension(pcols) :: snow ! snow from ZM convection - real(kind=c_real) :: ntprprd(pcols,pver) ! evap outfld: net precip production in layer - real(kind=c_real) :: ntsnprd(pcols,pver) ! evap outfld: net snow production in layer - real(kind=c_real), intent(out), dimension(pcols,pver) :: flxprec ! Convective-scale flux of precip at interfaces (kg/m2/s) - real(kind=c_real), intent(out), dimension(pcols,pver) :: flxsnow ! Convective-scale flux of snow at interfaces (kg/m2/s) - real(kind=c_real), intent(in) :: ztodt ! 2 delta t (model time increment) - real(kind=c_real) :: pguall(pcols, pver, 2) - real(kind=c_real) :: pgdall(pcols, pver, 2) - real(kind=c_real) :: icwu(pcols,pver, 2) - real(kind=c_real) :: icwd(pcols,pver, 2) - real(kind=c_real) :: seten(pcols, pver) - integer, intent(in) :: ncnst -! - logical :: domomtran(ncnst) - logical :: doconvtran(ncnst) - fake_dpdry(:,:) = 0._r8 - il1g = 1 - - call zm_convr(lchnk ,ncol , & - t ,qh ,prec ,jctop ,jcbot , & - pblh ,zm ,geos ,zi ,qtnd , & - heat ,pap ,paph ,dpp , & - delt ,mcon ,cme ,cape , & - tpert ,dlf ,pflx ,zdu ,rprd , & - mu ,md ,du ,eu ,ed , & - dp ,dsubcld ,jt ,maxg ,ideep , & - lengath ,ql ,rliq ,landfrac,hu_nm1 , & - cnv_nm1 ,tm1 ,qm1 ,t_star ,q_star, dcape) - - call zm_conv_evap(ncol ,lchnk, & - t ,pap ,dpp ,q , & - tend_s, tend_s_snwprd ,tend_s_snwevmlt , & - tend_q ,rprd, cld ,ztodt , & - prec, snow, ntprprd, ntsnprd, flxprec, flxsnow ) - - - - call momtran(lchnk, ncol, & - domomtran,q ,ncnst ,mu ,md , & - du ,eu ,ed ,dp ,dsubcld , & - jt ,mx ,ideep ,il1g ,lengath , & - nstep ,fake_dqdt ,pguall ,pgdall, icwu, icwd, ztodt, seten ) - - call convtran(lchnk , & - doconvtran,q ,ncnst ,mu ,md , & - du ,eu ,ed ,dp ,dsubcld , & - jt ,mx ,ideep ,il1g , lengath , & - 1 ,fracis ,fake_dqdt, fake_dpdry ) - - end subroutine zm_main_f90 - !====================================================================! - subroutine zm_finalize_f90 () bind(c) - - - end subroutine zm_finalize_f90 - !====================================================================! - -end module scream_zm_interface_mod diff --git a/components/eamxx/src/physics/zm/scream_zm_interface.hpp b/components/eamxx/src/physics/zm/scream_zm_interface.hpp deleted file mode 100644 index ffeb15707cd..00000000000 --- a/components/eamxx/src/physics/zm/scream_zm_interface.hpp +++ /dev/null @@ -1,39 +0,0 @@ -#ifndef SCREAM_ZM_INTERFACE_HPP -#define SCREAM_ZM_INTERFACE_HPP - -#include "ekat/ekat_assert.hpp" -//#include "ekat/util/scream_utils.hpp" -#include "ekat/util/ekat_file_utils.hpp" - -#include "share/scream_types.hpp" - -// Put everything into a scream namespace - - -namespace scream { - -extern "C" -{ - -// Fortran routines to be called from C -void zm_init_f90 (const Real& limcnv_in, const bool& no_deep_pbl_in); -void zm_main_f90(const Real& lchnk, const Real& ncol, Real* t, Real* qh, Real* prec, - Real* jctop, Real* jcbot, Real* pblh, Real *zm, Real* geos, Real* zi, - Real* qtnd, Real* heat, Real* pap, Real* paph, Real* dpp, const Real &delt, - Real* mcon, Real* cme, Real* cape, Real* tpert, Real* dlf, Real* plfx, - Real* zdu, Real* rprd, Real* mu, Real* md, Real* du, Real* eu, - Real* ed, Real* dp, Real* dsubcld, Real* jt, Real* maxg, Real* ideep, - const Real& lengath, Real* ql, Real* rliq, Real* landfrac, Real* hu_nm1, - Real* cnv_nm1, Real* tm1, Real* qm1, Real** t_star, Real** q_star, - Real* dcape, Real* q, Real** tend_s, - Real** tend_q, Real** cld, Real* snow, Real* ntprprd, Real* ntsnprd, - Real** flxprec, Real** flxsnow, const Real& ztodt, Real* pguall, Real* pgdall, - Real* icwu, const Real& ncnst, Real*** fracis - ); -void zm_finalize_f90 (); - -} // extern "C" - -} // namespace scream - -#endif // SCREAM_SHOC_INTERFACE_HPP diff --git a/components/eamxx/src/physics/zm/zm_conv.F90 b/components/eamxx/src/physics/zm/zm_conv.F90 deleted file mode 100644 index 507d6f43d95..00000000000 --- a/components/eamxx/src/physics/zm/zm_conv.F90 +++ /dev/null @@ -1,4348 +0,0 @@ - -module zm_conv - -!--------------------------------------------------------------------------------- -! Purpose: -! -! Interface from Zhang-McFarlane convection scheme, includes evaporation of convective -! precip from the ZM scheme -! -! Apr 2006: RBN: Code added to perform a dilute ascent for closure of the CM mass flux -! based on an entraining plume a la Raymond and Blythe (1992) -! -! Author: Byron Boville, from code in tphysbc -! -!--------------------------------------------------------------------------------- - use iso_c_binding, only: c_ptr, c_f_pointer, c_int, c_double, c_bool,C_NULL_CHAR, c_float - - use physics_utils, only: r8 => rtype, rtype8, itype, btype -! use shr_kind_mod, only: r8 => shr_kind_r8 -! use spmd_utils, only: masterproc - ! use ppgrid, only: pcols, pver, pverp -! use cloud_fraction, only: cldfrc_fice - ! use physconst, only: cpair, epsilo, gravit, latice, latvap, tmelt, rair, & - ! cpwv, cpliq, rh2o -! use cam_abortutils, only: endrun -! use cam_logfile, only: iulog - use scream_abortutils, only: endscreamrun - implicit none - -#include "scream_config.f" -#ifdef SCREAM_DOUBLE_PRECISION -# define c_real c_double -#else -# define c_real c_float -#endif - - save - private ! Make default type private to the module -! -! PUBLIC: interfaces -! - public zmconv_readnl ! read zmconv_nl namelist - public zm_convi ! ZM schemea - public zm_convr ! ZM schemea - public zm_conv_evap ! evaporation of precip from ZM schemea -! AaronDonahue - TODO Question: Are convtran and momtran going to be used? Do -! not appear to be used by zm_convr - public convtran ! convective transport - public momtran ! convective momentum transport - public trigmem ! true if convective memory - public trigdcape_ull ! true if to use dcape-ULL trigger - public is_first_step -! -! Private data -! - - integer(kind=c_int) :: pcols = 32 - integer(kind=c_int) :: pver = 72 - integer(kind=c_int) :: pverp = 73 - - logical :: masterproc = .false. - real(kind=c_real) :: cpair = 1004.64000000000 - real(kind=c_real) :: rh2o = 461.504639820160 - real(kind=c_real) :: gravit = 9.80616000000000 - real(kind=c_real) :: latvap = 2501000.00000000 - real(kind=c_real) :: latice = 333700.000000000 - real(kind=c_real) :: tmelt = 273.150000000000 - real(kind=c_real) :: rair = 287.042311365049 - real(kind=c_real) :: cpliq = 4188.00000000000 - real(kind=c_real) :: cpwv = 1.810e3 - real(kind=c_real) :: epsilo = 18.0160000000000/28.9660000000000 !mwh20/mwdry - - real(r8), parameter :: unset_r8 = huge(1.0_r8) - integer , parameter :: unset_int = huge(1) - real(r8) :: zmconv_alfa = unset_r8 - - real(r8) rl ! wg latent heat of vaporization. - real(r8) cpres ! specific heat at constant pressure in j/kg-degk. - real(r8), parameter :: capelmt = 70._r8 ! threshold value for cape for deep convection. -!songxl 2014-05-20------------------ - -!DCAPE-ULL - real(r8), parameter :: trigdcapelmt = 0._r8 ! threshold value of dcape for deep convection - logical :: trigdcape_ull = .false. !true to use DCAPE trigger and -ULL - integer, allocatable :: dcapemx(:) ! save maxi from 1st call for CAPE calculation and used in 2nd call when DCAPE-ULL active -! May need to change to use local variable ! as passed via dummy argument. For now, making it threadprivate as follows, -!$omp threadprivate (dcapemx) - - real(r8) :: ke ! Tunable evaporation efficiency set from namelist input zmconv_ke - real(r8) :: c0_lnd ! set from namelist input zmconv_c0_lnd - real(r8) :: c0_ocn ! set from namelist input zmconv_c0_ocn - real(r8) :: dmpdz = unset_r8 ! Parcel fractional mass entrainment rate (/m) - real(r8) :: alfa_scalar ! maximum downdraft mass flux fraction - real(r8) :: tiedke_add = unset_r8 - logical :: trigmem ! set from namelist input zmconv_trigmem - integer :: num_cin = unset_int !number of negative buoyancy regions that are allowed before the conv. top and CAPE calc are completed - integer :: mx_bot_lyr_adj = unset_int !bottom layer adjustment for setting "launching" level(mx) (to be at maximum moist static energy). - real(r8) tau ! convective time scale - real(r8),parameter :: c1 = 6.112_r8 - real(r8),parameter :: c2 = 17.67_r8 - real(r8),parameter :: c3 = 243.5_r8 - real(r8) :: tfreez - real(r8) :: eps1 - - - logical :: no_deep_pbl ! default = .false. - ! no_deep_pbl = .true. eliminates deep convection entirely within PBL - - -!moved from moistconvection.F90 - real(r8) :: rgrav ! reciprocal of grav - real(r8) :: rgas ! gas constant for dry air - real(r8) :: grav ! = gravit - real(r8) :: cp ! = cpres = cpair - - integer limcnv ! top interface level limit for convection - - real(r8) :: tp_fac = unset_r8 ! PMA tunes tpert - - logical :: is_first_step_local = .true. ! AaronDonahue - TODO, actually check if this is the first step given input from the SCREAM-AD - -contains - -subroutine zmconv_readnl()!nlfile) -! AaronDonahue - we do not need/have namelist capabilities yet. That being -! said, we probably still need the constants that are defined here. Switching -! to hard-coded values, no namelist. -! use namelist_utils, only: find_group_name -! use units, only: getunit, freeunit -! use mpishorthand - -! character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input - - ! Local variables -! integer :: unitn, ierr -! character(len=*), parameter :: subname = 'zmconv_readnl' - -! namelist /zmconv_nl/ zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke, zmconv_tau, & -! zmconv_dmpdz, zmconv_alfa, zmconv_trigmem, zmconv_tiedke_add, & -! zmconv_cape_cin, zmconv_mx_bot_lyr_adj, zmconv_tp_fac, zmconv_trigdcape_ull -! !----------------------------------------------------------------------------- - -! zmconv_tau = 3600._r8 -! if (masterproc) then -! unitn = getunit() -! open( unitn, file=trim(nlfile), status='old' ) -! call find_group_name(unitn, 'zmconv_nl', status=ierr) -! if (ierr == 0) then -! read(unitn, zmconv_nl, iostat=ierr) -! if (ierr /= 0) then -! call endrun(subname // ':: ERROR reading namelist') -! end if -! end if -! close(unitn) -! call freeunit(unitn) - - ! set local variables - c0_lnd = 0.007_r8 ! zmconv_c0_lnd - c0_ocn = 0.007_r8 ! zmconv_c0_ocn - ke = 1.5e-6_r8 ! zmconv_ke - tau = 3600 ! zmconv_tau - trigmem = .false. ! zmconv_trigmem - trigdcape_ull = .false. ! zmconv_trigdcape_ull - tiedke_add = 0.8_r8 ! zmconv_tiedke_add - num_cin = 1 ! zmconv_cape_cin - mx_bot_lyr_adj = 2 ! zmconv_mx_bot_lyr_adj - dmpdz = -0.7e-3_r8 ! zmconv_dmpdz - tp_fac = 0.0_r8 ! zmconv_tp_fac - - if ( zmconv_alfa /= unset_r8 ) then - alfa_scalar = zmconv_alfa - else - alfa_scalar = 0.1_r8 - end if - - if (trigdcape_ull) then - write(*,*)'**** ZMCONV-DCAPE trigger with unrestricted launch level:', trigdcape_ull - endif - -! end if -! -!#ifdef SPMD -! ! Broadcast namelist variables -! call mpibcast(c0_lnd, 1, mpir8, 0, mpicom) -! call mpibcast(c0_ocn, 1, mpir8, 0, mpicom) -! call mpibcast(ke, 1, mpir8, 0, mpicom) -! call mpibcast(tau, 1, mpir8, 0, mpicom) -! call mpibcast(dmpdz, 1, mpir8, 0, mpicom) -! call mpibcast(alfa_scalar, 1, mpir8, 0, mpicom) -! call mpibcast(trigmem, 1, mpilog, 0, mpicom) -! call mpibcast(trigdcape_ull, 1, mpilog, 0, mpicom) -! call mpibcast(tiedke_add, 1, mpir8, 0, mpicom) -! call mpibcast(num_cin, 1, mpiint, 0, mpicom) -! call mpibcast(mx_bot_lyr_adj, 1, mpiint, 0, mpicom) -! call mpibcast(tp_fac, 1, mpir8, 0, mpicom) -!#endif - -end subroutine zmconv_readnl - - -subroutine zm_convi(limcnv_in, no_deep_pbl_in) - - - integer, intent(in) :: limcnv_in ! top interface level limit for convection - logical, intent(in), optional :: no_deep_pbl_in ! no_deep_pbl = .true. eliminates ZM convection entirely within PBL - - ! local variables -! character(len=32) :: hgrid ! horizontal grid specifier ! AaronDonahue - Doesn't appear to be used - - ! Initialization of ZM constants - limcnv = limcnv_in - tfreez = tmelt - eps1 = epsilo - rl = latvap - cpres = cpair - rgrav = 1.0_r8/gravit - rgas = rair - grav = gravit - cp = cpres - - if ( present(no_deep_pbl_in) ) then - no_deep_pbl = no_deep_pbl_in - else - no_deep_pbl = .false. - endif - - ! tau=4800. were used in canadian climate center. however, in echam3 t42, - ! convection is too weak, thus adjusted to 2400. - -! hgrid = get_resolution() ! AaronDonahue - Doesn't appear to be used. - if(trigmem)tau = 3600._r8 - - if ( masterproc ) then - write(*,*) 'tuning parameters zm_convi: tau',tau - write(*,*) 'tuning parameters zm_convi: c0_lnd',c0_lnd, ', c0_ocn', c0_ocn - write(*,*) 'tuning parameters zm_convi: ke',ke - write(*,*) 'tuning parameters zm_convi: dmpdz',dmpdz - write(*,*) 'tuning parameters zm_convi: alfa',alfa_scalar - write(*,*) 'tuning parameters zm_convi: no_deep_pbl',no_deep_pbl - endif - - if (masterproc) write(*,*)'**** ZM: DILUTE Buoyancy Calculation ****' - -end subroutine zm_convi - - - -subroutine zm_convr(lchnk ,ncol , & - t ,qh ,prec ,jctop ,jcbot , & - pblh ,zm ,geos ,zi ,qtnd , & - heat ,pap ,paph ,dpp , & - delt ,mcon ,cme ,cape , & - tpert ,dlf ,pflx ,zdu ,rprd , & - mu ,md ,du ,eu ,ed , & - dp ,dsubcld ,jt ,maxg ,ideep , & - lengath ,ql ,rliq ,landfrac,hu_nm1 , & - cnv_nm1 ,tm1 ,qm1 ,t_star ,q_star, dcape) -!----------------------------------------------------------------------- -! -! Purpose: -! Main driver for zhang-mcfarlane convection scheme -! -! Method: -! performs deep convective adjustment based on mass-flux closure -! algorithm. -! -! Author:guang jun zhang, m.lazare, n.mcfarlane. CAM Contact: P. Rasch -! -! This is contributed code not fully standardized by the CAM core group. -! All variables have been typed, where most are identified in comments -! The current procedure will be reimplemented in a subsequent version -! of the CAM where it will include a more straightforward formulation -! and will make use of the standard CAM nomenclature -! -!----------------------------------------------------------------------- -! use time_manager, only: is_first_step, is_first_restart_step !songxl 2014-05-20 -! -! ************************ index of variables ********************** -! -! wg * alpha array of vertical differencing used (=1. for upstream). -! w * cape convective available potential energy. -! wg * capeg gathered convective available potential energy. -! c * capelmt threshold value for cape for deep convection. -! ic * cpres specific heat at constant pressure in j/kg-degk. -! i * dpp -! ic * delt length of model time-step in seconds. -! wg * dp layer thickness in mbs (between upper/lower interface). -! wg * dqdt mixing ratio tendency at gathered points. -! wg * dsdt dry static energy ("temp") tendency at gathered points. -! wg * dudt u-wind tendency at gathered points. -! wg * dvdt v-wind tendency at gathered points. -! wg * dsubcld layer thickness in mbs between lcl and maxi. -! ic * grav acceleration due to gravity in m/sec2. -! wg * du detrainment in updraft. specified in mid-layer -! wg * ed entrainment in downdraft. -! wg * eu entrainment in updraft. -! wg * hmn moist static energy. -! wg * hsat saturated moist static energy. -! w * ideep holds position of gathered points vs longitude index. -! ic * pver number of model levels. -! wg * j0 detrainment initiation level index. -! wg * jd downdraft initiation level index. -! ic * jlatpr gaussian latitude index for printing grids (if needed). -! wg * jt top level index of deep cumulus convection. -! w * lcl base level index of deep cumulus convection. -! wg * lclg gathered values of lcl. -! w * lel index of highest theoretical convective plume. -! wg * lelg gathered values of lel. -! w * lon index of onset level for deep convection. -! w * maxi index of level with largest moist static energy. -! wg * maxg gathered values of maxi. -! wg * mb cloud base mass flux. -! wg * mc net upward (scaled by mb) cloud mass flux. -! wg * md downward cloud mass flux (positive up). -! wg * mu upward cloud mass flux (positive up). specified -! at interface -! ic * msg number of missing moisture levels at the top of model. -! w * p grid slice of ambient mid-layer pressure in mbs. -! i * pblt row of pbl top indices. -! w * pcpdh scaled surface pressure. -! w * pf grid slice of ambient interface pressure in mbs. -! wg * pg grid slice of gathered values of p. -! w * q grid slice of mixing ratio. -! wg * qd grid slice of mixing ratio in downdraft. -! wg * qg grid slice of gathered values of q. -! i/o * qh grid slice of specific humidity. -! w * qh0 grid slice of initial specific humidity. -! wg * qhat grid slice of upper interface mixing ratio. -! wg * ql grid slice of cloud liquid water. -! wg * qs grid slice of saturation mixing ratio. -! w * qstp grid slice of parcel temp. saturation mixing ratio. -! wg * qstpg grid slice of gathered values of qstp. -! wg * qu grid slice of mixing ratio in updraft. -! ic * rgas dry air gas constant. -! wg * rl latent heat of vaporization. -! w * s grid slice of scaled dry static energy (t+gz/cp). -! wg * sd grid slice of dry static energy in downdraft. -! wg * sg grid slice of gathered values of s. -! wg * shat grid slice of upper interface dry static energy. -! wg * su grid slice of dry static energy in updraft. -! i/o * t -! o * jctop row of top-of-deep-convection indices passed out. -! O * jcbot row of base of cloud indices passed out. -! wg * tg grid slice of gathered values of t. -! w * tl row of parcel temperature at lcl. -! wg * tlg grid slice of gathered values of tl. -! w * tp grid slice of parcel temperatures. -! wg * tpg grid slice of gathered values of tp. -! i/o * u grid slice of u-wind (real). -! wg * ug grid slice of gathered values of u. -! i/o * utg grid slice of u-wind tendency (real). -! i/o * v grid slice of v-wind (real). -! w * va work array re-used by called subroutines. -! wg * vg grid slice of gathered values of v. -! i/o * vtg grid slice of v-wind tendency (real). -! i * w grid slice of diagnosed large-scale vertical velocity. -! w * z grid slice of ambient mid-layer height in metres. -! w * zf grid slice of ambient interface height in metres. -! wg * zfg grid slice of gathered values of zf. -! wg * zg grid slice of gathered values of z. -! -!----------------------------------------------------------------------- -! -! multi-level i/o fields: -! i => input arrays. -! i/o => input/output arrays. -! w => work arrays. -! wg => work arrays operating only on gathered points. -! ic => input data constants. -! c => data constants pertaining to subroutine itself. -! -! input arguments -! - integer, intent(in) :: lchnk ! chunk identifier - integer, intent(in) :: ncol ! number of atmospheric columns - - real(r8), intent(in) :: t(pcols,pver) ! grid slice of temperature at mid-layer. - real(r8), intent(in) :: qh(pcols,pver) ! grid slice of specific humidity. - real(r8), intent(in) :: pap(pcols,pver) - real(r8), intent(in) :: paph(pcols,pver+1) - real(r8), intent(in) :: dpp(pcols,pver) ! local sigma half-level thickness (i.e. dshj). - real(r8), intent(in) :: zm(pcols,pver) - real(r8), intent(in) :: geos(pcols) - real(r8), intent(in) :: zi(pcols,pver+1) - real(r8), intent(in) :: pblh(pcols) - real(r8), intent(in) :: tpert(pcols) - real(r8), intent(in) :: landfrac(pcols) ! RBN Landfrac -!songxl 2014-05-20------------------ - -!DCAPE-ULL - real(r8), intent(in), dimension(pcols,pver) :: t_star ! intermediate T between n and n-1 time step - real(r8), intent(in), dimension(pcols,pver) :: q_star ! intermediate q between n and n-1 time step - - -! -! output arguments -! - real(r8), intent(out) :: qtnd(pcols,pver) ! specific humidity tendency (kg/kg/s) - real(r8), intent(out) :: heat(pcols,pver) ! heating rate (dry static energy tendency, W/kg) - real(r8), intent(out) :: mcon(pcols,pverp) - real(r8), intent(out) :: dlf(pcols,pver) ! scattrd version of the detraining cld h2o tend - real(r8), intent(out) :: pflx(pcols,pverp) ! scattered precip flux at each level - real(r8), intent(out) :: cme(pcols,pver) - real(r8), intent(out) :: cape(pcols) ! w convective available potential energy. - real(r8), intent(out) :: zdu(pcols,pver) - real(r8), intent(out) :: rprd(pcols,pver) ! rain production rate -! move these vars from local storage to output so that convective -! transports can be done in outside of conv_cam. - real(r8), intent(out) :: mu(pcols,pver) - real(r8), intent(out) :: eu(pcols,pver) - real(r8), intent(out) :: du(pcols,pver) - real(r8), intent(out) :: md(pcols,pver) - real(r8), intent(out) :: ed(pcols,pver) - real(r8), intent(out) :: dp(pcols,pver) ! wg layer thickness in mbs (between upper/lower interface). - real(r8), intent(out) :: dsubcld(pcols) ! wg layer thickness in mbs between lcl and maxi. - real(r8), intent(out) :: jctop(pcols) ! o row of top-of-deep-convection indices passed out. - real(r8), intent(out) :: jcbot(pcols) ! o row of base of cloud indices passed out. - real(r8), intent(out) :: prec(pcols) - real(r8), intent(out) :: rliq(pcols) ! reserved liquid (not yet in cldliq) for energy integrals - real(r8), intent(out) :: dcape(pcols) ! output dynamical CAPE - -!songxl 2014-05-20------------------ - - real(r8) zs(pcols) - real(r8) dlg(pcols,pver) ! gathrd version of the detraining cld h2o tend - real(r8) pflxg(pcols,pverp) ! gather precip flux at each level - real(r8) cug(pcols,pver) ! gathered condensation rate - real(r8) evpg(pcols,pver) ! gathered evap rate of rain in downdraft - real(r8) mumax(pcols) - integer jt(pcols) ! wg top level index of deep cumulus convection. - integer maxg(pcols) ! wg gathered values of maxi. - integer ideep(pcols) ! w holds position of gathered points vs longitude index. - integer lengath -! diagnostic field used by chem/wetdep codes - real(r8) ql(pcols,pver) ! wg grid slice of cloud liquid water. -! - real(r8) pblt(pcols) ! i row of pbl top indices. - - - - -! -!----------------------------------------------------------------------- -! -! general work fields (local variables): -! - real(r8) q(pcols,pver) ! w grid slice of mixing ratio. - real(r8) p(pcols,pver) ! w grid slice of ambient mid-layer pressure in mbs. - real(r8) z(pcols,pver) ! w grid slice of ambient mid-layer height in metres. - real(r8) s(pcols,pver) ! w grid slice of scaled dry static energy (t+gz/cp). - real(r8) tp(pcols,pver) ! w grid slice of parcel temperatures. - real(r8) zf(pcols,pver+1) ! w grid slice of ambient interface height in metres. - real(r8) pf(pcols,pver+1) ! w grid slice of ambient interface pressure in mbs. - real(r8) qstp(pcols,pver) ! w grid slice of parcel temp. saturation mixing ratio. - - real(r8) tl(pcols) ! w row of parcel temperature at lcl. -!songxl 2014-05-20----------------- - - logical iclosure ! switch on sequence of call to buoyan_dilute to derive DCAPE - real(r8) capelmt_wk ! work capelmt to allow diff values passed to closure with trigdcape - - integer lcl(pcols) ! w base level index of deep cumulus convection. - integer lel(pcols) ! w index of highest theoretical convective plume. - - integer lon(pcols) ! w index of onset level for deep convection. - integer maxi(pcols) ! w index of level with largest moist static energy. - integer index(pcols) -! -! gathered work fields: -! - real(r8) qg(pcols,pver) ! wg grid slice of gathered values of q. - real(r8) tg(pcols,pver) ! w grid slice of temperature at interface. - real(r8) pg(pcols,pver) ! wg grid slice of gathered values of p. - real(r8) zg(pcols,pver) ! wg grid slice of gathered values of z. - real(r8) sg(pcols,pver) ! wg grid slice of gathered values of s. - real(r8) tpg(pcols,pver) ! wg grid slice of gathered values of tp. - real(r8) zfg(pcols,pver+1) ! wg grid slice of gathered values of zf. - real(r8) qstpg(pcols,pver) ! wg grid slice of gathered values of qstp. - real(r8) ug(pcols,pver) ! wg grid slice of gathered values of u. - real(r8) vg(pcols,pver) ! wg grid slice of gathered values of v. - real(r8) cmeg(pcols,pver) - - real(r8) hu_nm1g(pcols,pver) !songxl 2014-05-20 - - real(r8) rprdg(pcols,pver) ! wg gathered rain production rate - real(r8) capeg(pcols) ! wg gathered convective available potential energy. - real(r8) tlg(pcols) ! wg grid slice of gathered values of tl. - real(r8) landfracg(pcols) ! wg grid slice of landfrac - - integer lclg(pcols) ! wg gathered values of lcl. - integer lelg(pcols) -! -! work fields arising from gathered calculations. -! - real(r8) dqdt(pcols,pver) ! wg mixing ratio tendency at gathered points. - real(r8) dsdt(pcols,pver) ! wg dry static energy ("temp") tendency at gathered points. -! real(r8) alpha(pcols,pver) ! array of vertical differencing used (=1. for upstream). - real(r8) sd(pcols,pver) ! wg grid slice of dry static energy in downdraft. - real(r8) qd(pcols,pver) ! wg grid slice of mixing ratio in downdraft. - real(r8) mc(pcols,pver) ! wg net upward (scaled by mb) cloud mass flux. - real(r8) qhat(pcols,pver) ! wg grid slice of upper interface mixing ratio. - real(r8) qu(pcols,pver) ! wg grid slice of mixing ratio in updraft. - real(r8) su(pcols,pver) ! wg grid slice of dry static energy in updraft. - real(r8) qs(pcols,pver) ! wg grid slice of saturation mixing ratio. - real(r8) shat(pcols,pver) ! wg grid slice of upper interface dry static energy. - real(r8) hmn(pcols,pver) ! wg moist static energy. - real(r8) hsat(pcols,pver) ! wg saturated moist static energy. - real(r8) qlg(pcols,pver) - real(r8) dudt(pcols,pver) ! wg u-wind tendency at gathered points. - real(r8) dvdt(pcols,pver) ! wg v-wind tendency at gathered points. -! real(r8) ud(pcols,pver) -! real(r8) vd(pcols,pver) - - real(r8) mb(pcols) ! wg cloud base mass flux. - - integer jlcl(pcols) - integer j0(pcols) ! wg detrainment initiation level index. - integer jd(pcols) ! wg downdraft initiation level index. - - real(r8) delt ! length of model time-step in seconds. - - integer i - integer ii - integer k - integer msg ! ic number of missing moisture levels at the top of model. - - real(r8) qdifr - real(r8) sdifr - - logical :: is_first_step = .false. - -! -!--------------------------Data statements------------------------------ -! -! -! Set internal variable "msg" (convection limit) to "limcnv-1" -! - msg = limcnv - 1 -! -! initialize necessary arrays. -! zero out variables not used in cam -! - qtnd(:,:) = 0.0_r8 - heat(:,:) = 0.0_r8 - !real(r8), intent(out) :: mcon(pcols,pverp) -! mcon( :, :(pver+1) ) = !Meredith - This is the line that when uncommented causes issues. - rliq(:ncol) = 0.0_r8 - -! -! initialize convective tendencies -! - prec(:ncol) = 0._r8 - do k = 1,pver - do i = 1,ncol - dqdt(i,k) = 0._r8 - dsdt(i,k) = 0._r8 - dudt(i,k) = 0._r8 - dvdt(i,k) = 0._r8 - pflx(i,k) = 0._r8 - pflxg(i,k) = 0._r8 - cme(i,k) = 0._r8 - rprd(i,k) = 0._r8 - zdu(i,k) = 0._r8 - ql(i,k) = 0._r8 - qlg(i,k) = 0._r8 - dlf(i,k) = 0._r8 - dlg(i,k) = 0._r8 - end do - end do - do i = 1,ncol - pflx(i,pverp) = 0 - pflxg(i,pverp) = 0 - mcon(i, pverp) = 0 - end do -! - do i = 1,ncol - pblt(i) = pver - dsubcld(i) = 0._r8 - - jctop(i) = pver - jcbot(i) = 1 - if(trigmem)dcape(i) = 0._r8 !songxl 2014-05-20 - end do - - -! -! calculate local pressure (mbs) and height (m) for both interface -! and mid-layer locations. -! - do i = 1,ncol - zs(i) = geos(i)*rgrav - pf(i,pver+1) = paph(i,pver+1)*0.01_r8 - zf(i,pver+1) = zi(i,pver+1) + zs(i) - end do - do k = 1,pver - do i = 1,ncol - p(i,k) = pap(i,k)*0.01_r8 - pf(i,k) = paph(i,k)*0.01_r8 - z(i,k) = zm(i,k) + zs(i) - zf(i,k) = zi(i,k) + zs(i) - end do - end do -! - do k = pver - 1,msg + 1,-1 - do i = 1,ncol - if (abs(z(i,k)-zs(i)-pblh(i)) < (zf(i,k)-zf(i,k+1))*0.5_r8) pblt(i) = k - end do - end do -! -! store incoming specific humidity field for subsequent calculation -! of precipitation (through change in storage). -! define dry static energy (normalized by cp). -! - do k = 1,pver - do i = 1,ncol - q(i,k) = qh(i,k) - s(i,k) = t(i,k) + (grav/cpres)*z(i,k) - tp(i,k)=0.0_r8 - if(trigmem)tpm1(i,k) = 0.0_r8 !songxl 2014-05-20 - shat(i,k) = s(i,k) - qhat(i,k) = q(i,k) - end do - end do - -!songxl 2014-05-20--------------------- - - do i = 1,ncol - capeg(i) = 0._r8 - lclg(i) = 1 - lelg(i) = pver - maxg(i) = 1 - tlg(i) = 400._r8 - dsubcld(i) = 0._r8 - end do - - - ! Evaluate Tparcel, qs(Tparcel), buoyancy and CAPE, - ! lcl, lel, parcel launch level at index maxi()=hmax - ! - - ! 1. First call, iclosure = .true., standard calculation, scanning for launching level up to 600 hPa - ! 2. Second call, iclosure = .faklse. pass the launch level from 1st call to determine CAPE at previous step - ! The differewnce of CAPE values from the two calls is DCAPE, based on the same launch level - - iclosure = .true. - call buoyan_dilute(lchnk ,ncol , & - q ,t ,p ,z ,pf , & - tp ,qstp ,tl ,rl ,cape , & - pblt ,lcl ,lel ,lon ,maxi , & - rgas ,grav ,cpres ,msg , & - tpert ,iclosure) - - ! if (trigdcape_ull) then - ! if (.not. allocated(dcapemx)) then - ! allocate (dcapemx(pcols), stat=ierror) - ! if ( ierror /= 0 ) call endrun('ZM_CONVR error: allocation error dcapemx') - ! endif - ! dcapemx(:ncol) = maxi(:ncol) - ! endif - - if(trigmem)then - call buoyan_dilute(lchnk ,ncol , & - qm1 ,tm1 ,p ,z ,pf , & - tpm1 ,qstpm1 ,tlm1 ,rl ,capem1 , & - pblt ,lclm1 ,lelm1 ,lonm1 ,maxim1 , & - rgas ,grav ,cpres ,msg , & - tpert ,iclosure) - - do i=1,ncol - dcape(i) = (cape(i)-capem1(i))/(delt*2._r8) - end do - endif - - !DCAPE-ULL - if (.not. is_first_step .and. trigdcape_ull) then - iclosure = .false. - call buoyan_dilute(lchnk ,ncol , & - q_star ,t_star ,p ,z ,pf , & - tpm1 ,qstpm1 ,tlm1 ,rl ,capem1 , & - pblt ,lclm1 ,lelm1 ,lonm1 ,maxim1 , & - rgas ,grav ,cpres ,msg , & - tpert ,iclosure) - - dcape(:ncol) = (cape(:ncol)-capem1(:ncol))/(delt*2._r8) - endif - -! -! determine whether grid points will undergo some deep convection -! (ideep=1) or not (ideep=0), based on values of cape,lcl,lel -! (require cape.gt. 0 and lel capelmt) then - lengath = lengath + 1 - index(lengath) = i - end if - else - if (dcape(i) > dcapelmt) then - lengath = lengath + 1 - index(lengath) = i - end if - end if - else if (trigdcape_ull) then - ! DCAPE-ULL - if (is_first_step) then - !Will this cause restart to be non-BFB - if (cape(i) > capelmt) then - lengath = lengath + 1 - index(lengath) = i - end if - else if (cape(i) > 0.0_r8 .and. dcape(i) > trigdcapelmt) then - ! use constant 0 or a separate threshold for capt because capelmt is for default trigger - lengath = lengath + 1 - index(lengath) = i - endif - else - if (cape(i) > capelmt) then - lengath = lengath + 1 - index(lengath) = i - end if - end if -!>songxl 2014-05-20---------------- - end do - - if (lengath.eq.0) return - do ii=1,lengath - i=index(ii) - ideep(ii)=i - end do -! -! obtain gathered arrays necessary for ensuing calculations. -! - do k = 1,pver - do i = 1,lengath - dp(i,k) = 0.01_r8*dpp(ideep(i),k) - qg(i,k) = q(ideep(i),k) - tg(i,k) = t(ideep(i),k) - pg(i,k) = p(ideep(i),k) - zg(i,k) = z(ideep(i),k) - sg(i,k) = s(ideep(i),k) - tpg(i,k) = tp(ideep(i),k) - zfg(i,k) = zf(ideep(i),k) - qstpg(i,k) = qstp(ideep(i),k) - ug(i,k) = 0._r8 - vg(i,k) = 0._r8 - if(trigmem)hu_nm1g(i,k) = hu_nm1(ideep(i),k) !songxl 2014-05-20 - end do - end do -! - do i = 1,lengath - zfg(i,pver+1) = zf(ideep(i),pver+1) - end do - do i = 1,lengath - capeg(i) = cape(ideep(i)) - lclg(i) = lcl(ideep(i)) - lelg(i) = lel(ideep(i)) - maxg(i) = maxi(ideep(i)) - tlg(i) = tl(ideep(i)) - landfracg(i) = landfrac(ideep(i)) - end do -! -! calculate sub-cloud layer pressure "thickness" for use in -! closure and tendency routines. -! - do k = msg + 1,pver - do i = 1,lengath - if (k >= maxg(i)) then - dsubcld(i) = dsubcld(i) + dp(i,k) - end if - end do - end do -! -! define array of factors (alpha) which defines interfacial -! values, as well as interfacial values for (q,s) used in -! subsequent routines. -! - do k = msg + 2,pver - do i = 1,lengath -! alpha(i,k) = 0.5 - sdifr = 0._r8 - qdifr = 0._r8 - if (sg(i,k) > 0._r8 .or. sg(i,k-1) > 0._r8) & - sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k))) - if (qg(i,k) > 0._r8 .or. qg(i,k-1) > 0._r8) & - qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k))) - if (sdifr > 1.E-6_r8) then - shat(i,k) = log(sg(i,k-1)/sg(i,k))*sg(i,k-1)*sg(i,k)/(sg(i,k-1)-sg(i,k)) - else - shat(i,k) = 0.5_r8* (sg(i,k)+sg(i,k-1)) - end if - if (qdifr > 1.E-6_r8) then - qhat(i,k) = log(qg(i,k-1)/qg(i,k))*qg(i,k-1)*qg(i,k)/(qg(i,k-1)-qg(i,k)) - else - qhat(i,k) = 0.5_r8* (qg(i,k)+qg(i,k-1)) - end if - end do - end do -! -! obtain cloud properties. -! - - call cldprp(lchnk , & - qg ,tg ,ug ,vg ,pg , & - zg ,sg ,mu ,eu ,du , & - md ,ed ,sd ,qd ,mc , & - qu ,su ,zfg ,qs ,hmn , & - hsat ,shat ,qlg , & - cmeg ,maxg ,lelg ,jt ,jlcl , & - maxg ,j0 ,jd ,rl ,lengath , & - rgas ,grav ,cpres ,msg , & - pflxg ,evpg ,cug ,rprdg ,limcnv , & - landfracg, hu_nm1g, tpert) !songxl 2014-05-20 - !PMA adds tpert to the calculation -! -! convert detrainment from units of "1/m" to "1/mb". -! - do k = msg + 1,pver - do i = 1,lengath - du (i,k) = du (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) - eu (i,k) = eu (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) - ed (i,k) = ed (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) - cug (i,k) = cug (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) - cmeg (i,k) = cmeg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) - rprdg(i,k) = rprdg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) - evpg (i,k) = evpg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) - end do - end do - - call closure(lchnk , & - qg ,tg ,pg ,zg ,sg , & - tpg ,qs ,qu ,su ,mc , & - du ,mu ,md ,qd ,sd , & - qhat ,shat ,dp ,qstpg ,zfg , & - qlg ,dsubcld ,mb ,capeg ,tlg , & - lclg ,lelg ,jt ,maxg ,1 , & - lengath ,rgas ,grav ,cpres ,rl , & - msg ,capelmt_wk ) -! -! limit cloud base mass flux to theoretical upper bound. -! - do i=1,lengath - mumax(i) = 0 - end do - do k=msg + 2,pver - do i=1,lengath - mumax(i) = max(mumax(i), mu(i,k)/dp(i,k)) - end do - end do - - do i=1,lengath - if (mumax(i) > 0._r8) then - mb(i) = min(mb(i),0.5_r8/(delt*mumax(i))) - else - mb(i) = 0._r8 - endif - end do - ! If no_deep_pbl = .true., don't allow convection entirely - ! within PBL (suggestion of Bjorn Stevens, 8-2000) - - if (no_deep_pbl) then - do i=1,lengath - if (zm(ideep(i),jt(i)) < pblh(ideep(i))) mb(i) = 0 - end do - end if - - - do k=msg+1,pver - do i=1,lengath - mu (i,k) = mu (i,k)*mb(i) - md (i,k) = md (i,k)*mb(i) - mc (i,k) = mc (i,k)*mb(i) - du (i,k) = du (i,k)*mb(i) - eu (i,k) = eu (i,k)*mb(i) - ed (i,k) = ed (i,k)*mb(i) - cmeg (i,k) = cmeg (i,k)*mb(i) - rprdg(i,k) = rprdg(i,k)*mb(i) - cug (i,k) = cug (i,k)*mb(i) - evpg (i,k) = evpg (i,k)*mb(i) - pflxg(i,k+1)= pflxg(i,k+1)*mb(i)*100._r8/grav - end do - end do -! -! compute temperature and moisture changes due to convection. -! - call q1q2_pjr(lchnk , & - dqdt ,dsdt ,qg ,qs ,qu , & - su ,du ,qhat ,shat ,dp , & - mu ,md ,sd ,qd ,qlg , & - dsubcld ,jt ,maxg ,1 ,lengath , & - cpres ,rl ,msg , & - dlg ,evpg ,cug ) -! -!songxl 2014-05-20----------------- - -! gather back temperature and mixing ratio. -! - do k = msg + 1,pver -!DIR$ CONCURRENT - do i = 1,lengath -! -! q is updated to compute net precip. -! - q(ideep(i),k) = qh(ideep(i),k) + 2._r8*delt*dqdt(i,k) - qtnd(ideep(i),k) = dqdt (i,k) - cme (ideep(i),k) = cmeg (i,k) - rprd(ideep(i),k) = rprdg(i,k) - zdu (ideep(i),k) = du (i,k) - mcon(ideep(i),k) = mc (i,k) - heat(ideep(i),k) = dsdt (i,k)*cpres - dlf (ideep(i),k) = dlg (i,k) - pflx(ideep(i),k) = pflxg(i,k) - ql (ideep(i),k) = qlg (i,k) -!songxl 2014-05-20 - end do - end do -! -!DIR$ CONCURRENT - do i = 1,lengath - jctop(ideep(i)) = jt(i) -!++bee - jcbot(ideep(i)) = maxg(i) -!--bee - pflx(ideep(i),pverp) = pflxg(i,pverp) - end do - -! Compute precip by integrating change in water vapor minus detrained cloud water - do k = pver,msg + 1,-1 - do i = 1,ncol - prec(i) = prec(i) - dpp(i,k)* (q(i,k)-qh(i,k)) - dpp(i,k)*dlf(i,k)*2*delt - end do - end do - -! obtain final precipitation rate in m/s. - do i = 1,ncol - prec(i) = rgrav*max(prec(i),0._r8)/ (2._r8*delt)/1000._r8 - end do - -! Compute reserved liquid (not yet in cldliq) for energy integrals. -! Treat rliq as flux out bottom, to be added back later. - do k = 1, pver - do i = 1, ncol - rliq(i) = rliq(i) + dlf(i,k)*dpp(i,k)/gravit - end do - end do - rliq(:ncol) = rliq(:ncol) /1000._r8 - - is_first_step_local = .false. - return -end subroutine zm_convr - -!=============================================================================== -subroutine zm_conv_evap(ncol,lchnk, & - t,pmid,pdel,q, & - tend_s, tend_s_snwprd, tend_s_snwevmlt, tend_q, & - prdprec, cldfrc, deltat, & - prec, snow, ntprprd, ntsnprd, flxprec, flxsnow ) - -!----------------------------------------------------------------------- -! Compute tendencies due to evaporation of rain from ZM scheme -!-- -! Compute the total precipitation and snow fluxes at the surface. -! Add in the latent heat of fusion for snow formation and melt, since it not dealt with -! in the Zhang-MacFarlane parameterization. -! Evaporate some of the precip directly into the environment using a Sundqvist type algorithm -!----------------------------------------------------------------------- - -! use wv_saturation, only: qsat - !use phys_grid, only: get_rlat_all_p - -!------------------------------Arguments-------------------------------- - integer,intent(in) :: ncol, lchnk ! number of columns and chunk index - real(r8),intent(in), dimension(pcols,pver) :: t ! temperature (K) - real(r8),intent(in), dimension(pcols,pver) :: pmid ! midpoint pressure (Pa) - real(r8),intent(in), dimension(pcols,pver) :: pdel ! layer thickness (Pa) - real(r8),intent(in), dimension(pcols,pver) :: q ! water vapor (kg/kg) - real(r8),intent(inout), dimension(pcols,pver) :: tend_s ! heating rate (J/kg/s) - real(r8),intent(inout), dimension(pcols,pver) :: tend_q ! water vapor tendency (kg/kg/s) - real(r8),intent(out ), dimension(pcols,pver) :: tend_s_snwprd ! Heating rate of snow production - real(r8),intent(out ), dimension(pcols,pver) :: tend_s_snwevmlt ! Heating rate of evap/melting of snow - - - - real(r8), intent(in ) :: prdprec(pcols,pver)! precipitation production (kg/ks/s) - real(r8), intent(in ) :: cldfrc(pcols,pver) ! cloud fraction - real(r8), intent(in ) :: deltat ! time step - - real(r8), intent(inout) :: prec(pcols) ! Convective-scale preciptn rate - real(r8), intent(out) :: snow(pcols) ! Convective-scale snowfall rate -! -!---------------------------Local storage------------------------------- - - real(r8) :: es (pcols,pver) ! Saturation vapor pressure - real(r8) :: fice (pcols,pver) ! ice fraction in precip production - real(r8) :: fsnow_conv(pcols,pver) ! snow fraction in precip production - real(r8) :: qs (pcols,pver) ! saturation specific humidity - real(r8),intent(out) :: flxprec(pcols,pverp) ! Convective-scale flux of precip at interfaces (kg/m2/s) - real(r8),intent(out) :: flxsnow(pcols,pverp) ! Convective-scale flux of snow at interfaces (kg/m2/s) - real(r8),intent(out) :: ntprprd(pcols,pver) ! net precip production in layer - real(r8),intent(out) :: ntsnprd(pcols,pver) ! net snow production in layer - real(r8) :: work1 ! temp variable (pjr) - real(r8) :: work2 ! temp variable (pjr) - - real(r8) :: evpvint(pcols) ! vertical integral of evaporation - real(r8) :: evpprec(pcols) ! evaporation of precipitation (kg/kg/s) - real(r8) :: evpsnow(pcols) ! evaporation of snowfall (kg/kg/s) - real(r8) :: snowmlt(pcols) ! snow melt tendency in layer - real(r8) :: flxsntm(pcols) ! flux of snow into layer, after melting - - real(r8) :: evplimit ! temp variable for evaporation limits - - integer :: i,k ! longitude,level indices - - -!----------------------------------------------------------------------- - -! convert input precip to kg/m2/s - prec(:ncol) = prec(:ncol)*1000._r8 - -! determine saturation vapor pressure - do i = 1,ncol - do k = 1,pver - call qsat(t(i,k), pmid(i,k), & - es(i,k), qs(i,k), 0) - end do - end do - -! determine ice fraction in rain production (use cloud water parameterization fraction at present) - call cldfrc_fice(ncol, t, fice, fsnow_conv) - -! zero the flux integrals on the top boundary - flxprec(:ncol,1) = 0._r8 - flxsnow(:ncol,1) = 0._r8 - evpvint(:ncol) = 0._r8 - - do k = 1, pver - do i = 1, ncol - -! Melt snow falling into layer, if necessary. - if (t(i,k) > tmelt) then - flxsntm(i) = 0._r8 - snowmlt(i) = flxsnow(i,k) * gravit/ pdel(i,k) - else - flxsntm(i) = flxsnow(i,k) - snowmlt(i) = 0._r8 - end if - -! relative humidity depression must be > 0 for evaporation - evplimit = max(1._r8 - q(i,k)/qs(i,k), 0._r8) - -! total evaporation depends on flux in the top of the layer -! flux prec is the net production above layer minus evaporation into environmet - evpprec(i) = ke * (1._r8 - cldfrc(i,k)) * evplimit * sqrt(flxprec(i,k)) -!********************************************************** -!! evpprec(i) = 0. ! turn off evaporation for now -!********************************************************** - -! Don't let evaporation supersaturate layer (approx). Layer may already be saturated. -! Currently does not include heating/cooling change to qs - evplimit = max(0._r8, (qs(i,k)-q(i,k)) / deltat) - -! Don't evaporate more than is falling into the layer - do not evaporate rain formed -! in this layer but if precip production is negative, remove from the available precip -! Negative precip production occurs because of evaporation in downdrafts. -!!$ evplimit = flxprec(i,k) * gravit / pdel(i,k) + min(prdprec(i,k), 0.) - evplimit = min(evplimit, flxprec(i,k) * gravit / pdel(i,k)) - -! Total evaporation cannot exceed input precipitation - evplimit = min(evplimit, (prec(i) - evpvint(i)) * gravit / pdel(i,k)) - - evpprec(i) = min(evplimit, evpprec(i)) - -! evaporation of snow depends on snow fraction of total precipitation in the top after melting - if (flxprec(i,k) > 0._r8) then -! evpsnow(i) = evpprec(i) * flxsntm(i) / flxprec(i,k) -! prevent roundoff problems - work1 = min(max(0._r8,flxsntm(i)/flxprec(i,k)),1._r8) - evpsnow(i) = evpprec(i) * work1 - else - evpsnow(i) = 0._r8 - end if - -! vertically integrated evaporation - evpvint(i) = evpvint(i) + evpprec(i) * pdel(i,k)/gravit - -! net precip production is production - evaporation - ntprprd(i,k) = prdprec(i,k) - evpprec(i) -! net snow production is precip production * ice fraction - evaporation - melting -!pjrworks ntsnprd(i,k) = prdprec(i,k)*fice(i,k) - evpsnow(i) - snowmlt(i) -!pjrwrks2 ntsnprd(i,k) = prdprec(i,k)*fsnow_conv(i,k) - evpsnow(i) - snowmlt(i) -! the small amount added to flxprec in the work1 expression has been increased from -! 1e-36 to 8.64e-11 (1e-5 mm/day). This causes the temperature based partitioning -! scheme to be used for small flxprec amounts. This is to address error growth problems. -#ifdef PERGRO - work1 = min(max(0._r8,flxsnow(i,k)/(flxprec(i,k)+8.64e-11_r8)),1._r8) -#else - if (flxprec(i,k).gt.0._r8) then - work1 = min(max(0._r8,flxsnow(i,k)/flxprec(i,k)),1._r8) - else - work1 = 0._r8 - endif -#endif - work2 = max(fsnow_conv(i,k), work1) - if (snowmlt(i).gt.0._r8) work2 = 0._r8 -! work2 = fsnow_conv(i,k) - ntsnprd(i,k) = prdprec(i,k)*work2 - evpsnow(i) - snowmlt(i) - tend_s_snwprd (i,k) = prdprec(i,k)*work2*latice - tend_s_snwevmlt(i,k) = - ( evpsnow(i) + snowmlt(i) )*latice - -! precipitation fluxes - flxprec(i,k+1) = flxprec(i,k) + ntprprd(i,k) * pdel(i,k)/gravit - flxsnow(i,k+1) = flxsnow(i,k) + ntsnprd(i,k) * pdel(i,k)/gravit - -! protect against rounding error - flxprec(i,k+1) = max(flxprec(i,k+1), 0._r8) - flxsnow(i,k+1) = max(flxsnow(i,k+1), 0._r8) -! more protection (pjr) -! flxsnow(i,k+1) = min(flxsnow(i,k+1), flxprec(i,k+1)) - -! heating (cooling) and moistening due to evaporation -! - latent heat of vaporization for precip production has already been accounted for -! - snow is contained in prec - tend_s(i,k) =-evpprec(i)*latvap + ntsnprd(i,k)*latice - tend_q(i,k) = evpprec(i) - end do - end do - -! set output precipitation rates (m/s) - prec(:ncol) = flxprec(:ncol,pver+1) / 1000._r8 - snow(:ncol) = flxsnow(:ncol,pver+1) / 1000._r8 - -!********************************************************** -!!$ tend_s(:ncol,:) = 0. ! turn heating off -!********************************************************** - - end subroutine zm_conv_evap - - -! AaronDonahue - skip convtran and momtran for now -subroutine convtran(lchnk , & - doconvtran,q ,ncnst ,mu ,md , & - du ,eu ,ed ,dp ,dsubcld , & - jt ,mx ,ideep ,il1g ,il2g , & - nstep ,fracis ,dqdt ,dpdry ) -!----------------------------------------------------------------------- -! -! Purpose: -! Convective transport of trace species -! -! Mixing ratios may be with respect to either dry or moist air -! -! Method: -! -! -! -! Author: P. Rasch -! -!----------------------------------------------------------------------- -! use constituents, only: cnst_get_type_byind -! use ppgrid - use scream_abortutils, only: endscreamrun - - implicit none -!----------------------------------------------------------------------- -! -! Input arguments -! - integer, intent(in) :: lchnk ! chunk identifier - integer, intent(in) :: ncnst ! number of tracers to transport - logical, intent(in) :: doconvtran(ncnst) ! flag for doing convective transport - real(r8), intent(in) :: q(pcols,pver,ncnst) ! Tracer array including moisture - real(r8), intent(in) :: mu(pcols,pver) ! Mass flux up - real(r8), intent(in) :: md(pcols,pver) ! Mass flux down - real(r8), intent(in) :: du(pcols,pver) ! Mass detraining from updraft - real(r8), intent(in) :: eu(pcols,pver) ! Mass entraining from updraft - real(r8), intent(in) :: ed(pcols,pver) ! Mass entraining from downdraft - real(r8), intent(in) :: dp(pcols,pver) ! Delta pressure between interfaces - real(r8), intent(in) :: dsubcld(pcols) ! Delta pressure from cloud base to sfc - real(r8), intent(in) :: fracis(pcols,pver,ncnst) ! fraction of tracer that is insoluble - - integer, intent(in) :: jt(pcols) ! Index of cloud top for each column - integer, intent(in) :: mx(pcols) ! Index of cloud top for each column - integer, intent(in) :: ideep(pcols) ! Gathering array - integer, intent(in) :: il1g ! Gathered min lon indices over which to operate - integer, intent(in) :: il2g ! Gathered max lon indices over which to operate - integer, intent(in) :: nstep ! Time step index - - real(r8), intent(in) :: dpdry(pcols,pver) ! Delta pressure between interfaces - - -! input/output - - real(r8), intent(out) :: dqdt(pcols,pver,ncnst) ! Tracer tendency array - -!--------------------------Local Variables------------------------------ - - integer i ! Work index - integer k ! Work index - integer kbm ! Highest altitude index of cloud base - integer kk ! Work index - integer kkp1 ! Work index - integer km1 ! Work index - integer kp1 ! Work index - integer ktm ! Highest altitude index of cloud top - integer m ! Work index - - real(r8) cabv ! Mix ratio of constituent above - real(r8) cbel ! Mix ratio of constituent below - real(r8) cdifr ! Normalized diff between cabv and cbel - real(r8) chat(pcols,pver) ! Mix ratio in env at interfaces - real(r8) cond(pcols,pver) ! Mix ratio in downdraft at interfaces - real(r8) const(pcols,pver) ! Gathered tracer array - real(r8) fisg(pcols,pver) ! gathered insoluble fraction of tracer - real(r8) conu(pcols,pver) ! Mix ratio in updraft at interfaces - real(r8) dcondt(pcols,pver) ! Gathered tend array - real(r8) small ! A small number - real(r8) mbsth ! Threshold for mass fluxes - real(r8) mupdudp ! A work variable - real(r8) minc ! A work variable - real(r8) maxc ! A work variable - real(r8) fluxin ! A work variable - real(r8) fluxout ! A work variable - real(r8) netflux ! A work variable - - real(r8) dutmp(pcols,pver) ! Mass detraining from updraft - real(r8) eutmp(pcols,pver) ! Mass entraining from updraft - real(r8) edtmp(pcols,pver) ! Mass entraining from downdraft - real(r8) dptmp(pcols,pver) ! Delta pressure between interfaces -!----------------------------------------------------------------------- -! - small = 1.e-36_r8 -! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s) - mbsth = 1.e-15_r8 - -! Find the highest level top and bottom levels of convection - ktm = pver - kbm = pver - do i = il1g, il2g - ktm = min(ktm,jt(i)) - kbm = min(kbm,mx(i)) - end do - -! Loop ever each constituent - do m = 2, ncnst - if (doconvtran(m)) then - -! TODO: Meredith Fossitt, find a way to sub dependency for -! cnst_get_type_blind -! if (cnst_get_type_byind(m).eq.'dry') then - if (.true.) then - do k = 1,pver - do i =il1g,il2g - dptmp(i,k) = dpdry(i,k) - dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k) - eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k) - edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k) - end do - end do - else - do k = 1,pver - do i =il1g,il2g - dptmp(i,k) = dp(i,k) - dutmp(i,k) = du(i,k) - eutmp(i,k) = eu(i,k) - edtmp(i,k) = ed(i,k) - end do - end do - endif -! dptmp = dp - -! Gather up the constituent and set tend to zero - do k = 1,pver - do i =il1g,il2g - const(i,k) = q(ideep(i),k,m) - fisg(i,k) = fracis(ideep(i),k,m) - end do - end do - -! From now on work only with gathered data - -! Interpolate environment tracer values to interfaces - do k = 1,pver - km1 = max(1,k-1) - do i = il1g, il2g - minc = min(const(i,km1),const(i,k)) - maxc = max(const(i,km1),const(i,k)) - if (minc < 0) then - cdifr = 0._r8 - else - cdifr = abs(const(i,k)-const(i,km1))/max(maxc,small) - endif - -! If the two layers differ significantly use a geometric averaging -! procedure - if (cdifr > 1.E-6_r8) then - cabv = max(const(i,km1),maxc*1.e-12_r8) - cbel = max(const(i,k),maxc*1.e-12_r8) - chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel - - else ! Small diff, so just arithmetic mean - chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1)) - end if - -! Provisional up and down draft values - conu(i,k) = chat(i,k) - cond(i,k) = chat(i,k) - -! provisional tends - dcondt(i,k) = 0._r8 - - end do - end do - -! Do levels adjacent to top and bottom - k = 2 - km1 = 1 - kk = pver - do i = il1g,il2g - mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk) - if (mupdudp > mbsth) then - conu(i,kk) = (+eutmp(i,kk)*fisg(i,kk)*const(i,kk)*dptmp(i,kk))/mupdudp - endif - if (md(i,k) < -mbsth) then - cond(i,k) = (-edtmp(i,km1)*fisg(i,km1)*const(i,km1)*dptmp(i,km1))/md(i,k) - endif - end do - -! Updraft from bottom to top - do kk = pver-1,1,-1 - kkp1 = min(pver,kk+1) - do i = il1g,il2g - mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk) - if (mupdudp > mbsth) then - conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* & - const(i,kk)*dptmp(i,kk) )/mupdudp - endif - end do - end do - -! Downdraft from top to bottom - do k = 3,pver - km1 = max(1,k-1) - do i = il1g,il2g - if (md(i,k) < -mbsth) then - cond(i,k) = ( md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) & - *dptmp(i,km1) )/md(i,k) - endif - end do - end do - - - do k = ktm,pver - km1 = max(1,k-1) - kp1 = min(pver,k+1) - do i = il1g,il2g - -! version 1 hard to check for roundoff errors -! dcondt(i,k) = -! $ +(+mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) -! $ -mu(i,k)* (conu(i,k)-chat(i,k)) -! $ +md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) -! $ -md(i,k)* (cond(i,k)-chat(i,k)) -! $ )/dp(i,k) - -! version 2 hard to limit fluxes -! fluxin = mu(i,kp1)*conu(i,kp1) + mu(i,k)*chat(i,k) -! $ -(md(i,k) *cond(i,k) + md(i,kp1)*chat(i,kp1)) -! fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*chat(i,kp1) -! $ -(md(i,kp1)*cond(i,kp1) + md(i,k)*chat(i,k)) - -! version 3 limit fluxes outside convection to mass in appropriate layer -! these limiters are probably only safe for positive definite quantitities -! it assumes that mu and md already satify a courant number limit of 1 - fluxin = mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) & - -(md(i,k) *cond(i,k) + md(i,kp1)*min(chat(i,kp1),const(i,kp1))) - fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*min(chat(i,kp1),const(i,k)) & - -(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k))) - - netflux = fluxin - fluxout - if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then - netflux = 0._r8 - endif - dcondt(i,k) = netflux/dptmp(i,k) - end do - end do -! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! -!DIR$ NOINTERCHANGE - do k = kbm,pver - km1 = max(1,k-1) - do i = il1g,il2g - if (k == mx(i)) then - -! version 1 -! dcondt(i,k) = (1./dsubcld(i))* -! $ (-mu(i,k)*(conu(i,k)-chat(i,k)) -! $ -md(i,k)*(cond(i,k)-chat(i,k)) -! $ ) - -! version 2 -! fluxin = mu(i,k)*chat(i,k) - md(i,k)*cond(i,k) -! fluxout = mu(i,k)*conu(i,k) - md(i,k)*chat(i,k) -! version 3 - fluxin = mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k) - fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k)) - - netflux = fluxin - fluxout - if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then - netflux = 0._r8 - endif -! dcondt(i,k) = netflux/dsubcld(i) - dcondt(i,k) = netflux/dptmp(i,k) - else if (k > mx(i)) then -! dcondt(i,k) = dcondt(i,k-1) - dcondt(i,k) = 0._r8 - end if - end do - end do - -! Initialize to zero everywhere, then scatter tendency back to full array - dqdt(:,:,m) = 0._r8 - do k = 1,pver - kp1 = min(pver,k+1) -!DIR$ CONCURRENT - do i = il1g,il2g - dqdt(ideep(i),k,m) = dcondt(i,k) - end do - end do - - end if ! for doconvtran - - end do - - return -end subroutine convtran - -!========================================================================================= - -subroutine momtran(lchnk, ncol, & - domomtran,q ,ncnst ,mu ,md , & - du ,eu ,ed ,dp ,dsubcld , & - jt ,mx ,ideep ,il1g ,il2g , & - nstep ,dqdt ,pguall ,pgdall, icwu, icwd, dt, seten ) -!----------------------------------------------------------------------- -! -! Purpose: -! Convective transport of momentum -! -! Mixing ratios may be with respect to either dry or moist air -! -! Method: -! Based on the convtran subroutine by P. Rasch -! -! -! Author: J. Richter and P. Rasch -! -!----------------------------------------------------------------------- - use physics_utils, only: r8 => rtype, rtype8, itype, btype -! use constituents, only: cnst_get_type_byind -! use ppgrid - use scream_abortutils, only: endscreamrun - implicit none -!----------------------------------------------------------------------- -! -! Input arguments -! - integer, intent(in) :: lchnk ! chunk identifier - integer, intent(in) :: ncol ! number of atmospheric columns - integer, intent(in) :: ncnst ! number of tracers to transport - logical, intent(in) :: domomtran(ncnst) ! flag for doing convective transport - real(r8), intent(in) :: q(pcols,pver,ncnst) ! Wind array - real(r8), intent(in) :: mu(pcols,pver) ! Mass flux up - real(r8), intent(in) :: md(pcols,pver) ! Mass flux down - real(r8), intent(in) :: du(pcols,pver) ! Mass detraining from updraft - real(r8), intent(in) :: eu(pcols,pver) ! Mass entraining from updraft - real(r8), intent(in) :: ed(pcols,pver) ! Mass entraining from downdraft - real(r8), intent(in) :: dp(pcols,pver) ! Delta pressure between interfaces - real(r8), intent(in) :: dsubcld(pcols) ! Delta pressure from cloud base to sfc - real(r8), intent(in) :: dt ! time step in seconds : 2*delta_t - - integer, intent(in) :: jt(pcols) ! Index of cloud top for each column - integer, intent(in) :: mx(pcols) ! Index of cloud top for each column - integer, intent(in) :: ideep(pcols) ! Gathering array - integer, intent(in) :: il1g ! Gathered min lon indices over which to operate - integer, intent(in) :: il2g ! Gathered max lon indices over which to operate - integer, intent(in) :: nstep ! Time step index - - - -! input/output - - real(r8), intent(out) :: dqdt(pcols,pver,ncnst) ! Tracer tendency array - -!--------------------------Local Variables------------------------------ - - integer i ! Work index - integer k ! Work index - integer kbm ! Highest altitude index of cloud base - integer kk ! Work index - integer kkp1 ! Work index - integer kkm1 ! Work index - integer km1 ! Work index - integer kp1 ! Work index - integer ktm ! Highest altitude index of cloud top - integer m ! Work index - integer ii ! Work index - - real(r8) chat(pcols,pver) ! Mix ratio in env at interfaces - real(r8) cond(pcols,pver) ! Mix ratio in downdraft at interfaces - real(r8) const(pcols,pver) ! Gathered wind array - real(r8) conu(pcols,pver) ! Mix ratio in updraft at interfaces - real(r8) dcondt(pcols,pver) ! Gathered tend array - real(r8) small ! A small number - real(r8) mbsth ! Threshold for mass fluxes - real(r8) mupdudp ! A work variable - - real(r8) momcu ! constant for updraft pressure gradient term - real(r8) momcd ! constant for downdraft pressure gradient term - real(r8) sum ! sum - real(r8) sum2 ! sum2 - - real(r8) mududp(pcols,pver) ! working variable - real(r8) mddudp(pcols,pver) ! working variable - - real(r8) pgu(pcols,pver) ! Pressure gradient term for updraft - real(r8) pgd(pcols,pver) ! Pressure gradient term for downdraft - - real(r8),intent(out) :: pguall(pcols,pver,ncnst) ! Apparent force from updraft PG - real(r8),intent(out) :: pgdall(pcols,pver,ncnst) ! Apparent force from downdraft PG - - real(r8),intent(out) :: icwu(pcols,pver,ncnst) ! In-cloud winds in updraft - real(r8),intent(out) :: icwd(pcols,pver,ncnst) ! In-cloud winds in downdraft - - real(r8),intent(out) :: seten(pcols,pver) ! Dry static energy tendency - real(r8) gseten(pcols,pver) ! Gathered dry static energy tendency - - real(r8) mflux(pcols,pverp,ncnst) ! Gathered momentum flux - - real(r8) wind0(pcols,pver,ncnst) ! gathered wind before time step - real(r8) windf(pcols,pver,ncnst) ! gathered wind after time step - real(r8) fkeb, fket, ketend_cons, ketend, utop, ubot, vtop, vbot, gset2 - - -!----------------------------------------------------------------------- -! - -! Initialize outgoing fields - pguall(:,:,:) = 0.0_r8 - pgdall(:,:,:) = 0.0_r8 -! Initialize in-cloud winds to environmental wind - icwu(:ncol,:,:) = q(:ncol,:,:) - icwd(:ncol,:,:) = q(:ncol,:,:) - -! Initialize momentum flux and final winds - mflux(:,:,:) = 0.0_r8 - wind0(:,:,:) = 0.0_r8 - windf(:,:,:) = 0.0_r8 - -! Initialize dry static energy - - seten(:,:) = 0.0_r8 - gseten(:,:) = 0.0_r8 - -! Define constants for parameterization - - momcu = 0.4_r8 - momcd = 0.4_r8 - - small = 1.e-36_r8 -! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s) - mbsth = 1.e-15_r8 - -! Find the highest level top and bottom levels of convection - ktm = pver - kbm = pver - do i = il1g, il2g - ktm = min(ktm,jt(i)) - kbm = min(kbm,mx(i)) - end do - -! Loop ever each wind component - do m = 1, ncnst !start at m = 1 to transport momentum - if (domomtran(m)) then - -! Gather up the winds and set tend to zero - do k = 1,pver - do i =il1g,il2g - const(i,k) = q(ideep(i),k,m) - wind0(i,k,m) = const(i,k) - end do - end do - - -! From now on work only with gathered data - -! Interpolate winds to interfaces - - do k = 1,pver - km1 = max(1,k-1) - do i = il1g, il2g - - ! use arithmetic mean - chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1)) - -! Provisional up and down draft values - conu(i,k) = chat(i,k) - cond(i,k) = chat(i,k) - -! provisional tends - dcondt(i,k) = 0._r8 - - end do - end do - - -! -! Pressure Perturbation Term -! - - !Top boundary: assume mu is zero - - k=1 - pgu(:il2g,k) = 0.0_r8 - pgd(:il2g,k) = 0.0_r8 - - do k=2,pver-1 - km1 = max(1,k-1) - kp1 = min(pver,k+1) - do i = il1g,il2g - - !interior points - - mududp(i,k) = ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) & - + mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k)) - - pgu(i,k) = - momcu * 0.5_r8 * mududp(i,k) - - - mddudp(i,k) = ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) & - + md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k)) - - pgd(i,k) = - momcd * 0.5_r8 * mddudp(i,k) - - - end do - end do - - ! bottom boundary - k = pver - km1 = max(1,k-1) - do i=il1g,il2g - - mududp(i,k) = mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) - pgu(i,k) = - momcu * mududp(i,k) - - mddudp(i,k) = md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) - - pgd(i,k) = - momcd * mddudp(i,k) - - end do - - -! -! In-cloud velocity calculations -! - -! Do levels adjacent to top and bottom - k = 2 - km1 = 1 - kk = pver - kkm1 = max(1,kk-1) - do i = il1g,il2g - mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk) - if (mupdudp > mbsth) then - - conu(i,kk) = (+eu(i,kk)*const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp - endif - if (md(i,k) < -mbsth) then - cond(i,k) = (-ed(i,km1)*const(i,km1)*dp(i,km1))-pgd(i,km1)*dp(i,km1)/md(i,k) - endif - - - end do - - - -! Updraft from bottom to top - do kk = pver-1,1,-1 - kkm1 = max(1,kk-1) - kkp1 = min(pver,kk+1) - do i = il1g,il2g - mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk) - if (mupdudp > mbsth) then - - conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eu(i,kk)* & - const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp - endif - end do - - end do - - -! Downdraft from top to bottom - do k = 3,pver - km1 = max(1,k-1) - do i = il1g,il2g - if (md(i,k) < -mbsth) then - - cond(i,k) = ( md(i,km1)*cond(i,km1)-ed(i,km1)*const(i,km1) & - *dp(i,km1)-pgd(i,km1)*dp(i,km1) )/md(i,k) - - endif - end do - end do - - - sum = 0._r8 - sum2 = 0._r8 - - - do k = ktm,pver - km1 = max(1,k-1) - kp1 = min(pver,k+1) - do i = il1g,il2g - ii = ideep(i) - -! version 1 hard to check for roundoff errors - dcondt(i,k) = & - +(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) & - -mu(i,k)* (conu(i,k)-chat(i,k)) & - +md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) & - -md(i,k)* (cond(i,k)-chat(i,k)) & - )/dp(i,k) - - end do - end do - - ! dcont for bottom layer - ! - !DIR$ NOINTERCHANGE - do k = kbm,pver - km1 = max(1,k-1) - do i = il1g,il2g - if (k == mx(i)) then - - ! version 1 - dcondt(i,k) = (1._r8/dp(i,k))* & - (-mu(i,k)*(conu(i,k)-chat(i,k)) & - -md(i,k)*(cond(i,k)-chat(i,k)) & - ) - end if - end do - end do - -! Initialize to zero everywhere, then scatter tendency back to full array - dqdt(:,:,m) = 0._r8 - - do k = 1,pver - do i = il1g,il2g - ii = ideep(i) - dqdt(ii,k,m) = dcondt(i,k) - ! Output apparent force on the mean flow from pressure gradient - pguall(ii,k,m) = -pgu(i,k) - pgdall(ii,k,m) = -pgd(i,k) - icwu(ii,k,m) = conu(i,k) - icwd(ii,k,m) = cond(i,k) - end do - end do - - ! Calculate momentum flux in units of mb*m/s2 - - do k = ktm,pver - do i = il1g,il2g - ii = ideep(i) - mflux(i,k,m) = & - -mu(i,k)* (conu(i,k)-chat(i,k)) & - -md(i,k)* (cond(i,k)-chat(i,k)) - end do - end do - - - ! Calculate winds at the end of the time step - - do k = ktm,pver - do i = il1g,il2g - ii = ideep(i) - km1 = max(1,k-1) - kp1 = k+1 - windf(i,k,m) = const(i,k) - (mflux(i,kp1,m) - mflux(i,k,m)) * dt /dp(i,k) - - end do - end do - - end if ! for domomtran - end do - - ! Need to add an energy fix to account for the dissipation of kinetic energy - ! Formulation follows from Boville and Bretherton (2003) - ! formulation by PJR - - do k = ktm,pver - km1 = max(1,k-1) - kp1 = min(pver,k+1) - do i = il1g,il2g - - ii = ideep(i) - - ! calculate the KE fluxes at top and bot of layer - ! based on a discrete approximation to b&b eq(35) F_KE = u*F_u + v*F_v at interface - utop = (wind0(i,k,1)+wind0(i,km1,1))/2._r8 - vtop = (wind0(i,k,2)+wind0(i,km1,2))/2._r8 - ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2._r8 - vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2._r8 - fket = utop*mflux(i,k,1) + vtop*mflux(i,k,2) ! top of layer - fkeb = ubot*mflux(i,k+1,1) + vbot*mflux(i,k+1,2) ! bot of layer - - ! divergence of these fluxes should give a conservative redistribution of KE - ketend_cons = (fket-fkeb)/dp(i,k) - - ! tendency in kinetic energy resulting from the momentum transport - ketend = ((windf(i,k,1)**2 + windf(i,k,2)**2) - (wind0(i,k,1)**2 + wind0(i,k,2)**2))*0.5_r8/dt - - ! the difference should be the dissipation - gset2 = ketend_cons - ketend - gseten(i,k) = gset2 - - end do - - end do - - ! Scatter dry static energy to full array - do k = 1,pver - do i = il1g,il2g - ii = ideep(i) - seten(ii,k) = gseten(i,k) - - end do - end do - - return -end subroutine momtran -!endif -! AaronDonahue - skip convtran and momtran for now -!========================================================================================= - -subroutine buoyan(lchnk ,ncol , & - q ,t ,p ,z ,pf , & - tp ,qstp ,tl ,rl ,cape , & - pblt ,lcl ,lel ,lon ,mx , & - rd ,grav ,cp ,msg , & - tpert ) -!----------------------------------------------------------------------- -! -! Purpose: -! -! -! Method: -! -! -! -! Author: -! This is contributed code not fully standardized by the CCM core group. -! The documentation has been enhanced to the degree that we are able. -! Reviewed: P. Rasch, April 1996 -! -!----------------------------------------------------------------------- - implicit none -!----------------------------------------------------------------------- -! -! input arguments -! - integer, intent(in) :: lchnk ! chunk identifier - integer, intent(in) :: ncol ! number of atmospheric columns - - real(r8), intent(in) :: q(pcols,pver) ! spec. humidity - real(r8), intent(in) :: t(pcols,pver) ! temperature - real(r8), intent(in) :: p(pcols,pver) ! pressure - real(r8), intent(in) :: z(pcols,pver) ! height - real(r8), intent(in) :: pf(pcols,pver+1) ! pressure at interfaces - real(r8), intent(in) :: pblt(pcols) ! index of pbl depth - real(r8), intent(in) :: tpert(pcols) ! perturbation temperature by pbl processes - -! -! output arguments -! - real(r8), intent(out) :: tp(pcols,pver) ! parcel temperature - real(r8), intent(out) :: qstp(pcols,pver) ! saturation mixing ratio of parcel - real(r8), intent(out) :: tl(pcols) ! parcel temperature at lcl - real(r8), intent(out) :: cape(pcols) ! convective aval. pot. energy. - integer lcl(pcols) ! - integer lel(pcols) ! - integer lon(pcols) ! level of onset of deep convection - integer mx(pcols) ! level of max moist static energy -! -!--------------------------Local Variables------------------------------ -! - real(r8) capeten(pcols,5) ! provisional value of cape - real(r8) tv(pcols,pver) ! - real(r8) tpv(pcols,pver) ! - real(r8) buoy(pcols,pver) - - real(r8) a1(pcols) - real(r8) a2(pcols) - real(r8) estp(pcols) - real(r8) pl(pcols) - real(r8) plexp(pcols) - real(r8) hmax(pcols) - real(r8) hmn(pcols) - real(r8) y(pcols) - - logical plge600(pcols) - integer knt(pcols) - integer lelten(pcols,5) - - real(r8) cp - real(r8) e - real(r8) grav - - integer i - integer k - integer msg - integer n - - real(r8) rd - real(r8) rl -#ifdef PERGRO - real(r8) rhd -#endif -! -!----------------------------------------------------------------------- -! - do n = 1,5 - do i = 1,ncol - lelten(i,n) = pver - capeten(i,n) = 0._r8 - end do - end do -! - do i = 1,ncol - lon(i) = pver - knt(i) = 0 - lel(i) = pver - mx(i) = lon(i) - cape(i) = 0._r8 - hmax(i) = 0._r8 - end do - - tp(:ncol,:) = t(:ncol,:) - qstp(:ncol,:) = q(:ncol,:) - -!!! RBN - Initialize tv and buoy for output. -!!! tv=tv : tpv=tpv : qstp=q : buoy=0. - tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608_r8*q(:ncol,:))/ (1._r8+q(:ncol,:)) - tpv(:ncol,:) = tv(:ncol,:) - buoy(:ncol,:) = 0._r8 - -! -! set "launching" level(mx) to be at maximum moist static energy. -! search for this level stops at planetary boundary layer top. -! -#ifdef PERGRO - do k = pver,msg + 1,-1 - do i = 1,ncol - hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k) -! -! Reset max moist static energy level when relative difference exceeds 1.e-4 -! - rhd = (hmn(i) - hmax(i))/(hmn(i) + hmax(i)) - if (k >= nint(pblt(i)) .and. k <= lon(i) .and. rhd > -1.e-4_r8) then - hmax(i) = hmn(i) - mx(i) = k - end if - end do - end do -#else - do k = pver,msg + 1,-1 - do i = 1,ncol - hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k) - if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then - hmax(i) = hmn(i) - mx(i) = k - end if - end do - end do -#endif -! - do i = 1,ncol - lcl(i) = mx(i) - e = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i))) - tl(i) = 2840._r8/ (3.5_r8*log(t(i,mx(i)))-log(e)-4.805_r8) + 55._r8 - if (tl(i) < t(i,mx(i))) then - plexp(i) = (1._r8/ (0.2854_r8* (1._r8-0.28_r8*q(i,mx(i))))) - pl(i) = p(i,mx(i))* (tl(i)/t(i,mx(i)))**plexp(i) - else - tl(i) = t(i,mx(i)) - pl(i) = p(i,mx(i)) - end if - end do - -! -! calculate lifting condensation level (lcl). -! - do k = pver,msg + 2,-1 - do i = 1,ncol - if (k <= mx(i) .and. (p(i,k) > pl(i) .and. p(i,k-1) <= pl(i))) then - lcl(i) = k - 1 - end if - end do - end do -! -! if lcl is above the nominal level of non-divergence (600 mbs), -! no deep convection is permitted (ensuing calculations -! skipped and cape retains initialized value of zero). -! - do i = 1,ncol - plge600(i) = pl(i).ge.600._r8 - end do -! -! initialize parcel properties in sub-cloud layer below lcl. -! - do k = pver,msg + 1,-1 - do i=1,ncol - if (k > lcl(i) .and. k <= mx(i) .and. plge600(i)) then - tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k)) - qstp(i,k) = q(i,mx(i)) - tp(i,k) = t(i,mx(i))* (p(i,k)/p(i,mx(i)))**(0.2854_r8* (1._r8-0.28_r8*q(i,mx(i)))) -! -! buoyancy is increased by 0.5 k as in tiedtke -! -!-jjh tpv (i,k)=tp(i,k)*(1.+1.608*q(i,mx(i)))/ -!-jjh 1 (1.+q(i,mx(i))) - tpv(i,k) = (tp(i,k)+tpert(i))*(1._r8+1.608_r8*q(i,mx(i)))/ (1._r8+q(i,mx(i))) - buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add - end if - end do - end do - -! -! define parcel properties at lcl (i.e. level immediately above pl). -! - do k = pver,msg + 1,-1 - do i=1,ncol - if (k == lcl(i) .and. plge600(i)) then - tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k)) - qstp(i,k) = q(i,mx(i)) - tp(i,k) = tl(i)* (p(i,k)/pl(i))**(0.2854_r8* (1._r8-0.28_r8*qstp(i,k))) -! estp(i) =exp(21.656_r8 - 5418._r8/tp(i,k)) -! use of different formulas for es has about 1 g/kg difference -! in qs at t= 300k, and 0.02 g/kg at t=263k, with the formula -! above giving larger qs. - call qsat_hPa(tp(i,k), p(i,k), estp(i), qstp(i,k)) - a1(i) = cp / rl + qstp(i,k) * (1._r8+ qstp(i,k) / eps1) * rl * eps1 / & - (rd * tp(i,k) ** 2) - a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1*qstp(i,k))* & - (1._r8+qstp(i,k)/eps1)*eps1**2*rl*rl/ & - (rd**2*tp(i,k)**4)-qstp(i,k)* & - (1._r8+qstp(i,k)/eps1)*2._r8*eps1*rl/ & - (rd*tp(i,k)**3)) - a1(i) = 1._r8/a1(i) - a2(i) = -a2(i)*a1(i)**3 - y(i) = q(i,mx(i)) - qstp(i,k) - tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2 - call qsat_hPa(tp(i,k), p(i,k), estp(i), qstp(i,k)) -! -! buoyancy is increased by 0.5 k in cape calculation. -! dec. 9, 1994 -!-jjh tpv(i,k) =tp(i,k)*(1.+1.608*qstp(i,k))/(1.+q(i,mx(i))) -! - tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+q(i,mx(i))) - buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add - end if - end do - end do -! -! main buoyancy calculation. -! - do k = pver - 1,msg + 1,-1 - do i=1,ncol - if (k < lcl(i) .and. plge600(i)) then - tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k)) - qstp(i,k) = qstp(i,k+1) - tp(i,k) = tp(i,k+1)* (p(i,k)/p(i,k+1))**(0.2854_r8* (1._r8-0.28_r8*qstp(i,k))) - call qsat_hPa(tp(i,k), p(i,k), estp(i), qstp(i,k)) - a1(i) = cp/rl + qstp(i,k)* (1._r8+qstp(i,k)/eps1)*rl*eps1/ (rd*tp(i,k)**2) - a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1*qstp(i,k))* & - (1._r8+qstp(i,k)/eps1)*eps1**2*rl*rl/ & - (rd**2*tp(i,k)**4)-qstp(i,k)* & - (1._r8+qstp(i,k)/eps1)*2._r8*eps1*rl/ & - (rd*tp(i,k)**3)) - a1(i) = 1._r8/a1(i) - a2(i) = -a2(i)*a1(i)**3 - y(i) = qstp(i,k+1) - qstp(i,k) - tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2 - call qsat_hPa(tp(i,k), p(i,k), estp(i), qstp(i,k)) -!-jjh tpv(i,k) =tp(i,k)*(1.+1.608*qstp(i,k))/ -!jt (1.+q(i,mx(i))) - tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k))/(1._r8+q(i,mx(i))) - buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add - end if - end do - end do - -! - do k = msg + 2,pver - do i = 1,ncol - if (k < lcl(i) .and. plge600(i)) then - if (buoy(i,k+1) > 0._r8 .and. buoy(i,k) <= 0._r8) then - knt(i) = min(5,knt(i) + 1) - lelten(i,knt(i)) = k - end if - end if - end do - end do -! -! calculate convective available potential energy (cape). -! - do n = 1,5 - do k = msg + 1,pver - do i = 1,ncol - if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then - capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k)) - end if - end do - end do - end do -! -! find maximum cape from all possible tentative capes from -! one sounding, -! and use it as the final cape, april 26, 1995 -! - do n = 1,5 - do i = 1,ncol - if (capeten(i,n) > cape(i)) then - cape(i) = capeten(i,n) - lel(i) = lelten(i,n) - end if - end do - end do -! -! put lower bound on cape for diagnostic purposes. -! - do i = 1,ncol - cape(i) = max(cape(i), 0._r8) - end do -! - return -end subroutine buoyan - -subroutine cldprp(lchnk , & - q ,t ,u ,v ,p , & - z ,s ,mu ,eu ,du , & - md ,ed ,sd ,qd ,mc , & - qu ,su ,zf ,qst ,hmn , & - hsat ,shat ,ql , & - cmeg ,jb ,lel ,jt ,jlcl , & - mx ,j0 ,jd ,rl ,il2g , & - rd ,grav ,cp ,msg , & -!songxl 2014-05-20------- -!----------------------------------------------------------------------- -! -! Purpose: -! -! -! Method: -! may 09/91 - guang jun zhang, m.lazare, n.mcfarlane. -! original version cldprop. -! -! Author: See above, modified by P. Rasch -! This is contributed code not fully standardized by the CCM core group. -! -! this code is very much rougher than virtually anything else in the CCM -! there are debug statements left strewn about and code segments disabled -! these are to facilitate future development. We expect to release a -! cleaner code in a future release -! -! the documentation has been enhanced to the degree that we are able -! -!----------------------------------------------------------------------- - - implicit none - -!------------------------------------------------------------------------------ -! -! Input arguments -! - integer, intent(in) :: lchnk ! chunk identifier - - real(r8), intent(in) :: q(pcols,pver) ! spec. humidity of env - real(r8), intent(in) :: t(pcols,pver) ! temp of env - real(r8), intent(in) :: p(pcols,pver) ! pressure of env - real(r8), intent(in) :: z(pcols,pver) ! height of env - real(r8), intent(in) :: s(pcols,pver) ! normalized dry static energy of env - real(r8), intent(in) :: zf(pcols,pverp) ! height of interfaces - real(r8), intent(in) :: u(pcols,pver) ! zonal velocity of env - real(r8), intent(in) :: v(pcols,pver) ! merid. velocity of env - - real(r8), intent(in) :: landfrac(pcols) ! RBN Landfrac - - integer, intent(in) :: jb(pcols) ! updraft base level - integer, intent(in) :: lel(pcols) ! updraft launch level - integer, intent(out) :: jt(pcols) ! updraft plume top - integer, intent(out) :: jlcl(pcols) ! updraft lifting cond level - integer, intent(in) :: mx(pcols) ! updraft base level (same is jb) - integer, intent(out) :: j0(pcols) ! level where updraft begins detraining - integer, intent(out) :: jd(pcols) ! level of downdraft - integer, intent(in) :: limcnv ! convection limiting level - integer, intent(in) :: il2g !CORE GROUP REMOVE - integer, intent(in) :: msg ! missing moisture vals (always 0) - real(r8), intent(in) :: rl ! latent heat of vap - real(r8), intent(in) :: shat(pcols,pver) ! interface values of dry stat energy - real(r8), intent(in) :: tpert(pcols) -! -! output -! - real(r8), intent(out) :: rprd(pcols,pver) ! rate of production of precip at that layer - real(r8), intent(out) :: du(pcols,pver) ! detrainement rate of updraft - real(r8), intent(out) :: ed(pcols,pver) ! entrainment rate of downdraft - real(r8), intent(out) :: eu(pcols,pver) ! entrainment rate of updraft - real(r8), intent(out) :: hmn(pcols,pver) ! moist stat energy of env - real(r8), intent(out) :: hsat(pcols,pver) ! sat moist stat energy of env - real(r8), intent(out) :: mc(pcols,pver) ! net mass flux - real(r8), intent(out) :: md(pcols,pver) ! downdraft mass flux - real(r8), intent(out) :: mu(pcols,pver) ! updraft mass flux - real(r8), intent(out) :: pflx(pcols,pverp) ! precipitation flux thru layer - real(r8), intent(out) :: qd(pcols,pver) ! spec humidity of downdraft - real(r8), intent(out) :: ql(pcols,pver) ! liq water of updraft - real(r8), intent(out) :: qst(pcols,pver) ! saturation mixing ratio of env. - real(r8), intent(out) :: qu(pcols,pver) ! spec hum of updraft - real(r8), intent(out) :: sd(pcols,pver) ! normalized dry stat energy of downdraft - real(r8), intent(out) :: su(pcols,pver) ! normalized dry stat energy of updraft - -!songxl 2014-05-20---------------- - - real(r8) rd ! gas constant for dry air - real(r8) grav ! gravity - real(r8) cp ! heat capacity of dry air - -! -! Local workspace -! - real(r8) gamma(pcols,pver) - real(r8) dz(pcols,pver) - real(r8) iprm(pcols,pver) - real(r8) hu(pcols,pver) - real(r8) hd(pcols,pver) - real(r8) eps(pcols,pver) - real(r8) f(pcols,pver) - real(r8) k1(pcols,pver) - real(r8) i2(pcols,pver) - real(r8) ihat(pcols,pver) - real(r8) i3(pcols,pver) - real(r8) idag(pcols,pver) - real(r8) i4(pcols,pver) - real(r8) qsthat(pcols,pver) - real(r8) hsthat(pcols,pver) - real(r8) gamhat(pcols,pver) - real(r8) cu(pcols,pver) - real(r8) evp(pcols,pver) - real(r8) cmeg(pcols,pver) - real(r8) qds(pcols,pver) -! RBN For c0mask - real(r8) c0mask(pcols) - - real(r8) hmin(pcols) - real(r8) expdif(pcols) - real(r8) expnum(pcols) - real(r8) ftemp(pcols) - real(r8) eps0(pcols) - real(r8) rmue(pcols) - real(r8) zuef(pcols) - real(r8) zdef(pcols) - real(r8) epsm(pcols) - real(r8) ratmjb(pcols) - real(r8) est(pcols) - real(r8) totpcp(pcols) - real(r8) totevp(pcols) - real(r8) alfa(pcols) - real(r8) ql1 - real(r8) tu - real(r8) estu - real(r8) qstu - - real(r8) small - real(r8) mdt - - integer khighest - integer klowest - integer kount - integer i,k - - logical doit(pcols) - logical done(pcols) - -!songxl 2014-05-20---------------- -! -!------------------------------------------------------------------------------ -! - do i = 1,il2g - ftemp(i) = 0._r8 - expnum(i) = 0._r8 - expdif(i) = 0._r8 - c0mask(i) = c0_ocn * (1._r8-landfrac(i)) + c0_lnd * landfrac(i) - end do -! -!jr Change from msg+1 to 1 to prevent blowup -! - do k = 1,pver - do i = 1,il2g - dz(i,k) = zf(i,k) - zf(i,k+1) - end do - end do - -! -! initialize many output and work variables to zero -! - pflx(:il2g,1) = 0 - - do k = 1,pver - do i = 1,il2g - k1(i,k) = 0._r8 - i2(i,k) = 0._r8 - i3(i,k) = 0._r8 - i4(i,k) = 0._r8 - mu(i,k) = 0._r8 - f(i,k) = 0._r8 - eps(i,k) = 0._r8 - eu(i,k) = 0._r8 - du(i,k) = 0._r8 - ql(i,k) = 0._r8 - cu(i,k) = 0._r8 - evp(i,k) = 0._r8 - cmeg(i,k) = 0._r8 - qds(i,k) = q(i,k) - md(i,k) = 0._r8 - ed(i,k) = 0._r8 - sd(i,k) = s(i,k) - qd(i,k) = q(i,k) - mc(i,k) = 0._r8 - qu(i,k) = q(i,k) - su(i,k) = s(i,k) - call qsat_hPa(t(i,k), p(i,k), est(i), qst(i,k)) -!++bee - if ( p(i,k)-est(i) <= 0._r8 ) then - qst(i,k) = 1.0_r8 - end if -!--bee - gamma(i,k) = qst(i,k)*(1._r8 + qst(i,k)/eps1)*eps1*rl/(rd*t(i,k)**2)*rl/cp - hmn(i,k) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k) - hsat(i,k) = cp*t(i,k) + grav*z(i,k) + rl*qst(i,k) - hu(i,k) = hmn(i,k) - hd(i,k) = hmn(i,k) - rprd(i,k) = 0._r8 - end do - end do -! -!jr Set to zero things which make this routine blow up -! - do k=1,msg - do i=1,il2g - rprd(i,k) = 0._r8 - end do - end do -! -! interpolate the layer values of qst, hsat and gamma to -! layer interfaces -! - do k = 1, msg+1 - do i = 1,il2g - hsthat(i,k) = hsat(i,k) - qsthat(i,k) = qst(i,k) - gamhat(i,k) = gamma(i,k) - end do - end do - do i = 1,il2g - totpcp(i) = 0._r8 - totevp(i) = 0._r8 - end do - do k = msg + 2,pver - do i = 1,il2g - if (abs(qst(i,k-1)-qst(i,k)) > 1.E-6_r8) then - qsthat(i,k) = log(qst(i,k-1)/qst(i,k))*qst(i,k-1)*qst(i,k)/ (qst(i,k-1)-qst(i,k)) - else - qsthat(i,k) = qst(i,k) - end if - hsthat(i,k) = cp*shat(i,k) + rl*qsthat(i,k) - if (abs(gamma(i,k-1)-gamma(i,k)) > 1.E-6_r8) then - gamhat(i,k) = log(gamma(i,k-1)/gamma(i,k))*gamma(i,k-1)*gamma(i,k)/ & - (gamma(i,k-1)-gamma(i,k)) - else - gamhat(i,k) = gamma(i,k) - end if - end do - end do -! -! initialize cloud top to highest plume top. -!jr changed hard-wired 4 to limcnv+1 (not to exceed pver) -! - jt(:) = pver - do i = 1,il2g - jt(i) = max(lel(i),limcnv+1) - jt(i) = min(jt(i),pver) - jd(i) = pver - jlcl(i) = lel(i) - hmin(i) = 1.E6_r8 - end do -! -! find the level of minimum hsat, where detrainment starts -! - - do k = msg + 1,pver - do i = 1,il2g - if (hsat(i,k) <= hmin(i) .and. k >= jt(i) .and. k <= jb(i)) then - hmin(i) = hsat(i,k) - j0(i) = k - end if - end do - end do - do i = 1,il2g - j0(i) = min(j0(i),jb(i)-2) - j0(i) = max(j0(i),jt(i)+2) -! -! Fix from Guang Zhang to address out of bounds array reference -! - j0(i) = min(j0(i),pver) - end do -! -! Initialize certain arrays inside cloud -! - do k = msg + 1,pver - do i = 1,il2g - if (k >= jt(i) .and. k <= jb(i)) then - hu(i,k) = hmn(i,mx(i)) + cp*(tiedke_add+tp_fac*tpert(i)) !PMA - su(i,k) = s(i,mx(i)) + tiedke_add+tp_fac*tpert(i) - end if - end do - end do -! -! ********************************************************* -! compute taylor series for approximate eps(z) below -! ********************************************************* -! - do k = pver - 1,msg + 1,-1 - do i = 1,il2g - if (k < jb(i) .and. k >= jt(i)) then - k1(i,k) = k1(i,k+1) + (hmn(i,mx(i))-hmn(i,k))*dz(i,k) - ihat(i,k) = 0.5_r8* (k1(i,k+1)+k1(i,k)) - i2(i,k) = i2(i,k+1) + ihat(i,k)*dz(i,k) - idag(i,k) = 0.5_r8* (i2(i,k+1)+i2(i,k)) - i3(i,k) = i3(i,k+1) + idag(i,k)*dz(i,k) - iprm(i,k) = 0.5_r8* (i3(i,k+1)+i3(i,k)) - i4(i,k) = i4(i,k+1) + iprm(i,k)*dz(i,k) - end if - end do - end do -! -! re-initialize hmin array for ensuing calculation. -! - do i = 1,il2g - hmin(i) = 1.E6_r8 - end do - do k = msg + 1,pver - do i = 1,il2g - if (k >= j0(i) .and. k <= jb(i) .and. hmn(i,k) <= hmin(i)) then - hmin(i) = hmn(i,k) - expdif(i) = hmn(i,mx(i)) - hmin(i) - end if - end do - end do -! -! ********************************************************* -! compute approximate eps(z) using above taylor series -! ********************************************************* -! - do k = msg + 2,pver - do i = 1,il2g - expnum(i) = 0._r8 - ftemp(i) = 0._r8 - if (k < jt(i) .or. k >= jb(i)) then - k1(i,k) = 0._r8 - expnum(i) = 0._r8 - else - expnum(i) = hmn(i,mx(i)) - (hsat(i,k-1)*(zf(i,k)-z(i,k)) + & - hsat(i,k)* (z(i,k-1)-zf(i,k)))/(z(i,k-1)-z(i,k)) - end if - if ((expdif(i) > 100._r8 .and. expnum(i) > 0._r8) .and. & - k1(i,k) > expnum(i)*dz(i,k)) then - ftemp(i) = expnum(i)/k1(i,k) - f(i,k) = ftemp(i) + i2(i,k)/k1(i,k)*ftemp(i)**2 + & - (2._r8*i2(i,k)**2-k1(i,k)*i3(i,k))/k1(i,k)**2* & - ftemp(i)**3 + (-5._r8*k1(i,k)*i2(i,k)*i3(i,k)+ & - 5._r8*i2(i,k)**3+k1(i,k)**2*i4(i,k))/ & - k1(i,k)**3*ftemp(i)**4 - f(i,k) = max(f(i,k),0._r8) - f(i,k) = min(f(i,k),0.0002_r8) - end if - end do - end do - do i = 1,il2g - if (j0(i) < jb(i)) then - if (f(i,j0(i)) < 1.E-6_r8 .and. f(i,j0(i)+1) > f(i,j0(i))) j0(i) = j0(i) + 1 - end if - end do - do k = msg + 2,pver - do i = 1,il2g - if (k >= jt(i) .and. k <= j0(i)) then - f(i,k) = max(f(i,k),f(i,k-1)) - end if - end do - end do - do i = 1,il2g - eps0(i) = f(i,j0(i)) - eps(i,jb(i)) = eps0(i) - end do -! -! This is set to match the Rasch and Kristjansson paper -! - do k = pver,msg + 1,-1 - do i = 1,il2g - if (k >= j0(i) .and. k <= jb(i)) then - eps(i,k) = f(i,j0(i)) - end if - end do - end do - do k = pver,msg + 1,-1 - do i = 1,il2g - if (k < j0(i) .and. k >= jt(i)) eps(i,k) = f(i,k) - end do - end do -! -! specify the updraft mass flux mu, entrainment eu, detrainment du -! and moist static energy hu. -! here and below mu, eu,du, md and ed are all normalized by mb -! - do i = 1,il2g - if (eps0(i) > 0._r8) then - mu(i,jb(i)) = 1._r8 - eu(i,jb(i)) = mu(i,jb(i))/dz(i,jb(i)) - end if - end do - do k = pver,msg + 1,-1 - do i = 1,il2g - if (eps0(i) > 0._r8 .and. (k >= jt(i) .and. k < jb(i))) then - zuef(i) = zf(i,k) - zf(i,jb(i)) - rmue(i) = (1._r8/eps0(i))* (exp(eps(i,k+1)*zuef(i))-1._r8)/zuef(i) - mu(i,k) = (1._r8/eps0(i))* (exp(eps(i,k )*zuef(i))-1._r8)/zuef(i) - eu(i,k) = (rmue(i)-mu(i,k+1))/dz(i,k) - du(i,k) = (rmue(i)-mu(i,k))/dz(i,k) - end if - end do - end do -! - khighest = pverp - klowest = 1 - do i=1,il2g - khighest = min(khighest,lel(i)) - klowest = max(klowest,jb(i)) - end do - -!songxl 2014-05-20-------------- - - do k = klowest-1,khighest,-1 - do i = 1,il2g - if (k <= jb(i)-1 .and. k >= lel(i) .and. eps0(i) > 0._r8) then - if (mu(i,k) < 0.02_r8) then - hu(i,k) = hmn(i,k) - mu(i,k) = 0._r8 - eu(i,k) = 0._r8 - du(i,k) = mu(i,k+1)/dz(i,k) - else -!songxl 2014-05-20-------------- - end if - end if - end do - end do -! -! -! reset cloud top index beginning from two layers above the -! cloud base (i.e. if cloud is only one layer thick, top is not reset -! - do i=1,il2g - doit(i) = .true. - end do - do k=klowest-2,khighest-1,-1 - do i=1,il2g - if (doit(i) .and. k <= jb(i)-2 .and. k >= lel(i)-1) then - if(trigmem)then - if (hu(i,k) <= hsthat(i,k) .and. hu(i,k+1) > hsthat(i,k+1) & - .and. mu(i,k) >= 0.02_r8) then - if (hu(i,k)-hsthat(i,k) < -2000._r8) then - jt(i) = k + 1 - doit(i) = .false. - else - jt(i) = k - doit(i) = .false. - end if - else - if (hu_tot(i).eq.hmn_tot(i)) then - if (hu(i,k) > hu(i,jb(i)) .or. mu(i,k) < 0.02_r8) then - jt(i) = k + 1 - doit(i) = .false. - end if - else - if ( mu(i,k) < 0.02_r8) then - jt(i) = k + 1 - doit(i) = .false. - end if - end if - end if - else ! not trigmem - if (hu(i,k) <= hsthat(i,k) .and. hu(i,k+1) > hsthat(i,k+1) & - .and. mu(i,k) >= 0.02_r8) then - if (hu(i,k)-hsthat(i,k) < -2000._r8) then - jt(i) = k + 1 - doit(i) = .false. - else - jt(i) = k - doit(i) = .false. - end if - else if (hu(i,k) > hu(i,jb(i)) .or. mu(i,k) < 0.02_r8) then - jt(i) = k + 1 - doit(i) = .false. - end if - end if ! end trigmem -!>songxl 2014-06-20------------------------ - end if - end do - end do - do k = pver,msg + 1,-1 - do i = 1,il2g - if (k >= lel(i) .and. k <= jt(i) .and. eps0(i) > 0._r8) then - mu(i,k) = 0._r8 - eu(i,k) = 0._r8 - du(i,k) = 0._r8 - hu(i,k) = hmn(i,k) - end if - if (k == jt(i) .and. eps0(i) > 0._r8) then - du(i,k) = mu(i,k+1)/dz(i,k) - eu(i,k) = 0._r8 - mu(i,k) = 0._r8 - end if - end do - end do -! -!songxl 2014-05-20----------------- - -! specify downdraft properties (no downdrafts if jd.ge.jb). -! scale down downward mass flux profile so that net flux -! (up-down) at cloud base in not negative. -! - do i = 1,il2g -! -! in normal downdraft strength run alfa=0.2. In test4 alfa=0.1 -! - alfa(i) = alfa_scalar - jt(i) = min(jt(i),jb(i)-1) - jd(i) = max(j0(i),jt(i)+1) - jd(i) = min(jd(i),jb(i)) - hd(i,jd(i)) = hmn(i,jd(i)-1) - if (jd(i) < jb(i) .and. eps0(i) > 0._r8) then - epsm(i) = eps0(i) - md(i,jd(i)) = -alfa(i)*epsm(i)/eps0(i) - end if - end do - do k = msg + 1,pver - do i = 1,il2g - if ((k > jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8) then - zdef(i) = zf(i,jd(i)) - zf(i,k) - md(i,k) = -alfa(i)/ (2._r8*eps0(i))*(exp(2._r8*epsm(i)*zdef(i))-1._r8)/zdef(i) - end if - end do - end do - do k = msg + 1,pver - do i = 1,il2g - if ((k >= jt(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8 .and. jd(i) < jb(i)) then - ratmjb(i) = min(abs(mu(i,jb(i))/md(i,jb(i))),1._r8) - md(i,k) = md(i,k)*ratmjb(i) - end if - end do - end do - - small = 1.e-20_r8 - do k = msg + 1,pver - do i = 1,il2g - if ((k >= jt(i) .and. k <= pver) .and. eps0(i) > 0._r8) then - ed(i,k-1) = (md(i,k-1)-md(i,k))/dz(i,k-1) - mdt = min(md(i,k),-small) - hd(i,k) = (md(i,k-1)*hd(i,k-1) - dz(i,k-1)*ed(i,k-1)*hmn(i,k-1))/mdt - end if - end do - end do -! -! calculate updraft and downdraft properties. -! - do k = msg + 2,pver - do i = 1,il2g - if ((k >= jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8 .and. jd(i) < jb(i)) then - qds(i,k) = qsthat(i,k) + gamhat(i,k)*(hd(i,k)-hsthat(i,k))/ & - (rl*(1._r8 + gamhat(i,k))) - end if - end do - end do -! - do i = 1,il2g - done(i) = .false. - end do - kount = 0 - do k = pver,msg + 2,-1 - do i = 1,il2g - if (k == jb(i) .and. eps0(i) > 0._r8) then - qu(i,k) = q(i,mx(i)) - su(i,k) = (hu(i,k)-rl*qu(i,k))/cp - end if - if (( .not. done(i) .and. k > jt(i) .and. k < jb(i)) .and. eps0(i) > 0._r8) then - su(i,k) = mu(i,k+1)/mu(i,k)*su(i,k+1) + & - dz(i,k)/mu(i,k)* (eu(i,k)-du(i,k))*s(i,k) - qu(i,k) = mu(i,k+1)/mu(i,k)*qu(i,k+1) + dz(i,k)/mu(i,k)* (eu(i,k)*q(i,k)- & - du(i,k)*qst(i,k)) - tu = su(i,k) - grav/cp*zf(i,k) - call qsat_hPa(tu, (p(i,k)+p(i,k-1))/2._r8, estu, qstu) - if (qu(i,k) >= qstu) then - jlcl(i) = k - kount = kount + 1 - done(i) = .true. - end if - end if - end do - if (kount >= il2g) goto 690 - end do -690 continue - do k = msg + 2,pver - do i = 1,il2g - if ((k > jt(i) .and. k <= jlcl(i)) .and. eps0(i) > 0._r8) then - su(i,k) = shat(i,k) + (hu(i,k)-hsthat(i,k))/(cp* (1._r8+gamhat(i,k))) - qu(i,k) = qsthat(i,k) + gamhat(i,k)*(hu(i,k)-hsthat(i,k))/ & - (rl* (1._r8+gamhat(i,k))) - end if - end do - end do - -! compute condensation in updraft - do k = pver,msg + 2,-1 - do i = 1,il2g - if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8) then - cu(i,k) = ((mu(i,k)*su(i,k)-mu(i,k+1)*su(i,k+1))/ & - dz(i,k)- (eu(i,k)-du(i,k))*s(i,k))/(rl/cp) - if (k == jt(i)) cu(i,k) = 0._r8 - cu(i,k) = max(0._r8,cu(i,k)) - end if - end do - end do - -! compute condensed liquid, rain production rate -! accumulate total precipitation (condensation - detrainment of liquid) -! Note ql1 = ql(k) + rprd(k)*dz(k)/mu(k) -! The differencing is somewhat strange (e.g. du(i,k)*ql(i,k+1)) but is -! consistently applied. -! mu, ql are interface quantities -! cu, du, eu, rprd are midpoint quantites - do k = pver,msg + 2,-1 - do i = 1,il2g - rprd(i,k) = 0._r8 - if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8 .and. mu(i,k) >= 0.0_r8) then - if (mu(i,k) > 0._r8) then - ql1 = 1._r8/mu(i,k)* (mu(i,k+1)*ql(i,k+1)- & - dz(i,k)*du(i,k)*ql(i,k+1)+dz(i,k)*cu(i,k)) - ql(i,k) = ql1/ (1._r8+dz(i,k)*c0mask(i)) - else - ql(i,k) = 0._r8 - end if - totpcp(i) = totpcp(i) + dz(i,k)*(cu(i,k)-du(i,k)*ql(i,k+1)) - rprd(i,k) = c0mask(i)*mu(i,k)*ql(i,k) - end if - end do - end do -! - do i = 1,il2g - qd(i,jd(i)) = qds(i,jd(i)) - sd(i,jd(i)) = (hd(i,jd(i)) - rl*qd(i,jd(i)))/cp - end do -! - do k = msg + 2,pver - do i = 1,il2g - if (k >= jd(i) .and. k < jb(i) .and. eps0(i) > 0._r8) then - qd(i,k+1) = qds(i,k+1) - evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k)-md(i,k+1)*qd(i,k+1))/dz(i,k) - evp(i,k) = max(evp(i,k),0._r8) - mdt = min(md(i,k+1),-small) - sd(i,k+1) = ((rl/cp*evp(i,k)-ed(i,k)*s(i,k))*dz(i,k) + md(i,k)*sd(i,k))/mdt - totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k) - end if - end do - end do - do i = 1,il2g -!*guang totevp(i) = totevp(i) + md(i,jd(i))*q(i,jd(i)-1) - - totevp(i) = totevp(i) + md(i,jd(i))*qd(i,jd(i)) - md(i,jb(i))*qd(i,jb(i)) - end do -!!$ if (.true.) then - if (.false.) then - do i = 1,il2g - k = jb(i) - if (eps0(i) > 0._r8) then - evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k))/dz(i,k) - evp(i,k) = max(evp(i,k),0._r8) - totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k) - end if - end do - endif - - do i = 1,il2g - totpcp(i) = max(totpcp(i),0._r8) - totevp(i) = max(totevp(i),0._r8) - end do -! - do k = msg + 2,pver - do i = 1,il2g - if (totevp(i) > 0._r8 .and. totpcp(i) > 0._r8) then - md(i,k) = md (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i))) - ed(i,k) = ed (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i))) - evp(i,k) = evp(i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i))) - else - md(i,k) = 0._r8 - ed(i,k) = 0._r8 - evp(i,k) = 0._r8 - end if -! cmeg is the cloud water condensed - rain water evaporated -! rprd is the cloud water converted to rain - (rain evaporated) - cmeg(i,k) = cu(i,k) - evp(i,k) - rprd(i,k) = rprd(i,k)-evp(i,k) - end do - end do - -! compute the net precipitation flux across interfaces - pflx(:il2g,1) = 0._r8 - do k = 2,pverp - do i = 1,il2g - pflx(i,k) = pflx(i,k-1) + rprd(i,k-1)*dz(i,k-1) - end do - end do -! - do k = msg + 1,pver - do i = 1,il2g - mc(i,k) = mu(i,k) + md(i,k) - end do - end do -! - return -end subroutine cldprp - -subroutine closure(lchnk , & - q ,t ,p ,z ,s , & - tp ,qs ,qu ,su ,mc , & - du ,mu ,md ,qd ,sd , & - qhat ,shat ,dp ,qstp ,zf , & - ql ,dsubcld ,mb ,cape ,tl , & - lcl ,lel ,jt ,mx ,il1g , & - il2g ,rd ,grav ,cp ,rl , & - msg ,capelmt ) -!----------------------------------------------------------------------- -! -! Purpose: -! -! -! Method: -! -! -! -! Author: G. Zhang and collaborators. CCM contact:P. Rasch -! This is contributed code not fully standardized by the CCM core group. -! -! this code is very much rougher than virtually anything else in the CCM -! We expect to release cleaner code in a future release -! -! the documentation has been enhanced to the degree that we are able -! -!----------------------------------------------------------------------- - - implicit none - -! -!-----------------------------Arguments--------------------------------- -! - integer, intent(in) :: lchnk ! chunk identifier - - real(r8), intent(inout) :: q(pcols,pver) ! spec humidity - real(r8), intent(inout) :: t(pcols,pver) ! temperature - real(r8), intent(inout) :: p(pcols,pver) ! pressure (mb) - real(r8), intent(inout) :: mb(pcols) ! cloud base mass flux - real(r8), intent(in) :: z(pcols,pver) ! height (m) - real(r8), intent(in) :: s(pcols,pver) ! normalized dry static energy - real(r8), intent(in) :: tp(pcols,pver) ! parcel temp - real(r8), intent(in) :: qs(pcols,pver) ! sat spec humidity - real(r8), intent(in) :: qu(pcols,pver) ! updraft spec. humidity - real(r8), intent(in) :: su(pcols,pver) ! normalized dry stat energy of updraft - real(r8), intent(in) :: mc(pcols,pver) ! net convective mass flux - real(r8), intent(in) :: du(pcols,pver) ! detrainment from updraft - real(r8), intent(in) :: mu(pcols,pver) ! mass flux of updraft - real(r8), intent(in) :: md(pcols,pver) ! mass flux of downdraft - real(r8), intent(in) :: qd(pcols,pver) ! spec. humidity of downdraft - real(r8), intent(in) :: sd(pcols,pver) ! dry static energy of downdraft - real(r8), intent(in) :: qhat(pcols,pver) ! environment spec humidity at interfaces - real(r8), intent(in) :: shat(pcols,pver) ! env. normalized dry static energy at intrfcs - real(r8), intent(in) :: dp(pcols,pver) ! pressure thickness of layers - real(r8), intent(in) :: qstp(pcols,pver) ! spec humidity of parcel - real(r8), intent(in) :: zf(pcols,pver+1) ! height of interface levels - real(r8), intent(in) :: ql(pcols,pver) ! liquid water mixing ratio - - real(r8), intent(in) :: cape(pcols) ! available pot. energy of column - real(r8), intent(in) :: tl(pcols) - real(r8), intent(in) :: dsubcld(pcols) ! thickness of subcloud layer - - integer, intent(in) :: lcl(pcols) ! index of lcl - integer, intent(in) :: lel(pcols) ! index of launch leve - integer, intent(in) :: jt(pcols) ! top of updraft - integer, intent(in) :: mx(pcols) ! base of updraft -! -!--------------------------Local variables------------------------------ -! - real(r8) dtpdt(pcols,pver) - real(r8) dqsdtp(pcols,pver) - real(r8) dtmdt(pcols,pver) - real(r8) dqmdt(pcols,pver) - real(r8) dboydt(pcols,pver) - real(r8) thetavp(pcols,pver) - real(r8) thetavm(pcols,pver) - - real(r8) dtbdt(pcols),dqbdt(pcols),dtldt(pcols) - real(r8) beta - real(r8) capelmt - real(r8) cp - real(r8) dadt(pcols) - real(r8) debdt - real(r8) dltaa - real(r8) eb - real(r8) grav - - integer i - integer il1g - integer il2g - integer k, kmin, kmax - integer msg - - real(r8) rd - real(r8) rl -! change of subcloud layer properties due to convection is -! related to cumulus updrafts and downdrafts. -! mc(z)=f(z)*mb, mub=betau*mb, mdb=betad*mb are used -! to define betau, betad and f(z). -! note that this implies all time derivatives are in effect -! time derivatives per unit cloud-base mass flux, i.e. they -! have units of 1/mb instead of 1/sec. -! - do i = il1g,il2g - mb(i) = 0._r8 - eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i))) - dtbdt(i) = (1._r8/dsubcld(i))* (mu(i,mx(i))*(shat(i,mx(i))-su(i,mx(i)))+ & - md(i,mx(i))* (shat(i,mx(i))-sd(i,mx(i)))) - dqbdt(i) = (1._r8/dsubcld(i))* (mu(i,mx(i))*(qhat(i,mx(i))-qu(i,mx(i)))+ & - md(i,mx(i))* (qhat(i,mx(i))-qd(i,mx(i)))) - debdt = eps1*p(i,mx(i))/ (eps1+q(i,mx(i)))**2*dqbdt(i) - dtldt(i) = -2840._r8* (3.5_r8/t(i,mx(i))*dtbdt(i)-debdt/eb)/ & - (3.5_r8*log(t(i,mx(i)))-log(eb)-4.805_r8)**2 - end do -! -! dtmdt and dqmdt are cumulus heating and drying. -! - do k = msg + 1,pver - do i = il1g,il2g - dtmdt(i,k) = 0._r8 - dqmdt(i,k) = 0._r8 - end do - end do -! - do k = msg + 1,pver - 1 - do i = il1g,il2g - if (k == jt(i)) then - dtmdt(i,k) = (1._r8/dp(i,k))*(mu(i,k+1)* (su(i,k+1)-shat(i,k+1)- & - rl/cp*ql(i,k+1))+md(i,k+1)* (sd(i,k+1)-shat(i,k+1))) - dqmdt(i,k) = (1._r8/dp(i,k))*(mu(i,k+1)* (qu(i,k+1)- & - qhat(i,k+1)+ql(i,k+1))+md(i,k+1)*(qd(i,k+1)-qhat(i,k+1))) - end if - end do - end do -! - beta = 0._r8 - do k = msg + 1,pver - 1 - do i = il1g,il2g - if (k > jt(i) .and. k < mx(i)) then - dtmdt(i,k) = (mc(i,k)* (shat(i,k)-s(i,k))+mc(i,k+1)* (s(i,k)-shat(i,k+1)))/ & - dp(i,k) - rl/cp*du(i,k)*(beta*ql(i,k)+ (1-beta)*ql(i,k+1)) -! dqmdt(i,k)=(mc(i,k)*(qhat(i,k)-q(i,k)) -! 1 +mc(i,k+1)*(q(i,k)-qhat(i,k+1)))/dp(i,k) -! 2 +du(i,k)*(qs(i,k)-q(i,k)) -! 3 +du(i,k)*(beta*ql(i,k)+(1-beta)*ql(i,k+1)) - - dqmdt(i,k) = (mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)+cp/rl* (su(i,k+1)-s(i,k)))- & - mu(i,k)* (qu(i,k)-qhat(i,k)+cp/rl*(su(i,k)-s(i,k)))+md(i,k+1)* & - (qd(i,k+1)-qhat(i,k+1)+cp/rl*(sd(i,k+1)-s(i,k)))-md(i,k)* & - (qd(i,k)-qhat(i,k)+cp/rl*(sd(i,k)-s(i,k))))/dp(i,k) + & - du(i,k)* (beta*ql(i,k)+(1-beta)*ql(i,k+1)) - end if - end do - end do -! - do k = msg + 1,pver - do i = il1g,il2g - if (k >= lel(i) .and. k <= lcl(i)) then - thetavp(i,k) = tp(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i))) - thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,k)) - dqsdtp(i,k) = qstp(i,k)* (1._r8+qstp(i,k)/eps1)*eps1*rl/(rd*tp(i,k)**2) -! -! dtpdt is the parcel temperature change due to change of -! subcloud layer properties during convection. -! - dtpdt(i,k) = tp(i,k)/ (1._r8+rl/cp* (dqsdtp(i,k)-qstp(i,k)/tp(i,k)))* & - (dtbdt(i)/t(i,mx(i))+rl/cp* (dqbdt(i)/tl(i)-q(i,mx(i))/ & - tl(i)**2*dtldt(i))) -! -! dboydt is the integrand of cape change. -! - dboydt(i,k) = ((dtpdt(i,k)/tp(i,k)+1._r8/(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i)))* & - (1.608_r8 * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608_r8/ & - (1._r8+0.608_r8*q(i,k))*dqmdt(i,k)))*grav*thetavp(i,k)/thetavm(i,k) - end if - end do - end do -! - do k = msg + 1,pver - do i = il1g,il2g - if (k > lcl(i) .and. k < mx(i)) then - thetavp(i,k) = tp(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,mx(i))) - thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,k)) -! -! dboydt is the integrand of cape change. -! - dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608_r8/ (1._r8+0.608_r8*q(i,mx(i)))*dqbdt(i)- & - dtmdt(i,k)/t(i,k)-0.608_r8/ (1._r8+0.608_r8*q(i,k))*dqmdt(i,k))* & - grav*thetavp(i,k)/thetavm(i,k) - end if - end do - end do - -! -! buoyant energy change is set to 2/3*excess cape per 3 hours -! - dadt(il1g:il2g) = 0._r8 - kmin = minval(lel(il1g:il2g)) - kmax = maxval(mx(il1g:il2g)) - 1 - do k = kmin, kmax - do i = il1g,il2g - if ( k >= lel(i) .and. k <= mx(i) - 1) then - dadt(i) = dadt(i) + dboydt(i,k)* (zf(i,k)-zf(i,k+1)) - endif - end do - end do - do i = il1g,il2g - dltaa = -1._r8* (cape(i)-capelmt) - if (dadt(i) /= 0._r8) mb(i) = max(dltaa/tau/dadt(i),0._r8) - end do -! - return -end subroutine closure - -subroutine q1q2_pjr(lchnk , & - dqdt ,dsdt ,q ,qs ,qu , & - su ,du ,qhat ,shat ,dp , & - mu ,md ,sd ,qd ,ql , & - dsubcld ,jt ,mx ,il1g ,il2g , & - cp ,rl ,msg , & - dl ,evp ,cu ) - - - implicit none - -!----------------------------------------------------------------------- -! -! Purpose: -! -! -! Method: -! -! -! -! Author: phil rasch dec 19 1995 -! -!----------------------------------------------------------------------- - - - real(r8), intent(in) :: cp - - integer, intent(in) :: lchnk ! chunk identifier - integer, intent(in) :: il1g - integer, intent(in) :: il2g - integer, intent(in) :: msg - - real(r8), intent(in) :: q(pcols,pver) - real(r8), intent(in) :: qs(pcols,pver) - real(r8), intent(in) :: qu(pcols,pver) - real(r8), intent(in) :: su(pcols,pver) - real(r8), intent(in) :: du(pcols,pver) - real(r8), intent(in) :: qhat(pcols,pver) - real(r8), intent(in) :: shat(pcols,pver) - real(r8), intent(in) :: dp(pcols,pver) - real(r8), intent(in) :: mu(pcols,pver) - real(r8), intent(in) :: md(pcols,pver) - real(r8), intent(in) :: sd(pcols,pver) - real(r8), intent(in) :: qd(pcols,pver) - real(r8), intent(in) :: ql(pcols,pver) - real(r8), intent(in) :: evp(pcols,pver) - real(r8), intent(in) :: cu(pcols,pver) - real(r8), intent(in) :: dsubcld(pcols) - - real(r8),intent(out) :: dqdt(pcols,pver),dsdt(pcols,pver) - real(r8),intent(out) :: dl(pcols,pver) - integer kbm - integer ktm - integer jt(pcols) - integer mx(pcols) -! -! work fields: -! - integer i - integer k - - real(r8) emc - real(r8) rl -!------------------------------------------------------------------- - do k = msg + 1,pver - do i = il1g,il2g - dsdt(i,k) = 0._r8 - dqdt(i,k) = 0._r8 - dl(i,k) = 0._r8 - end do - end do -! -! find the highest level top and bottom levels of convection -! - ktm = pver - kbm = pver - do i = il1g, il2g - ktm = min(ktm,jt(i)) - kbm = min(kbm,mx(i)) - end do - - do k = ktm,pver-1 - do i = il1g,il2g - emc = -cu (i,k) & ! condensation in updraft - +evp(i,k) ! evaporating rain in downdraft - - dsdt(i,k) = -rl/cp*emc & - + (+mu(i,k+1)* (su(i,k+1)-shat(i,k+1)) & - -mu(i,k)* (su(i,k)-shat(i,k)) & - +md(i,k+1)* (sd(i,k+1)-shat(i,k+1)) & - -md(i,k)* (sd(i,k)-shat(i,k)) & - )/dp(i,k) - - dqdt(i,k) = emc + & - (+mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)) & - -mu(i,k)* (qu(i,k)-qhat(i,k)) & - +md(i,k+1)* (qd(i,k+1)-qhat(i,k+1)) & - -md(i,k)* (qd(i,k)-qhat(i,k)) & - )/dp(i,k) - - dl(i,k) = du(i,k)*ql(i,k+1) - - end do - end do - -! -!DIR$ NOINTERCHANGE! - do k = kbm,pver - do i = il1g,il2g - if (k == mx(i)) then - dsdt(i,k) = (1._r8/dsubcld(i))* & - (-mu(i,k)* (su(i,k)-shat(i,k)) & - -md(i,k)* (sd(i,k)-shat(i,k)) & - ) - dqdt(i,k) = (1._r8/dsubcld(i))* & - (-mu(i,k)*(qu(i,k)-qhat(i,k)) & - -md(i,k)*(qd(i,k)-qhat(i,k)) & - ) - else if (k > mx(i)) then - dsdt(i,k) = dsdt(i,k-1) - dqdt(i,k) = dqdt(i,k-1) - end if - end do - end do -! - return -end subroutine q1q2_pjr - -subroutine buoyan_dilute(lchnk ,ncol , & - q ,t ,p ,z ,pf , & - tp ,qstp ,tl ,rl ,cape , & - pblt ,lcl ,lel ,lon ,mx , & - rd ,grav ,cp ,msg , & - tpert ,iclosure) -!----------------------------------------------------------------------- -! -! Purpose: -! Calculates CAPE the lifting condensation level and the convective top -! where buoyancy is first -ve. -! -! Method: Calculates the parcel temperature based on a simple constant -! entraining plume model. CAPE is integrated from buoyancy. -! 09/09/04 - Simplest approach using an assumed entrainment rate for -! testing (dmpdp). -! 08/04/05 - Swap to convert dmpdz to dmpdp -! -! SCAM Logical Switches - DILUTE:RBN - Now Disabled -! --------------------- -! switch(1) = .T. - Uses the dilute parcel calculation to obtain tendencies. -! switch(2) = .T. - Includes entropy/q changes due to condensate loss and freezing. -! switch(3) = .T. - Adds the PBL Tpert for the parcel temperature at all levels. -! -! References: -! Raymond and Blythe (1992) JAS -! -! Author: -! Richard Neale - September 2004 -! -!----------------------------------------------------------------------- - implicit none -!----------------------------------------------------------------------- -! -! input arguments -! - integer, intent(in) :: lchnk ! chunk identifier - integer, intent(in) :: ncol ! number of atmospheric columns - - real(r8), intent(in) :: q(pcols,pver) ! spec. humidity - real(r8), intent(in) :: t(pcols,pver) ! temperature - real(r8), intent(in) :: p(pcols,pver) ! pressure - real(r8), intent(in) :: z(pcols,pver) ! height - real(r8), intent(in) :: pf(pcols,pver+1) ! pressure at interfaces - real(r8), intent(in) :: pblt(pcols) ! index of pbl depth - real(r8), intent(in) :: tpert(pcols) ! perturbation temperature by pbl processes - logical, intent(in) :: iclosure ! true for normal procedure, otherwise use dcapemx from 1st call -! -! output arguments -! - real(r8), intent(out) :: tp(pcols,pver) ! parcel temperature - real(r8), intent(out) :: qstp(pcols,pver) ! saturation mixing ratio of parcel (only above lcl, just q below). - real(r8), intent(out) :: tl(pcols) ! parcel temperature at lcl - real(r8), intent(out) :: cape(pcols) ! convective aval. pot. energy. - integer lcl(pcols) ! - integer lel(pcols) ! - integer lon(pcols) ! level of onset of deep convection - integer mx(pcols) ! level of max moist static energy -! -!--------------------------Local Variables------------------------------ -! - real(r8) capeten(pcols,num_cin) ! provisional value of cape - real(r8) tv(pcols,pver) ! - real(r8) tpv(pcols,pver) ! - real(r8) buoy(pcols,pver) - - real(r8) pl(pcols) - real(r8) hmax(pcols) - real(r8) hmn(pcols) - - logical plge600(pcols) - integer knt(pcols) - integer lelten(pcols,num_cin) - -! DCAPE-ULL - real(r8) pblt600(pcols) - - real(r8) cp - real(r8) grav - - integer i - integer k - integer msg - integer n - integer bot_layer - - real(r8) rd - real(r8) rl -#ifdef PERGRO - real(r8) rhd -#endif -! -!----------------------------------------------------------------------- -! - do n = 1,num_cin - do i = 1,ncol - lelten(i,n) = pver - capeten(i,n) = 0._r8 - end do - end do -! - do i = 1,ncol - lon(i) = pver - knt(i) = 0 - lel(i) = pver - mx(i) = lon(i) - cape(i) = 0._r8 - hmax(i) = 0._r8 - end do - - tp(:ncol,:) = t(:ncol,:) - qstp(:ncol,:) = q(:ncol,:) - -!DCAPE-ULL - if (trigdcape_ull) then - do k = pver - 1,msg + 1,-1 - do i = 1,ncol - if ((p(i,k).le.600._r8) .and. (p(i,k+1).gt.600._r8)) pblt600(i) = dble(k) - end do - end do - endif - -!!! RBN - Initialize tv and buoy for output. -!!! tv=tv : tpv=tpv : qstp=q : buoy=0. - tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608_r8*q(:ncol,:))/ (1._r8+q(:ncol,:)) - tpv(:ncol,:) = tv(:ncol,:) - buoy(:ncol,:) = 0._r8 - -! -! set "launching" level(mx) to be at maximum moist static energy. -! search for this level stops at planetary boundary layer top. -! - bot_layer = pver - mx_bot_lyr_adj - -! DCAPE-ULL - if (trigdcape_ull .and. (.not. iclosure)) then - mx(:ncol) = dcapemx(:ncol) - else -#ifdef PERGRO - do k = bot_layer,msg + 1,-1 - do i = 1,ncol - hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k) -! -! Reset max moist static energy level when relative difference exceeds 1.e-4 -! - rhd = (hmn(i) - hmax(i))/(hmn(i) + hmax(i)) - - !DCAPE-ULL - if (trigdcape_ull) then - if (k >= nint(pblt600(i)) .and. k <= lon(i) .and. rhd > -1.e-4_r8) then - hmax(i) = hmn(i) - mx(i) = k - end if - else - if (k >= nint(pblt(i)) .and. k <= lon(i) .and. rhd > -1.e-4_r8) then - hmax(i) = hmn(i) - mx(i) = k - end if - end if - end do - end do -#else - do k = bot_layer,msg + 1,-1 - do i = 1,ncol - hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k) - - !DCAPE-ULL - if (trigdcape_ull) then - if (k >= nint(pblt600(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then - hmax(i) = hmn(i) - mx(i) = k - end if - else - if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then - hmax(i) = hmn(i) - mx(i) = k - end if - end if - end do - end do -#endif - - end if - -! LCL dilute calculation - initialize to mx(i) -! Determine lcl in parcel_dilute and get pl,tl after parcel_dilute -! Original code actually sets LCL as level above wher condensate forms. -! Therefore in parcel_dilute lcl(i) will be at first level where qsmix < qtmix. - - do i = 1,ncol ! Initialise LCL variables. - lcl(i) = mx(i) - tl(i) = t(i,mx(i)) - pl(i) = p(i,mx(i)) - end do - -! -! main buoyancy calculation. -! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! DILUTE PLUME CALCULATION USING ENTRAINING PLUME !!! -!!! RBN 9/9/04 !!! - - call parcel_dilute(lchnk, ncol, msg, mx, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl) - -! If lcl is above the nominal level of non-divergence (600 mbs), -! no deep convection is permitted (ensuing calculations -! skipped and cape retains initialized value of zero). -! - do i = 1,ncol - plge600(i) = pl(i).ge.600._r8 ! Just change to always allow buoy calculation. - end do - -! -! Main buoyancy calculation. -! - do k = pver,msg + 1,-1 - do i=1,ncol - if (k <= mx(i) .and. plge600(i)) then ! Define buoy from launch level to cloud top. - tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k)) - buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add ! +0.5K or not? - else - qstp(i,k) = q(i,k) - tp(i,k) = t(i,k) - tpv(i,k) = tv(i,k) - endif - end do - end do - - - -!------------------------------------------------------------------------------- - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - -! - do k = msg + 2,pver - do i = 1,ncol - if (k < lcl(i) .and. plge600(i)) then - if (buoy(i,k+1) > 0._r8 .and. buoy(i,k) <= 0._r8) then - knt(i) = min(num_cin,knt(i) + 1) - lelten(i,knt(i)) = k - end if - end if - end do - end do -! -! calculate convective available potential energy (cape). -! - do n = 1,num_cin - do k = msg + 1,pver - do i = 1,ncol - if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then - capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k)) - end if - end do - end do - end do -! -! find maximum cape from all possible tentative capes from -! one sounding, -! and use it as the final cape, april 26, 1995 -! - do n = 1,num_cin - do i = 1,ncol - if (capeten(i,n) > cape(i)) then - cape(i) = capeten(i,n) - lel(i) = lelten(i,n) - end if - end do - end do -! -! put lower bound on cape for diagnostic purposes. -! - do i = 1,ncol - cape(i) = max(cape(i), 0._r8) - end do -! - return -end subroutine buoyan_dilute - -subroutine parcel_dilute (lchnk, ncol, msg, klaunch, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl) -! Routine to determine -! 1. Tp - Parcel temperature -! 2. qstp - Saturated mixing ratio at the parcel temperature. - -!-------------------- -implicit none -!-------------------- - -integer, intent(in) :: lchnk -integer, intent(in) :: ncol -integer, intent(in) :: msg - -integer, intent(in), dimension(pcols) :: klaunch(pcols) - -real(r8), intent(in), dimension(pcols,pver) :: p -real(r8), intent(in), dimension(pcols,pver) :: t -real(r8), intent(in), dimension(pcols,pver) :: q -real(r8), intent(in), dimension(pcols) :: tpert ! PBL temperature perturbation. - -real(r8), intent(inout), dimension(pcols,pver) :: tp ! Parcel temp. -real(r8), intent(inout), dimension(pcols,pver) :: qstp ! Parcel water vapour (sat value above lcl). -real(r8), intent(inout), dimension(pcols) :: tl ! Actual temp of LCL. -real(r8), intent(inout), dimension(pcols) :: pl ! Actual pressure of LCL. - -integer, intent(inout), dimension(pcols) :: lcl ! Lifting condesation level (first model level with saturation). - -real(r8), intent(out), dimension(pcols,pver) :: tpv ! Define tpv within this routine. - -!-------------------- - -! Have to be careful as s is also dry static energy. - - -! If we are to retain the fact that CAM loops over grid-points in the internal -! loop then we need to dimension sp,atp,mp,xsh2o with ncol. - - -real(r8) tmix(pcols,pver) ! Tempertaure of the entraining parcel. -real(r8) qtmix(pcols,pver) ! Total water of the entraining parcel. -real(r8) qsmix(pcols,pver) ! Saturated mixing ratio at the tmix. -real(r8) smix(pcols,pver) ! Entropy of the entraining parcel. -real(r8) xsh2o(pcols,pver) ! Precipitate lost from parcel. -real(r8) ds_xsh2o(pcols,pver) ! Entropy change due to loss of condensate. -real(r8) ds_freeze(pcols,pver) ! Entropy change sue to freezing of precip. - -real(r8) mp(pcols) ! Parcel mass flux. -real(r8) qtp(pcols) ! Parcel total water. -real(r8) sp(pcols) ! Parcel entropy. - -real(r8) sp0(pcols) ! Parcel launch entropy. -real(r8) qtp0(pcols) ! Parcel launch total water. -real(r8) mp0(pcols) ! Parcel launch relative mass flux. - -real(r8) lwmax ! Maximum condesate that can be held in cloud before rainout. -real(r8) dmpdp ! Parcel fractional mass entrainment rate (/mb). -!real(r8) dmpdpc ! In cloud parcel mass entrainment rate (/mb). -!real(r8) dmpdz ! Parcel fractional mass entrainment rate (/m) -real(r8) dpdz,dzdp ! Hydrstatic relation and inverse of. -real(r8) senv ! Environmental entropy at each grid point. -real(r8) qtenv ! Environmental total water " " ". -real(r8) penv ! Environmental total pressure " " ". -real(r8) tenv ! Environmental total temperature " " ". -real(r8) new_s ! Hold value for entropy after condensation/freezing adjustments. -real(r8) new_q ! Hold value for total water after condensation/freezing adjustments. -real(r8) dp ! Layer thickness (center to center) -real(r8) tfguess ! First guess for entropy inversion - crucial for efficiency! -real(r8) tscool ! Super cooled temperature offset (in degC) (eg -35). - -real(r8) qxsk, qxskp1 ! LCL excess water (k, k+1) -real(r8) dsdp, dqtdp, dqxsdp ! LCL s, qt, p gradients (k, k+1) -real(r8) slcl,qtlcl,qslcl ! LCL s, qt, qs values. - -integer rcall ! Number of ientropy call for errors recording -integer nit_lheat ! Number of iterations for condensation/freezing loop. -integer i,k,ii ! Loop counters. - -!====================================================================== -! SUMMARY -! -! 9/9/04 - Assumes parcel is initiated from level of maxh (klaunch) -! and entrains at each level with a specified entrainment rate. -! -! 15/9/04 - Calculates lcl(i) based on k where qsmix is first < qtmix. -! -!====================================================================== -! -! Set some values that may be changed frequently. -! - -nit_lheat = 2 ! iterations for ds,dq changes from condensation freezing. - - -!dmpdpc = 3.e-2_r8 ! In cloud entrainment rate (/mb). -lwmax = 1.e-3_r8 ! Need to put formula in for this. -tscool = 0.0_r8 ! Temp at which water loading freezes in the cloud. - -qtmix=0._r8 -smix=0._r8 - -qtenv = 0._r8 -senv = 0._r8 -tenv = 0._r8 -penv = 0._r8 - -qtp0 = 0._r8 -sp0 = 0._r8 -mp0 = 0._r8 - -qtp = 0._r8 -sp = 0._r8 -mp = 0._r8 - -new_q = 0._r8 -new_s = 0._r8 - -! **** Begin loops **** - -do k = pver, msg+1, -1 - do i=1,ncol - -! Initialize parcel values at launch level. - - if (k == klaunch(i)) then - qtp0(i) = q(i,k) ! Parcel launch total water (assuming subsaturated) - OK????. - sp0(i) = entropy(t(i,k),p(i,k),qtp0(i)) ! Parcel launch entropy. - mp0(i) = 1._r8 ! Parcel launch relative mass (i.e. 1 parcel stays 1 parcel for dmpdp=0, undilute). - smix(i,k) = sp0(i) - qtmix(i,k) = qtp0(i) - tfguess = t(i,k) - rcall = 1 - call ientropy (rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess) - end if - -! Entraining levels - - if (k < klaunch(i)) then - -! Set environmental values for this level. - - dp = (p(i,k)-p(i,k+1)) ! In -ve mb as p decreasing with height - difference between center of layers. - qtenv = 0.5_r8*(q(i,k)+q(i,k+1)) ! Total water of environment. - tenv = 0.5_r8*(t(i,k)+t(i,k+1)) - penv = 0.5_r8*(p(i,k)+p(i,k+1)) - - senv = entropy(tenv,penv,qtenv) ! Entropy of environment. - -! Determine fractional entrainment rate /pa given value /m. - - dpdz = -(penv*grav)/(rgas*tenv) ! in mb/m since p in mb. - dzdp = 1._r8/dpdz ! in m/mb - dmpdp = dmpdz*dzdp ! /mb Fractional entrainment - -! Sum entrainment to current level -! entrains q,s out of intervening dp layers, in which linear variation is assumed -! so really it entrains the mean of the 2 stored values. - - sp(i) = sp(i) - dmpdp*dp*senv - qtp(i) = qtp(i) - dmpdp*dp*qtenv - mp(i) = mp(i) - dmpdp*dp - -! Entrain s and qt to next level. - - smix(i,k) = (sp0(i) + sp(i)) / (mp0(i) + mp(i)) - qtmix(i,k) = (qtp0(i) + qtp(i)) / (mp0(i) + mp(i)) - -! Invert entropy from s and q to determine T and saturation-capped q of mixture. -! t(i,k) used as a first guess so that it converges faster. - - tfguess = tmix(i,k+1) - rcall = 2 - call ientropy(rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess) - -! -! Determine if this is lcl of this column if qsmix <= qtmix. -! FIRST LEVEL where this happens on ascending. - - if (qsmix(i,k) <= qtmix(i,k) .and. qsmix(i,k+1) > qtmix(i,k+1)) then - lcl(i) = k - qxsk = qtmix(i,k) - qsmix(i,k) - qxskp1 = qtmix(i,k+1) - qsmix(i,k+1) - dqxsdp = (qxsk - qxskp1)/dp - pl(i) = p(i,k+1) - qxskp1/dqxsdp ! pressure level of actual lcl. - dsdp = (smix(i,k) - smix(i,k+1))/dp - dqtdp = (qtmix(i,k) - qtmix(i,k+1))/dp - slcl = smix(i,k+1) + dsdp* (pl(i)-p(i,k+1)) - qtlcl = qtmix(i,k+1) + dqtdp*(pl(i)-p(i,k+1)) - - tfguess = tmix(i,k) - rcall = 3 - call ientropy (rcall,i,lchnk,slcl,pl(i),qtlcl,tl(i),qslcl,tfguess) - -! write(iulog,*)' ' -! write(iulog,*)' p',p(i,k+1),pl(i),p(i,lcl(i)) -! write(iulog,*)' t',tmix(i,k+1),tl(i),tmix(i,lcl(i)) -! write(iulog,*)' s',smix(i,k+1),slcl,smix(i,lcl(i)) -! write(iulog,*)'qt',qtmix(i,k+1),qtlcl,qtmix(i,lcl(i)) -! write(iulog,*)'qs',qsmix(i,k+1),qslcl,qsmix(i,lcl(i)) - - endif -! - end if ! k < klaunch - - - end do ! Levels loop -end do ! Columns loop - -!!!!!!!!!!!!!!!!!!!!!!!!!!END ENTRAINMENT LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -!! Could stop now and test with this as it will provide some estimate of buoyancy -!! without the effects of freezing/condensation taken into account for tmix. - -!! So we now have a profile of entropy and total water of the entraining parcel -!! Varying with height from the launch level klaunch parcel=environment. To the -!! top allowed level for the existence of convection. - -!! Now we have to adjust these values such that the water held in vaopor is < or -!! = to qsmix. Therefore, we assume that the cloud holds a certain amount of -!! condensate (lwmax) and the rest is rained out (xsh2o). This, obviously -!! provides latent heating to the mixed parcel and so this has to be added back -!! to it. But does this also increase qsmix as well? Also freezing processes - - -xsh2o = 0._r8 -ds_xsh2o = 0._r8 -ds_freeze = 0._r8 - -!!!!!!!!!!!!!!!!!!!!!!!!!PRECIPITATION/FREEZING LOOP!!!!!!!!!!!!!!!!!!!!!!!!!! -!! Iterate solution twice for accuracy - - - -do k = pver, msg+1, -1 - do i=1,ncol - -! Initialize variables at k=klaunch - - if (k == klaunch(i)) then - -! Set parcel values at launch level assume no liquid water. - - tp(i,k) = tmix(i,k) - qstp(i,k) = q(i,k) - tpv(i,k) = (tp(i,k) + tp_fac*tpert(i)) * (1._r8+1.608_r8*qstp(i,k)) / (1._r8+qstp(i,k)) - - end if - - if (k < klaunch(i)) then - -! Initiaite loop if switch(2) = .T. - RBN:DILUTE - TAKEN OUT BUT COULD BE RETURNED LATER. - -! Iterate nit_lheat times for s,qt changes. - - do ii=0,nit_lheat-1 - -! Rain (xsh2o) is excess condensate, bar LWMAX (Accumulated loss from qtmix). - - xsh2o(i,k) = max (0._r8, qtmix(i,k) - qsmix(i,k) - lwmax) - -! Contribution to ds from precip loss of condensate (Accumulated change from smix).(-ve) - - ds_xsh2o(i,k) = ds_xsh2o(i,k+1) - cpliq * log (tmix(i,k)/tfreez) * max(0._r8,(xsh2o(i,k)-xsh2o(i,k+1))) -! -! Entropy of freezing: latice times amount of water involved divided by T. -! - - if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) == 0._r8) then ! One off freezing of condensate. - ds_freeze(i,k) = (latice/tmix(i,k)) * max(0._r8,qtmix(i,k)-qsmix(i,k)-xsh2o(i,k)) ! Gain of LH - end if - - if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) /= 0._r8) then ! Continual freezing of additional condensate. - ds_freeze(i,k) = ds_freeze(i,k+1)+(latice/tmix(i,k)) * max(0._r8,(qsmix(i,k+1)-qsmix(i,k))) - end if - -! Adjust entropy and accordingly to sum of ds (be careful of signs). - - new_s = smix(i,k) + ds_xsh2o(i,k) + ds_freeze(i,k) - -! Adjust liquid water and accordingly to xsh2o. - - new_q = qtmix(i,k) - xsh2o(i,k) - -! Invert entropy to get updated Tmix and qsmix of parcel. - - tfguess = tmix(i,k) - rcall =4 - call ientropy (rcall,i,lchnk,new_s, p(i,k), new_q, tmix(i,k), qsmix(i,k), tfguess) - - end do ! Iteration loop for freezing processes. - -! tp - Parcel temp is temp of mixture. -! tpv - Parcel v. temp should be density temp with new_q total water. - - tp(i,k) = tmix(i,k) - -! tpv = tprho in the presence of condensate (i.e. when new_q > qsmix) - - if (new_q > qsmix(i,k)) then ! Super-saturated so condensate present - reduces buoyancy. - qstp(i,k) = qsmix(i,k) - else ! Just saturated/sub-saturated - no condensate virtual effects. - qstp(i,k) = new_q - end if - - tpv(i,k) = (tp(i,k)+tp_fac*tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+ new_q) - - end if ! k < klaunch - - end do ! Loop for columns - -end do ! Loop for vertical levels. - - -return -end subroutine parcel_dilute - -!----------------------------------------------------------------------------------------- -real(r8) function entropy(TK,p,qtot) -!----------------------------------------------------------------------------------------- -! -! TK(K),p(mb),qtot(kg/kg) -! from Raymond and Blyth 1992 -! - real(r8), intent(in) :: p,qtot,TK - real(r8) :: qv,qst,e,est,L - real(r8), parameter :: pref = 1000._r8 - -L = rl - (cpliq - cpwv)*(TK-tfreez) ! T IN CENTIGRADE - -call qsat_hPa(TK, p, est, qst) - -qv = min(qtot,qst) ! Partition qtot into vapor part only. -e = qv*p / (eps1 +qv) - -entropy = (cpres + qtot*cpliq)*log( TK/tfreez) - rgas*log( (p-e)/pref ) + & - L*qv/TK - qv*rh2o*log(qv/qst) - -end FUNCTION entropy - -! -!----------------------------------------------------------------------------------------- -SUBROUTINE ientropy (rcall,icol,lchnk,s,p,qt,T,qst,Tfg) -!----------------------------------------------------------------------------------------- -! -! p(mb), Tfg/T(K), qt/qv(kg/kg), s(J/kg). -! Inverts entropy, pressure and total water qt -! for T and saturated vapor mixing ratio -! - -! use phys_grid, only: get_rlon_p, get_rlat_p - - integer, intent(in) :: icol, lchnk, rcall - real(r8), intent(in) :: s, p, Tfg, qt - real(r8), intent(out) :: qst, T - real(r8) :: est - real(r8) :: a,b,c,d,ebr,fa,fb,fc,pbr,qbr,rbr,sbr,tol1,xm,tol - integer :: i - - logical :: converged - - ! Max number of iteration loops. - integer, parameter :: LOOPMAX = 100 - real(r8), parameter :: EPS = 3.e-8_r8 - - converged = .false. - - ! Invert the entropy equation -- use Brent's method - ! Brent, R. P. Ch. 3-4 in Algorithms for Minimization Without Derivatives. Englewood Cliffs, NJ: Prentice-Hall, 1973. - - T = Tfg ! Better first guess based on Tprofile from conv. - - a = Tfg-10 !low bracket - b = Tfg+10 !high bracket - - fa = entropy(a, p, qt) - s - fb = entropy(b, p, qt) - s - - c=b - fc=fb - tol=0.001_r8 - - converge: do i=0, LOOPMAX - if ((fb > 0.0_r8 .and. fc > 0.0_r8) .or. & - (fb < 0.0_r8 .and. fc < 0.0_r8)) then - c=a - fc=fa - d=b-a - ebr=d - end if - if (abs(fc) < abs(fb)) then - a=b - b=c - c=a - fa=fb - fb=fc - fc=fa - end if - - tol1=2.0_r8*EPS*abs(b)+0.5_r8*tol - xm=0.5_r8*(c-b) - converged = (abs(xm) <= tol1 .or. fb == 0.0_r8) - if (converged) exit converge - - if (abs(ebr) >= tol1 .and. abs(fa) > abs(fb)) then - sbr=fb/fa - if (a == c) then - pbr=2.0_r8*xm*sbr - qbr=1.0_r8-sbr - else - qbr=fa/fc - rbr=fb/fc - pbr=sbr*(2.0_r8*xm*qbr*(qbr-rbr)-(b-a)*(rbr-1.0_r8)) - qbr=(qbr-1.0_r8)*(rbr-1.0_r8)*(sbr-1.0_r8) - end if - if (pbr > 0.0_r8) qbr=-qbr - pbr=abs(pbr) - if (2.0_r8*pbr < min(3.0_r8*xm*qbr-abs(tol1*qbr),abs(ebr*qbr))) then - ebr=d - d=pbr/qbr - else - d=xm - ebr=d - end if - else - d=xm - ebr=d - end if - a=b - fa=fb - b=b+merge(d,sign(tol1,xm), abs(d) > tol1 ) - - fb = entropy(b, p, qt) - s - - end do converge - - T = b - call qsat_hPa(T, p, est, qst) - -! if (.not. converged) then -! this_lat = get_rlat_p(lchnk, icol)*57.296_r8 -! this_lon = get_rlon_p(lchnk, icol)*57.296_r8 -! write(*,*) '*** ZM_CONV: IENTROPY: Failed and about to exit, info follows ****' -! write(*,*) 'ZM_CONV: IENTROPY. Details: call#,lchnk,icol= ',rcall,lchnk,icol, & -! ' lat: ',this_lat,' lon: ',this_lon, & -! ' P(mb)= ', p, ' Tfg(K)= ', Tfg, ' qt(g/kg) = ', 1000._r8*qt, & -! ' qst(g/kg) = ', 1000._r8*qst,', s(J/kg) = ',s -! call endrun('**** ZM_CONV IENTROPY: Tmix did not converge ****') -! end if - -100 format (A,I1,I4,I4,7(A,F6.2)) - -end SUBROUTINE ientropy - -! Wrapper for qsat_water that does translation between Pa and hPa -! qsat_water uses Pa internally, so get it right, need to pass in Pa. -! Afterward, set es back to hPa. -subroutine qsat_hPa(t, p, es, qm) -! use wv_saturation, only: qsat_water - - ! Inputs - real(r8), intent(in) :: t ! Temperature (K) - real(r8), intent(in) :: p ! Pressure (hPa) - ! Outputs - real(r8), intent(out) :: es ! Saturation vapor pressure (hPa) - real(r8), intent(out) :: qm ! Saturation mass mixing ratio - ! (vapor mass over dry mass, kg/kg) - - call qsat(t, p*100._r8, es, qm, 0) - - es = es*0.01_r8 - -end subroutine qsat_hPa - -! AaronDonahue - Dummy function to handle whether or not this is the first -! step, until AD can provide this information directly. -function is_first_step() result(val) - logical :: val - val = is_first_step_local -end function is_first_step - -function is_first_restart_step() result(val) - logical :: val - val = .false. -end function is_first_restart_step -! AaronDonahue - Taken from cloud_fraction.F90 and adapted for local use. -! TODO, we should probably use the ice cloud fraction from macrophysics - subroutine cldfrc_fice(ncol, t, fice, fsnow) -! -! Compute the fraction of the total cloud water which is in ice phase. -! The fraction depends on temperature only. -! This is the form that was used for radiation, the code came from cldefr originally -! -! Author: B. A. Boville Sept 10, 2002 -! modified: PJR 3/13/03 (added fsnow to ascribe snow production for convection ) -! -! AaronDonahue - taken from original location as a stopgap to get ZM working -! locally. -!----------------------------------------------------------------------- - -! Arguments - integer, intent(in) :: ncol ! number of active columns - real(r8), intent(in) :: t(pcols,pver) ! temperature - - real(r8), intent(out) :: fice(pcols,pver) ! Fractional ice content within cloud - real(r8), intent(out) :: fsnow(pcols,pver) ! Fractional snow content for convection - -! Local variables - real(r8) :: tmax_fice ! max temperature for cloud ice formation - real(r8) :: tmin_fice ! min temperature for cloud ice formation - real(r8) :: tmax_fsnow ! max temperature for transition to convective snow - real(r8) :: tmin_fsnow ! min temperature for transition to convective snow - - integer :: i,k ! loop indexes - ! Top level - integer :: top_lev = 1 - -!----------------------------------------------------------------------- - - tmax_fice = tmelt - 10._r8 ! max temperature for cloud ice formation - tmin_fice = tmax_fice - 30._r8 ! min temperature for cloud ice formation - tmax_fsnow = tmelt ! max temperature for transition to convective snow - tmin_fsnow = tmelt - 5._r8 ! min temperature for transition to convective snow - - fice(:,:top_lev-1) = 0._r8 - fsnow(:,:top_lev-1) = 0._r8 - -! Define fractional amount of cloud that is ice - do k=top_lev,pver - do i=1,ncol - -! If warmer than tmax then water phase - if (t(i,k) > tmax_fice) then - fice(i,k) = 0.0_r8 - -! If colder than tmin then ice phase - else if (t(i,k) < tmin_fice) then - fice(i,k) = 1.0_r8 - -! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax - else - fice(i,k) =(tmax_fice - t(i,k)) / (tmax_fice - tmin_fice) - end if - -! snow fraction partitioning - -! If warmer than tmax then water phase - if (t(i,k) > tmax_fsnow) then - fsnow(i,k) = 0.0_r8 - -! If colder than tmin then ice phase - else if (t(i,k) < tmin_fsnow) then - fsnow(i,k) = 1.0_r8 - -! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax - else - fsnow(i,k) =(tmax_fsnow - t(i,k)) / (tmax_fsnow - tmin_fsnow) - end if - - end do - end do - - end subroutine cldfrc_fice -!--------------------------------------------------------------------- -! QSAT (SPECIFIC HUMIDITY) PROCEDURES -!--------------------------------------------------------------------- -! AaronDonahue - Taken from micro_p3.F90 to get ZM working locally. -! TODO - A unified qsat calculation that can be used by multiple physics -! processes? -!=========================================================================================== - - subroutine qsat(t_atm,p_atm,e_pres,qv_sat, i_wrt) - - !------------------------------------------------------------------------------------ - ! Calls polysvp1 to obtain the saturation vapor pressure, and then computes - ! and returns the saturation mixing ratio, with respect to either liquid or ice, - ! depending on value of 'i_wrt' - !------------------------------------------------------------------------------------ - - implicit none - - !Calling parameters: - real(r8), intent(IN) :: t_atm !temperature [K] - real(r8), intent(IN) :: p_atm !pressure [Pa] - real(r8), intent(OUT) :: e_pres !saturation vapor pressure [Pa] - real(r8), intent(OUT) :: qv_sat !saturation mixing ratio - integer, intent(IN) :: i_wrt !index, 0 = w.r.t. liquid, 1 = w.r.t. ice - - !Local variables: - real(r8) :: ep_2 - - !------------------ - ep_2 = epsilo - e_pres = polysvp1(t_atm,i_wrt) - qv_sat = ep_2*e_pres/max(1.e-3_r8,(p_atm-e_pres)) - - return - - end subroutine qsat -!==========================================================================================! - !_rtype - real(r8) function polysvp1(t,i_type) - - !------------------------------------------- - ! COMPUTE SATURATION VAPOR PRESSURE - ! POLYSVP1 RETURNED IN UNITS OF PA. - ! T IS INPUT IN UNITS OF K. - ! i_type REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1) - !------------------------------------------- - - use scream_abortutils, only : endscreamrun - - implicit none - - real(r8), intent(in) :: t - integer, intent(in) :: i_type - - ! REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM FLATAU ET AL. 1992, TABLE 4 (RIGHT-HAND COLUMN) - - !local variables - character(len=1000) :: err_msg - - ! ice - real(r8) a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i - data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /& - 6.11147274_r8, 0.503160820_r8, 0.188439774e-1_r8, & - 0.420895665e-3_r8, 0.615021634e-5_r8, 0.602588177e-7_r8, & - 0.385852041e-9_r8, 0.146898966e-11_r8, 0.252751365e-14_r8/ - - ! liquid - real(r8) a0,a1,a2,a3,a4,a5,a6,a7,a8 - - ! V1.7 - data a0,a1,a2,a3,a4,a5,a6,a7,a8 /& - 6.11239921_r8, 0.443987641_r8, 0.142986287e-1_r8, & - 0.264847430e-3_r8, 0.302950461e-5_r8, 0.206739458e-7_r8, & - 0.640689451e-10_r8,-0.952447341e-13_r8,-0.976195544e-15_r8/ - real(r8) dt - real(r8) zerodegc - !------------------------------------------- - zerodegc = tmelt - if (i_type.eq.1 .and. t.lt.zerodegc) then - ! ICE - - ! Flatau formulation: - dt = max(-80._r8,t-273.15_r8) - polysvp1 = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+ & - a8i*dt))))))) - polysvp1 = polysvp1*100._r8 - - ! Goff-Gratch formulation: - ! POLYSVP1 = 10.**(-9.09718*(273.16/T-1.)-3.56654* & - ! log10(273.16/T)+0.876793*(1.-T/273.16)+ & - ! log10(6.1071))*100. - - - elseif (i_type.eq.0 .or. t.ge.zerodegc) then - ! LIQUID - - ! Flatau formulation: - dt = max(-80._r8,t-273.15_r8) - polysvp1 = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) - polysvp1 = polysvp1*100._r8 - - ! Goff-Gratch formulation: - ! POLYSVP1 = 10.**(-7.90298*(373.16/T-1.)+ & - ! 5.02808*log10(373.16/T)- & - ! 1.3816E-7*(10**(11.344*(1.-T/373.16))-1.)+ & - ! 8.1328E-3*(10**(-3.49149*(373.16/T-1.))-1.)+ & - ! log10(1013.246))*100. - - !PMC added error checking - else - - write(err_msg,*)'** polysvp1 i_type must be 0 or 1 but is: ', & - i_type,' temperature is:',t,' in file: ',__FILE__,' at line:',__LINE__ - - call endscreamrun(err_msg) - endif - - return - - end function polysvp1 -! Below is what was originally used by ZM and is taken from wv_saturation -!elemental subroutine qsat(t, p, es, qs, gam, dqsdt, enthalpy) -! !------------------------------------------------------------------! -! ! Purpose: ! -! ! Look up and return saturation vapor pressure from precomputed ! -! ! table, then calculate and return saturation specific humidity. ! -! ! Optionally return various temperature derivatives or enthalpy ! -! ! at saturation. ! -! !------------------------------------------------------------------! -! -! ! Inputs -! real(r8), intent(in) :: t ! Temperature -! real(r8), intent(in) :: p ! Pressure -! ! Outputs -! real(r8), intent(out) :: es ! Saturation vapor pressure -! real(r8), intent(out) :: qs ! Saturation specific humidity -! -! real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) -! real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) -! real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q -! -! ! Local variables -! real(r8) :: hltalt ! Modified latent heat for T derivatives -! real(r8) :: tterm ! Account for d(es)/dT in transition region -! -! es = estblf(t) -! -! qs = svp_to_qsat(es, p) -! -! ! Ensures returned es is consistent with limiters on qs. -! es = min(es, p) -! -! ! Calculate optional arguments. -! if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then -! -! ! "generalized" analytic expression for t derivative of es -! ! accurate to within 1 percent for 173.16 < t < 373.16 -! call calc_hltalt(t, hltalt, tterm) -! -! if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt) -! -! call deriv_outputs(t, p, es, qs, hltalt, tterm, & -! gam=gam, dqsdt=dqsdt) -! -! end if -! -!end subroutine qsat -! -!elemental subroutine qsat_water(t, p, es, qs, gam, dqsdt, enthalpy) -! !------------------------------------------------------------------! -! ! Purpose: ! -! ! Calculate SVP over water at a given temperature, and then ! -! ! calculate and return saturation specific humidity. ! -! ! Optionally return various temperature derivatives or enthalpy ! -! ! at saturation. ! -! !------------------------------------------------------------------! -! -!! use wv_sat_methods, only: wv_sat_qsat_water -! -! ! Inputs -! real(r8), intent(in) :: t ! Temperature -! real(r8), intent(in) :: p ! Pressure -! ! Outputs -! real(r8), intent(out) :: es ! Saturation vapor pressure -! real(r8), intent(out) :: qs ! Saturation specific humidity -! -! real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) -! real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) -! real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q -! -! ! Local variables -! real(r8) :: hltalt ! Modified latent heat for T derivatives -! -! call wv_sat_qsat_water(t, p, es, qs) -! -! if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then -! -! ! "generalized" analytic expression for t derivative of es -! ! accurate to within 1 percent for 173.16 < t < 373.16 -! call no_ip_hltalt(t, hltalt) -! -! if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt) -! -! ! For pure water/ice transition term is 0. -! call deriv_outputs(t, p, es, qs, hltalt, 0._r8, & -! gam=gam, dqsdt=dqsdt) -! -! end if -! -!end subroutine qsat_water -! -!elemental subroutine wv_sat_qsat_water(t, p, es, qs, idx) -! !------------------------------------------------------------------! -! ! Purpose: ! -! ! Calculate SVP over water at a given temperature, and then ! -! ! calculate and return saturation specific humidity. ! -! !------------------------------------------------------------------! -! -! ! Inputs -! real(r8), intent(in) :: t ! Temperature -! real(r8), intent(in) :: p ! Pressure -! ! Outputs -! real(r8), intent(out) :: es ! Saturation vapor pressure -! real(r8), intent(out) :: qs ! Saturation specific humidity -! -! integer, intent(in), optional :: idx ! Scheme index -! -! es = wv_sat_svp_water(t, idx) -! -! qs = wv_sat_svp_to_qsat(es, p) -! -! ! Ensures returned es is consistent with limiters on qs. -! es = min(es, p) -! -!end subroutine wv_sat_qsat_water -!--------------------------------------------------------------------- - -end module zm_conv