Skip to content

[WIP] Sets up column-to-column connections in ELM #7778

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 8 commits into
base: master
Choose a base branch
from
Draft

[WIP] Sets up column-to-column connections in ELM #7778

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
20cee55
Changes the error check because the field may include data for ghost …
bishtgautam Aug 11, 2025
1ff566e
Initialize only the locally owned cells
bishtgautam Aug 11, 2025
b46458c
Adds code to do MOAB-based halo exchange
bishtgautam Aug 13, 2025
12b4e34
Adds code to initialize naturally vegetated ghost columns
bishtgautam Aug 13, 2025
86d586c
Perform halo exhcnage for soil properties
bishtgautam Aug 13, 2025
bc67638
Initialize ghost grid cells
bishtgautam Sep 3, 2025
53c63ee
Modifies computation of beg/end indices for ghost subgrid units
bishtgautam Sep 3, 2025
697bd7c
Adds column-to-column connections when using MOAB
bishtgautam Oct 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 113 additions & 2 deletions components/elm/src/biogeophys/SoilStateType.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

!------------------------------------------------------------------------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

!------------------------------------------------------------------------
Expand All @@ -1096,6 +1206,7 @@ subroutine InitColdGhost(this, bounds_proc)

end subroutine InitColdGhost

#endif
#endif
!------------------------------------------------------------------------

Expand Down
137 changes: 137 additions & 0 deletions components/elm/src/data_types/ColumnConnectionSetType.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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

18 changes: 9 additions & 9 deletions components/elm/src/data_types/ColumnDataType.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
8 changes: 4 additions & 4 deletions components/elm/src/main/accumulMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading