From 64c7c9a7ea3164533d2ba17a9077ad79161bb637 Mon Sep 17 00:00:00 2001 From: Benjamin Hillman Date: Tue, 21 May 2024 12:19:25 -0700 Subject: [PATCH] Enable MISR simulator and update cosp dimension names --- .../eamxx/src/physics/cosp/cosp_c2f.F90 | 26 ++++---- .../eamxx/src/physics/cosp/cosp_functions.hpp | 25 +++++--- .../eamxx/src/physics/cosp/eamxx_cosp.cpp | 64 ++++++++++++++++--- .../eamxx/src/physics/cosp/eamxx_cosp.hpp | 14 +++- 4 files changed, 97 insertions(+), 32 deletions(-) diff --git a/components/eamxx/src/physics/cosp/cosp_c2f.F90 b/components/eamxx/src/physics/cosp/cosp_c2f.F90 index febef82c6b0..77db4858ccd 100644 --- a/components/eamxx/src/physics/cosp/cosp_c2f.F90 +++ b/components/eamxx/src/physics/cosp/cosp_c2f.F90 @@ -40,11 +40,11 @@ module cosp_c2f ! Local variables; control what runs and what does not logical :: & - lsingle = .false., & ! True if using MMF_v3_single_moment CLOUDSAT microphysical scheme (default) - ldouble = .true., & ! True if using MMF_v3.5_two_moment CLOUDSAT microphysical scheme + lsingle = .false., & ! True if using MMF_v3_single_moment CLOUDSAT microphysical scheme (default) + ldouble = .true. , & ! True if using MMF_v3.5_two_moment CLOUDSAT microphysical scheme lisccp = .true. , & ! Local on/off switch for simulators (used by initialization) - lmodis = .true., & ! - lmisr = .false., & ! + lmodis = .true. , & ! + lmisr = .true. , & ! lcalipso = .false., & ! lgrLidar532 = .false., & ! latlid = .false., & ! @@ -64,7 +64,7 @@ module cosp_c2f Lmeantbisccp = .false., & ! ISCCP mean all-sky 10.5micron brightness temperature Lmeantbclrisccp = .false., & ! ISCCP mean clear-sky 10.5micron brightness temperature Lalbisccp = .false., & ! ISCCP mean cloud albedo - LclMISR = .false., & ! MISR cloud fraction + LclMISR = .true. , & ! MISR cloud fraction Lcltmodis = .false., & ! MODIS total cloud fraction Lclwmodis = .false., & ! MODIS liquid cloud fraction Lclimodis = .false., & ! MODIS ice cloud fraction @@ -270,17 +270,18 @@ subroutine cosp_c2f_init(npoints, ncolumns, nlevels) bind(c, name='cosp_c2f_init end subroutine cosp_c2f_init - subroutine cosp_c2f_run(npoints, ncolumns, nlevels, ntau, nctp, & - emsfc_lw, sunlit, skt, T_mid, p_mid, p_int, qv, & - cldfrac, reff_qc, reff_qi, dtau067, dtau105, isccp_cldtot, isccp_ctptau, modis_ctptau & + subroutine cosp_c2f_run(npoints, ncolumns, nlevels, ntau, nctp, ncth, & + emsfc_lw, sunlit, skt, T_mid, p_mid, p_int, z_mid, qv, & + cldfrac, reff_qc, reff_qi, dtau067, dtau105, isccp_cldtot, isccp_ctptau, modis_ctptau, misr_cthtau & ) bind(C, name='cosp_c2f_run') - integer(kind=c_int), value, intent(in) :: npoints, ncolumns, nlevels, ntau, nctp + integer(kind=c_int), value, intent(in) :: npoints, ncolumns, nlevels, ntau, nctp, ncth real(kind=c_double), value, intent(in) :: emsfc_lw real(kind=c_double), intent(in), dimension(npoints) :: sunlit, skt - real(kind=c_double), intent(in), dimension(npoints,nlevels) :: T_mid, p_mid, qv, cldfrac, reff_qc, reff_qi, dtau067, dtau105 + real(kind=c_double), intent(in), dimension(npoints,nlevels) :: T_mid, p_mid, z_mid, qv, cldfrac, reff_qc, reff_qi, dtau067, dtau105 real(kind=c_double), intent(in), dimension(npoints,nlevels+1) :: p_int real(kind=c_double), intent(inout), dimension(npoints) :: isccp_cldtot real(kind=c_double), intent(inout), dimension(npoints,ntau,nctp) :: isccp_ctptau, modis_ctptau + real(kind=c_double), intent(inout), dimension(npoints,ntau,ncth) :: misr_cthtau ! Takes normal arrays as input and populates COSP derived types character(len=256),dimension(100) :: cosp_status integer :: nptsperit @@ -326,7 +327,7 @@ subroutine cosp_c2f_run(npoints, ncolumns, nlevels, ntau, nctp, & ! Translate arrays to derived types cospIN%emsfc_lw = emsfc_lw cospIN%rcfg_cloudsat = rcfg_cloudsat -! cospstateIN%hgt_matrix = zlev(start_idx:end_idx,Nlevels:1:-1) ! km + cospstateIN%hgt_matrix = z_mid(start_idx:end_idx,1:Nlevels) ! m cospstateIN%sunlit = sunlit(start_idx:end_idx) ! 0-1 cospstateIN%skt = skt(start_idx:end_idx) ! K ! cospstateIN%surfelev = surfelev(start_idx:end_idx) ! m @@ -358,6 +359,9 @@ subroutine cosp_c2f_run(npoints, ncolumns, nlevels, ntau, nctp, & ! Modis modis_ctptau(:npoints,:,:) = cospOUT%modis_Optical_Thickness_vs_Cloud_Top_Pressure(:npoints,:,:) + ! MISR + misr_cthtau(:npoints,:,:) = cospOUT%misr_fq(:npoints,:,:) + end subroutine cosp_c2f_run subroutine cosp_c2f_final() bind(C, name='cosp_c2f_final') diff --git a/components/eamxx/src/physics/cosp/cosp_functions.hpp b/components/eamxx/src/physics/cosp/cosp_functions.hpp index ef6e5493d46..37377ca1475 100644 --- a/components/eamxx/src/physics/cosp/cosp_functions.hpp +++ b/components/eamxx/src/physics/cosp/cosp_functions.hpp @@ -4,11 +4,11 @@ using scream::Real; extern "C" void cosp_c2f_init(int ncol, int nsubcol, int nlay); extern "C" void cosp_c2f_final(); -extern "C" void cosp_c2f_run(const int ncol, const int nsubcol, const int nlay, const int ntau, const int nctp, +extern "C" void cosp_c2f_run(const int ncol, const int nsubcol, const int nlay, const int ntau, const int nctp, const int ncth, const Real emsfc_lw, const Real* sunlit, const Real* skt, - const Real* T_mid, const Real* p_mid, const Real* p_int, const Real* qv, + const Real* T_mid, const Real* p_mid, const Real* p_int, const Real* z_mid, const Real* qv, const Real* cldfrac, const Real* reff_qc, const Real* reff_qi, const Real* dtau067, const Real* dtau105, - Real* isccp_cldtot, Real* isccp_ctptau, Real* modis_ctptau); + Real* isccp_cldtot, Real* isccp_ctptau, Real* modis_ctptau, Real* misr_cthtau); namespace scream { @@ -30,28 +30,30 @@ namespace scream { cosp_c2f_final(); }; inline void main( - const Int ncol, const Int nsubcol, const Int nlay, const Int ntau, const Int nctp, const Real emsfc_lw, + const Int ncol, const Int nsubcol, const Int nlay, const Int ntau, const Int nctp, const Int ncth, const Real emsfc_lw, view_1d& sunlit , view_1d& skt, view_2d& T_mid , view_2d& p_mid , view_2d& p_int, - view_2d& qv , view_2d& cldfrac, + view_2d& z_mid , view_2d& qv , view_2d& cldfrac, view_2d& reff_qc, view_2d& reff_qi, view_2d& dtau067, view_2d& dtau105, - view_1d& isccp_cldtot , view_3d& isccp_ctptau, view_3d& modis_ctptau) { + view_1d& isccp_cldtot , view_3d& isccp_ctptau, view_3d& modis_ctptau, view_3d& misr_cthtau) { // Make host copies and permute data as needed lview_host_2d T_mid_h("T_mid_h", ncol, nlay), p_mid_h("p_mid_h", ncol, nlay), p_int_h("p_int_h", ncol, nlay+1), - qv_h("qv_h", ncol, nlay), cldfrac_h("cldfrac_h", ncol, nlay), + z_mid_h("z_mid_h", ncol, nlay), qv_h("qv_h", ncol, nlay), cldfrac_h("cldfrac_h", ncol, nlay), reff_qc_h("reff_qc_h", ncol, nlay), reff_qi_h("reff_qi_h", ncol, nlay), dtau067_h("dtau_067_h", ncol, nlay), dtau105_h("dtau105_h", ncol, nlay); lview_host_3d isccp_ctptau_h("isccp_ctptau_h", ncol, ntau, nctp); lview_host_3d modis_ctptau_h("modis_ctptau_h", ncol, ntau, nctp); + lview_host_3d misr_cthtau_h("misr_cthtau_h", ncol, ntau, ncth); // Copy to layoutLeft host views for (int i = 0; i < ncol; i++) { for (int j = 0; j < nlay; j++) { T_mid_h(i,j) = T_mid(i,j); p_mid_h(i,j) = p_mid(i,j); + z_mid_h(i,j) = z_mid(i,j); qv_h(i,j) = qv(i,j); cldfrac_h(i,j) = cldfrac(i,j); reff_qc_h(i,j) = reff_qc(i,j); @@ -69,11 +71,11 @@ namespace scream { // Subsample here? // Call COSP wrapper - cosp_c2f_run(ncol, nsubcol, nlay, ntau, nctp, + cosp_c2f_run(ncol, nsubcol, nlay, ntau, nctp, ncth, emsfc_lw, sunlit.data(), skt.data(), T_mid_h.data(), p_mid_h.data(), p_int_h.data(), - qv_h.data(), + z_mid_h.data(), qv_h.data(), cldfrac_h.data(), reff_qc_h.data(), reff_qi_h.data(), dtau067_h.data(), dtau105_h.data(), - isccp_cldtot.data(), isccp_ctptau_h.data(), modis_ctptau_h.data()); + isccp_cldtot.data(), isccp_ctptau_h.data(), modis_ctptau_h.data(), misr_cthtau_h.data()); // Copy outputs back to layoutRight views for (int i = 0; i < ncol; i++) { @@ -82,6 +84,9 @@ namespace scream { isccp_ctptau(i,j,k) = isccp_ctptau_h(i,j,k); modis_ctptau(i,j,k) = modis_ctptau_h(i,j,k); } + for (int k = 0; k < ncth; k++) { + misr_cthtau(i,j,k) = misr_cthtau_h(i,j,k); + } } } } diff --git a/components/eamxx/src/physics/cosp/eamxx_cosp.cpp b/components/eamxx/src/physics/cosp/eamxx_cosp.cpp index a86d1926e88..f01b619797e 100644 --- a/components/eamxx/src/physics/cosp/eamxx_cosp.cpp +++ b/components/eamxx/src/physics/cosp/eamxx_cosp.cpp @@ -41,6 +41,8 @@ void Cosp::set_grids(const std::shared_ptr grids_manager) auto percent = Units::nondimensional(); percent.set_string("%"); auto micron = m / 1000000; + auto m2 = m * m; + auto s2 = s * s; m_grid = grids_manager->get_grid("Physics"); const auto& grid_name = m_grid->name(); @@ -54,8 +56,11 @@ void Cosp::set_grids(const std::shared_ptr grids_manager) FieldLayout scalar3d_mid = m_grid->get_3d_scalar_layout(true); FieldLayout scalar3d_int = m_grid->get_3d_scalar_layout(false); FieldLayout scalar4d_ctptau ( {COL,CMP,CMP}, - {m_num_cols,m_num_isccptau,m_num_isccpctp}, - {e2str(COL), "ISCCPTAU", "ISCCPPRS"}); + {m_num_cols,m_num_tau,m_num_ctp}, + {e2str(COL), "cosp_tau", "cosp_prs"}); + FieldLayout scalar4d_cthtau ( {COL,CMP,CMP}, + {m_num_cols,m_num_tau,m_num_cth}, + {e2str(COL), "cosp_tau", "cosp_cth"}); // Set of fields used strictly as input // Name in AD Layout Units Grid Group @@ -68,6 +73,8 @@ void Cosp::set_grids(const std::shared_ptr grids_manager) //add_field("height_mid", scalar3d_mid, m, grid_name); //add_field("height_int", scalar3d_int, m, grid_name); add_field("T_mid", scalar3d_mid, K, grid_name); + add_field("phis", scalar2d , m2/s2, grid_name); + add_field("pseudo_density", scalar3d_mid, Pa, grid_name); add_field("qv", scalar3d_mid, Q, grid_name, "tracers"); add_field("qc", scalar3d_mid, Q, grid_name, "tracers"); add_field("qi", scalar3d_mid, Q, grid_name, "tracers"); @@ -86,6 +93,7 @@ void Cosp::set_grids(const std::shared_ptr grids_manager) add_field("isccp_cldtot", scalar2d, percent, grid_name); add_field("isccp_ctptau", scalar4d_ctptau, percent, grid_name, 1); add_field("modis_ctptau", scalar4d_ctptau, percent, grid_name, 1); + add_field("misr_cthtau" , scalar4d_cthtau, percent, grid_name, 1); add_field("isccp_mask" , scalar2d, nondim, grid_name); } @@ -99,8 +107,8 @@ void Cosp::initialize_impl (const RunType /* run_type */) // Add note to output files about processing ISCCP fields that are only valid during // daytime. This can go away once I/O can handle masked time averages. using stratts_t = std::map; - std::list vnames = {"isccp_cldtot", "isccp_ctptau", "modis_ctptau"}; - for (const auto field_name : {"isccp_cldtot", "isccp_ctptau", "modis_ctptau"}) { + std::list vnames = {"isccp_cldtot", "isccp_ctptau", "modis_ctptau", "misr_cthtau"}; + for (const auto field_name : {"isccp_cldtot", "isccp_ctptau", "modis_ctptau", "misr_cthtau"}) { auto& f = get_field_out(field_name); auto& atts = f.get_header().get_extra_data("io: string attributes"); atts["note"] = "Night values are zero; divide by isccp_mask to get daytime mean"; @@ -147,6 +155,8 @@ void Cosp::run_impl (const double dt) get_field_in("eff_radius_qi").sync_to_host(); get_field_in("dtau067").sync_to_host(); get_field_in("dtau105").sync_to_host(); + get_field_in("phis").sync_to_host(); + get_field_in("pseudo_density").sync_to_host(); auto qv = get_field_in("qv").get_view(); auto qc = get_field_in("qc").get_view(); @@ -156,6 +166,8 @@ void Cosp::run_impl (const double dt) auto T_mid = get_field_in("T_mid").get_view(); auto p_mid = get_field_in("p_mid").get_view(); auto p_int = get_field_in("p_int").get_view(); + auto phis = get_field_in("phis").get_view(); + auto pseudo_density = get_field_in("pseudo_density").get_view(); auto cldfrac = get_field_in("cldfrac_rad").get_view(); auto reff_qc = get_field_in("eff_radius_qc").get_view(); auto reff_qi = get_field_in("eff_radius_qi").get_view(); @@ -164,28 +176,60 @@ void Cosp::run_impl (const double dt) auto isccp_cldtot = get_field_out("isccp_cldtot").get_view(); auto isccp_ctptau = get_field_out("isccp_ctptau").get_view(); auto modis_ctptau = get_field_out("modis_ctptau").get_view(); + auto misr_cthtau = get_field_out("misr_cthtau" ).get_view(); auto isccp_mask = get_field_out("isccp_mask" ).get_view(); // Copy of sunlit flag with COSP frequency for proper averaging + // Compute heights + const auto z_mid = CospFunc::view_2d("z_mid", m_num_cols, m_num_levs); + const auto z_int = CospFunc::view_2d("z_int", m_num_cols, m_num_levs); + const auto dz = z_mid; // reuse tmp memory for dz + const auto ncol = m_num_cols; + const auto nlev = m_num_levs; + // calculate_z_int contains a team-level parallel_scan, which requires a special policy + // TODO: do this on device? + const auto scan_policy = ekat::ExeSpaceUtils::get_thread_range_parallel_scan_team_policy(ncol, nlev); + Kokkos::parallel_for(scan_policy, KOKKOS_LAMBDA (const KTH::MemberType& team) { + const int i = team.league_rank(); + const auto dz_s = ekat::subview(dz, i); + const auto p_mid_s = ekat::subview(p_mid, i); + const auto T_mid_s = ekat::subview(T_mid, i); + const auto qv_s = ekat::subview(qv, i); + const auto z_int_s = ekat::subview(z_int, i); + const auto z_mid_s = ekat::subview(z_mid, i); + const Real z_surf = phis(i) / 9.81; + const auto pseudo_density_s = ekat::subview(pseudo_density, i); + PF::calculate_dz(team, pseudo_density_s, p_mid_s, T_mid_s, qv_s, dz_s); + team.team_barrier(); + PF::calculate_z_int(team,nlev,dz_s,z_surf,z_int_s); + team.team_barrier(); + PF::calculate_z_mid(team,nlev,z_int_s,z_mid_s); + team.team_barrier(); + }); + // Call COSP wrapper routines if (update_cosp) { Real emsfc_lw = 0.99; Kokkos::deep_copy(isccp_mask, sunlit); + CospFunc::view_2d z_mid_c = z_mid; // Need a const version of z_mid for call to CospFunc::main CospFunc::main( - m_num_cols, m_num_subcols, m_num_levs, m_num_isccptau, m_num_isccpctp, - emsfc_lw, sunlit, skt, T_mid, p_mid, p_int, qv, + m_num_cols, m_num_subcols, m_num_levs, m_num_tau, m_num_ctp, m_num_cth, + emsfc_lw, sunlit, skt, T_mid, p_mid, p_int, z_mid_c, qv, cldfrac, reff_qc, reff_qi, dtau067, dtau105, - isccp_cldtot, isccp_ctptau, modis_ctptau + isccp_cldtot, isccp_ctptau, modis_ctptau, misr_cthtau ); // Remask night values to ZERO since our I/O does not know how to handle masked/missing values // in temporal averages; this is all host data, so we can just use host loops like its the 1980s for (int i = 0; i < m_num_cols; i++) { if (sunlit(i) == 0) { isccp_cldtot(i) = 0; - for (int j = 0; j < m_num_isccptau; j++) { - for (int k = 0; k < m_num_isccpctp; k++) { + for (int j = 0; j < m_num_tau; j++) { + for (int k = 0; k < m_num_ctp; k++) { isccp_ctptau(i,j,k) = 0; modis_ctptau(i,j,k) = 0; } + for (int k = 0; k < m_num_cth; k++) { + misr_cthtau (i,j,k) = 0; + } } } } @@ -202,11 +246,13 @@ void Cosp::run_impl (const double dt) Kokkos::deep_copy(isccp_cldtot, 0.0); Kokkos::deep_copy(isccp_ctptau, 0.0); Kokkos::deep_copy(modis_ctptau, 0.0); + Kokkos::deep_copy(misr_cthtau , 0.0); Kokkos::deep_copy(isccp_mask , 0.0); } get_field_out("isccp_cldtot").sync_to_dev(); get_field_out("isccp_ctptau").sync_to_dev(); get_field_out("modis_ctptau").sync_to_dev(); + get_field_out("misr_cthtau" ).sync_to_dev(); get_field_out("isccp_mask" ).sync_to_dev(); } diff --git a/components/eamxx/src/physics/cosp/eamxx_cosp.hpp b/components/eamxx/src/physics/cosp/eamxx_cosp.hpp index 99f3c179a36..f7ff41d5a2d 100644 --- a/components/eamxx/src/physics/cosp/eamxx_cosp.hpp +++ b/components/eamxx/src/physics/cosp/eamxx_cosp.hpp @@ -2,6 +2,7 @@ #define SCREAM_COSP_HPP #include "share/atm_process/atmosphere_process.hpp" +#include "share/util/scream_common_physics_functions.hpp" #include "ekat/ekat_parameter_list.hpp" #include @@ -19,6 +20,10 @@ class Cosp : public AtmosphereProcess { public: + //using Pack = ekat::Pack; + using PF = scream::PhysicsFunctions; + using KT = KokkosTypes; + using KTH = KokkosTypes; // Constructors Cosp (const ekat::Comm& comm, const ekat::ParameterList& params); @@ -49,7 +54,12 @@ class Cosp : public AtmosphereProcess // The three main overrides for the subcomponent void initialize_impl (const RunType run_type); +#ifdef KOKKOS_ENABLE_CUDA + // Cuda requires methods enclosing __device__ lambda's to be public +public: +#endif void run_impl (const double dt); +protected: void finalize_impl (); // cosp frequency; positive is interpreted as number of steps, negative as number of hours @@ -60,8 +70,8 @@ class Cosp : public AtmosphereProcess Int m_num_cols; Int m_num_subcols; Int m_num_levs; - Int m_num_isccptau = 7; - Int m_num_isccpctp = 7; + Int m_num_tau = 7; + Int m_num_ctp = 7; Int m_num_cth = 16; std::shared_ptr m_grid;