From 20cee55304a04983ae38e3649761e09d280ac804 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Mon, 11 Aug 2025 13:31:32 -0700 Subject: [PATCH 1/8] Changes the error check because the field may include data for ghost elements --- components/elm/src/main/accumulMod.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/components/elm/src/main/accumulMod.F90 b/components/elm/src/main/accumulMod.F90 index 4b261c898378..dd8fd41dbe76 100644 --- a/components/elm/src/main/accumulMod.F90 +++ b/components/elm/src/main/accumulMod.F90 @@ -278,7 +278,7 @@ subroutine extract_accum_field_sl (name, field, nstep) beg = accum(nf)%beg1d end = accum(nf)%end1d - if (size(field,dim=1) /= end-beg+1) then + if (size(field,dim=1) < end-beg+1) then write(iulog,*)'ERROR in extract_accum_field for field ',accum(nf)%name write(iulog,*)'size of first dimension of field is ',& size(field,dim=1),' and should be ',end-beg+1 @@ -341,7 +341,7 @@ subroutine extract_accum_field_ml (name, field, nstep) numlev = accum(nf)%numlev beg = accum(nf)%beg1d end = accum(nf)%end1d - if (size(field,dim=1) /= end-beg+1) then + if (size(field,dim=1) < end-beg+1) then write(iulog,*)'ERROR in extract_accum_field for field ',accum(nf)%name write(iulog,*)'size of first dimension of field is ',& size(field,dim=1),' and should be ',end-beg+1 @@ -406,7 +406,7 @@ subroutine update_accum_field_sl (name, field, nstep) beg = accum(nf)%beg1d end = accum(nf)%end1d - if (size(field,dim=1) /= end-beg+1) then + if (size(field,dim=1) < end-beg+1) then write(iulog,*)'ERROR in UPDATE_ACCUM_FIELD_SL for field ',accum(nf)%name write(iulog,*)'size of first dimension of field is ',size(field,dim=1),& ' and should be ',end-beg+1 @@ -500,7 +500,7 @@ subroutine update_accum_field_ml (name, field, nstep) numlev = accum(nf)%numlev beg = accum(nf)%beg1d end = accum(nf)%end1d - if (size(field,dim=1) /= end-beg+1) then + if (size(field,dim=1) < end-beg+1) then write(iulog,*)'ERROR in UPDATE_ACCUM_FIELD_ML for field ',accum(nf)%name write(iulog,*)'size of first dimension of field is ',size(field,dim=1),& ' and should be ',end-beg+1 From 1ff566ec8a6c7ca6e9722176b78432c70def6749 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Mon, 11 Aug 2025 13:33:36 -0700 Subject: [PATCH 2/8] Initialize only the locally owned cells --- .../elm/src/data_types/ColumnDataType.F90 | 18 +++++++++--------- components/elm/src/main/elm_instMod.F90 | 6 +++--- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90 index 3833186b2429..bc0a16563aa2 100644 --- a/components/elm/src/data_types/ColumnDataType.F90 +++ b/components/elm/src/data_types/ColumnDataType.F90 @@ -1097,7 +1097,7 @@ module ColumnDataType !------------------------------------------------------------------------ ! Subroutines to initialize and clean column energy state data structure !------------------------------------------------------------------------ - subroutine col_es_init(this, begc, endc) + subroutine col_es_init(this, begc, endc, endc_owned) ! ! !USES: use landunit_varcon, only : istice, istwet, istsoil, istdlak, istice_mec @@ -1107,7 +1107,7 @@ subroutine col_es_init(this, begc, endc) ! ! !ARGUMENTS: class(column_energy_state) :: this - integer, intent(in) :: begc,endc + integer, intent(in) :: begc,endc, endc_owned !------------------------------------------------------------------------ ! ! !LOCAL VARIABLES: @@ -1223,7 +1223,7 @@ subroutine col_es_init(this, begc, endc) !----------------------------------------------------------------------- ! Initialize soil+snow temperatures - do c = begc,endc + do c = begc,endc_owned l = col_pp%landunit(c) ! Snow level temperatures - all land points @@ -1381,12 +1381,12 @@ end subroutine col_es_clean !------------------------------------------------------------------------ ! Subroutines to initialize and clean column water state data structure !------------------------------------------------------------------------ - subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_input) + subroutine col_ws_init(this, begc, endc, endc_owned, h2osno_input, snow_depth_input, watsat_input) ! use elm_varctl , only : use_lake_wat_storage, use_arctic_init ! !ARGUMENTS: class(column_water_state) :: this - integer , intent(in) :: begc,endc + integer , intent(in) :: begc,endc, endc_owned real(r8), intent(in) :: h2osno_input(begc:) real(r8), intent(in) :: snow_depth_input(begc:) real(r8), intent(in) :: watsat_input(begc:, 1:) ! volumetric soil water at saturation (porosity) @@ -1675,7 +1675,7 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ ! Arrays that are initialized from input arguments this%wslake_col(begc:endc) = 0._r8 - do c = begc,endc + do c = begc,endc_owned l = col_pp%landunit(c) this%h2osno(c) = h2osno_input(c) this%int_snow(c) = h2osno_input(c) @@ -5768,11 +5768,11 @@ end subroutine col_ef_clean !------------------------------------------------------------------------ ! Subroutines to initialize and clean column water flux data structure !------------------------------------------------------------------------ - subroutine col_wf_init(this, begc, endc) + subroutine col_wf_init(this, begc, endc, endc_owned) ! ! !ARGUMENTS: class(column_water_flux) :: this - integer, intent(in) :: begc,endc + integer, intent(in) :: begc,endc, endc_owned ! !LOCAL VARIABLES: integer :: l,c integer :: ncells @@ -6042,7 +6042,7 @@ subroutine col_wf_init(this, begc, endc) this%qflx_from_uphill(begc:endc) = 0._r8 this%qflx_to_downhill(begc:endc) = 0._r8 ! needed for CNNLeaching - do c = begc, endc + do c = begc, endc_owned l = col_pp%landunit(c) if (lun_pp%itype(l) == istsoil .or. lun_pp%itype(l) == istcrop) then this%qflx_drain(c) = 0._r8 diff --git a/components/elm/src/main/elm_instMod.F90 b/components/elm/src/main/elm_instMod.F90 index 2ce43f15f4a9..8fa4ecaae667 100644 --- a/components/elm/src/main/elm_instMod.F90 +++ b/components/elm/src/main/elm_instMod.F90 @@ -417,7 +417,7 @@ subroutine elm_inst_biogeophys(bounds_proc) call grc_es%Init(bounds_proc%begg_all, bounds_proc%endg_all) call lun_es%Init(bounds_proc%begl_all, bounds_proc%endl_all) - call col_es%Init(bounds_proc%begc_all, bounds_proc%endc_all) + call col_es%Init(bounds_proc%begc_all, bounds_proc%endc_all, bounds_proc%endc) call veg_es%Init(bounds_proc%begp_all, bounds_proc%endp_all) call canopystate_vars%init(bounds_proc) @@ -432,7 +432,7 @@ subroutine elm_inst_biogeophys(bounds_proc) call grc_ws%Init(bounds_proc%begg_all, bounds_proc%endg_all) call lun_ws%Init(bounds_proc%begl_all, bounds_proc%endl_all) - call col_ws%Init(bounds_proc%begc_all, bounds_proc%endc_all, & + call col_ws%Init(bounds_proc%begc_all, bounds_proc%endc_all, bounds_proc%endc, & h2osno_col(begc:endc), & snow_depth_col(begc:endc), & soilstate_vars%watsat_col(begc:endc, 1:)) @@ -441,7 +441,7 @@ subroutine elm_inst_biogeophys(bounds_proc) call waterflux_vars%init(bounds_proc) call grc_wf%Init(bounds_proc%begg_all, bounds_proc%endg_all, bounds_proc) - call col_wf%Init(bounds_proc%begc_all, bounds_proc%endc_all) + call col_wf%Init(bounds_proc%begc_all, bounds_proc%endc_all, bounds_proc%endc) call veg_wf%Init(bounds_proc%begp_all, bounds_proc%endp_all) call chemstate_vars%Init(bounds_proc) From b46458cc6af45deb33bccfd4ddf9276d9386924a Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 13 Aug 2025 10:19:33 -0700 Subject: [PATCH 3/8] Adds code to do MOAB-based halo exchange --- components/elm/src/utils/domainLateralMod.F90 | 299 ++++++++++++++++++ 1 file changed, 299 insertions(+) diff --git a/components/elm/src/utils/domainLateralMod.F90 b/components/elm/src/utils/domainLateralMod.F90 index b46299a13661..a8535152551a 100644 --- a/components/elm/src/utils/domainLateralMod.F90 +++ b/components/elm/src/utils/domainLateralMod.F90 @@ -545,6 +545,303 @@ subroutine ExchangeColumnLevelGhostData(bounds_proc, nvals_per_col, & end subroutine ExchangeColumnLevelGhostData +#else + +#ifdef HAVE_MOAB + + !----------------------------------------------------------------------- + ! This is a stub for the case when PETSc is unavailable + ! + use shr_kind_mod, only : r8 => shr_kind_r8 + use shr_sys_mod , only : shr_sys_abort + use spmdMod , only : masterproc + use elm_varctl , only : iulog + use spmdMod , only : masterproc, iam, npes, mpicom, comp_id + use abortutils , only : endrun + use MOABGridType, only : moab_gcell, mlndghostid + ! + ! !PUBLIC TYPES: + implicit none + private + ! + + type, public :: oneD_int_data_for_moab + integer :: moab_app_id ! ID of MAOB app + character(len=1024) :: tag_name ! MOAB tag name + integer :: tag_type ! type of MOAB tag: 0 = dense, int; 1 = dense, double + integer :: num_tags ! Number of tags + integer :: entity_type(1) ! vertex or element based type + integer :: tag_index(1) ! Index of tag after it is registered in MOAB + integer :: num_comp ! number of components + integer :: ngcells ! number of grid cells + integer, allocatable :: values(:) ! data + end type oneD_int_data_for_moab + + type, public :: twoD_real_data_for_moab + integer :: moab_app_id ! ID of MAOB app + character(len=1024) :: tag_name ! MOAB tag name + integer :: tag_type ! type of MOAB tag: 0 = dense, int; 1 = dense, double + integer :: num_tags ! Number of tags + integer :: entity_type(1) ! vertex or element based type + integer :: tag_index(1) ! Index of tag after it is registered in MOAB + integer :: num_comp ! number of components + integer :: ngcells ! number of grid cells + real(r8), allocatable :: values(:,:) ! data + end type twoD_real_data_for_moab + + type, public :: domainlateral_type + type(oneD_int_data_for_moab) :: grid_level_count + type(twoD_real_data_for_moab) :: soil_lyr_data_real + end type domainlateral_type + + type(domainlateral_type) , public :: ldomain_lateral + ! + ! !PUBLIC MEMBER FUNCTIONS: + public domainlateral_init ! initializes + public GridLevelIntegerDataHaloExchange + public GridLevelSoilLayerDataHaloExchange ! + ! + !EOP + !------------------------------------------------------------------------------ + +contains + + !------------------------------------------------------------------------------ + subroutine setup_oneD_int_data_for_moab(moab_app_id, tag_name, num_cells_ghosted, data) + ! + ! DESCRIPTION: + ! Sets up 1D integer-type data structure for performing MOAB-based halo exchange + ! of data. This supports a single value per grid cell to be exchanged. + ! + use iso_c_binding + use iMOAB, only : iMOAB_DefineTagStorage + ! + implicit none + ! + ! ARGUMENTS: + integer , intent(in) :: moab_app_id + character(len=*) , intent(in) :: tag_name + integer , intent(in) :: num_cells_ghosted + type(oneD_int_data_for_moab) , intent(out) :: data + ! + ! LOCAL VARIABLES: + integer :: ierr + + data%moab_app_id = moab_app_id + data%tag_name = trim(tag_name) // C_NULL_CHAR ! name + data%tag_type = 0 ! 0 = dense, int + data%num_tags = 1 ! a single tag + data%entity_type(1) = 1 ! element (== cell) based data + data%num_comp = 1 ! number components in the tag + data%ngcells = num_cells_ghosted ! number of grid cells + + ! allocate memory + allocate(data%values(data%ngcells)) + + ! define the tag in MOAB + ierr = iMOAB_DefineTagStorage(data%moab_app_id, data%tag_name, data%tag_type, data%num_comp, data%tag_index(1)) + + end subroutine setup_oneD_int_data_for_moab + + !------------------------------------------------------------------------------ + subroutine setup_twoD_real_data_for_moab(moab_app_id, tag_name, num_comp, num_cells_ghosted, data) + ! + ! DESCRIPTION: + ! Sets up 2D real-type data structure for performing MOAB-based halo exchange + ! of data. This supports 'num_comp' values per grid cell to be exchanged. + ! + use iso_c_binding + use iMOAB, only : iMOAB_DefineTagStorage + ! + ! ARGUMENT: + integer , intent(in) :: moab_app_id + character(len=*) , intent(in) :: tag_name + integer , intent(in) :: num_comp + integer , intent(in) :: num_cells_ghosted + type(twoD_real_data_for_moab) , intent(out) :: data + ! + ! LOCAL VARIABLES: + integer :: ierr + + data%moab_app_id = moab_app_id + data%tag_name = trim(tag_name) // C_NULL_CHAR ! name + data%tag_type = 1 ! 1 = dense, double + data%num_tags = 1 ! a single tag + data%entity_type(1) = 1 ! element (== cell) based data + data%num_comp = num_comp ! number components in the tag + data%ngcells = num_cells_ghosted ! number of grid cells + + ! allocate memory + allocate(data%values(data%num_comp, data%ngcells)) + + ! define the tag in MOAB + ierr = iMOAB_DefineTagStorage(data%moab_app_id, data%tag_name, data%tag_type, data%num_comp, data%tag_index(1)) + + end subroutine setup_twoD_real_data_for_moab + + !------------------------------------------------------------------------------ + subroutine domainlateral_init(domain_l) + ! + ! DESCRIPTION: + ! Creates data structure for doing halo exchanges using MOAB + ! + use elm_varpar, only : nlevgrnd + ! + implicit none + ! + ! ARGUMENTS: + type(domainlateral_type) :: domain_l ! domain datatype + + ! creates the MOAB tag 1D data at grid level + call setup_oneD_int_data_for_moab(mlndghostid, 'grid_level_count', moab_gcell%num_ghosted, domain_l%grid_level_count) + + ! creates the MOAB tag for exchanging vertically distributed soil dataset + call setup_twoD_real_data_for_moab(mlndghostid, 'soil_data', nlevgrnd, moab_gcell%num_ghosted, domain_l%soil_lyr_data_real) + + end subroutine domainlateral_init + + !------------------------------------------------------------------------------ + subroutine do_haloexchange_oneD_integer_data_for_moab(data) + ! + ! DESCRIPTION: + ! Perform MOAB-based halo exchange + ! + use iMOAB, only : iMOAB_SetIntTagStorage, iMOAB_GetIntTagStorage, iMOAB_SynchronizeTags + ! + ! ARGUMENT: + type(oneD_int_data_for_moab) , intent(inout) :: data + ! + ! LOCAL VARIABLE: + integer :: ierr + + ! set the data in MOAB tag + ierr = iMOAB_SetIntTagStorage(data%moab_app_id, data%tag_name, data%ngcells * data%num_comp, data%entity_type(1), data%values) + if (ierr > 0) call endrun('Error: setting values in MOAB tag failed.') + + ! do the halo-exchange + ierr = iMOAB_SynchronizeTags(data%moab_app_id, data%num_tags, data%tag_index(1), data%entity_type(1)) + if (ierr > 0) call endrun('Error: synchronization of MOAB tag failed.') + + ! get the data from MOAB tag + ierr = iMOAB_GetIntTagStorage(data%moab_app_id, data%tag_name, data%ngcells * data%num_comp, data%entity_type(1), data%values) + if (ierr > 0) call endrun('Error: setting values in MOAB tag failed.') + + end subroutine do_haloexchange_oneD_integer_data_for_moab + + !------------------------------------------------------------------------------ + subroutine do_haloexchange_twoD_real_data_for_moab(data) + ! + ! DESCRIPTION: + ! Perform MOAB-based halo exchange + ! + use iMOAB, only : iMOAB_SetDoubleTagStorage, iMOAB_GetDoubleTagStorage, iMOAB_SynchronizeTags + ! + ! INPUT ARGUMENT: + type(twoD_real_data_for_moab) , intent(inout) :: data + ! + ! LOCAL VARIABLE: + integer :: ierr + + ! set the data in MOAB tag + ierr = iMOAB_SetDoubleTagStorage(data%moab_app_id, data%tag_name, data%ngcells * data%num_comp, data%entity_type(1), data%values) + if (ierr > 0) call endrun('Error: setting values in MOAB tag failed.') + + ! do the halo-exchange + ierr = iMOAB_SynchronizeTags(data%moab_app_id, data%num_tags, data%tag_index(1), data%entity_type(1)) + if (ierr > 0) call endrun('Error: synchronization of MOAB tag failed.') + + ! get the data from MOAB tag + ierr = iMOAB_GetDoubleTagStorage(data%moab_app_id, data%tag_name, data%ngcells * data%num_comp, data%entity_type(1), data%values) + if (ierr > 0) call endrun('Error: setting values in MOAB tag failed.') + + end subroutine do_haloexchange_twoD_real_data_for_moab + + !------------------------------------------------------------------------------ + subroutine GridLevelIntegerDataHaloExchange(domain_l, begg, endg_owned, endg_all, elm_data) + ! + ! DESCRIPTION: + ! Performs halo exchange of integer data. It is assumed that there is only one value + ! per grid cell. elm_data has data in ELM-format such that owned grid cells at the beginning + ! followed by ghost grid cells. After MOAB-based halo exchange values are filled in + ! elm_data corresponding to ghost cells. + ! + implicit none + ! + ! ARGUMENTS: + type(domainlateral_type) :: domain_l ! domain datatype + integer, intent(in) :: begg ! beginning index of grid cell + integer, intent(in) :: endg_owned ! ending index for owned grid cells + integer, intent(in) :: endg_all ! ending index for all (owned + ghost) grid cells + integer, intent(inout) , pointer :: elm_data(:) ! data packed in ELM's format + ! + ! LOCAL VARAIBLES: + integer :: g, j, idx + integer :: ierr + + ! convert data from ELM format to MOAB format + do g = begg, endg_owned + idx = moab_gcell%elm2moab(g) + domain_l%grid_level_count%values(idx) = elm_data(g) + end do + + ! perform halo exchange + call do_haloexchange_oneD_integer_data_for_moab(domain_l%grid_level_count) + + ! convert data from MOAB format to ELM format + do idx = 1, moab_gcell%num_ghosted + if (.not.moab_gcell%is_owned(idx)) then + g = moab_gcell%moab2elm(idx) + elm_data(g) = domain_l%grid_level_count%values(idx) + end if + end do + + end subroutine GridLevelIntegerDataHaloExchange + + !------------------------------------------------------------------------------ + subroutine GridLevelSoilLayerDataHaloExchange(domain_l, begg, endg_owned, endg_all, elm_data) + ! + ! DESCRIPTION: + ! Performs halo exchange of real data. It is assumed that there are nlevgrnd values + ! per grid cell. elm_data has data in ELM-format such that owned grid cells at the beginning + ! followed by ghost grid cells. After MOAB-based halo exchange values are filled in + ! elm_data corresponding to ghost cells. + ! + ! + ! !ARGUMENTS: + implicit none + ! + type(domainlateral_type) :: domain_l ! domain datatype + integer, intent(in) :: begg ! beginning index of grid cell + integer, intent(in) :: endg_owned ! ending index for owned grid cells + integer, intent(in) :: endg_all ! ending index for all (owned + ghost) grid cells + real(r8), intent(inout) , pointer :: elm_data(:,:) ! data packed in ELM's format + ! + integer :: g, j, idx + integer :: ierr + + ! convert data from ELM format to MOAB format + do g = begg, endg_owned + idx = moab_gcell%elm2moab(g) + do j = 1, domain_l%soil_lyr_data_real%num_comp + domain_l%soil_lyr_data_real%values(j, idx) = elm_data(g, j) + end do + end do + + ! perform halo exchange + call do_haloexchange_twoD_real_data_for_moab(domain_l%soil_lyr_data_real) + + ! convert data from MOAB format to ELM format + do idx = 1, moab_gcell%num_ghosted + if (.not.moab_gcell%is_owned(idx)) then + g = moab_gcell%moab2elm(idx) + do j = 1, domain_l%soil_lyr_data_real%num_comp + elm_data(g, j) = domain_l%soil_lyr_data_real%values(j, idx) + end do + end if + end do + + end subroutine GridLevelSoilLayerDataHaloExchange + #else !----------------------------------------------------------------------- @@ -610,4 +907,6 @@ end subroutine domainlateral_init #endif +#endif + end module domainLateralMod From 12b4e345c3924031fad7dc94093b8965a4236433 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 13 Aug 2025 11:40:55 -0700 Subject: [PATCH 4/8] Adds code to initialize naturally vegetated ghost columns --- components/elm/src/main/initGridCellsMod.F90 | 55 ++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/components/elm/src/main/initGridCellsMod.F90 b/components/elm/src/main/initGridCellsMod.F90 index 2d8547334eea..ec90c0a05cad 100644 --- a/components/elm/src/main/initGridCellsMod.F90 +++ b/components/elm/src/main/initGridCellsMod.F90 @@ -715,6 +715,9 @@ subroutine initGhostGridcells() call CheckGhostSubgridHierarchy() #endif +#ifdef HAVE_MOAB + call initGhostColumnsMOAB() +#endif end subroutine initGhostGridcells #ifdef USE_PETSC_LIB @@ -1326,4 +1329,56 @@ end subroutine CheckGhostSubgridHierarchy #endif !^ifdef USE_PETSC_LIB + +#if HAVE_MOAB + !------------------------------------------------------------------------ + subroutine initGhostColumnsMOAB() + ! + ! DESCRIPTION + ! + use decompMod , only : get_proc_bounds + use domainLateralMod , only : ldomain_lateral, GridLevelIntegerDataHaloExchange + ! + implicit none + ! + ! LOCAL VARIABLES + type(bounds_type) :: bounds_proc ! temporary + integer, pointer :: num_nat_veg_columns(:) + integer :: g, c + integer, parameter :: icol_nat_veg = 1 + + call get_proc_bounds(bounds_proc) + + allocate(num_nat_veg_columns(bounds_proc%begg:bounds_proc%endg_all)) + num_nat_veg_columns(:) = 0 + + ! for owned grid cells, save the number of naturally vegetated soil columns + do c = bounds_proc%begc, bounds_proc%endc + g = col_pp%gridcell(c) + + if (col_pp%itype(c) == icol_nat_veg) then + if (num_nat_veg_columns(g) /= 0) then + call endrun(msg='ERROR: initGhostColumnsMOAB only supports the one natural soil column per grid cell') + else + num_nat_veg_columns(g) = num_nat_veg_columns(g) + 1 + end if + end if + + end do + + ! do the halo exchange + call GridLevelIntegerDataHaloExchange(ldomain_lateral, bounds_proc%begg, bounds_proc%endg, bounds_proc%endg_all, num_nat_veg_columns) + + ! for ghost columns, set mapping to ghost grid cell + c = bounds_proc%endc + do g = bounds_proc%endg + 1, bounds_proc%endg_all + if (num_nat_veg_columns(g) > 0) then + c = c + 1 + col_pp%gridcell(c) = g + col_pp%itype(c) = icol_nat_veg + endif + end do + + end subroutine initGhostColumnsMOAB +#endif end module initGridCellsMod From 86d586c7c61016a33f6e8565cdc2102d7099ed54 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 13 Aug 2025 13:49:18 -0700 Subject: [PATCH 5/8] Perform halo exhcnage for soil properties --- .../elm/src/biogeophys/SoilStateType.F90 | 115 +++++++++++++++++- 1 file changed, 113 insertions(+), 2 deletions(-) diff --git a/components/elm/src/biogeophys/SoilStateType.F90 b/components/elm/src/biogeophys/SoilStateType.F90 index 79ab33fcb8ea..46bf0cd94336 100644 --- a/components/elm/src/biogeophys/SoilStateType.F90 +++ b/components/elm/src/biogeophys/SoilStateType.F90 @@ -6,7 +6,7 @@ module SoilStateType use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type use abortutils , only : endrun - use spmdMod , only : mpicom, MPI_INTEGER, masterproc + use spmdMod , only : mpicom, MPI_INTEGER, masterproc, iam use ncdio_pio , only : file_desc_t, ncd_defvar, ncd_io, ncd_double, ncd_int, ncd_inqvdlen use ncdio_pio , only : ncd_pio_openfile, ncd_inqfdims, ncd_pio_closefile, ncd_inqdid, ncd_inqdlen use elm_varpar , only : more_vertlayers, numpft, numrad @@ -117,6 +117,10 @@ subroutine Init(this, bounds) call this%InitHistory(bounds) call this%InitCold(bounds) +#ifdef HAVE_MOAB + call this%InitColdGhost(bounds) +#endif + end subroutine Init !------------------------------------------------------------------------ @@ -163,7 +167,7 @@ subroutine InitAllocate(this, bounds) allocate(this%watopt_col (begc:endc,nlevgrnd)) ; this%watopt_col (:,:) = spval allocate(this%watfc_col (begc:endc,nlevgrnd)) ; this%watfc_col (:,:) = spval allocate(this%watmin_col (begc:endc,nlevgrnd)) ; this%watmin_col (:,:) = spval - allocate(this%sucsat_col (begc:endc,nlevgrnd)) ; this%sucsat_col (:,:) = spval + allocate(this%sucsat_col (begc_all:endc_all,nlevgrnd)) ; this%sucsat_col (:,:) = spval allocate(this%sucmin_col (begc:endc,nlevgrnd)) ; this%sucmin_col (:,:) = spval allocate(this%soilbeta_col (begc:endc)) ; this%soilbeta_col (:) = spval allocate(this%soilalpha_col (begc:endc)) ; this%soilalpha_col (:) = spval @@ -1074,6 +1078,112 @@ subroutine InitColdGhost(this, bounds_proc) end subroutine InitColdGhost +#else + +#ifdef HAVE_MOAB + + !------------------------------------------------------------------------ + subroutine PackOwnedGridLevelDataForMOAB(bounds_proc, col_itype, data_c_in, data_g_out) + ! + implicit none + ! + type(bounds_type) , intent(in) :: bounds_proc + integer , intent(in) :: col_itype + real(r8), pointer , intent(in) :: data_c_in(:,:) + real(r8), pointer , intent(inout) :: data_g_out(:,:) + ! + integer :: c, g, j + + data_g_out(:,:) = 0._r8 + + do c = bounds_proc%begc, bounds_proc%endc + if (col_pp%itype(c) == col_itype) then + g = col_pp%gridcell(c) + do j = 1, nlevgrnd + data_g_out(g, j) = data_c_in(c, j) + end do + end if + end do + + end subroutine PackOwnedGridLevelDataForMOAB + + !------------------------------------------------------------------------ + subroutine UnpackGhostGridLevelDataFromMOAB(bounds_proc, col_itype, data_g_in, data_c_out) + ! + implicit none + ! + type(bounds_type) , intent(in) :: bounds_proc + integer , intent(in) :: col_itype + real(r8) , pointer, intent(in) :: data_g_in(:,:) + real(r8) , pointer, intent(inout) :: data_c_out(:,:) + ! + integer :: c, g, j + + do c = bounds_proc%endc + 1, bounds_proc%endc_all + if (col_pp%itype(c) == col_itype) then + g = col_pp%gridcell(c) + do j = 1, nlevgrnd + data_c_out(c, j) = data_g_in(g, j) + end do + end if + end do + + end subroutine UnpackGhostGridLevelDataFromMOAB + + !------------------------------------------------------------------------ + subroutine ExchangeAFieldUsingMOAB(bounds_proc, col_itype, field_col) + ! + use domainLateralMod , only : ldomain_lateral, GridLevelSoilLayerDataHaloExchange + ! + implicit none + ! + ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds_proc + integer , intent(in) :: col_itype + real(r8), pointer , intent(inout) :: field_col(:,:) + ! + real(r8), pointer :: data(:,:) + + ! allocate memory for owned+ghost cells + allocate(data(bounds_proc%begg:bounds_proc%endg_all, nlevgrnd)) + + ! pack data + call PackOwnedGridLevelDataForMOAB(bounds_proc, col_itype, field_col, data) + + ! do the halo exchange + call GridLevelSoilLayerDataHaloExchange(ldomain_lateral, bounds_proc%begg, bounds_proc%endg, bounds_proc%endg_all, data) + + ! unpack data + call UnpackGhostGridLevelDataFromMOAB(bounds_proc, col_itype, data, field_col) + + ! free memory + deallocate(data) + + end subroutine ExchangeAFieldUsingMOAB + + !------------------------------------------------------------------------ + subroutine InitColdGhost(this, bounds_proc) + ! + ! !DESCRIPTION: + ! Assign soil properties for ghost/halo columns + ! + ! !USES: + ! + implicit none + ! + ! !ARGUMENTS: + class(soilstate_type) :: this + type(bounds_type), intent(in) :: bounds_proc + ! + integer, parameter :: nat_veg_col_itype = 1 + + call ExchangeAFieldUsingMOAB(bounds_proc, nat_veg_col_itype, this%watsat_col) + call ExchangeAFieldUsingMOAB(bounds_proc, nat_veg_col_itype, this%hksat_col ) + call ExchangeAFieldUsingMOAB(bounds_proc, nat_veg_col_itype, this%bsw_col ) + call ExchangeAFieldUsingMOAB(bounds_proc, nat_veg_col_itype, this%sucsat_col) + + end subroutine InitColdGhost + #else !------------------------------------------------------------------------ @@ -1096,6 +1206,7 @@ subroutine InitColdGhost(this, bounds_proc) end subroutine InitColdGhost +#endif #endif !------------------------------------------------------------------------ From bc676381a91a763e12697c43fb60f55401ebcc35 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 3 Sep 2025 12:46:50 -0700 Subject: [PATCH 6/8] Initialize ghost grid cells --- components/elm/src/main/elm_initializeMod.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/components/elm/src/main/elm_initializeMod.F90 b/components/elm/src/main/elm_initializeMod.F90 index 43e124372c83..264e02fc267b 100755 --- a/components/elm/src/main/elm_initializeMod.F90 +++ b/components/elm/src/main/elm_initializeMod.F90 @@ -229,11 +229,15 @@ subroutine initialize1( ) 'Unsupported domain_decomp_type = ' // trim(domain_decomp_type)) end select +#ifdef HAVE_MOAB + call domainlateral_init(ldomain_lateral) +#else if (lateral_connectivity) then call domainlateral_init(ldomain_lateral, cellsOnCell, edgesOnCell, & nEdgesOnCell, areaCell, dcEdge, dvEdge, & nCells_loc, nEdges_loc, maxEdges) endif +#endif ! *** Get JUST gridcell processor bounds *** ! Remaining bounds (landunits, columns, patches) will be determined @@ -421,6 +425,7 @@ subroutine initialize1( ) ! This is needed here for the following call to decompInit_glcp call initGridCells() + call initGhostGridCells() ! Set global seg maps for gridcells, topounits, landlunits, columns and patches !if(max_topounits > 1) then From 53c63ee1d609a2e47a1f99d322806c3d9bde4832 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 3 Sep 2025 13:48:12 -0700 Subject: [PATCH 7/8] Modifies computation of beg/end indices for ghost subgrid units --- components/elm/src/main/decompInitMod.F90 | 93 ++++++++++++----------- 1 file changed, 47 insertions(+), 46 deletions(-) diff --git a/components/elm/src/main/decompInitMod.F90 b/components/elm/src/main/decompInitMod.F90 index c03856404056..58241b805724 100644 --- a/components/elm/src/main/decompInitMod.F90 +++ b/components/elm/src/main/decompInitMod.F90 @@ -2380,52 +2380,7 @@ subroutine decompInit_ghosts(glcmask) if (.not.lateral_connectivity) then - ! No ghost cells - procinfo%ncells_ghost = 0 - procinfo%ntunits_ghost = 0 - procinfo%nlunits_ghost = 0 - procinfo%ncols_ghost = 0 - procinfo%npfts_ghost = 0 - procinfo%nCohorts_ghost = 0 - - procinfo%begg_ghost = 0 - procinfo%begt_ghost = 0 - procinfo%begl_ghost = 0 - procinfo%begc_ghost = 0 - procinfo%begp_ghost = 0 - procinfo%begCohort_ghost = 0 - procinfo%endg_ghost = 0 - procinfo%endt_ghost = 0 - procinfo%endl_ghost = 0 - procinfo%endc_ghost = 0 - procinfo%endp_ghost = 0 - procinfo%endCohort_ghost = 0 - - ! All = local (as no ghost cells) - procinfo%ncells_all = procinfo%ncells - procinfo%ntunits_all = procinfo%ntunits - procinfo%nlunits_all = procinfo%nlunits - procinfo%ncols_all = procinfo%ncols - procinfo%npfts_all = procinfo%npfts - procinfo%nCohorts_all = procinfo%nCohorts - - procinfo%begg_all = procinfo%begg - procinfo%begt_all = procinfo%begt - procinfo%begl_all = procinfo%begl - procinfo%begc_all = procinfo%begc - procinfo%begp_all = procinfo%begp - procinfo%begCohort_all = procinfo%begCohort - procinfo%endg_all = procinfo%endg - procinfo%endt_all = procinfo%endt - procinfo%endl_all = procinfo%endl - procinfo%endc_all = procinfo%endc - procinfo%endp_all = procinfo%endp - procinfo%endCohort_all = procinfo%endCohort - - else - #if defined(HAVE_MOAB) - call get_proc_bounds(begg, endg) tagname = 'nsubgrid'//C_NULL_CHAR @@ -2528,7 +2483,53 @@ subroutine decompInit_ghosts(glcmask) procinfo%endp_all = procinfo%endp + procinfo%npfts_ghost procinfo%endCohort_all = procinfo%endCohort + procinfo%nCohorts_ghost -#elif defined(USE_PETSC_LIB) +#else + ! No ghost cells + procinfo%ncells_ghost = 0 + procinfo%ntunits_ghost = 0 + procinfo%nlunits_ghost = 0 + procinfo%ncols_ghost = 0 + procinfo%npfts_ghost = 0 + procinfo%nCohorts_ghost = 0 + + procinfo%begg_ghost = 0 + procinfo%begt_ghost = 0 + procinfo%begl_ghost = 0 + procinfo%begc_ghost = 0 + procinfo%begp_ghost = 0 + procinfo%begCohort_ghost = 0 + procinfo%endg_ghost = 0 + procinfo%endt_ghost = 0 + procinfo%endl_ghost = 0 + procinfo%endc_ghost = 0 + procinfo%endp_ghost = 0 + procinfo%endCohort_ghost = 0 + + ! All = local (as no ghost cells) + procinfo%ncells_all = procinfo%ncells + procinfo%ntunits_all = procinfo%ntunits + procinfo%nlunits_all = procinfo%nlunits + procinfo%ncols_all = procinfo%ncols + procinfo%npfts_all = procinfo%npfts + procinfo%nCohorts_all = procinfo%nCohorts + + procinfo%begg_all = procinfo%begg + procinfo%begt_all = procinfo%begt + procinfo%begl_all = procinfo%begl + procinfo%begc_all = procinfo%begc + procinfo%begp_all = procinfo%begp + procinfo%begCohort_all = procinfo%begCohort + procinfo%endg_all = procinfo%endg + procinfo%endt_all = procinfo%endt + procinfo%endl_all = procinfo%endl + procinfo%endc_all = procinfo%endc + procinfo%endp_all = procinfo%endp + procinfo%endCohort_all = procinfo%endCohort +#endif + + else + +#if defined(USE_PETSC_LIB) call get_proc_bounds(begg, endg) From 697bd7cf90c2e46c2af30bc5e7671feb4c61b0ea Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Thu, 9 Oct 2025 09:53:13 -0700 Subject: [PATCH 8/8] Adds column-to-column connections when using MOAB --- .../data_types/ColumnConnectionSetType.F90 | 137 ++++++++++++++++++ components/elm/src/main/elm_initializeMod.F90 | 4 + 2 files changed, 141 insertions(+) create mode 100644 components/elm/src/data_types/ColumnConnectionSetType.F90 diff --git a/components/elm/src/data_types/ColumnConnectionSetType.F90 b/components/elm/src/data_types/ColumnConnectionSetType.F90 new file mode 100644 index 000000000000..e2a2fdf5fee2 --- /dev/null +++ b/components/elm/src/data_types/ColumnConnectionSetType.F90 @@ -0,0 +1,137 @@ +module ColumnConnectionSetType + + + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_infnan_mod , only : isnan => shr_infnan_isnan, nan => shr_infnan_nan, assignment(=) + use decompMod , only : bounds_type + use abortutils , only : endrun + use ColumnType , only : col_pp + implicit none + save + public + + type, public :: col_connection_set_type + Integer :: nconn ! number of connections + Integer, pointer :: col_id_up(:) => null() ! list of ids of upwind cells + Integer, pointer :: col_id_dn(:) => null() ! list of ids of downwind cells + Integer, pointer :: grid_id_up(:) => null() ! list of ids of upwind cells + Integer, pointer :: grid_id_dn(:) => null() ! list of ids of downwind cells + integer, pointer :: grid_id_up_norder(:) => null() ! list of ids of upwind cells in natural order + integer, pointer :: grid_id_dn_norder(:) => null() ! list of ids of downwind cells in natural order + integer, pointer :: col_up_forder(:) => null() ! the order in which the lateral flux should be added for upwind cells + integer, pointer :: col_dn_forder(:) => null() ! the order in which the lateral flux should be added for downwind cells + Real(r8), pointer :: dist(:) => null() ! list of distance vectors + Real(r8), pointer :: face_length(:) => null() ! list of edge of faces normal to distance vectors + Real(r8), pointer :: uparea(:) => null() ! list of up cell areas of horizaontal faces + Real(r8), pointer :: downarea(:) => null() ! list of down cell areas of horizaontal faces + Real(r8), pointer :: dzg(:) => null() ! list of areas of dz between downwind and upwind cells + Real(r8), pointer :: facecos(:) => null() ! dot product of the cell face normal vector and cell centroid vector + Real(r8), pointer :: vertcos(:) => null() ! dot product of the cell face normal vector and cell centroid vector for vertical flux, the rank for vertcos + ! is from 1 to column size which is different from rank of lateral faces + contains +#ifdef HAVE_MOAB + procedure, public :: Init => InitViaMOAB +#endif + end type col_connection_set_type + + type (col_connection_set_type), public, target :: c2c_connections ! connection type + +contains + +#ifdef HAVE_MOAB + !------------------------------------------------------------------------ + subroutine InitViaMOAB(this, bounds_proc) + ! + use MOABGridType, only : moab_edge_internal, moab_gcell + use decompMod , only : bounds_type + ! + implicit none + ! + class (col_connection_set_type) :: this + type(bounds_type), intent(in) :: bounds_proc ! bound information at processor level + ! + integer :: g, g_up_moab, g_dn_moab, g_up_elm, g_dn_elm + integer :: c, c_up, c_dn + integer :: iconn, nconn + integer, parameter :: nat_veg_col_itype = 1 + integer, pointer :: nat_col_id(:) + + + ! allocate memory and initialize + allocate(nat_col_id(bounds_proc%begg_all:bounds_proc%endg_all)) + nat_col_id(:) = -1 + + ! loop over columns to determine the naturally-vegetated column for each grid cell. + do c = bounds_proc%begc_all, bounds_proc%endc_all + if (col_pp%itype(c) == nat_veg_col_itype) then + g = col_pp%gridcell(c) + + if (nat_col_id(g) /= -1) then + call endrun('ERROR: More than one naturally vegetated column found.') + end if + + nat_col_id(g) = c + end if + end do + + ! loop over grid level connections and determine number of column level connections + nconn = 0 + do iconn = 1, moab_edge_internal%num + g_up_moab = moab_edge_internal%cell_ids(iconn, 1) + g_dn_moab = moab_edge_internal%cell_ids(iconn, 2) + + g_up_elm = moab_gcell%moab2elm(g_up_moab) + g_dn_elm = moab_gcell%moab2elm(g_dn_moab) + + if (nat_col_id(g_up_elm) /= -1 .and. nat_col_id(g_dn_elm) /= -1) then + nconn = nconn + 1 + end if + end do + + ! allocate and initialize data structure + this%nconn = nconn + allocate(this%col_id_up(nconn)) ; this%col_id_up(:) = 0 + allocate(this%col_id_dn(nconn)) ; this%col_id_dn(:) = 0 + allocate(this%grid_id_up(nconn)) ; this%grid_id_up(:) = 0 + allocate(this%grid_id_dn(nconn)) ; this%grid_id_dn(:) = 0 + allocate(this%grid_id_up_norder(nconn)) ; this%grid_id_up_norder(:) = 0 + allocate(this%grid_id_dn_norder(nconn)) ; this%grid_id_dn_norder(:) = 0 + allocate(this%col_up_forder(nconn)) ; this%col_up_forder(:) = 0 + allocate(this%col_dn_forder(nconn)) ; this%col_dn_forder(:) = 0 + allocate(this%face_length(nconn)) ; this%face_length(:) = 0 + allocate(this%uparea(nconn)) ; this%uparea(:) = 0 + allocate(this%downarea(nconn)) ; this%downarea(:) = 0 + allocate(this%dist(nconn)) ; this%dist(:) = 0 + allocate(this%dzg(nconn)) ; this%dzg(:) = 0 + allocate(this%facecos(nconn)) ; this%facecos(:) = 0 + + nconn = 0 + do iconn = 1, moab_edge_internal%num + g_up_moab = moab_edge_internal%cell_ids(iconn, 1) + g_dn_moab = moab_edge_internal%cell_ids(iconn, 2) + + g_up_elm = moab_gcell%moab2elm(g_up_moab) + g_dn_elm = moab_gcell%moab2elm(g_dn_moab) + + if (nat_col_id(g_up_elm) /= -1 .and. nat_col_id(g_dn_elm) /= -1) then + nconn = nconn + 1 + + this%col_id_up(nconn) = nat_col_id(g_up_elm) + this%col_id_dn(nconn) = nat_col_id(g_dn_elm) + + this%grid_id_up(nconn) = g_up_elm + this%grid_id_dn(nconn) = g_dn_elm + + this%grid_id_up_norder(nconn) = moab_gcell%natural_id(g_up_moab) + this%grid_id_dn_norder(nconn) = moab_gcell%natural_id(g_dn_moab) + end if + end do + + ! free up memory + deallocate(nat_col_id) + + end subroutine InitViaMOAB +#endif + +end module ColumnConnectionSetType + diff --git a/components/elm/src/main/elm_initializeMod.F90 b/components/elm/src/main/elm_initializeMod.F90 index 264e02fc267b..59725603b612 100755 --- a/components/elm/src/main/elm_initializeMod.F90 +++ b/components/elm/src/main/elm_initializeMod.F90 @@ -24,6 +24,7 @@ module elm_initializeMod use ELMFatesParamInterfaceMod, only: FatesReadPFTs use BeTRSimulationELM, only : create_betr_simulation_elm use SoilLittVertTranspMod, only : CreateLitterTransportList + use ColumnConnectionSetType, only : c2c_connections use iso_c_binding ! !----------------------------------------- @@ -426,6 +427,9 @@ subroutine initialize1( ) call initGridCells() call initGhostGridCells() +#ifdef HAVE_MOAB + call c2c_connections%Init(bounds_proc) +#endif ! Set global seg maps for gridcells, topounits, landlunits, columns and patches !if(max_topounits > 1) then