Skip to content

Commit

Permalink
Enable MISR simulator and update cosp dimension names
Browse files Browse the repository at this point in the history
  • Loading branch information
@brhillman
brhillman committed May 21, 2024
1 parent 4c11ca5 commit 64c7c9a
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 32 deletions.
26 changes: 15 additions & 11 deletions components/eamxx/src/physics/cosp/cosp_c2f.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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., & !
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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')
Expand Down
25 changes: 15 additions & 10 deletions components/eamxx/src/physics/cosp/cosp_functions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {

Expand All @@ -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<const Real>& sunlit , view_1d<const Real>& skt,
view_2d<const Real>& T_mid , view_2d<const Real>& p_mid , view_2d<const Real>& p_int,
view_2d<const Real>& qv , view_2d<const Real>& cldfrac,
view_2d<const Real>& z_mid , view_2d<const Real>& qv , view_2d<const Real>& cldfrac,
view_2d<const Real>& reff_qc, view_2d<const Real>& reff_qi,
view_2d<const Real>& dtau067, view_2d<const Real>& dtau105,
view_1d<Real>& isccp_cldtot , view_3d<Real>& isccp_ctptau, view_3d<Real>& modis_ctptau) {
view_1d<Real>& isccp_cldtot , view_3d<Real>& isccp_ctptau, view_3d<Real>& modis_ctptau, view_3d<Real>& 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);
Expand All @@ -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++) {
Expand All @@ -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);
}
}
}
}
Expand Down
64 changes: 55 additions & 9 deletions components/eamxx/src/physics/cosp/eamxx_cosp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@ void Cosp::set_grids(const std::shared_ptr<const GridsManager> 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();
Expand All @@ -54,8 +56,11 @@ void Cosp::set_grids(const std::shared_ptr<const GridsManager> 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
Expand All @@ -68,6 +73,8 @@ void Cosp::set_grids(const std::shared_ptr<const GridsManager> grids_manager)
//add_field<Required>("height_mid", scalar3d_mid, m, grid_name);
//add_field<Required>("height_int", scalar3d_int, m, grid_name);
add_field<Required>("T_mid", scalar3d_mid, K, grid_name);
add_field<Required>("phis", scalar2d , m2/s2, grid_name);
add_field<Required>("pseudo_density", scalar3d_mid, Pa, grid_name);
add_field<Required>("qv", scalar3d_mid, Q, grid_name, "tracers");
add_field<Required>("qc", scalar3d_mid, Q, grid_name, "tracers");
add_field<Required>("qi", scalar3d_mid, Q, grid_name, "tracers");
Expand All @@ -86,6 +93,7 @@ void Cosp::set_grids(const std::shared_ptr<const GridsManager> grids_manager)
add_field<Computed>("isccp_cldtot", scalar2d, percent, grid_name);
add_field<Computed>("isccp_ctptau", scalar4d_ctptau, percent, grid_name, 1);
add_field<Computed>("modis_ctptau", scalar4d_ctptau, percent, grid_name, 1);
add_field<Computed>("misr_cthtau" , scalar4d_cthtau, percent, grid_name, 1);
add_field<Computed>("isccp_mask" , scalar2d, nondim, grid_name);
}

Expand All @@ -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::string,std::string>;
std::list<std::string> vnames = {"isccp_cldtot", "isccp_ctptau", "modis_ctptau"};
for (const auto field_name : {"isccp_cldtot", "isccp_ctptau", "modis_ctptau"}) {
std::list<std::string> 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<stratts_t>("io: string attributes");
atts["note"] = "Night values are zero; divide by isccp_mask to get daytime mean";
Expand Down Expand Up @@ -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<const Real**, Host>();
auto qc = get_field_in("qc").get_view<const Real**, Host>();
Expand All @@ -156,6 +166,8 @@ void Cosp::run_impl (const double dt)
auto T_mid = get_field_in("T_mid").get_view<const Real**, Host>();
auto p_mid = get_field_in("p_mid").get_view<const Real**, Host>();
auto p_int = get_field_in("p_int").get_view<const Real**, Host>();
auto phis = get_field_in("phis").get_view<const Real*, Host>();
auto pseudo_density = get_field_in("pseudo_density").get_view<const Real**, Host>();
auto cldfrac = get_field_in("cldfrac_rad").get_view<const Real**, Host>();
auto reff_qc = get_field_in("eff_radius_qc").get_view<const Real**, Host>();
auto reff_qi = get_field_in("eff_radius_qi").get_view<const Real**, Host>();
Expand All @@ -164,28 +176,60 @@ void Cosp::run_impl (const double dt)
auto isccp_cldtot = get_field_out("isccp_cldtot").get_view<Real*, Host>();
auto isccp_ctptau = get_field_out("isccp_ctptau").get_view<Real***, Host>();
auto modis_ctptau = get_field_out("modis_ctptau").get_view<Real***, Host>();
auto misr_cthtau = get_field_out("misr_cthtau" ).get_view<Real***, Host>();
auto isccp_mask = get_field_out("isccp_mask" ).get_view<Real*, Host>(); // Copy of sunlit flag with COSP frequency for proper averaging

// Compute heights
const auto z_mid = CospFunc::view_2d<Real>("z_mid", m_num_cols, m_num_levs);
const auto z_int = CospFunc::view_2d<Real>("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<KTH::ExeSpace>::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<const Real> 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;
}
}
}
}
Expand All @@ -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();
}

Expand Down
14 changes: 12 additions & 2 deletions components/eamxx/src/physics/cosp/eamxx_cosp.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <string>
Expand All @@ -19,6 +20,10 @@ class Cosp : public AtmosphereProcess
{

public:
//using Pack = ekat::Pack<Real,SCREAM_PACK_SIZE>;
using PF = scream::PhysicsFunctions<HostDevice>;
using KT = KokkosTypes<DefaultDevice>;
using KTH = KokkosTypes<HostDevice>;

// Constructors
Cosp (const ekat::Comm& comm, const ekat::ParameterList& params);
Expand Down Expand Up @@ -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
Expand All @@ -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<const AbstractGrid> m_grid;
Expand Down

0 comments on commit 64c7c9a

Please sign in to comment.