diff --git a/cime_config/machines/config_machines.xml b/cime_config/machines/config_machines.xml
index f9e08fb444e..3a852f4fe6c 100644
--- a/cime_config/machines/config_machines.xml
+++ b/cime_config/machines/config_machines.xml
@@ -443,6 +443,9 @@
1
+
+ $SHELL{if [ -z "$MOAB_ROOT" ]; then echo /global/cfs/cdirs/e3sm/software/moab/gnugpu ; else echo "$MOAB_ROOT"; fi}
+
$SHELL{if [ -z "$ADIOS2_ROOT" ]; then echo /global/cfs/cdirs/e3sm/3rdparty/adios2/2.9.1/cray-mpich-8.1.25/gcc-11.2.0; else echo "$ADIOS2_ROOT"; fi}
diff --git a/components/eam/src/dynamics/se/dyn_comp.F90 b/components/eam/src/dynamics/se/dyn_comp.F90
index b420be4ee69..faff6ffd314 100644
--- a/components/eam/src/dynamics/se/dyn_comp.F90
+++ b/components/eam/src/dynamics/se/dyn_comp.F90
@@ -269,7 +269,7 @@ subroutine dyn_init1(fh, NLFileName, dyn_in, dyn_out)
end if
#ifdef HAVE_MOAB
- call create_moab_meshes(par, elem)
+ call create_moab_meshes(par, elem, fv_nphys)
#endif
! Define the CAM grids (this has to be after dycore spinup).
! Physics-grid will be defined later by phys_grid_init
diff --git a/components/eamxx/src/control/atmosphere_driver.cpp b/components/eamxx/src/control/atmosphere_driver.cpp
index c3da1bacd3f..211b7536898 100644
--- a/components/eamxx/src/control/atmosphere_driver.cpp
+++ b/components/eamxx/src/control/atmosphere_driver.cpp
@@ -300,6 +300,9 @@ void AtmosphereDriver::create_grids()
void AtmosphereDriver::setup_surface_coupling_data_manager(SurfaceCouplingTransferType transfer_type,
const int num_cpl_fields, const int num_scream_fields,
const int field_size, Real* data_ptr,
+#ifdef HAVE_MOAB
+ Real* data_ptr_moab,
+#endif
char* names_ptr, int* cpl_indices_ptr, int* vec_comps_ptr,
Real* constant_multiple_ptr, bool* do_transfer_during_init_ptr)
{
@@ -318,6 +321,9 @@ void AtmosphereDriver::setup_surface_coupling_data_manager(SurfaceCouplingTransf
} else EKAT_ERROR_MSG("Error! Unexpected SurfaceCouplingTransferType.");
sc_data_mgr->setup_internals(num_cpl_fields, num_scream_fields, field_size, data_ptr,
+#ifdef HAVE_MOAB
+ data_ptr_moab,
+#endif
names_ptr, cpl_indices_ptr, vec_comps_ptr,
constant_multiple_ptr, do_transfer_during_init_ptr);
}
diff --git a/components/eamxx/src/control/atmosphere_driver.hpp b/components/eamxx/src/control/atmosphere_driver.hpp
index d43d233f4e9..7e96225210d 100644
--- a/components/eamxx/src/control/atmosphere_driver.hpp
+++ b/components/eamxx/src/control/atmosphere_driver.hpp
@@ -87,6 +87,9 @@ class AtmosphereDriver
void setup_surface_coupling_data_manager(SurfaceCouplingTransferType transfer_type,
const int num_cpl_fields, const int num_scream_fields,
const int field_size, Real* data_ptr,
+#ifdef HAVE_MOAB
+ Real* data_ptr_moab,
+#endif
char* names_ptr, int* cpl_indices_ptr, int* vec_comps_ptr,
Real* constant_multiple_ptr, bool* do_transfer_during_init_ptr);
diff --git a/components/eamxx/src/control/atmosphere_surface_coupling_exporter.cpp b/components/eamxx/src/control/atmosphere_surface_coupling_exporter.cpp
index 7a0253ec908..9044df7e76c 100644
--- a/components/eamxx/src/control/atmosphere_surface_coupling_exporter.cpp
+++ b/components/eamxx/src/control/atmosphere_surface_coupling_exporter.cpp
@@ -137,6 +137,13 @@ void SurfaceCouplingExporter::setup_surface_coupling_data(const SCDataManager &s
m_num_cols, m_num_cpl_exports);
m_cpl_exports_view_d = Kokkos::create_mirror_view(DefaultDevice(), m_cpl_exports_view_h);
+#ifdef HAVE_MOAB
+ // The export data is of size num_cpl_exports,ncols. All other data is of size num_scream_exports
+ m_moab_cpl_exports_view_h = decltype(m_moab_cpl_exports_view_h) (sc_data_manager.get_field_data_moab_ptr(),
+ m_num_cpl_exports, m_num_cols);
+ m_moab_cpl_exports_view_d = Kokkos::create_mirror_view(DefaultDevice(), m_moab_cpl_exports_view_h);
+#endif
+
m_export_field_names = new name_t[m_num_scream_exports];
std::memcpy(m_export_field_names, sc_data_manager.get_field_name_ptr(), m_num_scream_exports*32*sizeof(char));
@@ -536,6 +543,12 @@ void SurfaceCouplingExporter::do_export_to_cpl(const bool called_during_initiali
// Any field not exported by scream, or not exported
// during initialization, is set to 0.0
Kokkos::deep_copy(m_cpl_exports_view_d, 0.0);
+#ifdef HAVE_MOAB
+ // Any field not exported by scream, or not exported
+ // during initialization, is set to 0.0
+ Kokkos::deep_copy(m_moab_cpl_exports_view_d, 0.0);
+ const auto moab_cpl_exports_view_d = m_moab_cpl_exports_view_d;
+#endif
const auto cpl_exports_view_d = m_cpl_exports_view_d;
const int num_exports = m_num_scream_exports;
const int num_cols = m_num_cols;
@@ -554,9 +567,27 @@ void SurfaceCouplingExporter::do_export_to_cpl(const bool called_during_initiali
cpl_exports_view_d(icol,info.cpl_indx) = info.constant_multiple*info.data[offset];
}
});
+#ifdef HAVE_MOAB
+ Kokkos::parallel_for(export_policy, KOKKOS_LAMBDA(const int& i) {
+ const int ifield = i / num_cols;
+ const int icol = i % num_cols;
+ const auto& info = col_info(ifield);
+ const auto offset = icol*info.col_stride + info.col_offset;
+ // if this is during initialization, check whether or not the field should be exported
+ bool do_export = (not called_during_initialization || info.transfer_during_initialization);
+ if (do_export) {
+ moab_cpl_exports_view_d(info.cpl_indx, icol) = info.constant_multiple*info.data[offset];
+ }
+ });
+#endif
// Deep copy fields from device to cpl host array
Kokkos::deep_copy(m_cpl_exports_view_h,m_cpl_exports_view_d);
+#ifdef HAVE_MOAB
+ // Deep copy fields from device to cpl host array
+ Kokkos::deep_copy(m_moab_cpl_exports_view_h,m_moab_cpl_exports_view_d);
+#endif
+
}
// =========================================================================================
void SurfaceCouplingExporter::finalize_impl()
diff --git a/components/eamxx/src/control/atmosphere_surface_coupling_exporter.hpp b/components/eamxx/src/control/atmosphere_surface_coupling_exporter.hpp
index d8ccf5862f3..ca7772fdf1b 100644
--- a/components/eamxx/src/control/atmosphere_surface_coupling_exporter.hpp
+++ b/components/eamxx/src/control/atmosphere_surface_coupling_exporter.hpp
@@ -141,6 +141,13 @@ class SurfaceCouplingExporter : public AtmosphereProcess
view_2d m_cpl_exports_view_d;
uview_2d m_cpl_exports_view_h;
+#ifdef HAVE_MOAB
+ // Views storing a 2d array with dims (num_fields, num_cols) for moab cpl export data.
+ // The field cols strides faster, since that's what moab does (so we can "view" the
+ // pointer to the whole a2x_am(:,:) array from Fortran)
+ view_2d m_moab_cpl_exports_view_d;
+ uview_2d m_moab_cpl_exports_view_h;
+#endif
// Array storing the field names for exports
name_t* m_export_field_names;
std::vector m_export_field_names_vector;
diff --git a/components/eamxx/src/control/atmosphere_surface_coupling_importer.cpp b/components/eamxx/src/control/atmosphere_surface_coupling_importer.cpp
index 7573993a318..ee3e21e7461 100644
--- a/components/eamxx/src/control/atmosphere_surface_coupling_importer.cpp
+++ b/components/eamxx/src/control/atmosphere_surface_coupling_importer.cpp
@@ -74,6 +74,13 @@ void SurfaceCouplingImporter::set_grids(const std::shared_ptrget_params().get("iop_srf_prop")) {
// Overwrite imports with data from IOP file
@@ -232,6 +263,9 @@ void SurfaceCouplingImporter::overwrite_iop_imports (const bool called_during_in
const auto& info_d = col_info_d(ifield);
const auto offset = icol*info_d.col_stride + info_d.col_offset;
info_d.data[offset] = col_val;
+#ifdef HAVE_MOAB
+ // TODO
+#endif
});
}
}
diff --git a/components/eamxx/src/control/atmosphere_surface_coupling_importer.hpp b/components/eamxx/src/control/atmosphere_surface_coupling_importer.hpp
index 3a34a8b2951..6cf04c4166d 100644
--- a/components/eamxx/src/control/atmosphere_surface_coupling_importer.hpp
+++ b/components/eamxx/src/control/atmosphere_surface_coupling_importer.hpp
@@ -83,6 +83,14 @@ class SurfaceCouplingImporter : public AtmosphereProcess
view_2d m_cpl_imports_view_d;
uview_2d m_cpl_imports_view_h;
+#ifdef HAVE_MOAB
+ // Views storing a 2d array with dims (num_fields,num_cols) for import data.
+ // The colums index strides faster, since that's what moab does (so we can "view" the
+ // pointer to the whole x2a_am(:,:) array from Fortran)
+ view_2d m_moab_cpl_imports_view_d;
+ uview_2d m_moab_cpl_imports_view_h;
+#endif
+
// Array storing the field names for imports
name_t* m_import_field_names;
diff --git a/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90 b/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90
index cc1194aed21..901cd18796e 100644
--- a/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90
+++ b/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90
@@ -549,6 +549,14 @@ subroutine phys_grid_init (pgN)
use gllfvremap_mod, only: gfr_init
use homme_context_mod, only: elem, par
use dimensions_mod, only: nelem, nelemd
+#ifdef HAVE_MOAB
+ use seq_comm_mct, only: MHID, MHFID ! id of homme moab coarse and fine applications
+ use seq_comm_mct, only: ATMID
+ use seq_comm_mct, only: mhpgid ! id of pgx moab application
+ use semoab_mod, only: create_moab_meshes
+ use iMOAB, only : iMOAB_RegisterApplication
+ use iso_c_binding
+#endif
!
! Input(s)
!
@@ -561,6 +569,10 @@ subroutine phys_grid_init (pgN)
character(2) :: str
type(pg_specs_t), pointer :: pg
+#ifdef HAVE_MOAB
+ integer :: ATM_ID1
+ character*32 appname
+#endif
pg => pg_specs(pgN)
if (pg%inited) then
@@ -612,6 +624,48 @@ subroutine phys_grid_init (pgN)
call compute_global_dofs (pg)
call compute_global_coords (pg)
call compute_global_area (pg)
+#ifdef HAVE_MOAB
+ if (pgN > 0) then
+ appname="HM_COARSE"//C_NULL_CHAR
+ ATM_ID1 = 120 !
+ ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, MHID)
+ if (ierr > 0 ) &
+ call abortmp('Error: cannot register moab app')
+ if(par%masterproc) then
+ write(iulog,*) " "
+ write(iulog,*) "register MOAB app:", trim(appname), " MHID=", MHID
+ write(iulog,*) " "
+ endif
+ appname="HM_FINE"//C_NULL_CHAR
+ ATM_ID1 = 119 ! this number should not conflict with other components IDs; how do we know?
+ ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, MHFID)
+ if (ierr > 0 ) &
+ call abortmp('Error: cannot register moab app for fine mesh')
+ if(par%masterproc) then
+ write(iulog,*) " "
+ write(iulog,*) "register MOAB app:", trim(appname), " MHFID=", MHFID
+ write(iulog,*) " "
+ endif
+ appname="HM_PGX"//C_NULL_CHAR
+ ATM_ID1 = ATMID(1) ! this number should not conflict with other components IDs; how do we know?
+ !
+ ! in this case, we reuse the main atm id, mhid will not be used for intersection anymore
+ ! still, need to be careful
+ ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, mhpgid)
+ if (ierr > 0 ) &
+ call abortmp('Error: cannot register moab app for fine mesh')
+ if(par%masterproc) then
+ write(iulog,*) " "
+ write(iulog,*) "register MOAB app:", trim(appname), " MHPGID=", mhpgid
+ write(iulog,*) " "
+ endif
+! instance distributed moab meshes from elem structures
+! 1 ) spectral coarse mesh
+! 2 ) GLL fine quad mesh (used mostly for visualization)
+! 3 ) pgN FV type mesh, (most of the time pg2 mesh), used for coupling with other components;
+ call create_moab_meshes(par, elem, pgN)
+ endif
+#endif
end subroutine phys_grid_init
diff --git a/components/eamxx/src/mct_coupling/atm_comp_mct.F90 b/components/eamxx/src/mct_coupling/atm_comp_mct.F90
index dc57d828bbd..5e1b762270b 100644
--- a/components/eamxx/src/mct_coupling/atm_comp_mct.F90
+++ b/components/eamxx/src/mct_coupling/atm_comp_mct.F90
@@ -17,6 +17,14 @@ module atm_comp_mct
! MPI
use mpi
+#ifdef HAVE_MOAB
+ use seq_comm_mct , only: mphaid ! atm physics grid id in MOAB, on atm pes
+ use iso_c_binding
+#ifdef MOABCOMP
+ use seq_comm_mct, only: seq_comm_compare_mb_mct
+#endif
+#endif
+
implicit none
save
private
@@ -45,6 +53,17 @@ module atm_comp_mct
integer :: atm_log_unit
+
+#ifdef HAVE_MOAB
+ ! to store all fields to be set in moab
+ integer , private :: mblsize, totalmbls, nsend, totalmbls_r, nrecv
+ real(r8) , allocatable, target, private :: a2x_am(:,:) ! atm to coupler, on atm mesh, on atm component pes
+ real(r8) , allocatable, target, private :: x2a_am(:,:) ! coupler to atm, on atm mesh, on atm component pes
+#ifdef MOABCOMP
+ integer :: mpicom_atm_moab ! used just for mpi-reducing the difference between moab tags and mct avs
+ integer :: rank2
+#endif
+#endif
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CONTAINS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -80,7 +99,14 @@ subroutine atm_init_mct( EClock, cdata, x2a, a2x, NLFilename )
use seq_comm_mct, only: seq_comm_inst, seq_comm_name, seq_comm_suffix
use shr_file_mod, only: shr_file_getunit, shr_file_setIO
use shr_sys_mod, only: shr_sys_abort
-
+#ifdef HAVE_MOAB
+ use shr_kind_mod , only: cxx=>shr_kind_cxx
+ use mct_mod, only: mct_avect_nRattr
+ use iMOAB, only: iMOAB_SetDoubleTagStorage
+#ifdef MOABDEBUG
+ use iMOAB, only: iMOAB_WriteMesh
+#endif
+#endif
! !INPUT/OUTPUT PARAMETERS:
type(ESMF_Clock) , intent(inout) :: EClock
type(seq_cdata) , intent(inout) :: cdata
@@ -100,14 +126,19 @@ subroutine atm_init_mct( EClock, cdata, x2a, a2x, NLFilename )
type(c_ptr) :: x2a_ptr, a2x_ptr
character(len=256) :: atm_log_fname ! name of ATM log file
character(CL) :: calendar ! calendar string
-
+#ifdef HAVE_MOAB
+ integer :: ent_type
+ character(CXX) :: tagname ! will store all seq_flds_a2x_fields
+#ifdef MOABDEBUG
+ character*100 outfile, wopts
+#endif
+#endif
! TODO: read this from the namelist?
character(len=256) :: yaml_fname = "./data/scream_input.yaml"
character(kind=c_char,len=256), target :: yaml_fname_c, atm_log_fname_c
character(len=256) :: caseid, username, hostname
character(kind=c_char,len=256), target :: caseid_c, username_c, hostname_c, calendar_c
logical (kind=c_bool) :: restarted_run
-
!-------------------------------------------------------------------------------
! Grab some data from the cdata structure (coming from the coupler)
@@ -178,6 +209,14 @@ subroutine atm_init_mct( EClock, cdata, x2a, a2x, NLFilename )
! Init MCT domain structure
call atm_domain_mct (lsize, gsMap_atm, dom_atm)
+#ifdef HAVE_MOAB
+ call moab_atm_phys_scream()
+#ifdef MOABCOMP
+ mpicom_atm_moab = mpicom_atm ! just to store it to be used later in MOABCOMP
+ rank2 = my_task ! again, just to use it for MOABCOMP output
+#endif
+#endif
+
! Init import/export mct attribute vectors
call mct_aVect_init(x2a, rList=seq_flds_x2a_fields, lsize=lsize)
call mct_aVect_init(a2x, rList=seq_flds_a2x_fields, lsize=lsize)
@@ -195,12 +234,38 @@ subroutine atm_init_mct( EClock, cdata, x2a, a2x, NLFilename )
! Init surface coupling stuff in the AD
call scream_set_cpl_indices (x2a, a2x)
+#ifdef HAVE_MOAB
+! for MOAB driver, need to setup surface coupling for MOAB arrays that will hold
+! fields data for import and export
+! setup is similar to setup for the attribute vectors real data arrays (x2a%rAttr
+! and a2x%rAttr)
+! the major difference is that moab arrays are transposed compared to AVs arrays
+ mblsize = lsize
+ nsend = mct_avect_nRattr(a2x)
+ totalmbls = mblsize * nsend ! size of the double array
+ allocate (a2x_am(mblsize, nsend) )
+ a2x_am = 0. ! initialize everything with 0
+
+ nrecv = mct_avect_nRattr(x2a)
+ totalmbls_r = mblsize * nrecv ! size of the double array used to receive
+ allocate (x2a_am(mblsize, nrecv) ) ! these will be received by moab tags, then used to set cam in surf data
+ x2a_am = 0. ! initialize everything with 0
+#endif
+
call scream_setup_surface_coupling (c_loc(import_field_names), c_loc(import_cpl_indices), &
- c_loc(x2a%rAttr), c_loc(import_vector_components), &
+ c_loc(x2a%rAttr), &
+#ifdef HAVE_MOAB
+ c_loc(x2a_am), &
+#endif
+ c_loc(import_vector_components), &
c_loc(import_constant_multiple), c_loc(do_import_during_init), &
num_cpl_imports, num_scream_imports, import_field_size, &
c_loc(export_field_names), c_loc(export_cpl_indices), &
- c_loc(a2x%rAttr), c_loc(export_vector_components), &
+ c_loc(a2x%rAttr), &
+#ifdef HAVE_MOAB
+ c_loc(a2x_am), &
+#endif
+ c_loc(export_vector_components), &
c_loc(export_constant_multiple), c_loc(do_export_during_init), &
num_cpl_exports, num_scream_exports, export_field_size)
@@ -208,6 +273,24 @@ subroutine atm_init_mct( EClock, cdata, x2a, a2x, NLFilename )
call string_f2c(trim(username),username_c)
call string_f2c(trim(hostname),hostname_c)
call scream_init_atm (caseid_c,hostname_c,username_c)
+#ifdef HAVE_MOAB
+ ! data should be set now inside moab from import and export fields
+ ! do we import and export or just export at init stage ?
+ tagname=trim(seq_flds_a2x_fields)//C_NULL_CHAR
+ ent_type = 0 ! vertices, point cloud
+ ierr = iMOAB_SetDoubleTagStorage ( mphaid, tagname, totalmbls , ent_type, a2x_am )
+ if (ierr /= 0) then
+ print *,'moab error in setting data in a2x_am '
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ end if
+#ifdef MOABDEBUG
+ outfile = 'AtmPhys_init.h5m'//C_NULL_CHAR
+ wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR
+ ierr = iMOAB_WriteMesh(mphaid, outfile, wopts)
+ if (ierr > 0 ) &
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+#endif
+#endif
end subroutine atm_init_mct
@@ -215,6 +298,27 @@ end subroutine atm_init_mct
subroutine atm_run_mct(EClock, cdata, x2a, a2x)
use iso_c_binding, only: c_double
use scream_f2c_mod, only: scream_run
+#ifdef HAVE_MOAB
+ use seq_flds_mod, only: seq_flds_a2x_fields, seq_flds_x2a_fields
+ use shr_kind_mod , only: cxx=>shr_kind_cxx
+ use iMOAB, only: iMOAB_SetDoubleTagStorage, iMOAB_GetDoubleTagStorage
+#ifdef MOABDEBUG
+ use iMOAB, only: iMOAB_WriteMesh
+#endif
+#ifdef MOABCOMP
+ use mct_mod
+ use seq_comm_mct, only : num_moab_exports
+#endif
+
+ integer :: ent_type
+#ifdef MOABCOMP
+ real(r8) :: difference
+ type(mct_list) :: temp_list
+ integer :: size_list, index_list
+ type(mct_string) :: mctOStr !
+ character(CXX) :: mct_field, modelStr
+#endif
+#endif
! !INPUT/OUTPUT PARAMETERS:
@@ -229,7 +333,14 @@ subroutine atm_run_mct(EClock, cdata, x2a, a2x)
type(mct_gGrid) , pointer :: ggrid
real(R8) :: nextsw_cday ! calendar of next atm sw
integer :: dt_scream
-
+#ifdef HAVE_MOAB
+ integer :: ierr
+ character(CXX) :: tagname ! will store all seq_flds_a2x_fields , seq_flds_x2a_fields
+#ifdef MOABDEBUG
+ integer :: cur_atm_stepno
+ character*100 :: outfile, wopts, lnum
+#endif
+#endif
!-------------------------------------------------------------------------------
call seq_cdata_setptrs(cdata, &
@@ -239,10 +350,47 @@ subroutine atm_run_mct(EClock, cdata, x2a, a2x)
! Get time step info
call seq_timemgr_EClockGetData (EClock, next_cday=nextsw_cday, dtime=dt_scream)
+#ifdef HAVE_MOAB
+
+#ifdef MOABCOMP
+ ! loop over all fields in seq_flds_x2a_fields
+ call mct_list_init(temp_list ,seq_flds_x2a_fields)
+ size_list=mct_list_nitem (temp_list)
+ ent_type = 0 ! entity type is vertex for phys atm
+ if (rank2 .eq. 0) print *, num_moab_exports, trim(seq_flds_x2a_fields)
+ modelStr ='atm run'
+ do index_list = 1, size_list
+ call mct_list_get(mctOStr,index_list,temp_list)
+ mct_field = mct_string_toChar(mctOStr)
+ tagname= trim(mct_field)//C_NULL_CHAR
+ call seq_comm_compare_mb_mct(modelStr, mpicom_atm_moab, x2a, mct_field, mphaid, tagname, ent_type, difference)
+ enddo
+ call mct_list_clean(temp_list)
+#endif
+
+
+ ! import data from moab data structures
+ tagname=trim(seq_flds_x2a_fields)//C_NULL_CHAR
+ ent_type = 0 ! vertices, point cloud
+ ierr = iMOAB_GetDoubleTagStorage ( mphaid, tagname, totalmbls_r , ent_type, x2a_am )
+#endif
! Run scream
call scream_run( dt_scream )
+#ifdef HAVE_MOAB
+ ! export data to moab data structures
+ tagname=trim(seq_flds_a2x_fields)//C_NULL_CHAR
+ ent_type = 0 ! vertices, point cloud
+ ierr = iMOAB_SetDoubleTagStorage ( mphaid, tagname, totalmbls , ent_type, a2x_am )
+#ifdef MOABDEBUG
+ call seq_timemgr_EClockGetData( EClock, stepno=cur_atm_stepno )
+ write(lnum,"(I0.2)")cur_atm_stepno
+ outfile = 'AtmPhys_'//trim(lnum)//'.h5m'//C_NULL_CHAR
+ wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR
+ ierr = iMOAB_WriteMesh(mphaid, outfile, wopts)
+#endif
+#endif
! Set time of next radiadtion computation
call seq_infodata_PutData(infodata, nextsw_cday=nextsw_cday)
@@ -356,4 +504,180 @@ subroutine atm_domain_mct( lsize, gsMap_atm, dom_atm )
call mct_gGrid_importRAttr(dom_atm,"aream",data1,lsize)
end subroutine atm_domain_mct
+#ifdef HAVE_MOAB
+! as part of initialization, moab version of the MCT type grid needs to be instanced too
+! it is what we call "point cloud mesh", it contains just the vertices corresponding to the
+! grid for surface coupling
+! its tags are corresponding to the fields from a2x and x2a attribute vectors, and its
+! GLOBAL_ID tag corresponds to the global degrees of freedom in global atmosphere
+ subroutine moab_atm_phys_scream()
+
+ use iMOAB, only : iMOAB_RegisterApplication, iMOAB_CreateVertices, iMOAB_WriteMesh, &
+ iMOAB_DefineTagStorage, iMOAB_SetIntTagStorage, iMOAB_SetDoubleTagStorage, &
+ iMOAB_ResolveSharedEntities, iMOAB_UpdateMeshInfo
+ use iso_c_binding, only: c_int, c_loc
+ use scream_f2c_mod, only: scream_get_num_global_cols, scream_get_num_local_cols, &
+ scream_get_local_cols_gids
+ use scream_f2c_mod, only: scream_get_cols_latlon, scream_get_cols_area
+ use seq_flds_mod, only: seq_flds_dom_fields
+ use shr_kind_mod , only: r8 => shr_kind_r8, cl=>shr_kind_cl, cxx=>shr_kind_cxx
+ use shr_const_mod, only: SHR_CONST_PI
+ use seq_flds_mod, only: seq_flds_a2x_fields, seq_flds_x2a_fields
+
+ ! Local variables
+ !
+ integer(kind=c_int) :: num_local_cols, num_global_cols
+#ifdef MOABDEBUG
+ character*100 outfile, wopts
+#endif
+ character(CXX) :: tagname ! will store all seq_flds_a2x_fields
+ character*32 appname
+ integer :: ATM_PHYS ! used as global identifier for iMOAB app on pphys grid atmosphere (200+ atm id)
+ !-------------------------------------------------------------------
+ !
+ ! Local Variables
+ !
+ integer :: ierr,mpi_ierr ! error codes
+ integer :: i ! for loops along dofs
+ integer :: tagtype, tagindex, numco, ent_type
+ real(r8), pointer :: data1(:)
+ real(r8), pointer :: data2(:) ! temporary
+ integer(kind=c_int), allocatable, target :: col_gids(:)
+ real(r8), allocatable, target :: moab_vert_coords(:)
+ real(r8) :: latv, lonv
+ !-------------------------------------------------------------------
+
+ appname="ATM_PHYS_SCREAM"//C_NULL_CHAR
+ ATM_PHYS = 200 + ATM_ID !
+ ierr = iMOAB_RegisterApplication(appname, mpicom_atm, ATM_PHYS, mphaid)
+ if (ierr > 0 ) then
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+ if (my_task == master_task) then
+ print *, " moab_atm_phys_scream:: register MOAB app:", trim(appname), " mphaid=", mphaid
+ endif
+ ! find global ids
+ num_local_cols = scream_get_num_local_cols()
+ num_global_cols = scream_get_num_global_cols()
+
+ allocate(col_gids(num_local_cols))
+
+ call scream_get_local_cols_gids(c_loc(col_gids))
+
+ ! Fill in correct values for domain components
+ allocate(moab_vert_coords(num_local_cols*3))
+ allocate( data1(num_local_cols), data2(num_local_cols) )
+ call scream_get_cols_latlon(c_loc(data1),c_loc(data2))
+ do i=1,num_local_cols
+ latv = data1(i) * SHR_CONST_PI/180
+ lonv = data2(i) * SHR_CONST_PI/180
+ moab_vert_coords(3*i-2)=COS(latv)*COS(lonv)
+ moab_vert_coords(3*i-1)=COS(latv)*SIN(lonv)
+ moab_vert_coords(3*i )=SIN(latv)
+ enddo
+ ierr = iMOAB_CreateVertices(mphaid, num_local_cols*3, 3, moab_vert_coords)
+ if (ierr > 0 ) then
+ print *, "Error: fail to create MOAB vertices in phys atm model"
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+
+ tagtype = 0 ! dense, integer
+ numco = 1
+ tagname='GLOBAL_ID'//C_NULL_CHAR
+ ierr = iMOAB_DefineTagStorage(mphaid, tagname, tagtype, numco, tagindex )
+ if (ierr > 0 ) then
+ print *, "Error: fail to define GLOBAL_ID tag in phys atm model"
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+ ent_type = 0 ! vertex type
+ ierr = iMOAB_SetIntTagStorage(mphaid, tagname, num_local_cols , ent_type, col_gids)
+ if (ierr > 0 ) then
+ print *, "Error: fail to set GLOBAL_ID tag in phys atm model"
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+ ierr = iMOAB_ResolveSharedEntities( mphaid, num_local_cols, col_gids)
+ if (ierr > 0 ) then
+ print *, "Error: fail to resolve shared ents in phys atm model"
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+ ierr = iMOAB_UpdateMeshInfo(mphaid)
+ if (ierr > 0 ) then
+ print *, "Error: fail to update mesh info in phys atm model"
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+ tagname=trim(seq_flds_dom_fields)//C_NULL_CHAR ! mask is double too lat:lon:hgt:area:aream:mask:frac
+ tagtype = 1 ! dense, double
+ ierr = iMOAB_DefineTagStorage(mphaid, tagname, tagtype, numco, tagindex )
+ if (ierr > 0 ) then
+ print *, 'Error: fail to create tags from seq_flds_dom_fields '
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+
+ call scream_get_cols_area(c_loc(data1))
+ tagname = 'area'//C_NULL_CHAR !
+ ierr = iMOAB_SetDoubleTagStorage(mphaid, tagname, num_local_cols , ent_type, data1)
+ if (ierr > 0 ) then
+ print *, 'Error: fail to set area '
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+
+ ! Mask and frac are both exactly 1
+ data1 = 1.0
+ tagname = 'mask'//C_NULL_CHAR !
+ ierr = iMOAB_SetDoubleTagStorage(mphaid, tagname, num_local_cols , ent_type, data1)
+ if (ierr > 0 ) then
+ print *, 'Error: fail to set mask '
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+ tagname = 'frac'//C_NULL_CHAR !
+ ierr = iMOAB_SetDoubleTagStorage(mphaid, tagname, num_local_cols , ent_type, data1)
+ if (ierr > 0 ) then
+ print *, 'Error: fail to set frac '
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+
+ !call mct_gGrid_importRAttr(dom_atm,"mask",data1,lsize)
+ !call mct_gGrid_importRAttr(dom_atm,"frac",data1,lsize)
+
+ ! Aream is computed by mct, so give invalid initial value
+ data1 = -9999.0_R8
+ tagname = 'aream'//C_NULL_CHAR !
+ ierr = iMOAB_SetDoubleTagStorage(mphaid, tagname, num_local_cols , ent_type, data1)
+ if (ierr > 0 ) then
+ print *, 'Error: fail to set aream '
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+
+#ifdef MOABDEBUG
+ outfile = 'AtmPhys.h5m'//C_NULL_CHAR
+ wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR
+ ierr = iMOAB_WriteMesh(mphaid, outfile, wopts)
+ if (ierr > 0 ) then
+ print *, "Error: fail to write PhysAtm mesh "
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+#endif
+ ! define fields seq_flds_a2x_fields
+ tagtype = 1 ! dense, double
+ numco = 1 ! one value per vertex / entity
+ tagname = trim(seq_flds_a2x_fields)//C_NULL_CHAR
+ ierr = iMOAB_DefineTagStorage(mphaid, tagname, tagtype, numco, tagindex )
+ if (ierr > 0 ) then
+ print *, "Error: fail to define seq_flds_a2x_fields "
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+ ! make sure this is defined too; it could have the same fields, but in different order, or really different
+ ! fields; need to make sure we have them
+ tagname = trim(seq_flds_x2a_fields)//C_NULL_CHAR
+ ierr = iMOAB_DefineTagStorage(mphaid, tagname, tagtype, numco, tagindex )
+ if (ierr > 0 ) then
+ print *, "Error: fail to define seq_flds_x2a_fields "
+ call mpi_abort(mpicom_atm,ierr,mpi_ierr)
+ endif
+ deallocate(col_gids)
+ deallocate(data1)
+ deallocate(data2)
+ end subroutine moab_atm_phys_scream
+#endif
+
end module atm_comp_mct
diff --git a/components/eamxx/src/mct_coupling/scream_cxx_f90_interface.cpp b/components/eamxx/src/mct_coupling/scream_cxx_f90_interface.cpp
index 3e1514b3865..8bf027e10bf 100644
--- a/components/eamxx/src/mct_coupling/scream_cxx_f90_interface.cpp
+++ b/components/eamxx/src/mct_coupling/scream_cxx_f90_interface.cpp
@@ -178,11 +178,19 @@ void scream_create_atm_instance (const MPI_Fint f_comm, const int atm_id,
}
void scream_setup_surface_coupling (const char*& import_field_names, int*& import_cpl_indices,
- double*& x2a_ptr, int*& import_vector_components,
+ double*& x2a_ptr,
+#ifdef HAVE_MOAB
+ double*& x2a_moab_ptr,
+#endif
+ int*& import_vector_components,
double*& import_constant_multiple, bool*& do_import_during_init,
const int& num_cpl_imports, const int& num_scream_imports, const int& import_field_size,
char*& export_field_names, int*& export_cpl_indices,
- double*& a2x_ptr, int*& export_vector_components,
+ double*& a2x_ptr,
+#ifdef HAVE_MOAB
+ double*& a2x_moab_ptr,
+#endif
+ int*& export_vector_components,
double*& export_constant_multiple, bool*& do_export_during_init,
const int& num_cpl_exports, const int& num_scream_exports, const int& export_field_size)
{
@@ -210,10 +218,16 @@ void scream_setup_surface_coupling (const char*& import_field_names, int*& impor
ad.setup_surface_coupling_data_manager(scream::SurfaceCouplingTransferType::Import,
num_cpl_imports, num_scream_imports, import_field_size, x2a_ptr,
+#ifdef HAVE_MOAB
+ x2a_moab_ptr,
+#endif
names_in[0], import_cpl_indices, import_vector_components,
import_constant_multiple, do_import_during_init);
ad.setup_surface_coupling_data_manager(scream::SurfaceCouplingTransferType::Export,
num_cpl_exports, num_scream_exports, export_field_size, a2x_ptr,
+#ifdef HAVE_MOAB
+ a2x_moab_ptr,
+#endif
names_out[0], export_cpl_indices, export_vector_components,
export_constant_multiple, do_export_during_init);
});
diff --git a/components/eamxx/src/mct_coupling/scream_f2c_mod.F90 b/components/eamxx/src/mct_coupling/scream_f2c_mod.F90
index f7a0bc9c379..74b45f8330b 100644
--- a/components/eamxx/src/mct_coupling/scream_f2c_mod.F90
+++ b/components/eamxx/src/mct_coupling/scream_f2c_mod.F90
@@ -47,11 +47,19 @@ end subroutine scream_get_cols_area
! responsible to handle import/export operation from/into the component
! coupler surface fluxes/state structures
subroutine scream_setup_surface_coupling (import_field_names, import_cpl_indices, &
- x2a_ptr, import_vector_components, &
+ x2a_ptr, &
+#ifdef HAVE_MOAB
+ x2a_moab_ptr, &
+#endif
+ import_vector_components, &
import_constant_multiple, do_import_during_init, &
num_cpl_imports, num_scream_imports, import_field_size, &
export_field_names, export_cpl_indices, &
- a2x_ptr, export_vector_components, &
+ a2x_ptr, &
+#ifdef HAVE_MOAB
+ a2x_moab_ptr, &
+#endif
+ export_vector_components, &
export_constant_multiple, do_export_during_init, &
num_cpl_exports, num_scream_exports, export_field_size) bind(c)
use iso_c_binding, only: c_ptr, c_int
@@ -64,6 +72,9 @@ subroutine scream_setup_surface_coupling (import_field_names, import_cpl_indices
type(c_ptr), intent(in) :: export_field_names, export_cpl_indices, &
a2x_ptr, export_vector_components, &
export_constant_multiple, do_export_during_init
+#ifdef HAVE_MOAB
+ type(c_ptr), intent(in) :: x2a_moab_ptr, a2x_moab_ptr
+#endif
integer(kind=c_int), intent(in) :: num_cpl_imports, num_scream_imports, &
num_cpl_exports, num_scream_exports
integer(kind=c_int), intent(in) :: import_field_size, export_field_size
diff --git a/components/eamxx/src/share/atm_process/SCDataManager.hpp b/components/eamxx/src/share/atm_process/SCDataManager.hpp
index 02af15423e6..76bc9e5bd84 100644
--- a/components/eamxx/src/share/atm_process/SCDataManager.hpp
+++ b/components/eamxx/src/share/atm_process/SCDataManager.hpp
@@ -22,7 +22,11 @@ struct SCDataManager {
~SCDataManager() = default;
void setup_internals (const int num_cpl_fields, const int num_scream_fields, const int field_size,
- Real* field_data_ptr, char* field_names, int* field_cpl_indices_ptr,
+ Real* field_data_ptr,
+#ifdef HAVE_MOAB
+ Real* field_data_moab_ptr,
+#endif
+ char* field_names, int* field_cpl_indices_ptr,
int* field_vector_components_ptr, Real* field_constant_multiple_ptr,
bool* transfer_during_init_ptr)
{
@@ -38,6 +42,9 @@ struct SCDataManager {
EKAT_ASSERT_MSG(field_constant_multiple_ptr !=nullptr, "Error! Ptr for constant multiple is null.");
EKAT_ASSERT_MSG(transfer_during_init_ptr !=nullptr, "Error! Ptr for initial transfer boolean is null.");
m_field_data = decltype(m_field_data) (field_data_ptr, m_field_size, m_num_cpl_fields);
+#ifdef HAVE_MOAB
+ m_field_data_moab = decltype(m_field_data_moab) (field_data_moab_ptr, m_num_cpl_fields, m_field_size);
+#endif
m_field_cpl_indices = decltype(m_field_cpl_indices) (field_cpl_indices_ptr, m_num_scream_fields);
m_field_vector_components = decltype(m_field_vector_components) (field_vector_components_ptr, m_num_scream_fields);
m_field_constant_multiple = decltype(m_field_constant_multiple) (field_constant_multiple_ptr, m_num_scream_fields);
@@ -65,7 +72,11 @@ struct SCDataManager {
Real* get_field_data_ptr () const {
return m_field_data.data();
}
-
+#ifdef HAVE_MOAB
+ Real* get_field_data_moab_ptr () const {
+ return m_field_data_moab.data();
+ }
+#endif
Real get_field_data_view_entry(const int i, const int f) {
return m_field_data(i, f);
}
@@ -105,6 +116,9 @@ struct SCDataManager {
int m_num_scream_fields;
view_2d m_field_data;
+#ifdef HAVE_MOAB
+ view_2d m_field_data_moab;
+#endif
name_t* m_field_names;
view_1d m_field_cpl_indices;
view_1d m_field_vector_components;
diff --git a/components/homme/src/share/semoab_mod.F90 b/components/homme/src/share/semoab_mod.F90
index 2d8bb099ade..4841cb84b54 100644
--- a/components/homme/src/share/semoab_mod.F90
+++ b/components/homme/src/share/semoab_mod.F90
@@ -23,7 +23,7 @@ module semoab_mod
use seq_comm_mct, only: MHPGID ! app id on moab side, for PGx style mesh, uniform from se
use seq_comm_mct, only: atm_pg_active ! turn it on when PG style mesh active
- use dyn_grid, only: fv_nphys ! phys grid mesh will be replicated too
+ !use dyn_grid, only: fv_nphys ! phys grid mesh will be replicated too
use gllfvremap_mod, only: gfr_f_get_corner_latlon
use control_mod, only : west, east, south, north ! 1, 2, 3, 4
@@ -76,7 +76,7 @@ integer function search_in(intarr, leng, value)
end function search_in
- subroutine create_moab_meshes(par, elem)
+ subroutine create_moab_meshes(par, elem, fv_nphys)
use ISO_C_BINDING
use iMOAB, only: iMOAB_CreateVertices, iMOAB_WriteMesh, iMOAB_CreateElements, &
@@ -86,6 +86,7 @@ subroutine create_moab_meshes(par, elem)
type (element_t), intent(inout) :: elem(:)
type (parallel_t) , intent(in) :: par
+ integer , intent(in) :: fv_nphys
integer ierr, i, j, ie, iv, block_ID, k, numvals
integer icol, irow, je, linx ! local indices in fine el connect
diff --git a/components/homme/src/theta-l_kokkos/CMakeLists.txt b/components/homme/src/theta-l_kokkos/CMakeLists.txt
index 739866960e3..a585090ea76 100644
--- a/components/homme/src/theta-l_kokkos/CMakeLists.txt
+++ b/components/homme/src/theta-l_kokkos/CMakeLists.txt
@@ -106,6 +106,7 @@ MACRO(THETAL_KOKKOS_SETUP)
${SRC_SHARE_DIR}/planar_mod.F90
${SRC_SHARE_DIR}/geometry_mod.F90
${SRC_SHARE_DIR}/planar_mesh_mod.F90
+ ${SRC_SHARE_DIR}/semoab_mod.F90
)
IF (NOT HOMME_ENABLE_COMPOSE)