From 9aa8ae5d2cb3b2fc85be5971fd0c5bb3ebde7d75 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Fri, 26 Sep 2025 17:37:32 -0500 Subject: [PATCH 01/29] Remove setting of xao Av Remove setting of xao Av. Loop bound is no longer set by MCT and debugger was complaining. --- driver-moab/main/seq_flux_mct.F90 | 44 ------------------------------- 1 file changed, 44 deletions(-) diff --git a/driver-moab/main/seq_flux_mct.F90 b/driver-moab/main/seq_flux_mct.F90 index b5bcd6642ce9..1965b04db809 100644 --- a/driver-moab/main/seq_flux_mct.F90 +++ b/driver-moab/main/seq_flux_mct.F90 @@ -1937,50 +1937,6 @@ subroutine seq_flux_atmocn_moab(infodata, tod, dt, a2x, o2x, xao, mbid, mbfid) call mbSetCellTagVals(mbfid, 'So_u10withgusts', u10gust, nloc) call mbSetCellTagVals(mbfid, 'So_fswpen', fswpen, nloc) - do n = 1,nloc - if (mask(n) == 0) then - xao%rAttr(index_xao_Faox_sen ,n) = 0.0_r8 - xao%rAttr(index_xao_Faox_lat ,n) = 0.0_r8 - xao%rAttr(index_xao_Faox_taux,n) = 0.0_r8 - xao%rAttr(index_xao_Faox_tauy,n) = 0.0_r8 - xao%rAttr(index_xao_Faox_evap,n) = 0.0_r8 - if ( index_xao_Faox_evap_16O /= 0 ) xao%rAttr(index_xao_Faox_evap_16O,n) = evap_16O(n) - if ( index_xao_Faox_evap_HDO /= 0 ) xao%rAttr(index_xao_Faox_evap_HDO,n) = evap_HDO(n) - if ( index_xao_Faox_evap_18O /= 0 ) xao%rAttr(index_xao_Faox_evap_18O,n) = evap_18O(n) - xao%rAttr(index_xao_So_tref ,n) = 0.0_r8 - xao%rAttr(index_xao_So_qref ,n) = 0.0_r8 - xao%rAttr(index_xao_So_ustar ,n) = 0.0_r8 ! friction velocity - xao%rAttr(index_xao_So_re ,n) = 0.0_r8 ! reynolds number - xao%rAttr(index_xao_So_ssq ,n) = 0.0_r8 ! s.hum. saturation at Ts - xao%rAttr(index_xao_Faox_lwup,n) = 0.0_r8 - xao%rAttr(index_xao_So_duu10n,n) = 0.0_r8 - xao%rAttr(index_xao_So_u10 ,n) = 0.0_r8 - xao%rAttr(index_xao_So_u10withgusts ,n) = 0.0_r8 - if (flux_diurnal) then - xao%rAttr(index_xao_So_warm_diurn ,n) = warm(n) - xao%rAttr(index_xao_So_salt_diurn ,n) = salt(n) - xao%rAttr(index_xao_So_speed_diurn ,n) = speed(n) - xao%rAttr(index_xao_So_regime_diurn ,n) = regime(n) - xao%rAttr(index_xao_So_warmMax_diurn ,n) = warmMax(n) - xao%rAttr(index_xao_So_windMax_diurn ,n) = windMax(n) - xao%rAttr(index_xao_So_qSolAvg_diurn ,n) = qSolAvg(n) - xao%rAttr(index_xao_So_windAvg_diurn ,n) = windAvg(n) - xao%rAttr(index_xao_So_warmMaxInc_diurn ,n) = warmMaxInc(n) - xao%rAttr(index_xao_So_windMaxInc_diurn ,n) = windMaxInc(n) - xao%rAttr(index_xao_So_qSolInc_diurn ,n) = qSolInc(n) - xao%rAttr(index_xao_So_windInc_diurn ,n) = windInc(n) - xao%rAttr(index_xao_So_nInc_diurn ,n) = nInc(n) - xao%rAttr(index_xao_So_tbulk_diurn ,n) = tbulk(n) - xao%rAttr(index_xao_So_tskin_diurn ,n) = tskin(n) - xao%rAttr(index_xao_So_tskin_day_diurn ,n) = tskin_day(n) - xao%rAttr(index_xao_So_tskin_night_diurn,n) = tskin_night(n) - xao%rAttr(index_xao_So_cskin_diurn ,n) = cskin(n) - xao%rAttr(index_xao_So_cskin_night_diurn,n) = cskin_night(n) - xao%rAttr(index_xao_So_fswpen ,n) = fswpen(n) - endif - end if - enddo - #ifdef MOABDEBUG ! debug out file write(lnum,"(I0.2)")num_moab_exports From 3acde49d4861f1dc6d1a8a2bbcb822073718992a Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Fri, 26 Sep 2025 17:38:50 -0500 Subject: [PATCH 02/29] Use correct temporary array sizes Use correct temporary array sizes for ice and ocean instead of the maximum. Requires one more temp array. Also use a logical to allocate them just once. --- driver-moab/main/seq_frac_mct.F90 | 39 +++++++++++++++++-------------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/driver-moab/main/seq_frac_mct.F90 b/driver-moab/main/seq_frac_mct.F90 index 151f09773216..8c8587fc04b8 100644 --- a/driver-moab/main/seq_frac_mct.F90 +++ b/driver-moab/main/seq_frac_mct.F90 @@ -901,9 +901,10 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ integer :: n integer :: ki, kl, ko, kf real(r8),allocatable :: fcorr(:) - real(r8),allocatable :: tagValues(:) ! used for setting some default tags - real(r8),allocatable :: tagValues2(:) ! used for setting some default tags - real(r8),allocatable :: tagValues3(:) ! used for setting some default tags + real(r8),allocatable,save :: tagValues(:) ! used for setting some default tags + real(r8),allocatable,save :: tagValues2(:) ! used for setting some default tags + real(r8),allocatable,save :: tagValues3(:) ! used for setting some default tags + real(r8),allocatable,save :: tagValues4(:) ! used for setting some default tags logical, save :: first_time = .true. ! moab @@ -939,11 +940,9 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ dom_i => component_get_dom_cx(ice) i2x_i => component_get_c2x_cx(ice) - ! make local array big enough to old ocean or ice data. + ! make local array big enough to hold ocean or ice data. if (mbixid .ge. 0) arrSize_i = mbGetnCells(mbixid) if (mboxid .ge. 0) arrSize_o = mbGetnCells(mboxid) - arrSize = arrSize_i - if (arrSize_o .gt. arrSize) arrSize = arrSize_o ! entire update depends on if ice present. if (ice_present) then @@ -960,26 +959,26 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ call seq_frac_check(fractions_i,'ice set') ! MOAB - allocate(tagValues(arrSize) ) - allocate(tagValues2(arrSize) ) - allocate(tagValues3(arrSize) ) + if (first_time) then + allocate(tagValues(arrSize_i) ) + allocate(tagValues2(arrSize_i) ) + allocate(tagValues3(arrSize_i) ) + endif ! copy Si_ifrac to ifrac - call mbGetCellTagVals(mbixid, 'Si_ifrac',tagValues,arrSize) - call mbSetCellTagVals(mbixid, 'ifrac',tagValues,arrSize) + call mbGetCellTagVals(mbixid, 'Si_ifrac',tagValues,arrSize_i) + call mbSetCellTagVals(mbixid, 'ifrac',tagValues,arrSize_i) ! update ifrac and ofrac - call mbGetCellTagVals(mbixid, 'ifrac',tagValues,arrSize) - call mbGetCellTagVals(mbixid, 'frac',tagValues2,arrSize) - call mbGetCellTagVals(mbixid, 'ofrac',tagValues3,arrSize) + call mbGetCellTagVals(mbixid, 'frac',tagValues2,arrSize_i) !fractions_i%rAttr(ki,:) = fractions_i%rAttr(ki,:) * dom_i%data%rAttr(kf,:) tagValues(:) = tagValues(:)*tagValues2(:) - call mbSetCellTagVals(mbixid, 'ifrac',tagValues,arrSize) + call mbSetCellTagVals(mbixid, 'ifrac',tagValues,arrSize_i) !fractions_i%rAttr(ko,:) = dom_i%data%rAttr(kf,:) - fractions_i%rAttr(ki,:) tagValues3(:) = tagValues2(:) - tagValues(:) - call mbSetCellTagVals(mbixid, 'ofrac',tagValues3,arrSize) + call mbSetCellTagVals(mbixid, 'ofrac',tagValues3,arrSize_i) #ifdef MOABDEBUG write(lnum,"(I0.2)")num_moab_exports @@ -989,16 +988,20 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ #endif if (ocn_present) then + if (first_time) then + allocate(tagValues4(arrSize_o) ) + endif mapper_i2o => prep_ocn_get_mapper_SFi2o() call seq_map_map(mapper_i2o, fractions_i, fractions_o, & fldlist='ofrac:ifrac',norm=.false.) call seq_frac_check(fractions_o, 'ocn set') ! set the ofrac on mbofxid instance, because it is needed for prep_aoxflux - call mbGetCellTagVals(mboxid, 'ofrac',tagValues,arrSize) - call mbSetCellTagVals(mbofxid, 'ofrac',tagValues,arrSize) + call mbGetCellTagVals(mboxid, 'ofrac',tagValues4,arrSize_o) + call mbSetCellTagVals(mbofxid, 'ofrac',tagValues4,arrSize_o) endif + first_time=.false. if (atm_present) then mapper_i2a => prep_atm_get_mapper_Fi2a() From 014ede233628f2c45db0f439ab117d639c4a8cb2 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Mon, 1 Sep 2025 19:09:12 -0500 Subject: [PATCH 03/29] first pass at atm spectral separation need to add ext tags fix data atm too --- components/eam/src/cpl/atm_comp_mct.F90 | 9 +- components/eam/src/dynamics/se/dyn_comp.F90 | 4 +- .../homme/interface/phys_grid_mod.F90 | 49 +- components/homme/src/share/semoab_mod.F90 | 1189 +++++++++++++++++ 4 files changed, 1223 insertions(+), 28 deletions(-) diff --git a/components/eam/src/cpl/atm_comp_mct.F90 b/components/eam/src/cpl/atm_comp_mct.F90 index 302d6eedfb8b..1ebce43dfaf8 100644 --- a/components/eam/src/cpl/atm_comp_mct.F90 +++ b/components/eam/src/cpl/atm_comp_mct.F90 @@ -504,12 +504,15 @@ subroutine atm_init_mct( EClock, cdata_a, x2a_a, a2x_a, NLFilename ) call atm_export_moab(Eclock, cam_out) #endif else ! if (StepNo != 0) then + +#ifdef HAVE_MOAB + call t_startf('atm_read_srfrest_moab') + call atm_read_srfrest_moab ( EClock ) + call t_stopf('atm_read_srfrest_moab') +#else call t_startf('atm_read_srfrest_mct') call atm_read_srfrest_mct( EClock, x2a_a, a2x_a ) call t_stopf('atm_read_srfrest_mct') -#ifdef HAVE_MOAB - call atm_read_srfrest_moab ( EClock ) - #endif ! Sent .true. as an optional argument so that restart_init is set to .true. in atm_import diff --git a/components/eam/src/dynamics/se/dyn_comp.F90 b/components/eam/src/dynamics/se/dyn_comp.F90 index faff6ffd314b..b2e04cc21adf 100644 --- a/components/eam/src/dynamics/se/dyn_comp.F90 +++ b/components/eam/src/dynamics/se/dyn_comp.F90 @@ -106,7 +106,7 @@ subroutine dyn_init1(fh, NLFileName, dyn_in, dyn_out) 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 semoab_mod, only: create_moab_pg_mesh use iMOAB, only : iMOAB_RegisterApplication use iso_c_binding #endif @@ -269,7 +269,7 @@ subroutine dyn_init1(fh, NLFileName, dyn_in, dyn_out) end if #ifdef HAVE_MOAB - call create_moab_meshes(par, elem, fv_nphys) + call create_moab_pg_mesh(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/dynamics/homme/interface/phys_grid_mod.F90 b/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90 index 2f13cbc274d0..7c0abdc73f2c 100644 --- a/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90 +++ b/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90 @@ -552,7 +552,7 @@ subroutine phys_grid_init (pgN) 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 semoab_mod, only: create_moab_pg_mesh, create_spectral_moab_meshes use iMOAB, only : iMOAB_RegisterApplication use iso_c_binding #endif @@ -624,27 +624,28 @@ subroutine phys_grid_init (pgN) call compute_global_coords (pg) call compute_global_area (pg) #ifdef HAVE_MOAB + 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 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? ! @@ -662,8 +663,10 @@ subroutine phys_grid_init (pgN) ! 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 + call create_moab_pg_mesh(par, elem, pgN) + else + call create_spectral_moab_meshes(par, elem) + endif #endif end subroutine phys_grid_init diff --git a/components/homme/src/share/semoab_mod.F90 b/components/homme/src/share/semoab_mod.F90 index 4841cb84b543..8d46a441173f 100644 --- a/components/homme/src/share/semoab_mod.F90 +++ b/components/homme/src/share/semoab_mod.F90 @@ -76,6 +76,1195 @@ integer function search_in(intarr, leng, value) end function search_in + subroutine create_spectral_moab_meshes(par, elem) +use ISO_C_BINDING + use iMOAB, only: iMOAB_CreateVertices, iMOAB_WriteMesh, iMOAB_CreateElements, & + iMOAB_ResolveSharedEntities, iMOAB_UpdateMeshInfo, iMOAB_DefineTagStorage, & + iMOAB_SetIntTagStorage, iMOAB_ReduceTagsMax, iMOAB_GetIntTagStorage + use coordinate_systems_mod, only : cartesian3D_t, spherical_to_cart, spherical_polar_t + + type (element_t), intent(inout) :: elem(:) + type (parallel_t) , intent(in) :: par + + integer ierr, i, j, ie, iv, block_ID, k, numvals + integer icol, irow, je, linx ! local indices in fine el connect + + real(kind=real_kind), allocatable :: moab_vert_coords(:) + + integer moab_dim_cquads, ix, idx, nverts, nverts_c ! used for indexing in loops; nverts will have the number of local vertices + + integer nelemd2 ! do not confuse this with dimensions_mod::nelemd + + integer(kind=long_kind), dimension(:), allocatable :: gdofv + ! this will be moab vertex handle locally + integer, dimension(:), allocatable :: moabvh + integer, dimension(:), allocatable :: indx ! this will be ordered + + integer, dimension(:), allocatable :: vdone, elemids, vgids, gdofel + integer, dimension(:), allocatable :: vdone_c, moabconn_c, moabvh_c + integer currentval, dimcoord, dimen, num_el, mbtype, nve + + character*100 outfile, wopts, localmeshfile, lnum, tagname, newtagg + integer tagtype, numco, tag_sto_len, ent_type, tagindex + type (cartesian3D_t) :: cart + integer igcol, ii, neigh + + ! for np=4, + ! 28, 32, 36, 35 + ! 25, 29, 33, 34 + ! j | 13, 17, 21, 22 + ! 1, 5, 9, 10 + !(1,1) i-> + + ! character*100 outfile, wopts, localmeshfile, lnum, tagname + ! integer tagtype, numco, tag_sto_len, ent_type, tagindex + do j=1,np-1 + do i =1, np-1 + ix = (j-1)*(np-1)+i-1 + local_map(i,j) = ix*4 + 1 + enddo + enddo + do j=1, np-1 + i = j + local_map(np, j) = ((np-1)*j-1)*4 + 2 + local_map(i, np) = ( (np-1)*(np-2)+i-1)*4 + 4 + enddo + local_map(np, np) = ((np-1)*(np-1)-1)*4 + 3 + + nelemd2 = (nelemd)*(np-1)*(np-1) + moab_dim_cquads = (nelemd)*4*(np-1)*(np-1) + + if(par%masterproc) then + write (iulog, *) " MOAB: semoab_mod module: create_spectral_moab_meshes on processor " , par%rank ," nelemd: " , nelemd + endif + + if ( nelemd > 0 ) then + allocate(gdofv(moab_dim_cquads)) + allocate(elemids(nelemd2)) + endif + + k=0 ! will be the index for element global dofs + do ie=1,nelemd + do j=1,np-1 + do i=1,np-1 + ix = (ie-1)*(np-1)*(np-1)+(j-1)*(np-1)+i-1 + gdofv(ix*4+1) = elem(ie)%gdofP(i,j) + gdofv(ix*4+2) = elem(ie)%gdofP(i+1,j) + gdofv(ix*4+3) = elem(ie)%gdofP(i+1,j+1) + gdofv(ix*4+4) = elem(ie)%gdofP(i,j+1) + elemids(ix+1) = (elem(ie)%GlobalId-1)*(np-1)*(np-1)+(j-1)*(np-1)+i + enddo + enddo + enddo + +! order according to global dofs + + if ( nelemd > 0 ) then + allocate(moabvh(moab_dim_cquads)) + allocate(indx(moab_dim_cquads)) + + allocate(moabconn(moab_dim_cquads)) + call IndexSet(moab_dim_cquads, indx) + call IndexSort(moab_dim_cquads, indx, gdofv, descend=.false.) +! after sort, gdofv( indx(i)) < gdofv( indx(i+1) ) + endif + idx=0 + currentval = 0 + if ( nelemd > 0 ) then + currentval = gdofv( indx(1)) + idx = 1 + endif + + do ix=1,moab_dim_cquads + if (gdofv(indx(ix)) .ne. currentval ) then + idx=idx+1 + currentval = gdofv(indx(ix)) + endif + moabvh(ix) = idx ! the vertex in connectivity array will be at this local index + ! this will be the moab connectivity + moabconn(indx(ix)) = idx + enddo + + nverts = idx + if(par%masterproc) then + write (iulog, *) " MOAB: there are ", nverts, " local vertices on master task ", currentval, " is the max local gdof" + endif + if ( nelemd > 0 ) then + allocate(moab_vert_coords(3*nverts) ) + allocate(vdone(nverts)) + else + allocate(vdone(1)) ! will be passed to iMOAB_ResolveSharedEnts, and compilers are complaining about non-allocated arrays + endif + vdone = 0; + if ( nelemd > 0 ) currentval = gdofv( indx(1)) ! start over to identify coordinates of the vertices + + do ix=1,moab_dim_cquads + idx = indx(ix) ! index in initial array, vertices in all fine quads + k = (idx-1)/(4*(np-1)*(np-1)) ! index of coarse quad, locally, starts at 0 + ie = 1 + k ! this is the element number; starts at nets=1 + je = ( idx -1 -k*(np-1)*(np-1)*4 ) / 4 + 1 ! local fine quad in coarse, 1 to (np-1) ^ 2 + irow = (je-1)/(np-1)+1 + icol = je - (np-1)*(irow-1) + linx = idx - k*(np-1)*(np-1)*4 -(je-1)*4 ! this should be 1, 2, 3, 4 + if( linx == 1) then + j = irow + i = icol + else if (linx == 2) then + j = irow + i = icol + 1 + else if (linx == 3) then + j = irow + 1 + i = icol + 1 + else ! linx == 4 + j = irow + 1 + i = icol + endif + + iv = moabvh(ix) + if (vdone(iv) .eq. 0) then + cart = spherical_to_cart (elem(ie)%spherep(i,j) ) + moab_vert_coords ( 3*(iv-1)+1 ) = cart%x + moab_vert_coords ( 3*(iv-1)+2 ) = cart%y + moab_vert_coords ( 3*(iv-1)+3 ) = cart%z + vdone(iv) = gdofv(indx(ix)) ! this will be now our tag used for resolving shared entities ! convert to int, from long int + endif + + enddo + + dimcoord = nverts*3 + dimen = 3 + if ( nelemd > 0 ) then + ierr = iMOAB_CreateVertices(MHFID, dimcoord, dimen, moab_vert_coords) + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB vertices ') + endif + !!num_el = nelemd2 + mbtype = 3 ! quadrilateral + nve = 4; + block_ID = 200 ! this will be for coarse mesh + + if ( nelemd > 0 ) then + ierr = iMOAB_CreateElements( MHFID, nelemd2, mbtype, nve, moabconn, block_ID ); + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB elements') + endif + ! nverts: num vertices; vdone will store now the markers used in global resolve + ! for this particular problem, markers are the global dofs at corner nodes +! set the global id for vertices +! first, retrieve the tag + tagname='GDOF'//C_NULL_CHAR + tagtype = 0 ! dense, integer + numco = 1 + ierr = iMOAB_DefineTagStorage(MHFID, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve global id tag') + ! now set the values + ent_type = 0 ! vertex type + if ( nverts > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHFID, tagname, nverts , ent_type, vdone) + if (ierr > 0 ) & + call abortmp('Error: fail to set marker id tag for vertices') + endif + + ! we need to call this even when no mesh locally, it involves a collective + ierr = iMOAB_ResolveSharedEntities( MHFID, nverts, vdone ); + if (ierr > 0 ) & + call abortmp('Error: fail to resolve shared entities') + + if ( nelemd > 0) then + vdone = -1 ! reuse vdone for the new tag, GLOBAL_ID (actual tag that we want to store global dof ) + endif +! use element offset for actual global dofs + ! tagtype = 0 ! dense, integer + ! numco = 1 + newtagg='GLOBAL_ID'//C_NULL_CHAR + ierr = iMOAB_DefineTagStorage(MHFID, newtagg, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to create new GDOF tag') + do ie=1,nelemd + do ii=1,elem(ie)%idxp%NumUniquePts + i=elem(ie)%idxp%ia(ii) + j=elem(ie)%idxp%ja(ii) + igcol = elem(ie)%idxp%UniquePtoffset+ii-1 + ix = local_map(i,j) + idx = moabconn((ie-1)*(np-1)*(np-1)*4 + ix) ! should + vdone ( idx ) = igcol + end do + end do + ! now set the values + ent_type = 0 ! vertex type + if ( nverts > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHFID, newtagg, nverts , ent_type, vdone) + if (ierr > 0 ) & + call abortmp('Error: fail to set global dof tag for vertices') + endif + + ierr = iMOAB_ReduceTagsMax ( MHFID, tagindex, ent_type) + if (ierr > 0 ) & + call abortmp('Error: fail to reduce max tag') + + ! set global id tag for elements + ent_type = 1 ! now set the global id tag on elements + if ( nelemd2 > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHFID, newtagg, nelemd2 , ent_type, elemids) + if (ierr > 0 ) & + call abortmp('Error: fail to set global id tag for elements') + endif + +! now, after reduction, we can get the actual global ids for each vertex in the fine mesh +! before, some vertices that were owned in MOAB but not owned in CAM did not have the right global ID tag +! so vdone will be now correct on every task (no -1 anymore ) + ent_type = 0 ! vertex type + if ( nverts > 0 ) then + allocate(vgids(nverts)) + ierr = iMOAB_GetIntTagStorage ( MHFID, newtagg, nverts , ent_type, vgids) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve GLOBAL ID on each task') + endif + ierr = iMOAB_UpdateMeshInfo(MHFID) + if (ierr > 0 ) & + call abortmp('Error: fail to update mesh info') +#ifdef MOABDEBUG +! write out the mesh file to disk, in parallel + outfile = 'wholeFineATM.h5m'//C_NULL_CHAR + wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(MHFID, outfile, wopts) + if (ierr > 0 ) & + call abortmp('Error: fail to write the mesh file') +#endif + + +! now create the coarse mesh, but the global dofs will come from fine mesh, after solving + ! nelemd2 = nelemd + moab_dim_cquads = (nelemd)*4 + + if ( nelemd > 0 ) then + allocate(gdofel(nelemd*np*np)) + endif + k=0 ! will be the index for element global dofs + do ie=1,nelemd + ix = ie-1 + ! + gdofv(ix*4+1) = elem(ie)%gdofP(1,1) + gdofv(ix*4+2) = elem(ie)%gdofP(np,1) + gdofv(ix*4+3) = elem(ie)%gdofP(np,np) + gdofv(ix*4+4) = elem(ie)%gdofP(1,np) + elemids(ix+1) = elem(ie)%GlobalId + enddo +! now original order + +! order according to global dofs +! allocate(indx(moab_dim_cquads)) + if ( nelemd > 0 ) then + call IndexSet(moab_dim_cquads, indx) + call IndexSort(moab_dim_cquads, indx, gdofv, descend=.false.) +! after sort, gdofv( indx(i)) < gdofv( indx(i+1) ) + + allocate(moabvh_c(moab_dim_cquads)) + + allocate(moabconn_c(moab_dim_cquads)) + endif + idx = 0 + if ( nelemd > 0 ) then + idx=1 + currentval = gdofv( indx(1)) + endif + do ix=1,moab_dim_cquads + if (gdofv(indx(ix)) .ne. currentval ) then + idx=idx+1 + currentval = gdofv(indx(ix)) + endif + moabvh_c(ix) = idx ! the vertex in connectivity array will be at this local index + ! this will be the moab connectivity + moabconn_c(indx(ix)) = idx + enddo + nverts_c = idx + if(par%masterproc) then + write (iulog, *) " MOAB: there are ", nverts_c, " local vertices on master task, coarse mesh" + endif +! allocate(moab_vert_coords(3*idx) ) + if ( nelemd > 0 ) then + allocate(vdone_c(nverts_c)) + vdone_c = 0; + currentval = gdofv( indx(1)) ! start over to identify coordinates of the vertices + endif + + do ix=1,moab_dim_cquads + i = indx(ix) ! index in initial array + ie = 1+ (i-1)/4 ! this is the element number + j = i - ( i-1)/4*4 ! local index of vertex in element i + iv = moabvh_c(ix) + if (vdone_c(iv) .eq. 0) then + moab_vert_coords ( 3*(iv-1)+1 ) = elem(ie)%corners3d(j)%x + moab_vert_coords ( 3*(iv-1)+2 ) = elem(ie)%corners3d(j)%y + moab_vert_coords ( 3*(iv-1)+3 ) = elem(ie)%corners3d(j)%z + vdone_c(iv) = gdofv(indx(ix)) ! this will be now our tag used for resolving shared entities + endif + + enddo + + dimcoord = nverts_c*3 + dimen = 3 + if ( nverts_c > 0 ) then + ierr = iMOAB_CreateVertices(MHID, dimcoord, dimen, moab_vert_coords) + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB vertices ') + endif + ! num_el = nelemd + mbtype = 3 ! quadrilateral + nve = 4; + block_ID = 100 ! this will be for coarse mesh + + if ( nelemd > 0 ) then + ierr = iMOAB_CreateElements( MHID, nelemd, mbtype, nve, moabconn_c, block_ID ); + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB elements') + endif + ! idx: num vertices; vdone will store now the markers used in global resolve + ! for this particular problem, markers are the global dofs at corner nodes +! set the global id for vertices +! first, retrieve the tag + tagname='GDOFV'//C_NULL_CHAR + tagtype = 0 ! dense, integer + numco = 1 + ierr = iMOAB_DefineTagStorage(MHID, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve GDOFV id tag') + ierr = iMOAB_DefineTagStorage(MHID, newtagg, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve GLOBAL_ID tag on coarse mesh') + ! now set the values + ent_type = 0 ! vertex type + if ( nverts_c > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHID, tagname, nverts_c , ent_type, vdone_c) + if (ierr > 0 ) & + call abortmp('Error: fail to set GDOFV tag for vertices') + endif + ! set global id tag for coarse elements, too; they will start at nets=1, end at nete=nelemd + ent_type = 1 ! now set the global id tag on elements + if ( nelemd > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHID, newtagg, nelemd , ent_type, elemids) + if (ierr > 0 ) & + call abortmp('Error: fail to set global id tag for vertices') + endif + + ierr = iMOAB_ResolveSharedEntities( MHID, idx, vdone_c ); + if (ierr > 0 ) & + call abortmp('Error: fail to resolve shared entities') + +! global dofs are the GLL points are set for each element + tagname='GLOBAL_DOFS'//C_NULL_CHAR + tagtype = 0 ! dense, integer + numco = np*np ! usually, it is 16; each element will have the dofs in order + ierr = iMOAB_DefineTagStorage(MHID, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to create global DOFS tag') + ! now set the values + ! set global dofs tag for coarse elements, too; they will start at nets=1, end at nete=nelemd + ent_type = 1 ! now set the global id tag on elements + numvals = nelemd*np*np ! input is the total number of values + ! form gdofel from vgids + do ie=1, nelemd + ix = (ie-1)*np*np ! ie: index in coarse element + je = (ie-1) * 4 * (np-1) * (np -1) ! index in moabconn array + ! vgids are global ids for fine vertices (1,nverts) + iv = 1 + do j=1,np + do i=1,np + k = local_map(i,j) + gdofel(ix+iv) = vgids( moabconn( je + k ) ) + iv = iv + 1 + enddo + enddo + ! extract global ids + vdone_c( moabconn_c( (ie-1)*4+1) ) = vgids ( moabconn(je+1 )) + vdone_c( moabconn_c( (ie-1)*4+2) ) = vgids ( moabconn(je+ 4*(np-2)+2 )) ! valid for np = 4, 10 + vdone_c( moabconn_c( (ie-1)*4+3) ) = vgids ( moabconn(je+ 4*((np-1)*(np-1)-1) + 3 )) ! for np = 4, 35 + vdone_c( moabconn_c( (ie-1)*4+4) ) = vgids ( moabconn(je+ 4*(np-2)*(np-1) + 4 )) ! 28 for np = 4 + enddo + if ( nelemd > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHID, tagname, numvals, ent_type, gdofel) + if (ierr > 0 ) & + call abortmp('Error: fail to set globalDOFs tag for coarse elements') + endif + + ! set the global ids for coarse vertices the same as corresponding fine vertices + ent_type = 0 ! vertex type + if ( nverts_c > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHID, newtagg, nverts_c , ent_type, vdone_c) + if (ierr > 0 ) & + call abortmp('Error: fail to set GLOBAL_DOFS tag values') + endif + + ierr = iMOAB_UpdateMeshInfo(MHID) + if (ierr > 0 ) & + call abortmp('Error: fail to update mesh info') +#ifdef MOABDEBUG +! write out the mesh file to disk, in parallel + outfile = 'wholeATM.h5m'//C_NULL_CHAR + wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(MHID, outfile, wopts) + if (ierr > 0 ) & + call abortmp('Error: fail to write the mesh file') +#endif + ! deallocate + if ( nelemd > 0 ) then + deallocate(moabvh) + deallocate(moabconn) ! do not keep it anymore, we are not setting another tag on fine mesh + deallocate(gdofel) + deallocate(indx) + deallocate(elemids) + deallocate(gdofv) + deallocate(moabvh_c) + deallocate(moabconn_c) + deallocate(vdone_c) + endif + deallocate(vdone) ! we are always allocating this now +! end copy + end subroutine create_spectral_moab_meshes + + subroutine create_moab_pg_mesh(par, elem, fv_nphys) + + use ISO_C_BINDING + use iMOAB, only: iMOAB_CreateVertices, iMOAB_WriteMesh, iMOAB_CreateElements, & + iMOAB_ResolveSharedEntities, iMOAB_UpdateMeshInfo, iMOAB_DefineTagStorage, & + iMOAB_SetIntTagStorage, iMOAB_ReduceTagsMax, iMOAB_GetIntTagStorage + use coordinate_systems_mod, only : cartesian3D_t, spherical_to_cart, spherical_polar_t + + 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 + + real(kind=real_kind), allocatable :: moab_vert_coords(:) + + integer moab_dim_cquads, ix, idx, nverts, nverts_c ! used for indexing in loops; nverts will have the number of local vertices + + integer nelemd2 ! do not confuse this with dimensions_mod::nelemd + + integer(kind=long_kind), dimension(:), allocatable :: gdofv + ! this will be moab vertex handle locally + integer, dimension(:), allocatable :: moabvh + integer, dimension(:), allocatable :: indx ! this will be ordered + + integer, dimension(:), allocatable :: vdone, elemids, vgids, gdofel + integer, dimension(:), allocatable :: vdone_c, moabconn_c, moabvh_c + integer currentval, dimcoord, dimen, num_el, mbtype, nve + + character*100 outfile, wopts, localmeshfile, lnum, tagname, newtagg + integer tagtype, numco, tag_sto_len, ent_type, tagindex + type (cartesian3D_t) :: cart + integer igcol, ii, neigh + + integer nedges_c, nverts_pg, nelem_pg, edge_index, j1 + integer, dimension(:), allocatable :: local_cell_gids, indx_cell + integer, dimension(:,:), allocatable :: elem_edge, edge + integer, dimension(:), allocatable :: vdone_pg, moabconn_pg + + integer nat_edge_order(4) + integer internal_edges, boundary_edges, reverse_edges + integer edge_verts(4) ! local per coarse element ! nverts_c < edge_verts <= nverts_c + edge_index + integer middle_vertex ! nverts_c + edge_index < middle_vertex <= verts_pg + type (spherical_polar_t) :: current_2d_vertex + logical pos_edge ! when looping over edges , use gdof for marking !! + + nat_edge_order = (/south, east, north, west/) + + ! for np=4, + ! 28, 32, 36, 35 + ! 25, 29, 33, 34 + ! j | 13, 17, 21, 22 + ! 1, 5, 9, 10 + !(1,1) i-> + + ! character*100 outfile, wopts, localmeshfile, lnum, tagname + ! integer tagtype, numco, tag_sto_len, ent_type, tagindex + do j=1,np-1 + do i =1, np-1 + ix = (j-1)*(np-1)+i-1 + local_map(i,j) = ix*4 + 1 + enddo + enddo + do j=1, np-1 + i = j + local_map(np, j) = ((np-1)*j-1)*4 + 2 + local_map(i, np) = ( (np-1)*(np-2)+i-1)*4 + 4 + enddo + local_map(np, np) = ((np-1)*(np-1)-1)*4 + 3 + + nelemd2 = (nelemd)*(np-1)*(np-1) + moab_dim_cquads = (nelemd)*4*(np-1)*(np-1) + + if(par%masterproc) then + write (iulog, *) " MOAB: semoab_mod module: create_moab_mesh_fine; on processor " , par%rank ," nelemd: " , nelemd + endif + + if ( nelemd > 0 ) then + allocate(gdofv(moab_dim_cquads)) + allocate(elemids(nelemd2)) + endif + + k=0 ! will be the index for element global dofs + do ie=1,nelemd + do j=1,np-1 + do i=1,np-1 + ix = (ie-1)*(np-1)*(np-1)+(j-1)*(np-1)+i-1 + gdofv(ix*4+1) = elem(ie)%gdofP(i,j) + gdofv(ix*4+2) = elem(ie)%gdofP(i+1,j) + gdofv(ix*4+3) = elem(ie)%gdofP(i+1,j+1) + gdofv(ix*4+4) = elem(ie)%gdofP(i,j+1) + elemids(ix+1) = (elem(ie)%GlobalId-1)*(np-1)*(np-1)+(j-1)*(np-1)+i + enddo + enddo + enddo + +! order according to global dofs + + if ( nelemd > 0 ) then + allocate(moabvh(moab_dim_cquads)) + allocate(indx(moab_dim_cquads)) + + allocate(moabconn(moab_dim_cquads)) + call IndexSet(moab_dim_cquads, indx) + call IndexSort(moab_dim_cquads, indx, gdofv, descend=.false.) +! after sort, gdofv( indx(i)) < gdofv( indx(i+1) ) + endif + idx=0 + currentval = 0 + if ( nelemd > 0 ) then + currentval = gdofv( indx(1)) + idx = 1 + endif + + do ix=1,moab_dim_cquads + if (gdofv(indx(ix)) .ne. currentval ) then + idx=idx+1 + currentval = gdofv(indx(ix)) + endif + moabvh(ix) = idx ! the vertex in connectivity array will be at this local index + ! this will be the moab connectivity + moabconn(indx(ix)) = idx + enddo + + nverts = idx + if(par%masterproc) then + write (iulog, *) " MOAB: there are ", nverts, " local vertices on master task ", currentval, " is the max local gdof" + endif + if ( nelemd > 0 ) then + allocate(moab_vert_coords(3*nverts) ) + allocate(vdone(nverts)) + else + allocate(vdone(1)) ! will be passed to iMOAB_ResolveSharedEnts, and compilers are complaining about non-allocated arrays + endif + vdone = 0; + if ( nelemd > 0 ) currentval = gdofv( indx(1)) ! start over to identify coordinates of the vertices + + do ix=1,moab_dim_cquads + idx = indx(ix) ! index in initial array, vertices in all fine quads + k = (idx-1)/(4*(np-1)*(np-1)) ! index of coarse quad, locally, starts at 0 + ie = 1 + k ! this is the element number; starts at nets=1 + je = ( idx -1 -k*(np-1)*(np-1)*4 ) / 4 + 1 ! local fine quad in coarse, 1 to (np-1) ^ 2 + irow = (je-1)/(np-1)+1 + icol = je - (np-1)*(irow-1) + linx = idx - k*(np-1)*(np-1)*4 -(je-1)*4 ! this should be 1, 2, 3, 4 + if( linx == 1) then + j = irow + i = icol + else if (linx == 2) then + j = irow + i = icol + 1 + else if (linx == 3) then + j = irow + 1 + i = icol + 1 + else ! linx == 4 + j = irow + 1 + i = icol + endif + + iv = moabvh(ix) + if (vdone(iv) .eq. 0) then + cart = spherical_to_cart (elem(ie)%spherep(i,j) ) + moab_vert_coords ( 3*(iv-1)+1 ) = cart%x + moab_vert_coords ( 3*(iv-1)+2 ) = cart%y + moab_vert_coords ( 3*(iv-1)+3 ) = cart%z + vdone(iv) = gdofv(indx(ix)) ! this will be now our tag used for resolving shared entities ! convert to int, from long int + endif + + enddo + + dimcoord = nverts*3 + dimen = 3 + if ( nelemd > 0 ) then + ierr = iMOAB_CreateVertices(MHFID, dimcoord, dimen, moab_vert_coords) + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB vertices ') + endif + !!num_el = nelemd2 + mbtype = 3 ! quadrilateral + nve = 4; + block_ID = 200 ! this will be for coarse mesh + + if ( nelemd > 0 ) then + ierr = iMOAB_CreateElements( MHFID, nelemd2, mbtype, nve, moabconn, block_ID ); + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB elements') + endif + ! nverts: num vertices; vdone will store now the markers used in global resolve + ! for this particular problem, markers are the global dofs at corner nodes +! set the global id for vertices +! first, retrieve the tag + tagname='GDOF'//C_NULL_CHAR + tagtype = 0 ! dense, integer + numco = 1 + ierr = iMOAB_DefineTagStorage(MHFID, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve global id tag') + ! now set the values + ent_type = 0 ! vertex type + if ( nverts > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHFID, tagname, nverts , ent_type, vdone) + if (ierr > 0 ) & + call abortmp('Error: fail to set marker id tag for vertices') + endif + + ! we need to call this even when no mesh locally, it involves a collective + ierr = iMOAB_ResolveSharedEntities( MHFID, nverts, vdone ); + if (ierr > 0 ) & + call abortmp('Error: fail to resolve shared entities') + + if ( nelemd > 0) then + vdone = -1 ! reuse vdone for the new tag, GLOBAL_ID (actual tag that we want to store global dof ) + endif +! use element offset for actual global dofs + ! tagtype = 0 ! dense, integer + ! numco = 1 + newtagg='GLOBAL_ID'//C_NULL_CHAR + ierr = iMOAB_DefineTagStorage(MHFID, newtagg, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to create new GDOF tag') + do ie=1,nelemd + do ii=1,elem(ie)%idxp%NumUniquePts + i=elem(ie)%idxp%ia(ii) + j=elem(ie)%idxp%ja(ii) + igcol = elem(ie)%idxp%UniquePtoffset+ii-1 + ix = local_map(i,j) + idx = moabconn((ie-1)*(np-1)*(np-1)*4 + ix) ! should + vdone ( idx ) = igcol + end do + end do + ! now set the values + ent_type = 0 ! vertex type + if ( nverts > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHFID, newtagg, nverts , ent_type, vdone) + if (ierr > 0 ) & + call abortmp('Error: fail to set global dof tag for vertices') + endif + + ierr = iMOAB_ReduceTagsMax ( MHFID, tagindex, ent_type) + if (ierr > 0 ) & + call abortmp('Error: fail to reduce max tag') + + ! set global id tag for elements + ent_type = 1 ! now set the global id tag on elements + if ( nelemd2 > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHFID, newtagg, nelemd2 , ent_type, elemids) + if (ierr > 0 ) & + call abortmp('Error: fail to set global id tag for elements') + endif + +! now, after reduction, we can get the actual global ids for each vertex in the fine mesh +! before, some vertices that were owned in MOAB but not owned in CAM did not have the right global ID tag +! so vdone will be now correct on every task (no -1 anymore ) + ent_type = 0 ! vertex type + if ( nverts > 0 ) then + allocate(vgids(nverts)) + ierr = iMOAB_GetIntTagStorage ( MHFID, newtagg, nverts , ent_type, vgids) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve GLOBAL ID on each task') + endif + ierr = iMOAB_UpdateMeshInfo(MHFID) + if (ierr > 0 ) & + call abortmp('Error: fail to update mesh info') +#ifdef MOABDEBUG +! write out the mesh file to disk, in parallel + outfile = 'wholeFineATM.h5m'//C_NULL_CHAR + wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(MHFID, outfile, wopts) + if (ierr > 0 ) & + call abortmp('Error: fail to write the mesh file') +#endif + + +! now create the coarse mesh, but the global dofs will come from fine mesh, after solving + ! nelemd2 = nelemd + moab_dim_cquads = (nelemd)*4 + + if ( nelemd > 0 ) then + allocate(gdofel(nelemd*np*np)) + endif + k=0 ! will be the index for element global dofs + do ie=1,nelemd + ix = ie-1 + ! + gdofv(ix*4+1) = elem(ie)%gdofP(1,1) + gdofv(ix*4+2) = elem(ie)%gdofP(np,1) + gdofv(ix*4+3) = elem(ie)%gdofP(np,np) + gdofv(ix*4+4) = elem(ie)%gdofP(1,np) + elemids(ix+1) = elem(ie)%GlobalId + enddo +! now original order + +! order according to global dofs +! allocate(indx(moab_dim_cquads)) + if ( nelemd > 0 ) then + call IndexSet(moab_dim_cquads, indx) + call IndexSort(moab_dim_cquads, indx, gdofv, descend=.false.) +! after sort, gdofv( indx(i)) < gdofv( indx(i+1) ) + + allocate(moabvh_c(moab_dim_cquads)) + + allocate(moabconn_c(moab_dim_cquads)) + endif + idx = 0 + if ( nelemd > 0 ) then + idx=1 + currentval = gdofv( indx(1)) + endif + do ix=1,moab_dim_cquads + if (gdofv(indx(ix)) .ne. currentval ) then + idx=idx+1 + currentval = gdofv(indx(ix)) + endif + moabvh_c(ix) = idx ! the vertex in connectivity array will be at this local index + ! this will be the moab connectivity + moabconn_c(indx(ix)) = idx + enddo + nverts_c = idx + if(par%masterproc) then + write (iulog, *) " MOAB: there are ", nverts_c, " local vertices on master task, coarse mesh" + endif +! allocate(moab_vert_coords(3*idx) ) + if ( nelemd > 0 ) then + allocate(vdone_c(nverts_c)) + vdone_c = 0; + currentval = gdofv( indx(1)) ! start over to identify coordinates of the vertices + endif + + do ix=1,moab_dim_cquads + i = indx(ix) ! index in initial array + ie = 1+ (i-1)/4 ! this is the element number + j = i - ( i-1)/4*4 ! local index of vertex in element i + iv = moabvh_c(ix) + if (vdone_c(iv) .eq. 0) then + moab_vert_coords ( 3*(iv-1)+1 ) = elem(ie)%corners3d(j)%x + moab_vert_coords ( 3*(iv-1)+2 ) = elem(ie)%corners3d(j)%y + moab_vert_coords ( 3*(iv-1)+3 ) = elem(ie)%corners3d(j)%z + vdone_c(iv) = gdofv(indx(ix)) ! this will be now our tag used for resolving shared entities + endif + + enddo + + dimcoord = nverts_c*3 + dimen = 3 + if ( nverts_c > 0 ) then + ierr = iMOAB_CreateVertices(MHID, dimcoord, dimen, moab_vert_coords) + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB vertices ') + endif + ! num_el = nelemd + mbtype = 3 ! quadrilateral + nve = 4; + block_ID = 100 ! this will be for coarse mesh + + if ( nelemd > 0 ) then + ierr = iMOAB_CreateElements( MHID, nelemd, mbtype, nve, moabconn_c, block_ID ); + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB elements') + endif + ! idx: num vertices; vdone will store now the markers used in global resolve + ! for this particular problem, markers are the global dofs at corner nodes +! set the global id for vertices +! first, retrieve the tag + tagname='GDOFV'//C_NULL_CHAR + tagtype = 0 ! dense, integer + numco = 1 + ierr = iMOAB_DefineTagStorage(MHID, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve GDOFV id tag') + ierr = iMOAB_DefineTagStorage(MHID, newtagg, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve GLOBAL_ID tag on coarse mesh') + ! now set the values + ent_type = 0 ! vertex type + if ( nverts_c > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHID, tagname, nverts_c , ent_type, vdone_c) + if (ierr > 0 ) & + call abortmp('Error: fail to set GDOFV tag for vertices') + endif + ! set global id tag for coarse elements, too; they will start at nets=1, end at nete=nelemd + ent_type = 1 ! now set the global id tag on elements + if ( nelemd > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHID, newtagg, nelemd , ent_type, elemids) + if (ierr > 0 ) & + call abortmp('Error: fail to set global id tag for vertices') + endif + + ierr = iMOAB_ResolveSharedEntities( MHID, idx, vdone_c ); + if (ierr > 0 ) & + call abortmp('Error: fail to resolve shared entities') + +! global dofs are the GLL points are set for each element + tagname='GLOBAL_DOFS'//C_NULL_CHAR + tagtype = 0 ! dense, integer + numco = np*np ! usually, it is 16; each element will have the dofs in order + ierr = iMOAB_DefineTagStorage(MHID, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to create global DOFS tag') + ! now set the values + ! set global dofs tag for coarse elements, too; they will start at nets=1, end at nete=nelemd + ent_type = 1 ! now set the global id tag on elements + numvals = nelemd*np*np ! input is the total number of values + ! form gdofel from vgids + do ie=1, nelemd + ix = (ie-1)*np*np ! ie: index in coarse element + je = (ie-1) * 4 * (np-1) * (np -1) ! index in moabconn array + ! vgids are global ids for fine vertices (1,nverts) + iv = 1 + do j=1,np + do i=1,np + k = local_map(i,j) + gdofel(ix+iv) = vgids( moabconn( je + k ) ) + iv = iv + 1 + enddo + enddo + ! extract global ids + vdone_c( moabconn_c( (ie-1)*4+1) ) = vgids ( moabconn(je+1 )) + vdone_c( moabconn_c( (ie-1)*4+2) ) = vgids ( moabconn(je+ 4*(np-2)+2 )) ! valid for np = 4, 10 + vdone_c( moabconn_c( (ie-1)*4+3) ) = vgids ( moabconn(je+ 4*((np-1)*(np-1)-1) + 3 )) ! for np = 4, 35 + vdone_c( moabconn_c( (ie-1)*4+4) ) = vgids ( moabconn(je+ 4*(np-2)*(np-1) + 4 )) ! 28 for np = 4 + enddo + if ( nelemd > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHID, tagname, numvals, ent_type, gdofel) + if (ierr > 0 ) & + call abortmp('Error: fail to set globalDOFs tag for coarse elements') + endif + + ! set the global ids for coarse vertices the same as corresponding fine vertices + ent_type = 0 ! vertex type + if ( nverts_c > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHID, newtagg, nverts_c , ent_type, vdone_c) + if (ierr > 0 ) & + call abortmp('Error: fail to set GLOBAL_DOFS tag values') + endif + + ierr = iMOAB_UpdateMeshInfo(MHID) + if (ierr > 0 ) & + call abortmp('Error: fail to update mesh info') +#ifdef MOABDEBUG +! write out the mesh file to disk, in parallel + outfile = 'wholeATM.h5m'//C_NULL_CHAR + wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(MHID, outfile, wopts) + if (ierr > 0 ) & + call abortmp('Error: fail to write the mesh file') +#endif + + if (fv_nphys > 0 ) then + ! create FV mesh, base on PGx + atm_pg_active = .true. ! from now on, we will migrate / send tags for FV / ATM_PG2 mesh + ! first count the number of edges in the coarse mesh; + ! use euler: v-m+f = 2 => m = v + f - 2 + nedges_c = nverts_c + nelemd - 1 ! could be more, if unconnected regions ? + if ( nedges_c < 0 ) nedges_c = 0 ! it cannot be negative + internal_edges = 0 + boundary_edges = 0 + reverse_edges = 0 + nelem_pg = fv_nphys * fv_nphys * nelemd ! each coarse cell is divided in fv_nphys x fv_nphys subcells + ! + ! there are new vertices on each coarse edge (fv_phys - 1) , and (fv_nphys - 1) * (fv_nphys - 1) + ! new vertices on each coarse cell + if ( nelemd > 0 ) then + allocate (local_cell_gids(nelemd)) + allocate (indx_cell(nelemd)) + allocate (edge(2,nedges_c)) ! + endif + do ie=1, nelemd ! + local_cell_gids(ie) = elem(ie)%GlobalID + enddo + if ( nelemd > 0 ) then + call IndexSet(nelemd, indx_cell) + call IndexSort(nelemd, indx_cell, local_cell_gids, descend=.false.) + ! print *, ' local_cell_gids ', local_cell_gids + ! print *, ' indx_cell ', indx_cell + allocate( elem_edge (4, nelemd) ) + ! print *, '------------------------------- ' + ! print *, "RANK:", par%rank + endif + edge_index = 0 + do ie=1, nelemd ! + ! we need to check if neighbor is with id smaller; that means it was already created ? + ! print *, '-------------- ' + ! print *, ' elem ', ie, elem(ie)%desc%actual_neigh_edges, elem(ie)%vertex%number, elem(ie)%GlobalID + ! print *, ' nodes ', ( moabconn_c( (ie-1)*4+j), j=1,4 ) + ! print *, ' ids ', (vdone_c( moabconn_c( (ie-1)*4+j) ), j=1,4) + ! print *, ' neigh: ', (elem(ie)%desc%globalID(j), j=1,4) + ! print *, ' neigh order ', ( elem(ie)%desc%globalID(nat_edge_order(j)),j = 1,4 ) + k = elem(ie)%GlobalID ! current id + do j = 1,4 + ix = j+1 + if (ix .eq. 5) ix = 1 ! next vertex in connectivity array + neigh = elem(ie)%desc%globalID(nat_edge_order(j)) ! id neighbor + idx = search_in(local_cell_gids, nelemd, neigh) ! index in local cells + ! print *, ' ', j, 'neigh:', neigh, ' index' , idx + + if ( idx .gt. 0 ) then + ! a local edge is interior + + if (k < neigh) then ! form the edge, increment edge index + edge_index = edge_index + 1 + edge(1, edge_index) = moabconn_c(4*(ie-1) + j) ! first vertex + edge(2, edge_index) = moabconn_c(4*(ie-1) + ix) ! second vertex index + elem_edge(j, ie) = edge_index + internal_edges = internal_edges + 1 + ! print *, ' edge:', edge_index, edge(1, edge_index), edge(2, edge_index), 'verts:' , & + ! vdone_c(edge(1, edge_index)), vdone_c(edge(2, edge_index)), ' element ', ie, ' intedge:', internal_edges + + else + ! find the edge in the other list elem(idx)%globalID( nat_edge_order(j) ) + do j1 = 1,4 + if ( elem(idx)%desc%globalID( nat_edge_order(j1) ) .eq. k ) then + elem_edge(j, ie) = - elem_edge(j1, idx) ! inverse oriented + reverse_edges = reverse_edges + 1 + ! print *, ' negative edge: ', elem_edge(j, ie), edge(1, -elem_edge(j, ie)), edge(2, -elem_edge(j, ie)), & + ! 'verts:', vdone_c(edge(1, -elem_edge(j, ie))), vdone_c(edge(2, -elem_edge(j, ie))), 'indx neg', reverse_edges + endif + enddo + + endif + else ! idx is 0, so it means the edge is on the boundary, form it + edge_index = edge_index + 1 + edge(1, edge_index) = moabconn_c(4*(ie-1) + j) ! first vertex + edge(2, edge_index) = moabconn_c(4*(ie-1) + ix) ! second vertex index + elem_edge(j, ie) = edge_index + boundary_edges = boundary_edges + 1 + ! print *, ' edge:', edge_index, edge(1, edge_index), edge(2, edge_index), 'verts:' , & + ! vdone_c(edge(1, edge_index)), vdone_c(edge(2, edge_index)), ' element ', ie, & + ! ' bedge:', boundary_edges + endif + enddo + enddo + ! show off + nverts_pg = nverts_c + (fv_nphys - 1) * edge_index + (fv_nphys - 1) * (fv_nphys - 1) * nelemd + ! print *, " MOAB: there are ", nverts_pg, " local vertices on master task on pg mesh ", edge_index , " local coarse edges ", & + ! boundary_edges , ' boundary edges ' + if(par%masterproc) then + write (iulog, *) " MOAB: there are ", nverts_pg, " local vertices on master task on pg mesh ", edge_index , " local coarse edges " + endif + !print *, '\n ELEMENTS: ' + !do ie=1,nelemd + ! print *, ie, elem(ie)%GlobalID, ' local nodes:', ( moabconn_c( (ie-1)*4+j), j=1,4 ), ' edges:', (elem_edge(j, ie), j=1,4) + !enddo + !print *, '\n EDGES:' + !do ie=1,edge_index + ! print *, ie, (edge(j, ie), j=1,2) + !enddo + ! now generate phys grid, uniform FV type mesh; + ! 2 cases: fv_nphys is 1 or 2; when 2, we need new nodes; will use the same id as + ! the gdof on edge is used, with the smaller id chosen, among + if ( nelemd > 0 ) then + allocate(moabconn_pg(4*nelem_pg)) ! connectivity + ! reuse moab_vert_coords for coordinates of pg mesh + ! the first nverts_c coords are the same as coarse mesh; but we do create new + allocate(vdone_pg(nverts_pg)) + else + allocate(vdone_pg(1)) + endif + do iv = 1, nverts_c + vdone_pg(iv) = vdone_c(iv) ! also the coordinates will be the same !! + enddo + + ! copy the coordinates from the middle + j1 = 0 ! index in edge vertices; increase only for positive edges + ! still need some + if (fv_nphys .eq. 2) then + current_2d_vertex%r = 1. + do ie = 1,nelemd + ix = (ie-1)*np*np ! ie: index in coarse element + do j=1,4 + idx = elem_edge(j, ie) ! + if (idx .gt. 0) then ! increment edges, add vertex ! + j1 = j1 + 1 ! index in moab_vert_coords for edges ! nverts_c + j1 for vertex edges ! + ! current_2d_vertex%lat, current_2d_vertex%lon + pos_edge = .true. + iv = nverts_c + j1 + edge_verts(j) = iv ! to form the local connectivity array + if ( vdone_c(edge(1, idx)) .gt. vdone_c(edge(2, idx)) ) pos_edge = .false. + if (j .eq. 1) then + call gfr_f_get_corner_latlon(ie, 1, 1, 2, current_2d_vertex%lat, current_2d_vertex%lon) + if (pos_edge) then + vdone_pg (iv) = gdofel(ix + 2) ! elem(ie)%gdofP(2,1) ! gdofel(ix+ (j-1)*np + i) + else + vdone_pg (iv) = gdofel(ix + np - 1) !elem(ie)%gdofP(np-1,1) ! + endif + else if (j .eq. 2) then + call gfr_f_get_corner_latlon(ie, 2, 1, 3, current_2d_vertex%lat, current_2d_vertex%lon) + if (pos_edge) then + vdone_pg (iv) = gdofel(ix + (2 - 1) * np + np)!elem(ie)%gdofP(np,2) ! ! gdofel(ix+ (j-1)*np + i) + else + vdone_pg (iv) = gdofel(ix + (np - 2) * np + np)!elem(ie)%gdofP(np,np - 1) ! + endif + else if (j .eq. 3) then + call gfr_f_get_corner_latlon(ie, 2, 2, 4, current_2d_vertex%lat, current_2d_vertex%lon) + if (pos_edge) then + vdone_pg (iv) = gdofel(ix+ (np - 1) * np + np - 1)!elem(ie)%gdofP(np-1,np) ! + else + vdone_pg (iv) = gdofel(ix+ (np-1)*np + 2) !elem(ie)%gdofP(2,np) ! + endif + else ! if (j .eq. 4) + call gfr_f_get_corner_latlon(ie, 1, 2, 1, current_2d_vertex%lat, current_2d_vertex%lon) + if (pos_edge) then + vdone_pg (iv) = gdofel(ix+ (np - 2)*np + 1) !elem(ie)%gdofP(1,np-1) ! + else + vdone_pg (iv) = gdofel(ix+ ( 2 - 1 )*np + 1) ! elem(ie)%gdofP(1,2) ! + endif + endif + ! create the 3d vertex ! + cart = spherical_to_cart (current_2d_vertex ) + moab_vert_coords ( 3*(iv-1)+1 ) = cart%x + moab_vert_coords ( 3*(iv-1)+2 ) = cart%y + moab_vert_coords ( 3*(iv-1)+3 ) = cart%z + ! print *, 'ie, j, iv, vdone_pg(iv): ', ie, j, iv, vdone_pg(iv) + else ! the vertex was already created, but we need the index for connectivity of local fv cells + edge_verts(j) = nverts_c + ( -idx ) ! idx is index of edge (negative for already created) + endif + + enddo ! do j=1,4 + ! create the middle vertex too, in the center + call gfr_f_get_corner_latlon(ie, 1, 1, 3, current_2d_vertex%lat, current_2d_vertex%lon) + iv = nverts_c + edge_index + ie ! middle vertices are after corners, and edge vertices + middle_vertex = iv + vdone_pg (middle_vertex) = gdofel(ix+ np + 2)!elem(ie)%gdofP(2,2) ! first in the interior, not on edges! + ! print *, 'ie, middle = iv, vdone_pg(iv): ', ie, iv, vdone_pg(iv) + cart = spherical_to_cart (current_2d_vertex ) + moab_vert_coords ( 3*(iv-1)+1 ) = cart%x + moab_vert_coords ( 3*(iv-1)+2 ) = cart%y + moab_vert_coords ( 3*(iv-1)+3 ) = cart%z + + ! now form the local 2x2 cells, one by one; set the global id tag too! + idx = (ie-1)*4 + ix = idx * 4 ! + ! first + moabconn_pg(ix + 1) = moabconn_c(4*(ie-1)+1) + moabconn_pg(ix + 2) = edge_verts(1) + moabconn_pg(ix + 3) = middle_vertex + moabconn_pg(ix + 4) = edge_verts(4) + elemids(idx+1) = (elem(ie)%GlobalId-1)*4+1 + ! second + moabconn_pg(ix + 4 + 1) = edge_verts(1) + moabconn_pg(ix + 4 + 2) = moabconn_c(4*(ie-1)+2) + moabconn_pg(ix + 4 + 3) = edge_verts(2) + moabconn_pg(ix + 4 + 4) = middle_vertex + elemids(idx+2) = (elem(ie)%GlobalId-1)*4+2 + ! third + moabconn_pg(ix + 8 + 1) = edge_verts(4) + moabconn_pg(ix + 8 + 2) = middle_vertex + moabconn_pg(ix + 8 + 3) = edge_verts(3) + moabconn_pg(ix + 8 + 4) = moabconn_c(4*(ie-1)+4) + elemids(idx+3) = (elem(ie)%GlobalId-1)*4+3 + ! fourth + moabconn_pg(ix + 12 + 1) = middle_vertex + moabconn_pg(ix + 12 + 2) = edge_verts(2) + moabconn_pg(ix + 12 + 3) = moabconn_c(4*(ie-1)+3) + moabconn_pg(ix + 12 + 4) = edge_verts(3) + elemids(idx+4) = (elem(ie)%GlobalId-1)*4+4 + + enddo + ! now copy from coarse for pg mesh + + dimcoord = nverts_pg*3 + dimen = 3 + if ( nverts_pg > 0 ) then + ierr = iMOAB_CreateVertices(MHPGID, dimcoord, dimen, moab_vert_coords) + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB vertices ') + endif + ! num_el = nelem_pg * + mbtype = 3 ! quadrilateral + nve = 4; + block_ID = 300 ! this will be for pg mesh + + if ( nelem_pg > 0 ) then + ierr = iMOAB_CreateElements( MHPGID, nelem_pg, mbtype, nve, moabconn_pg, block_ID ); + if (ierr > 0 ) & + call abortmp('Error: fail to create MOAB elements') + endif + tagname='GLOBAL_ID'//C_NULL_CHAR + tagtype = 0 ! dense, integer + numco = 1 + ierr = iMOAB_DefineTagStorage(MHPGID, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call abortmp('Error: fail to retrieve GLOBAL id tag') + + ! now set the values + ent_type = 0 ! vertex type + if ( nverts_pg > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHPGID, tagname, nverts_pg , ent_type, vdone_pg) + if (ierr > 0 ) & + call abortmp('Error: fail to set global id tag for vertices') + endif + ! set global id tag for pg2 elements, too; they will start at nets=1, end at nete=nelemd*4 + ent_type = 1 ! now set the global id tag on elements + if ( nelem_pg > 0 ) then + ierr = iMOAB_SetIntTagStorage ( MHPGID, tagname, nelem_pg , ent_type, elemids) + if (ierr > 0 ) & + call abortmp('Error: fail to set global id tag for edges') + endif + + ! this involves a collective, vdone_pg can be empty + ierr = iMOAB_ResolveSharedEntities( MHPGID, nverts_pg, vdone_pg ); + if (ierr > 0 ) & + call abortmp('Error: fail to resolve shared ents for pg2 mesh') + + ierr = iMOAB_UpdateMeshInfo(MHPGID) + if (ierr > 0 ) & + call abortmp('Error: fail to update mesh info for pg2 mesh') +#ifdef MOABDEBUG + ! write out the mesh file to disk, in parallel + outfile = 'wholeATM_PG2.h5m'//C_NULL_CHAR + wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(MHPGID, outfile, wopts) + if (ierr > 0 ) & + call abortmp('Error: fail to write the mesh file') +#endif + endif ! only valid for pg == 2 + if ( nelemd > 0 ) then + deallocate (local_cell_gids) + deallocate (indx_cell) + deallocate (edge) ! + deallocate(moabconn_pg) ! connectivity + endif + deallocate(vdone_pg) ! this is now always allocated/deallocated, even for no mesh here + endif + + ! deallocate + if ( nelemd > 0 ) then + deallocate(moabvh) + deallocate(moabconn) ! do not keep it anymore, we are not setting another tag on fine mesh + deallocate(gdofel) + deallocate(indx) + deallocate(elemids) + deallocate(gdofv) + deallocate(moabvh_c) + deallocate(moabconn_c) + deallocate(vdone_c) + endif + deallocate(vdone) ! we are always allocating this now +! end copy + + end subroutine create_moab_pg_mesh + subroutine create_moab_meshes(par, elem, fv_nphys) use ISO_C_BINDING From 2f2054bdb4c50f8fe20d9455f02bbd8d9fdeab6d Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Thu, 18 Sep 2025 12:40:54 -0500 Subject: [PATCH 04/29] add more extension fields for _ext tags, we need more methods to retrieve data Or maybe we will go back to dual mesh, and then dual mesh will be like any other FV mesh on coupler side the _ext tags use about twice as much memory while the dual mesh would be more optimal --- driver-moab/main/cplcomp_exchange_mod.F90 | 37 +++++++---- driver-moab/shr/seq_flds_mod.F90 | 76 +++++++++++------------ 2 files changed, 61 insertions(+), 52 deletions(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 6d5c4bdc638c..1fa5273a3ce3 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -9,8 +9,8 @@ module cplcomp_exchange_mod use seq_map_type_mod use component_type_mod use seq_flds_mod, only: seq_flds_dom_coord, seq_flds_dom_other - use seq_flds_mod, only: seq_flds_dom_fields - use seq_flds_mod, only: seq_flds_a2x_ext_fields, seq_flds_a2x_fields, seq_flds_x2a_fields ! + use seq_flds_mod, only: seq_flds_dom_fields, seq_flds_dom_ext_fields + use seq_flds_mod, only: seq_flds_a2x_ext_fields, seq_flds_a2x_fields, seq_flds_x2a_ext_fields, seq_flds_x2a_fields ! use seq_flds_mod, only: seq_flds_o2x_fields ! needed for MOAB init of ocean fields o2x to be able to transfer to coupler use seq_flds_mod, only: seq_flds_x2o_fields ! needed for MOAB init of ocean fields x2o to be able to transfer from coupler use seq_flds_mod, only: seq_flds_i2x_fields, seq_flds_x2i_fields ! needed for MOAB init of ice fields x2o on coupler side, to save them @@ -1231,8 +1231,14 @@ subroutine cplcomp_moab_Init(infodata,comp) write(logunit,*) subname,' error in defining tags on atm on coupler ' call shr_sys_abort(subname//' ERROR in defining tags ') endif + if (atm_pg_active) then + tagname = trim(seq_flds_x2a_fields)//C_NULL_CHAR + numco = 1 ! usually 1 value per cell + else ! this is not supported now, but leave it here + tagname = trim(seq_flds_x2a_ext_fields)//C_NULL_CHAR ! MOAB versions of a2x for spectral + numco = 16 ! np*np ! usually 16 values per cell, GLL points; should be 4 x 4 = 16 + endif - tagname = trim(seq_flds_x2a_fields)//C_NULL_CHAR ! TODO should be also x2a_ext for spectral case ierr = iMOAB_DefineTagStorage(mbaxid, tagname, tagtype, numco, tagindex ) if (ierr .ne. 0) then write(logunit,*) subname,' error in defining tags seq_flds_x2a_fields on atm on coupler ' @@ -1240,19 +1246,28 @@ subroutine cplcomp_moab_Init(infodata,comp) endif !add the normalization tag - tagname = trim(seq_flds_dom_fields)//":norm8wt"//C_NULL_CHAR + if (atm_pg_active) then + tagname = trim(seq_flds_dom_fields)//C_NULL_CHAR + numco = 1 ! usually 1 value per cell + else ! this is not supported now, but leave it here + tagname = trim(seq_flds_dom_ext_fields)//C_NULL_CHAR ! MOAB versions of a2x for spectral + numco = 16 ! np*np ! usually 16 values per cell, GLL points; should be 4 x 4 = 16 + endif + ierr = iMOAB_DefineTagStorage(mbaxid, tagname, tagtype, numco, tagindex ) if (ierr .ne. 0) then write(logunit,*) subname,' error in defining tags seq_flds_dom_fields on atm on coupler ' call shr_sys_abort(subname//' ERROR in defining tags ') endif - ! also, frac, area, masks has to come from atm mphaid, not from domain file reader - ! this is hard to digest :( - tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR - ! TODO: this should be called on the joint procs, not coupler only. - call component_exch_moab(comp, mphaid, mbaxid, 0, tagname, context_exch='doma') - ! copy aream from area in case atm_mesh - call copy_aream_from_area(mbaxid) + if (atm_pg_active) then + ! also, frac, area, masks has to come from atm mphaid, not from domain file reader + ! this is hard to digest :( + tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR + ! TODO: this should be called on the joint procs, not coupler only. + call component_exch_moab(comp, mphaid, mbaxid, 0, tagname, context_exch='doma') + ! copy aream from area in case atm_mesh + call copy_aream_from_area(mbaxid) + endif endif ! coupler pes #ifdef MOABDEBUG diff --git a/driver-moab/shr/seq_flds_mod.F90 b/driver-moab/shr/seq_flds_mod.F90 index f8abf9e4c396..efb75f41d3d7 100644 --- a/driver-moab/shr/seq_flds_mod.F90 +++ b/driver-moab/shr/seq_flds_mod.F90 @@ -190,8 +190,6 @@ module seq_flds_mod character(CXX) :: seq_flds_a2x_states character(CXX) :: seq_flds_a2x_fluxes - character(CXX) :: seq_flds_a2x_ext_states - character(CXX) :: seq_flds_a2x_ext_fluxes character(CXX) :: seq_flds_a2x_states_to_rof character(CXX) :: seq_flds_a2x_fluxes_to_rof character(CXX) :: seq_flds_x2a_states @@ -259,10 +257,13 @@ module seq_flds_mod !---------------------------------------------------------------------------- character(CXX) :: seq_flds_dom_fields + character(CXX) :: seq_flds_dom_ext_fields character(CXX) :: seq_flds_a2x_fields character(CXX) :: seq_flds_a2x_ext_fields character(CXX) :: seq_flds_a2x_fields_to_rof + character(CXX) :: seq_flds_a2x_ext_fields_to_rof character(CXX) :: seq_flds_x2a_fields + character(CXX) :: seq_flds_x2a_ext_fields character(CXX) :: seq_flds_i2x_fields character(CXX) :: seq_flds_x2i_fields character(CXX) :: seq_flds_l2x_fields @@ -390,7 +391,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) !------ namelist ----- character(len=CSS) :: fldname, fldflow - character(len=CSS) :: fldname_ext ! use for moab extensions + type(mct_string) :: mctOStr ! mct string for output outfield logical :: is_state, is_flux integer :: i,n @@ -423,9 +424,6 @@ subroutine seq_flds_set(nmlfile, ID, infodata) character(len=*),parameter :: subname = '(seq_flds_set) ' - type(mct_list) :: temp_list - integer :: size_list, index_list - !------------------------------------------------------------------------------- call seq_comm_setptrs(ID,mpicom=mpicom) @@ -4099,45 +4097,14 @@ subroutine seq_flds_set(nmlfile, ID, infodata) call catFields(seq_flds_x2w_fields, seq_flds_x2w_states, seq_flds_x2w_fluxes) call catFields(seq_flds_o2x_fields_to_rof, seq_flds_o2x_states_to_rof, seq_flds_o2x_fluxes_to_rof) ! form character(CXX) :: seq_flds_a2x_ext_states from seq_flds_a2x_states by adding _ext in each field - ! first form a list - call mct_list_init(temp_list ,seq_flds_a2x_fields) - size_list=mct_list_nitem (temp_list) - seq_flds_a2x_ext_fields='' - do index_list = 1, size_list - call mct_list_get(mctOStr,index_list,temp_list) - fldname = mct_string_toChar(mctOStr) - fldname_ext = trim(fldname)//'_ext' - call seq_flds_add(seq_flds_a2x_ext_fields,trim(fldname_ext)) - enddo - seq_flds_a2x_ext_fluxes='' - call mct_list_clean(temp_list) - ! seq_flds_a2x_fluxes - call mct_list_init(temp_list ,seq_flds_a2x_fluxes) - size_list=mct_list_nitem (temp_list) - do index_list = 1, size_list - call mct_list_get(mctOStr,index_list,temp_list) - fldname = mct_string_toChar(mctOStr) - fldname_ext = trim(fldname)//'_ext' - call seq_flds_add(seq_flds_a2x_ext_fluxes,trim(fldname_ext)) - enddo - call mct_list_clean(temp_list) - call mct_list_init(temp_list ,seq_flds_a2x_states) - size_list=mct_list_nitem (temp_list) - seq_flds_a2x_ext_states='' - do index_list = 1, size_list - call mct_list_get(mctOStr,index_list,temp_list) - fldname = mct_string_toChar(mctOStr) - fldname_ext = trim(fldname)//'_ext' - call seq_flds_add(seq_flds_a2x_ext_states,trim(fldname_ext)) - enddo - call mct_list_clean(temp_list) + call create_moab_ext_fields(seq_flds_a2x_fields,seq_flds_a2x_ext_fields) + call create_moab_ext_fields(seq_flds_x2a_fields,seq_flds_x2a_ext_fields) + call create_moab_ext_fields(seq_flds_dom_fields,seq_flds_dom_ext_fields) + call create_moab_ext_fields(seq_flds_a2x_ext_fields_to_rof, seq_flds_a2x_fields_to_rof) - if (seq_comm_iamroot(ID)) then write(logunit,*) subname//': seq_flds_dom_fields= ',trim(seq_flds_dom_fields) - write(logunit,*) subname//': seq_flds_a2x_ext_states= ',trim(seq_flds_a2x_ext_states) - write(logunit,*) subname//': seq_flds_a2x_ext_fluxes= ',trim(seq_flds_a2x_ext_fluxes) write(logunit,*) subname//': seq_flds_a2x_ext_fields= ',trim(seq_flds_a2x_ext_fields) write(logunit,*) subname//': seq_flds_xao_fields= ',trim(seq_flds_xao_fields) endif @@ -4557,4 +4524,31 @@ subroutine moab_set_tag_from_av(tagname, avx, index, mbapid, dataarr, lsize) end subroutine moab_set_tag_from_av + subroutine create_moab_ext_fields(initial_fields, extented_fields) +! !USES: + implicit none + + ! !INPUT/OUTPUT PARAMETERS: + character(CXX), intent(in) :: initial_fields + character(CXX), intent(out) :: extented_fields + + type(mct_list) :: temp_list + character(len=CSS) :: fldname + character(len=CSS) :: fldname_ext ! use for moab extensions + type(mct_string) :: mctOStr ! mct string for output outfield + integer :: size_list, index_list + + ! first form a list + call mct_list_init(temp_list ,initial_fields) + size_list=mct_list_nitem (temp_list) + extented_fields='' + do index_list = 1, size_list + call mct_list_get(mctOStr,index_list,temp_list) + fldname = mct_string_toChar(mctOStr) + fldname_ext = trim(fldname)//'_ext' + call seq_flds_add(extented_fields,trim(fldname_ext)) + enddo + call mct_list_clean(temp_list) + + end subroutine create_moab_ext_fields end module seq_flds_mod From 9d603138704e33000c2acf50b9cf027caa9280fe Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Thu, 18 Sep 2025 18:48:22 -0500 Subject: [PATCH 05/29] get rid of _ext tags they are a relic from the past we were supporting atm spectral to ocean projection, using tempestremap moved to pg2 mesh, so spectral will be used only in monogrid scenarions there should be no mapping in spectral cases (ne4np4, for exanple) All grids are related to atm np4 gll points; either dual mesh, or culled dual mesh on coupler side (land) atm will be on coupler side just point cloud (migrated from comp side, directly) --- driver-moab/main/cplcomp_exchange_mod.F90 | 65 ++++++++++------------- driver-moab/shr/seq_flds_mod.F90 | 40 +------------- 2 files changed, 30 insertions(+), 75 deletions(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 1fa5273a3ce3..6796cd6d899b 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -9,8 +9,8 @@ module cplcomp_exchange_mod use seq_map_type_mod use component_type_mod use seq_flds_mod, only: seq_flds_dom_coord, seq_flds_dom_other - use seq_flds_mod, only: seq_flds_dom_fields, seq_flds_dom_ext_fields - use seq_flds_mod, only: seq_flds_a2x_ext_fields, seq_flds_a2x_fields, seq_flds_x2a_ext_fields, seq_flds_x2a_fields ! + use seq_flds_mod, only: seq_flds_dom_fields + use seq_flds_mod, only: seq_flds_a2x_fields, seq_flds_x2a_fields ! use seq_flds_mod, only: seq_flds_o2x_fields ! needed for MOAB init of ocean fields o2x to be able to transfer to coupler use seq_flds_mod, only: seq_flds_x2o_fields ! needed for MOAB init of ocean fields x2o to be able to transfer from coupler use seq_flds_mod, only: seq_flds_i2x_fields, seq_flds_x2i_fields ! needed for MOAB init of ice fields x2o on coupler side, to save them @@ -1126,7 +1126,7 @@ subroutine cplcomp_moab_Init(infodata,comp) ierr = iMOAB_SendMesh(mhpgid, mpicom_join, mpigrp_cplid, id_join, partMethod) else ! still use the mhid, original coarse mesh - ierr = iMOAB_SendMesh(mhid, mpicom_join, mpigrp_cplid, id_join, partMethod) + ierr = iMOAB_SendMesh(mphaid, mpicom_join, mpigrp_cplid, id_join, partMethod) endif if (ierr .ne. 0) then write(logunit,*) subname,' error in sending mesh from atm comp ' @@ -1145,6 +1145,9 @@ subroutine cplcomp_moab_Init(infodata,comp) endif !!!! FULL ATM if ( trim(atm_mesh) == 'none' ) then ! full atm + ! will receive either pg2 mesh, or point cloud mesh corresponding to GLL points + ! (mphaid app) for spectral case + ! this cannot be used for maps (either computed online or read) ierr = iMOAB_ReceiveMesh(mbaxid, mpicom_join, mpigrp_old, id_old) if (ierr .ne. 0) then write(logunit,*) subname,' error in receiving mesh on atm coupler ' @@ -1181,14 +1184,15 @@ subroutine cplcomp_moab_Init(infodata,comp) ! iMOAB_FreeSenderBuffers needs to be called after receiving the mesh !!!!!!!! ATM COMPONENT - if (mhid .ge. 0) then ! we are on component atm pes + if (mphaid .ge. 0) then ! we are on component atm pes !!!!! FULL ATM if ( trim(atm_mesh) == 'none' ) then ! full atmosphere context_id = id_join if (atm_pg_active) then! we send mesh from mhpgid app ierr = iMOAB_FreeSenderBuffers(mhpgid, context_id) else - ierr = iMOAB_FreeSenderBuffers(mhid, context_id) + ! we send mesh from point cloud data + ierr = iMOAB_FreeSenderBuffers(mphaid, context_id) endif if (ierr .ne. 0) then write(logunit,*) subname,' error in freeing send buffers ' @@ -1201,7 +1205,7 @@ subroutine cplcomp_moab_Init(infodata,comp) ! graph between atm phys, mphaid, and atm dyn on coupler, mbaxid ! phys atm group is mpigrp_old, coupler group is mpigrp_cplid typeA = 2 ! point cloud for mphaid - typeB = 1 ! spectral elements + typeB = 2 ! in spectral case, we will have just point cloud on coupler PEs if (atm_pg_active) then typeB = 3 ! in this case, we will have cells associated with DOFs as GLOBAL_ID tag endif @@ -1209,6 +1213,8 @@ subroutine cplcomp_moab_Init(infodata,comp) ! components/cam/src/cpl/atm_comp_mct.F90 ! components/data_comps/datm/src/atm_comp_mct.F90 ! line 177 !! + ! this is not needed for migrating point cloud to point cloud ! + ! it is needed only after migrating pg2 mesh to cpupler ierr = iMOAB_ComputeCommGraph( mphaid, mbaxid, mpicom_join, mpigrp_old, mpigrp_cplid, & typeA, typeB, ATM_PHYS_CID, id_join) ! ID_JOIN is now 6 @@ -1218,26 +1224,17 @@ subroutine cplcomp_moab_Init(infodata,comp) if (mbaxid .ge. 0 ) then ! coupler pes tagtype = 1 ! dense, double - if (atm_pg_active) then - tagname = trim(seq_flds_a2x_fields)//C_NULL_CHAR - numco = 1 ! usually 1 value per cell - else ! this is not supported now, but leave it here - tagname = trim(seq_flds_a2x_ext_fields)//C_NULL_CHAR ! MOAB versions of a2x for spectral - numco = 16 ! np*np ! usually 16 values per cell, GLL points; should be 4 x 4 = 16 - endif + tagname = trim(seq_flds_a2x_fields)//C_NULL_CHAR + numco = 1 ! usually 1 value per cell ierr = iMOAB_DefineTagStorage(mbaxid, tagname, tagtype, numco, tagindex ) if (ierr .ne. 0) then write(logunit,*) subname,' error in defining tags on atm on coupler ' call shr_sys_abort(subname//' ERROR in defining tags ') endif - if (atm_pg_active) then - tagname = trim(seq_flds_x2a_fields)//C_NULL_CHAR - numco = 1 ! usually 1 value per cell - else ! this is not supported now, but leave it here - tagname = trim(seq_flds_x2a_ext_fields)//C_NULL_CHAR ! MOAB versions of a2x for spectral - numco = 16 ! np*np ! usually 16 values per cell, GLL points; should be 4 x 4 = 16 - endif + + tagname = trim(seq_flds_x2a_fields)//C_NULL_CHAR + numco = 1 ! usually 1 value per cell ierr = iMOAB_DefineTagStorage(mbaxid, tagname, tagtype, numco, tagindex ) if (ierr .ne. 0) then @@ -1246,28 +1243,24 @@ subroutine cplcomp_moab_Init(infodata,comp) endif !add the normalization tag - if (atm_pg_active) then - tagname = trim(seq_flds_dom_fields)//C_NULL_CHAR - numco = 1 ! usually 1 value per cell - else ! this is not supported now, but leave it here - tagname = trim(seq_flds_dom_ext_fields)//C_NULL_CHAR ! MOAB versions of a2x for spectral - numco = 16 ! np*np ! usually 16 values per cell, GLL points; should be 4 x 4 = 16 - endif + + tagname = trim(seq_flds_dom_fields)//C_NULL_CHAR + numco = 1 ! usually 1 value per cell ierr = iMOAB_DefineTagStorage(mbaxid, tagname, tagtype, numco, tagindex ) if (ierr .ne. 0) then write(logunit,*) subname,' error in defining tags seq_flds_dom_fields on atm on coupler ' call shr_sys_abort(subname//' ERROR in defining tags ') endif - if (atm_pg_active) then - ! also, frac, area, masks has to come from atm mphaid, not from domain file reader - ! this is hard to digest :( - tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR - ! TODO: this should be called on the joint procs, not coupler only. - call component_exch_moab(comp, mphaid, mbaxid, 0, tagname, context_exch='doma') - ! copy aream from area in case atm_mesh - call copy_aream_from_area(mbaxid) - endif + + ! also, frac, area, masks has to come from atm mphaid, not from domain file reader + ! this is hard to digest :( + tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR + ! TODO: this should be called on the joint procs, not coupler only. + call component_exch_moab(comp, mphaid, mbaxid, 0, tagname, context_exch='doma') + ! copy aream from area in case atm_mesh + call copy_aream_from_area(mbaxid) + endif ! coupler pes #ifdef MOABDEBUG diff --git a/driver-moab/shr/seq_flds_mod.F90 b/driver-moab/shr/seq_flds_mod.F90 index efb75f41d3d7..9aa3be8d5b1a 100644 --- a/driver-moab/shr/seq_flds_mod.F90 +++ b/driver-moab/shr/seq_flds_mod.F90 @@ -257,13 +257,9 @@ module seq_flds_mod !---------------------------------------------------------------------------- character(CXX) :: seq_flds_dom_fields - character(CXX) :: seq_flds_dom_ext_fields character(CXX) :: seq_flds_a2x_fields - character(CXX) :: seq_flds_a2x_ext_fields character(CXX) :: seq_flds_a2x_fields_to_rof - character(CXX) :: seq_flds_a2x_ext_fields_to_rof character(CXX) :: seq_flds_x2a_fields - character(CXX) :: seq_flds_x2a_ext_fields character(CXX) :: seq_flds_i2x_fields character(CXX) :: seq_flds_x2i_fields character(CXX) :: seq_flds_l2x_fields @@ -4096,16 +4092,9 @@ subroutine seq_flds_set(nmlfile, ID, infodata) call catFields(seq_flds_w2x_fields, seq_flds_w2x_states, seq_flds_w2x_fluxes) call catFields(seq_flds_x2w_fields, seq_flds_x2w_states, seq_flds_x2w_fluxes) call catFields(seq_flds_o2x_fields_to_rof, seq_flds_o2x_states_to_rof, seq_flds_o2x_fluxes_to_rof) - ! form character(CXX) :: seq_flds_a2x_ext_states from seq_flds_a2x_states by adding _ext in each field - - call create_moab_ext_fields(seq_flds_a2x_fields,seq_flds_a2x_ext_fields) - call create_moab_ext_fields(seq_flds_x2a_fields,seq_flds_x2a_ext_fields) - call create_moab_ext_fields(seq_flds_dom_fields,seq_flds_dom_ext_fields) - call create_moab_ext_fields(seq_flds_a2x_ext_fields_to_rof, seq_flds_a2x_fields_to_rof) - + if (seq_comm_iamroot(ID)) then write(logunit,*) subname//': seq_flds_dom_fields= ',trim(seq_flds_dom_fields) - write(logunit,*) subname//': seq_flds_a2x_ext_fields= ',trim(seq_flds_a2x_ext_fields) write(logunit,*) subname//': seq_flds_xao_fields= ',trim(seq_flds_xao_fields) endif @@ -4524,31 +4513,4 @@ subroutine moab_set_tag_from_av(tagname, avx, index, mbapid, dataarr, lsize) end subroutine moab_set_tag_from_av - subroutine create_moab_ext_fields(initial_fields, extented_fields) -! !USES: - implicit none - - ! !INPUT/OUTPUT PARAMETERS: - character(CXX), intent(in) :: initial_fields - character(CXX), intent(out) :: extented_fields - - type(mct_list) :: temp_list - character(len=CSS) :: fldname - character(len=CSS) :: fldname_ext ! use for moab extensions - type(mct_string) :: mctOStr ! mct string for output outfield - integer :: size_list, index_list - - ! first form a list - call mct_list_init(temp_list ,initial_fields) - size_list=mct_list_nitem (temp_list) - extented_fields='' - do index_list = 1, size_list - call mct_list_get(mctOStr,index_list,temp_list) - fldname = mct_string_toChar(mctOStr) - fldname_ext = trim(fldname)//'_ext' - call seq_flds_add(extented_fields,trim(fldname_ext)) - enddo - call mct_list_clean(temp_list) - - end subroutine create_moab_ext_fields end module seq_flds_mod From c2cfd434f3593b4924779022aa4694e5388f31a2 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Fri, 19 Sep 2025 14:59:41 -0500 Subject: [PATCH 06/29] initialize aream for land too it may be left 0 for spectral case; no maps, so no aream to read from anywhere --- driver-moab/main/cplcomp_exchange_mod.F90 | 35 ++++++++++++++++------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 6796cd6d899b..83d411483fa7 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -995,12 +995,18 @@ subroutine copy_aream_from_area(mbappid) integer :: ierr, ent_type ! copy aream from area - if (mbappid >= 0) then ! coupler atm procs + if (mbappid >= 0) then ! coupler procs ierr = iMOAB_GetMeshInfo ( mbappid, nvert, nvise, nbl, nsurf, nvisBC ) - arrSize = nvise(1) ! cells + if ( (mbappid .eq. mbaxid) .and. (.not. atm_pg_active)) then + ! this is the spectral case, for atm only + arrSize = nvert(1) ! cells + ent_type = 0 ! vertices + else + arrSize = nvise(1) ! cells + ent_type = 1 ! cells + endif allocate(tagValues(arrSize)) tagname = 'area'//C_NULL_CHAR - ent_type = 1 ! cells ierr = iMOAB_GetDoubleTagStorage( mbappid, tagname, arrsize , ent_type, tagValues ) tagname = 'aream'//C_NULL_CHAR ierr = iMOAB_SetDoubleTagStorage( mbappid, tagname, arrsize , ent_type, tagValues ) @@ -1252,12 +1258,13 @@ subroutine cplcomp_moab_Init(infodata,comp) write(logunit,*) subname,' error in defining tags seq_flds_dom_fields on atm on coupler ' call shr_sys_abort(subname//' ERROR in defining tags ') endif - - ! also, frac, area, masks has to come from atm mphaid, not from domain file reader - ! this is hard to digest :( - tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR - ! TODO: this should be called on the joint procs, not coupler only. - call component_exch_moab(comp, mphaid, mbaxid, 0, tagname, context_exch='doma') + endif ! coupler pes + ! also, frac, area, masks has to come from atm mphaid, not from domain file reader + ! this is hard to digest :( + tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR + ! TODO: this should be called on the joint procs, not coupler only. + call component_exch_moab(comp, mphaid, mbaxid, 0, tagname, context_exch='doma') + if (mbaxid .ge. 0 ) then ! coupler pes only ! copy aream from area in case atm_mesh call copy_aream_from_area(mbaxid) @@ -1657,6 +1664,12 @@ subroutine cplcomp_moab_Init(infodata,comp) tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR call component_exch_moab(comp, mlnid, mblxid, 0, tagname, context_exch='doml') + ! initialize aream with area; in some cases lnd will not have any maps, so aream still need to + ! have some values (not only 0) + ! spectral case seems to need this, as everything is monogrid + if (MPI_COMM_NULL /= mpicom_new ) then ! we are on the coupler pes + call copy_aream_from_area(mblxid) + endif #ifdef MOABDEBUG outfile = 'recMeshLand.h5m'//C_NULL_CHAR wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! @@ -1906,7 +1919,9 @@ subroutine cplcomp_moab_Init(infodata,comp) call component_exch_moab(comp, mrofid, mbrxid, 0, tagname, context_exch='domr') ! copy aream from area in all cases ! initialize aream from area; it may have different values in the end, or reset again - call copy_aream_from_area(mbrxid) + if (mbrxid > 0) then ! on coupler pes only + call copy_aream_from_area(mbrxid) + endif #ifdef MOABDEBUG if (mbrxid >= 0) then outfile = 'recMeshRof.h5m'//C_NULL_CHAR From 23cfa9adbb64fcde0f384fdad467273b1616737c Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Fri, 19 Sep 2025 15:55:26 -0500 Subject: [PATCH 07/29] do not project aream for lnd for spectral case in general, aream for lnd is projected from atm in same grid cases, spectral (monogrid), do not need to project anymore --- driver-moab/main/component_mod.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/driver-moab/main/component_mod.F90 b/driver-moab/main/component_mod.F90 index 60f94f6992b3..fe02ddb032b6 100644 --- a/driver-moab/main/component_mod.F90 +++ b/driver-moab/main/component_mod.F90 @@ -467,6 +467,7 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, use seq_comm_mct, only: mblxid ! iMOAB id for lnd migrated mesh to coupler pes use seq_comm_mct, only: mbaxid ! iMOAB id for atm migrated mesh to coupler pes use seq_comm_mct, only: mbrxid ! iMOAB id for rof migrated mesh to coupler pes + use seq_comm_mct, only: atm_pg_active ! ! ! Arguments type (seq_infodata_type) , intent(inout) :: infodata @@ -594,7 +595,9 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, dom_s => component_get_dom_cx(atm(1)) !dom_ax dom_d => component_get_dom_cx(lnd(1)) !dom_lx - call seq_map_map(mapper_Sa2l, av_s=dom_s%data, av_d=dom_d%data, fldlist='aream') + if (atm_pg_active) then + call seq_map_map(mapper_Sa2l, av_s=dom_s%data, av_d=dom_d%data, fldlist='aream') + endif else gsmap_d => component_get_gsmap_cx(lnd(1)) ! gsmap_lx dom_d => component_get_dom_cx(lnd(1)) ! dom_lx From d2c84c4894251fa1d3f61928434ab345ef02a3d8 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Tue, 23 Sep 2025 10:52:13 -0500 Subject: [PATCH 08/29] spectral case fixes atm will be type 2 on coupler side, when atm_pg_active is false it cannot be used for any map read or compute, it can be only monogrid case there are still differences shown with MOABCOMP, so something is still not right --- driver-moab/main/component_mod.F90 | 13 ++++++++----- driver-moab/main/cplcomp_exchange_mod.F90 | 6 ------ driver-moab/main/prep_atm_mod.F90 | 18 +++++++++++++++--- driver-moab/main/prep_lnd_mod.F90 | 2 +- driver-moab/main/prep_ocn_mod.F90 | 2 +- 5 files changed, 25 insertions(+), 16 deletions(-) diff --git a/driver-moab/main/component_mod.F90 b/driver-moab/main/component_mod.F90 index fe02ddb032b6..68c378e92bbb 100644 --- a/driver-moab/main/component_mod.F90 +++ b/driver-moab/main/component_mod.F90 @@ -529,7 +529,12 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, nloc = mct_avect_lsize(dom_s%data) allocate(data1(nloc)) data1 = dom_s%data%rAttr(ka,:) - ent_type = 1 ! element dense double tags + if (atm_pg_active) then + ent_type = 1 ! element dense double tags + else ! this is true only for spectral atm now + ent_type = 0 ! for pure spectral case, the atm is PC on coupler side + endif + allocate(gids(nloc)) gids = dom_s%data%iAttr(mct_aVect_indexIA(dom_s%data,"GlobGridNum"),:) ! ! now set data on the coupler side too @@ -594,10 +599,8 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, if (samegrid_al) then dom_s => component_get_dom_cx(atm(1)) !dom_ax dom_d => component_get_dom_cx(lnd(1)) !dom_lx - - if (atm_pg_active) then - call seq_map_map(mapper_Sa2l, av_s=dom_s%data, av_d=dom_d%data, fldlist='aream') - endif + ! it should work for FV and spectral too + call seq_map_map(mapper_Sa2l, av_s=dom_s%data, av_d=dom_d%data, fldlist='aream') else gsmap_d => component_get_gsmap_cx(lnd(1)) ! gsmap_lx dom_d => component_get_dom_cx(lnd(1)) ! dom_lx diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 83d411483fa7..262237bb7912 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -1664,12 +1664,6 @@ subroutine cplcomp_moab_Init(infodata,comp) tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR call component_exch_moab(comp, mlnid, mblxid, 0, tagname, context_exch='doml') - ! initialize aream with area; in some cases lnd will not have any maps, so aream still need to - ! have some values (not only 0) - ! spectral case seems to need this, as everything is monogrid - if (MPI_COMM_NULL /= mpicom_new ) then ! we are on the coupler pes - call copy_aream_from_area(mblxid) - endif #ifdef MOABDEBUG outfile = 'recMeshLand.h5m'//C_NULL_CHAR wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! diff --git a/driver-moab/main/prep_atm_mod.F90 b/driver-moab/main/prep_atm_mod.F90 index d5a923722497..687b679bacf6 100644 --- a/driver-moab/main/prep_atm_mod.F90 +++ b/driver-moab/main/prep_atm_mod.F90 @@ -1088,7 +1088,11 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) write(logunit,*) subname,' error in getting info ' call shr_sys_abort(subname//' error in getting info ') endif - lsize = nvise(1) ! number of active cells + if (atm_pg_active) then + lsize = nvise(1) ! number of active cells + else + lsize = nvert(1) ! for the spectral case, everything is pc from now on + endif ! mct avs are used just for their fields metadata, not the actual reals ! (name of the fields) ! need these always, not only the first time @@ -1296,7 +1300,11 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) endif ! end first-time ! Get data from MOAB - ent_type = 1 ! cells + if (atm_pg_active) then + ent_type = 1 ! cells + else + ent_type = 0 ! vertices, it is PC + endif tagname = trim(seq_flds_x2a_fields)//C_NULL_CHAR arrsize = naflds * lsize ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, x2a_am) @@ -1536,7 +1544,11 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) ! 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 = 1 ! cell for atm, atm_pg_active + if (atm_pg_active) then + ent_type = 1 ! cell for atm, atm_pg_active + else + ent_type = 0 ! vertices is spectral case + endif if (iamroot) print *, subname, num_moab_exports, trim(seq_flds_x2a_fields) do index_list = 1, size_list call mct_list_get(mctOStr,index_list,temp_list) diff --git a/driver-moab/main/prep_lnd_mod.F90 b/driver-moab/main/prep_lnd_mod.F90 index 08d208d1c60e..bd737a1fccf5 100644 --- a/driver-moab/main/prep_lnd_mod.F90 +++ b/driver-moab/main/prep_lnd_mod.F90 @@ -592,7 +592,7 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln if (atm_pg_active) then type1 = 3; ! fv for atm; cgll does not work anyway else - type1 = 1 ! this projection works (cgll to fv), but reverse does not ( fv - cgll) + type1 = 2 ! in this case, atm is just PC endif type2 = 3; ! FV mesh on coupler land ierr = iMOAB_ComputeCommGraph( mbaxid, mblxid, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & diff --git a/driver-moab/main/prep_ocn_mod.F90 b/driver-moab/main/prep_ocn_mod.F90 index efa90630d3fd..806abdd3e696 100644 --- a/driver-moab/main/prep_ocn_mod.F90 +++ b/driver-moab/main/prep_ocn_mod.F90 @@ -586,7 +586,7 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc if (atm_pg_active) then type1 = 3; ! FV for ATM; CGLL does not work correctly in parallel at the moment else - type1 = 1 ! This projection works (CGLL to FV), but reverse does not (FV - CGLL) + type1 = 2 ! in the spectral case, the type on coupler will be point cloud endif type2 = 3; ! FV mesh on coupler OCN From 7c93be544f14e2cd4af05757c7af00e5246d1b7d Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Tue, 23 Sep 2025 15:06:03 -0500 Subject: [PATCH 09/29] more spectral fixes the utilities in shr_moab_mod.F90 need to know if spectral atm is involved; then, instead of cells, use vertices nb and tags --- driver-moab/main/prep_aoflux_mod.F90 | 11 ++++++++-- driver-moab/main/prep_atm_mod.F90 | 32 +++++++++++++++++++--------- driver-moab/shr/shr_moab_mod.F90 | 21 +++++++++++++++--- 3 files changed, 49 insertions(+), 15 deletions(-) diff --git a/driver-moab/main/prep_aoflux_mod.F90 b/driver-moab/main/prep_aoflux_mod.F90 index 958403da0e55..3e0b150632b9 100644 --- a/driver-moab/main/prep_aoflux_mod.F90 +++ b/driver-moab/main/prep_aoflux_mod.F90 @@ -12,6 +12,7 @@ module prep_aoflux_mod use seq_comm_mct, only : mbox2id ! used only for debugging ocn and mct #endif use seq_comm_mct, only : mbaxid ! iMOAB app id for atm on cpl pes + use seq_comm_mct, only: atm_pg_active ! whether the atm uses FV mesh or not ; made true if fv_nphys > 0 use seq_comm_mct, only: seq_comm_getData=>seq_comm_setptrs use seq_comm_mct, only : num_moab_exports use seq_infodata_mod, only: seq_infodata_getdata, seq_infodata_type @@ -227,9 +228,15 @@ subroutine prep_aoflux_init (infodata) ! find out the number of local elements in moab mesh ierr = iMOAB_GetMeshInfo ( mbaxid, nvert, nvise, nbl, nsurf, nvisBC ); ! could be different of lsize_o ! local size of vertices is different from lsize_o - arrSize = nvise(1) * size_list ! there are size_list tags that need to be zeroed out + if(atm_pg_active) then + arrSize = nvise(1) * size_list ! there are size_list tags that need to be zeroed out + ent_type = 1 ! cell type now, not a point cloud anymore + else + arrSize = nvert(1) * size_list + ent_type = 0 ! vertex type now, point cloud + endif allocate(tagValues(arrSize) ) - ent_type = 1 ! cell type now, not a point cloud anymore + tagValues = 0._r8 ierr = iMOAB_SetDoubleTagStorage ( mbaxid, tagname, arrSize , ent_type, tagValues) deallocate(tagValues) diff --git a/driver-moab/main/prep_atm_mod.F90 b/driver-moab/main/prep_atm_mod.F90 index 687b679bacf6..295efb8aeb87 100644 --- a/driver-moab/main/prep_atm_mod.F90 +++ b/driver-moab/main/prep_atm_mod.F90 @@ -393,7 +393,7 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at if (atm_pg_active) then type2 = 3; ! FV for ATM; CGLL does not work correctly in parallel at the moment else - type2 = 1 ! This projection works (CGLL to FV), but reverse does not (FV - CGLL) + type2 = 2 ! from now on, spectral is on PC on coupler side, too; no mapping allowed, just reorder? endif ierr = iMOAB_ComputeCommGraph( mboxid, mbaxid, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & ocn(1)%cplcompid, atm(1)%cplcompid ) @@ -515,14 +515,19 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at mapper_Fof2a%weight_identifier = wgtIdFo2a mapper_Fof2a%mbname = 'mapper_Fof2a' - type1 = 3; ! fv for ocean and atm; fv-cgll does not work anyway - type2 = 3; + type1 = 3; ! fv for ocean and atm; + if (atm_pg_active) then ! + type2 = 3 + else + type2 = 2 ! PC cloud + endif if (.not. samegrid_ao) then ! data-OCN case ! we use the same intx, because the mesh will be the same, between mbofxid and mboxid ierr = iMOAB_ComputeCommGraph( mbofxid, mbintxoa, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & context_id, idintx) else ! this is a case appearing in the data ocean case --res ne4pg2_ne4pg2 --compset FAQP + ! also in spectral case, monogrid --res ne4_ne4 --compset F2010-SCREAMv1 ( type2 is 2, point cloud ) ierr = iMOAB_ComputeCommGraph( mbofxid, mbaxid, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & context_id, atm(1)%cplcompid ) if (ierr .ne. 0) then @@ -748,8 +753,12 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at mapper_Fi2a%weight_identifier = wgtIdFi2a mapper_Fi2a%mbname = 'mapper_Fi2a' if ( samegrid_ao ) then ! this case can appear in cice case - type1 = 3; ! fv for ice and atm; - type2 = 3; + type1 = 3 ! fv for ice + if (atm_pg_active) then + type2 = 3 + else + type2 = 2 ! this is spectral case , PC cloud for atm + endif ierr = iMOAB_ComputeCommGraph( mbixid, mbaxid, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & ice(1)%cplcompid, atm(1)%cplcompid ) if (ierr .ne. 0) then @@ -905,8 +914,8 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at endif - type1 = 3; ! fv for lnd and atm; fv-cgll does not work anyway - type2 = 3; + type1 = 3 ! fv for lnd and atm; fv-cgll does not work anyway + type2 = 3 ierr = iMOAB_ComputeCommGraph( mblxid, mbintxla, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & lnd(1)%cplcompid, idintx) if (ierr .ne. 0) then @@ -919,8 +928,12 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at ! so we compute just a comm graph, between lnd and atm dofs, on the coupler; target is atm ! land is point cloud in this case, type1 = 2 call seq_comm_getinfo(CPLID, mpigrp=mpigrp_CPLID) ! make sure we have the right MPI group - type1 = 3; ! full mesh for land now - type2 = 3; ! fv for target atm + type1 = 3 ! full mesh for land now + if (atm_pg_active) then + type2 = 3 ! fv for target atm + else + type2 = 2 ! point cloud for spectral + endif ierr = iMOAB_ComputeCommGraph( mblxid, mbaxid, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & lnd(1)%cplcompid, atm(1)%cplcompid) if (ierr .ne. 0) then @@ -1561,7 +1574,6 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) ! loop over all fields in seq_flds_o2x_fields call mct_list_init(temp_list ,seq_flds_o2x_fields) size_list=mct_list_nitem (temp_list) - ent_type = 1 ! cell for atm, atm_pg_active if (iamroot) print *, subname, num_moab_exports, trim(seq_flds_o2x_fields) do index_list = 1, size_list call mct_list_get(mctOStr,index_list,temp_list) diff --git a/driver-moab/shr/shr_moab_mod.F90 b/driver-moab/shr/shr_moab_mod.F90 index f07d4da4ad6a..b51a0def95bd 100644 --- a/driver-moab/shr/shr_moab_mod.F90 +++ b/driver-moab/shr/shr_moab_mod.F90 @@ -7,6 +7,7 @@ module shr_moab_mod use shr_sys_mod, only: shr_sys_abort use seq_comm_mct, only: logunit + use seq_comm_mct, only: atm_pg_active, mbaxid use shr_kind_mod, only: CXX => shr_kind_CXX use shr_kind_mod , only: R8 => SHR_KIND_R8 use iso_c_binding @@ -53,7 +54,12 @@ integer function mbGetnCells(moabid) call shr_sys_abort(subname//' error in getting info ') endif - mbGetnCells = nvise(1) + ! only case on coupler side that we actually want number of vertices is when we use spectral atm + if ( .not. atm_pg_active .and. mbaxid .eq. moabid) then + mbGetnCells = nvert(1) + else + mbGetnCells = nvise(1) + endif end function mbGetnCells @@ -85,7 +91,12 @@ subroutine mbGetCellTagVals(mbid, intag,inarray,nMax) !----------------------------------------------------------------------- ! - ent_type=1 + if ( .not. atm_pg_active .and. mbaxid .eq. mbid) then + ent_type = 0 + else + ent_type = 1 + endif + tagname = trim(intag)//C_NULL_CHAR ierr = iMOAB_GetDoubleTagStorage (mbid, tagname, nMax , ent_type, inarray) if (ierr .ne. 0) then @@ -123,7 +134,11 @@ subroutine mbSetCellTagVals(mbid, intag,inarray,nMax) !----------------------------------------------------------------------- ! - ent_type=1 + if ( .not. atm_pg_active .and. mbaxid .eq. mbid) then + ent_type = 0 + else + ent_type = 1 + endif tagname = trim(intag)//C_NULL_CHAR ierr = iMOAB_SetDoubleTagStorage (mbid, tagname, nMax , ent_type, inarray) if (ierr .ne. 0) then From d20ec582ecd5fb3d9836ddd395acba6630cedd04 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Wed, 24 Sep 2025 11:59:40 -0500 Subject: [PATCH 10/29] io needs also spectral case changes the only coupler instance that uses vertices for primary data is the atm in spectral case anoter fix needed is for regular pg2 case; not working anymore --- driver-moab/main/seq_io_mod.F90 | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/driver-moab/main/seq_io_mod.F90 b/driver-moab/main/seq_io_mod.F90 index 6831e21dbce7..7af70c1e840f 100644 --- a/driver-moab/main/seq_io_mod.F90 +++ b/driver-moab/main/seq_io_mod.F90 @@ -33,6 +33,7 @@ module seq_io_mod use shr_sys_mod, only: shr_sys_abort use seq_comm_mct, only: logunit, CPLID, seq_comm_setptrs use seq_comm_mct, only: seq_comm_namelen, seq_comm_name + use seq_comm_mct, only: mbaxid, atm_pg_active use seq_flds_mod, only : seq_flds_lookup use mct_mod ! mct wrappers use pio @@ -1615,7 +1616,8 @@ end subroutine seq_io_write_time ! file_ind [in, optional] - File index for multi-file support ! ! !NOTES: - ! - Only cell-type entities are supported (ent_type=1). + ! - Only cell-type entities are supported (ent_type=1). + ! - not anymore: spectral case, atm needs vertices ent_type=0 ! - Handles reordering of local/global cell IDs for correct output. ! - If matrix is not present, data is read from MOAB tags; if present, matrix is written directly. ! - Skips writing the field "hgt" as a temporary exclusion. @@ -1623,6 +1625,7 @@ end subroutine seq_io_write_time ! ! !REVISION HISTORY: ! 2025-07-20 - Cursor - initial documentation + ! 2025-09-24 allow vertex type for entity, for spectral case for atmosphere ! ! !INTERFACE: ------------------------------------------------------------------ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, matrix, nx, ny, nt, file_ind, dims2din, dims2do, mask ) @@ -1691,7 +1694,6 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, ! should we write a warning? return endif - ent_type = 1 ! cells type lfile_ind = 0 if (present(file_ind)) lfile_ind=file_ind @@ -1700,7 +1702,6 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, call mct_list_init(temp_list ,trim(tag_list)) size_list=mct_list_nitem (temp_list) ! role of nf, number fields - ent_type = 1 ! cell for atm, atm_pg_active if (size_list < 1) then write(logunit,*) subname,' ERROR: size_list = ',size_list,trim(dname) @@ -1727,8 +1728,13 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, ! get the local size ns ierr = iMOAB_GetMeshInfo ( mbxid, nvert, nvise, nbl, nsurf, nvisBC ) - ns = nvise(1) ! local cells - + if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then + ent_type = 0 + ns = nvert(1) ! local vertices + else + ent_type = 1 + ns = nvise(1) ! local cells + endif if (lwhead) then if (present(dims2din)) then dimid2(1)=dims2din(1) @@ -2623,6 +2629,7 @@ end subroutine seq_io_read_char ! ! !NOTES: ! - Only cell-type entities are supported (ent_type=1). + ! - not anymore: spectral case, atm, uses vertices, ent_type = 0 ! - Handles reordering of local/global cell IDs for correct mapping. ! - If matrix is present, data is stored in the matrix; otherwise, data is set as MOAB tags. ! - Skips reading the field "hgt" as a temporary exclusion. @@ -2630,6 +2637,7 @@ end subroutine seq_io_read_char ! ! !REVISION HISTORY: ! 2025-07-20 - Cursor - initial documentation + ! 2025-09-24 spectral case atm ! ! !INTERFACE: ------------------------------------------------------------------ @@ -2687,7 +2695,6 @@ subroutine seq_io_read_moab_tags(filename, mbxid, dname, tag_list, matrix, nx) call mct_list_init(temp_list ,trim(tag_list)) size_list=mct_list_nitem (temp_list) ! role of nf, number fields - ent_type = 1 ! cell for atm, atm_pg_active if (size_list < 1) then write(logunit,*) subname,' ERROR: size_list = ',size_list,trim(dname) @@ -2723,7 +2730,13 @@ subroutine seq_io_read_moab_tags(filename, mbxid, dname, tag_list, matrix, nx) endif lny = 1 ! do we need 2 var, or just 1 ierr = iMOAB_GetMeshInfo ( mbxid, nvert, nvise, nbl, nsurf, nvisBC ) - ns = nvise(1) ! local cells + if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then + ent_type = 0 + ns = nvert(1) ! local vertices + else + ent_type = 1 + ns = nvise(1) ! local cells + endif allocate(data1(ns)) allocate(data_reorder(ns)) allocate(dof(ns)) From 3c12589f9cd9917aae2998ca630e37f1af75cc30 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Wed, 24 Sep 2025 12:57:43 -0500 Subject: [PATCH 11/29] field norm8wt removed by mistake not sure yet when --- driver-moab/main/cplcomp_exchange_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 262237bb7912..95c34f3e39a5 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -1250,7 +1250,7 @@ subroutine cplcomp_moab_Init(infodata,comp) !add the normalization tag - tagname = trim(seq_flds_dom_fields)//C_NULL_CHAR + tagname = trim(seq_flds_dom_fields)//":norm8wt"//C_NULL_CHAR numco = 1 ! usually 1 value per cell ierr = iMOAB_DefineTagStorage(mbaxid, tagname, tagtype, numco, tagindex ) From 8073f1ec7806815f2dcd419213c6916f45a3cfd5 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Thu, 25 Sep 2025 14:33:25 -0500 Subject: [PATCH 12/29] remove spectral mesh only method spectral mesh will not be instanced at all for pure spectral; it is not used at all right now only in the pgx case the spectral mesh is instanced and that is just because it is easier to instance all meshes, not just pgx FV mesh maybe it should be cleaned out completely --- .../homme/interface/phys_grid_mod.F90 | 73 ++- components/homme/src/share/semoab_mod.F90 | 447 ------------------ 2 files changed, 35 insertions(+), 485 deletions(-) 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 7c0abdc73f2c..b4989b0e40fc 100644 --- a/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90 +++ b/components/eamxx/src/dynamics/homme/interface/phys_grid_mod.F90 @@ -552,7 +552,7 @@ subroutine phys_grid_init (pgN) 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_pg_mesh, create_spectral_moab_meshes + use semoab_mod, only: create_moab_pg_mesh use iMOAB, only : iMOAB_RegisterApplication use iso_c_binding #endif @@ -624,48 +624,45 @@ subroutine phys_grid_init (pgN) call compute_global_coords (pg) call compute_global_area (pg) #ifdef HAVE_MOAB - 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 if (pgN > 0) then - - 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 + 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_pg_mesh(par, elem, pgN) - else - call create_spectral_moab_meshes(par, elem) + call create_moab_pg_mesh(par, elem, pgN) endif #endif end subroutine phys_grid_init diff --git a/components/homme/src/share/semoab_mod.F90 b/components/homme/src/share/semoab_mod.F90 index 8d46a441173f..4461c2fb39a7 100644 --- a/components/homme/src/share/semoab_mod.F90 +++ b/components/homme/src/share/semoab_mod.F90 @@ -76,453 +76,6 @@ integer function search_in(intarr, leng, value) end function search_in - subroutine create_spectral_moab_meshes(par, elem) -use ISO_C_BINDING - use iMOAB, only: iMOAB_CreateVertices, iMOAB_WriteMesh, iMOAB_CreateElements, & - iMOAB_ResolveSharedEntities, iMOAB_UpdateMeshInfo, iMOAB_DefineTagStorage, & - iMOAB_SetIntTagStorage, iMOAB_ReduceTagsMax, iMOAB_GetIntTagStorage - use coordinate_systems_mod, only : cartesian3D_t, spherical_to_cart, spherical_polar_t - - type (element_t), intent(inout) :: elem(:) - type (parallel_t) , intent(in) :: par - - integer ierr, i, j, ie, iv, block_ID, k, numvals - integer icol, irow, je, linx ! local indices in fine el connect - - real(kind=real_kind), allocatable :: moab_vert_coords(:) - - integer moab_dim_cquads, ix, idx, nverts, nverts_c ! used for indexing in loops; nverts will have the number of local vertices - - integer nelemd2 ! do not confuse this with dimensions_mod::nelemd - - integer(kind=long_kind), dimension(:), allocatable :: gdofv - ! this will be moab vertex handle locally - integer, dimension(:), allocatable :: moabvh - integer, dimension(:), allocatable :: indx ! this will be ordered - - integer, dimension(:), allocatable :: vdone, elemids, vgids, gdofel - integer, dimension(:), allocatable :: vdone_c, moabconn_c, moabvh_c - integer currentval, dimcoord, dimen, num_el, mbtype, nve - - character*100 outfile, wopts, localmeshfile, lnum, tagname, newtagg - integer tagtype, numco, tag_sto_len, ent_type, tagindex - type (cartesian3D_t) :: cart - integer igcol, ii, neigh - - ! for np=4, - ! 28, 32, 36, 35 - ! 25, 29, 33, 34 - ! j | 13, 17, 21, 22 - ! 1, 5, 9, 10 - !(1,1) i-> - - ! character*100 outfile, wopts, localmeshfile, lnum, tagname - ! integer tagtype, numco, tag_sto_len, ent_type, tagindex - do j=1,np-1 - do i =1, np-1 - ix = (j-1)*(np-1)+i-1 - local_map(i,j) = ix*4 + 1 - enddo - enddo - do j=1, np-1 - i = j - local_map(np, j) = ((np-1)*j-1)*4 + 2 - local_map(i, np) = ( (np-1)*(np-2)+i-1)*4 + 4 - enddo - local_map(np, np) = ((np-1)*(np-1)-1)*4 + 3 - - nelemd2 = (nelemd)*(np-1)*(np-1) - moab_dim_cquads = (nelemd)*4*(np-1)*(np-1) - - if(par%masterproc) then - write (iulog, *) " MOAB: semoab_mod module: create_spectral_moab_meshes on processor " , par%rank ," nelemd: " , nelemd - endif - - if ( nelemd > 0 ) then - allocate(gdofv(moab_dim_cquads)) - allocate(elemids(nelemd2)) - endif - - k=0 ! will be the index for element global dofs - do ie=1,nelemd - do j=1,np-1 - do i=1,np-1 - ix = (ie-1)*(np-1)*(np-1)+(j-1)*(np-1)+i-1 - gdofv(ix*4+1) = elem(ie)%gdofP(i,j) - gdofv(ix*4+2) = elem(ie)%gdofP(i+1,j) - gdofv(ix*4+3) = elem(ie)%gdofP(i+1,j+1) - gdofv(ix*4+4) = elem(ie)%gdofP(i,j+1) - elemids(ix+1) = (elem(ie)%GlobalId-1)*(np-1)*(np-1)+(j-1)*(np-1)+i - enddo - enddo - enddo - -! order according to global dofs - - if ( nelemd > 0 ) then - allocate(moabvh(moab_dim_cquads)) - allocate(indx(moab_dim_cquads)) - - allocate(moabconn(moab_dim_cquads)) - call IndexSet(moab_dim_cquads, indx) - call IndexSort(moab_dim_cquads, indx, gdofv, descend=.false.) -! after sort, gdofv( indx(i)) < gdofv( indx(i+1) ) - endif - idx=0 - currentval = 0 - if ( nelemd > 0 ) then - currentval = gdofv( indx(1)) - idx = 1 - endif - - do ix=1,moab_dim_cquads - if (gdofv(indx(ix)) .ne. currentval ) then - idx=idx+1 - currentval = gdofv(indx(ix)) - endif - moabvh(ix) = idx ! the vertex in connectivity array will be at this local index - ! this will be the moab connectivity - moabconn(indx(ix)) = idx - enddo - - nverts = idx - if(par%masterproc) then - write (iulog, *) " MOAB: there are ", nverts, " local vertices on master task ", currentval, " is the max local gdof" - endif - if ( nelemd > 0 ) then - allocate(moab_vert_coords(3*nverts) ) - allocate(vdone(nverts)) - else - allocate(vdone(1)) ! will be passed to iMOAB_ResolveSharedEnts, and compilers are complaining about non-allocated arrays - endif - vdone = 0; - if ( nelemd > 0 ) currentval = gdofv( indx(1)) ! start over to identify coordinates of the vertices - - do ix=1,moab_dim_cquads - idx = indx(ix) ! index in initial array, vertices in all fine quads - k = (idx-1)/(4*(np-1)*(np-1)) ! index of coarse quad, locally, starts at 0 - ie = 1 + k ! this is the element number; starts at nets=1 - je = ( idx -1 -k*(np-1)*(np-1)*4 ) / 4 + 1 ! local fine quad in coarse, 1 to (np-1) ^ 2 - irow = (je-1)/(np-1)+1 - icol = je - (np-1)*(irow-1) - linx = idx - k*(np-1)*(np-1)*4 -(je-1)*4 ! this should be 1, 2, 3, 4 - if( linx == 1) then - j = irow - i = icol - else if (linx == 2) then - j = irow - i = icol + 1 - else if (linx == 3) then - j = irow + 1 - i = icol + 1 - else ! linx == 4 - j = irow + 1 - i = icol - endif - - iv = moabvh(ix) - if (vdone(iv) .eq. 0) then - cart = spherical_to_cart (elem(ie)%spherep(i,j) ) - moab_vert_coords ( 3*(iv-1)+1 ) = cart%x - moab_vert_coords ( 3*(iv-1)+2 ) = cart%y - moab_vert_coords ( 3*(iv-1)+3 ) = cart%z - vdone(iv) = gdofv(indx(ix)) ! this will be now our tag used for resolving shared entities ! convert to int, from long int - endif - - enddo - - dimcoord = nverts*3 - dimen = 3 - if ( nelemd > 0 ) then - ierr = iMOAB_CreateVertices(MHFID, dimcoord, dimen, moab_vert_coords) - if (ierr > 0 ) & - call abortmp('Error: fail to create MOAB vertices ') - endif - !!num_el = nelemd2 - mbtype = 3 ! quadrilateral - nve = 4; - block_ID = 200 ! this will be for coarse mesh - - if ( nelemd > 0 ) then - ierr = iMOAB_CreateElements( MHFID, nelemd2, mbtype, nve, moabconn, block_ID ); - if (ierr > 0 ) & - call abortmp('Error: fail to create MOAB elements') - endif - ! nverts: num vertices; vdone will store now the markers used in global resolve - ! for this particular problem, markers are the global dofs at corner nodes -! set the global id for vertices -! first, retrieve the tag - tagname='GDOF'//C_NULL_CHAR - tagtype = 0 ! dense, integer - numco = 1 - ierr = iMOAB_DefineTagStorage(MHFID, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call abortmp('Error: fail to retrieve global id tag') - ! now set the values - ent_type = 0 ! vertex type - if ( nverts > 0 ) then - ierr = iMOAB_SetIntTagStorage ( MHFID, tagname, nverts , ent_type, vdone) - if (ierr > 0 ) & - call abortmp('Error: fail to set marker id tag for vertices') - endif - - ! we need to call this even when no mesh locally, it involves a collective - ierr = iMOAB_ResolveSharedEntities( MHFID, nverts, vdone ); - if (ierr > 0 ) & - call abortmp('Error: fail to resolve shared entities') - - if ( nelemd > 0) then - vdone = -1 ! reuse vdone for the new tag, GLOBAL_ID (actual tag that we want to store global dof ) - endif -! use element offset for actual global dofs - ! tagtype = 0 ! dense, integer - ! numco = 1 - newtagg='GLOBAL_ID'//C_NULL_CHAR - ierr = iMOAB_DefineTagStorage(MHFID, newtagg, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call abortmp('Error: fail to create new GDOF tag') - do ie=1,nelemd - do ii=1,elem(ie)%idxp%NumUniquePts - i=elem(ie)%idxp%ia(ii) - j=elem(ie)%idxp%ja(ii) - igcol = elem(ie)%idxp%UniquePtoffset+ii-1 - ix = local_map(i,j) - idx = moabconn((ie-1)*(np-1)*(np-1)*4 + ix) ! should - vdone ( idx ) = igcol - end do - end do - ! now set the values - ent_type = 0 ! vertex type - if ( nverts > 0 ) then - ierr = iMOAB_SetIntTagStorage ( MHFID, newtagg, nverts , ent_type, vdone) - if (ierr > 0 ) & - call abortmp('Error: fail to set global dof tag for vertices') - endif - - ierr = iMOAB_ReduceTagsMax ( MHFID, tagindex, ent_type) - if (ierr > 0 ) & - call abortmp('Error: fail to reduce max tag') - - ! set global id tag for elements - ent_type = 1 ! now set the global id tag on elements - if ( nelemd2 > 0 ) then - ierr = iMOAB_SetIntTagStorage ( MHFID, newtagg, nelemd2 , ent_type, elemids) - if (ierr > 0 ) & - call abortmp('Error: fail to set global id tag for elements') - endif - -! now, after reduction, we can get the actual global ids for each vertex in the fine mesh -! before, some vertices that were owned in MOAB but not owned in CAM did not have the right global ID tag -! so vdone will be now correct on every task (no -1 anymore ) - ent_type = 0 ! vertex type - if ( nverts > 0 ) then - allocate(vgids(nverts)) - ierr = iMOAB_GetIntTagStorage ( MHFID, newtagg, nverts , ent_type, vgids) - if (ierr > 0 ) & - call abortmp('Error: fail to retrieve GLOBAL ID on each task') - endif - ierr = iMOAB_UpdateMeshInfo(MHFID) - if (ierr > 0 ) & - call abortmp('Error: fail to update mesh info') -#ifdef MOABDEBUG -! write out the mesh file to disk, in parallel - outfile = 'wholeFineATM.h5m'//C_NULL_CHAR - wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR - ierr = iMOAB_WriteMesh(MHFID, outfile, wopts) - if (ierr > 0 ) & - call abortmp('Error: fail to write the mesh file') -#endif - - -! now create the coarse mesh, but the global dofs will come from fine mesh, after solving - ! nelemd2 = nelemd - moab_dim_cquads = (nelemd)*4 - - if ( nelemd > 0 ) then - allocate(gdofel(nelemd*np*np)) - endif - k=0 ! will be the index for element global dofs - do ie=1,nelemd - ix = ie-1 - ! - gdofv(ix*4+1) = elem(ie)%gdofP(1,1) - gdofv(ix*4+2) = elem(ie)%gdofP(np,1) - gdofv(ix*4+3) = elem(ie)%gdofP(np,np) - gdofv(ix*4+4) = elem(ie)%gdofP(1,np) - elemids(ix+1) = elem(ie)%GlobalId - enddo -! now original order - -! order according to global dofs -! allocate(indx(moab_dim_cquads)) - if ( nelemd > 0 ) then - call IndexSet(moab_dim_cquads, indx) - call IndexSort(moab_dim_cquads, indx, gdofv, descend=.false.) -! after sort, gdofv( indx(i)) < gdofv( indx(i+1) ) - - allocate(moabvh_c(moab_dim_cquads)) - - allocate(moabconn_c(moab_dim_cquads)) - endif - idx = 0 - if ( nelemd > 0 ) then - idx=1 - currentval = gdofv( indx(1)) - endif - do ix=1,moab_dim_cquads - if (gdofv(indx(ix)) .ne. currentval ) then - idx=idx+1 - currentval = gdofv(indx(ix)) - endif - moabvh_c(ix) = idx ! the vertex in connectivity array will be at this local index - ! this will be the moab connectivity - moabconn_c(indx(ix)) = idx - enddo - nverts_c = idx - if(par%masterproc) then - write (iulog, *) " MOAB: there are ", nverts_c, " local vertices on master task, coarse mesh" - endif -! allocate(moab_vert_coords(3*idx) ) - if ( nelemd > 0 ) then - allocate(vdone_c(nverts_c)) - vdone_c = 0; - currentval = gdofv( indx(1)) ! start over to identify coordinates of the vertices - endif - - do ix=1,moab_dim_cquads - i = indx(ix) ! index in initial array - ie = 1+ (i-1)/4 ! this is the element number - j = i - ( i-1)/4*4 ! local index of vertex in element i - iv = moabvh_c(ix) - if (vdone_c(iv) .eq. 0) then - moab_vert_coords ( 3*(iv-1)+1 ) = elem(ie)%corners3d(j)%x - moab_vert_coords ( 3*(iv-1)+2 ) = elem(ie)%corners3d(j)%y - moab_vert_coords ( 3*(iv-1)+3 ) = elem(ie)%corners3d(j)%z - vdone_c(iv) = gdofv(indx(ix)) ! this will be now our tag used for resolving shared entities - endif - - enddo - - dimcoord = nverts_c*3 - dimen = 3 - if ( nverts_c > 0 ) then - ierr = iMOAB_CreateVertices(MHID, dimcoord, dimen, moab_vert_coords) - if (ierr > 0 ) & - call abortmp('Error: fail to create MOAB vertices ') - endif - ! num_el = nelemd - mbtype = 3 ! quadrilateral - nve = 4; - block_ID = 100 ! this will be for coarse mesh - - if ( nelemd > 0 ) then - ierr = iMOAB_CreateElements( MHID, nelemd, mbtype, nve, moabconn_c, block_ID ); - if (ierr > 0 ) & - call abortmp('Error: fail to create MOAB elements') - endif - ! idx: num vertices; vdone will store now the markers used in global resolve - ! for this particular problem, markers are the global dofs at corner nodes -! set the global id for vertices -! first, retrieve the tag - tagname='GDOFV'//C_NULL_CHAR - tagtype = 0 ! dense, integer - numco = 1 - ierr = iMOAB_DefineTagStorage(MHID, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call abortmp('Error: fail to retrieve GDOFV id tag') - ierr = iMOAB_DefineTagStorage(MHID, newtagg, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call abortmp('Error: fail to retrieve GLOBAL_ID tag on coarse mesh') - ! now set the values - ent_type = 0 ! vertex type - if ( nverts_c > 0 ) then - ierr = iMOAB_SetIntTagStorage ( MHID, tagname, nverts_c , ent_type, vdone_c) - if (ierr > 0 ) & - call abortmp('Error: fail to set GDOFV tag for vertices') - endif - ! set global id tag for coarse elements, too; they will start at nets=1, end at nete=nelemd - ent_type = 1 ! now set the global id tag on elements - if ( nelemd > 0 ) then - ierr = iMOAB_SetIntTagStorage ( MHID, newtagg, nelemd , ent_type, elemids) - if (ierr > 0 ) & - call abortmp('Error: fail to set global id tag for vertices') - endif - - ierr = iMOAB_ResolveSharedEntities( MHID, idx, vdone_c ); - if (ierr > 0 ) & - call abortmp('Error: fail to resolve shared entities') - -! global dofs are the GLL points are set for each element - tagname='GLOBAL_DOFS'//C_NULL_CHAR - tagtype = 0 ! dense, integer - numco = np*np ! usually, it is 16; each element will have the dofs in order - ierr = iMOAB_DefineTagStorage(MHID, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call abortmp('Error: fail to create global DOFS tag') - ! now set the values - ! set global dofs tag for coarse elements, too; they will start at nets=1, end at nete=nelemd - ent_type = 1 ! now set the global id tag on elements - numvals = nelemd*np*np ! input is the total number of values - ! form gdofel from vgids - do ie=1, nelemd - ix = (ie-1)*np*np ! ie: index in coarse element - je = (ie-1) * 4 * (np-1) * (np -1) ! index in moabconn array - ! vgids are global ids for fine vertices (1,nverts) - iv = 1 - do j=1,np - do i=1,np - k = local_map(i,j) - gdofel(ix+iv) = vgids( moabconn( je + k ) ) - iv = iv + 1 - enddo - enddo - ! extract global ids - vdone_c( moabconn_c( (ie-1)*4+1) ) = vgids ( moabconn(je+1 )) - vdone_c( moabconn_c( (ie-1)*4+2) ) = vgids ( moabconn(je+ 4*(np-2)+2 )) ! valid for np = 4, 10 - vdone_c( moabconn_c( (ie-1)*4+3) ) = vgids ( moabconn(je+ 4*((np-1)*(np-1)-1) + 3 )) ! for np = 4, 35 - vdone_c( moabconn_c( (ie-1)*4+4) ) = vgids ( moabconn(je+ 4*(np-2)*(np-1) + 4 )) ! 28 for np = 4 - enddo - if ( nelemd > 0 ) then - ierr = iMOAB_SetIntTagStorage ( MHID, tagname, numvals, ent_type, gdofel) - if (ierr > 0 ) & - call abortmp('Error: fail to set globalDOFs tag for coarse elements') - endif - - ! set the global ids for coarse vertices the same as corresponding fine vertices - ent_type = 0 ! vertex type - if ( nverts_c > 0 ) then - ierr = iMOAB_SetIntTagStorage ( MHID, newtagg, nverts_c , ent_type, vdone_c) - if (ierr > 0 ) & - call abortmp('Error: fail to set GLOBAL_DOFS tag values') - endif - - ierr = iMOAB_UpdateMeshInfo(MHID) - if (ierr > 0 ) & - call abortmp('Error: fail to update mesh info') -#ifdef MOABDEBUG -! write out the mesh file to disk, in parallel - outfile = 'wholeATM.h5m'//C_NULL_CHAR - wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR - ierr = iMOAB_WriteMesh(MHID, outfile, wopts) - if (ierr > 0 ) & - call abortmp('Error: fail to write the mesh file') -#endif - ! deallocate - if ( nelemd > 0 ) then - deallocate(moabvh) - deallocate(moabconn) ! do not keep it anymore, we are not setting another tag on fine mesh - deallocate(gdofel) - deallocate(indx) - deallocate(elemids) - deallocate(gdofv) - deallocate(moabvh_c) - deallocate(moabconn_c) - deallocate(vdone_c) - endif - deallocate(vdone) ! we are always allocating this now -! end copy - end subroutine create_spectral_moab_meshes - subroutine create_moab_pg_mesh(par, elem, fv_nphys) use ISO_C_BINDING From 5622cb495e4c7c470320cf56087ae4a81649bb56 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Fri, 26 Sep 2025 15:11:01 -0500 Subject: [PATCH 13/29] for spectral case, global size changes when writing history files, the global size needs to be the number of global vertices, for spectral case ( atm_pg_active false) --- driver-moab/main/seq_io_mod.F90 | 49 +++++++++++++++++++-------------- 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/driver-moab/main/seq_io_mod.F90 b/driver-moab/main/seq_io_mod.F90 index 7af70c1e840f..28fd14053c09 100644 --- a/driver-moab/main/seq_io_mod.F90 +++ b/driver-moab/main/seq_io_mod.F90 @@ -1711,9 +1711,23 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, lpre = trim(dname) ! find out the number of global cells, needed for defining the variables length ierr = iMOAB_GetGlobalInfo( mbxid, dummy, ng) - lnx = ng lny = 1 +#ifdef MOABCOMP + write(logunit,*) subname, 'lnx, lny, mbxid ', lnx, lny, mbxid +#endif + + ! get the local size ns + ierr = iMOAB_GetMeshInfo ( mbxid, nvert, nvise, nbl, nsurf, nvisBC ) + if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then + ent_type = 0 + ns = nvert(1) ! local vertices + lnx = dummy ! number of global vertices + else + ent_type = 1 + ns = nvise(1) ! local cells + lnx = ng + endif ! it is needed to overwrite that for land, ng is too small ! ( for ne4pg2 it is 201 instead of 384) if (present(nx)) then @@ -1725,16 +1739,10 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, if (present(nt)) then frame = nt endif - - ! get the local size ns - ierr = iMOAB_GetMeshInfo ( mbxid, nvert, nvise, nbl, nsurf, nvisBC ) - if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then - ent_type = 0 - ns = nvert(1) ! local vertices - else - ent_type = 1 - ns = nvise(1) ! local cells - endif +#ifdef MOABCOMP + write(logunit,*) subname, ' ent_type, ns ', ent_type, ns + write(logunit,*) subname, ' tag_list ', trim(tag_list) +#endif if (lwhead) then if (present(dims2din)) then dimid2(1)=dims2din(1) @@ -2719,23 +2727,24 @@ subroutine seq_io_read_moab_tags(filename, mbxid, dname, tag_list, matrix, nx) ! find out the number of global cells, needed for defining the variables length ierr = iMOAB_GetGlobalInfo( mbxid, dummy, ng) - lnx = ng - ! it is needed to overwrite that for land, ng is too small - ! ( for ne4pg2 it is 201 instead of 384) - if (present(nx)) then -#ifdef MOABCOMP - if (iam==0) write(logunit,*) subname, ' nx present: ', nx -#endif - lnx = nx - endif lny = 1 ! do we need 2 var, or just 1 ierr = iMOAB_GetMeshInfo ( mbxid, nvert, nvise, nbl, nsurf, nvisBC ) if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then ent_type = 0 ns = nvert(1) ! local vertices + lnx = dummy ! number of global vertices else ent_type = 1 ns = nvise(1) ! local cells + lnx = ng + endif + ! it is needed to overwrite that for land, ng is too small + ! ( for ne4pg2 it is 201 instead of 384) + if (present(nx)) then +#ifdef MOABCOMP + if (iam==0) write(logunit,*) subname, ' nx present: ', nx +#endif + lnx = nx endif allocate(data1(ns)) allocate(data_reorder(ns)) From 789b94a6c798915ed9102cb675b5c03227b1a9d2 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Fri, 3 Oct 2025 13:12:30 -0500 Subject: [PATCH 14/29] replace dummy with ngv ngv is the number of global vertices for the imoab app it is used only for spectral case the other calls iMOAB_GetGlobalInfo need the number of cells, not the number of vertices --- driver-moab/main/seq_hist_mod.F90 | 1 - driver-moab/main/seq_io_mod.F90 | 12 ++++++------ driver-moab/main/seq_rest_mod.F90 | 20 ++++++++++---------- 3 files changed, 16 insertions(+), 17 deletions(-) diff --git a/driver-moab/main/seq_hist_mod.F90 b/driver-moab/main/seq_hist_mod.F90 index c79666ae7b75..b7093df89d90 100644 --- a/driver-moab/main/seq_hist_mod.F90 +++ b/driver-moab/main/seq_hist_mod.F90 @@ -46,7 +46,6 @@ module seq_hist_mod use seq_flds_mod, only: seq_flds_i2x_fields, seq_flds_r2x_fields,seq_flds_dom_fields use seq_flds_mod, only: seq_flds_l2x_fields, seq_flds_x2a_fields, seq_flds_x2i_fields use seq_flds_mod, only: seq_flds_x2l_fields, seq_flds_x2r_fields - use iMOAB, only: iMOAB_GetGlobalInfo use shr_moab_mod, only: mbGetnCells,mbGetCellTagVals use component_type_mod diff --git a/driver-moab/main/seq_io_mod.F90 b/driver-moab/main/seq_io_mod.F90 index 28fd14053c09..a50726a02eb3 100644 --- a/driver-moab/main/seq_io_mod.F90 +++ b/driver-moab/main/seq_io_mod.F90 @@ -1673,7 +1673,7 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, integer(in),target :: dimid2(2) integer(in),target :: dimid3(3) integer(in),pointer :: dimid(:) - integer(in) :: dummy, ent_type, ierr + integer(in) :: ngv, ent_type, ierr real(r8) :: lfillvalue integer, allocatable :: indx(:) ! this will be ordered integer, allocatable :: Dof(:) ! will be filled with global ids from cells @@ -1710,7 +1710,7 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, lpre = trim(dname) ! find out the number of global cells, needed for defining the variables length - ierr = iMOAB_GetGlobalInfo( mbxid, dummy, ng) + ierr = iMOAB_GetGlobalInfo( mbxid, ngv, ng) lny = 1 #ifdef MOABCOMP write(logunit,*) subname, 'lnx, lny, mbxid ', lnx, lny, mbxid @@ -1722,7 +1722,7 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then ent_type = 0 ns = nvert(1) ! local vertices - lnx = dummy ! number of global vertices + lnx = ngv ! number of global vertices else ent_type = 1 ns = nvise(1) ! local cells @@ -2693,7 +2693,7 @@ subroutine seq_io_read_moab_tags(filename, mbxid, dname, tag_list, matrix, nx) type(mct_string) :: mctOStr ! character(CXX) ::tagname, field - integer(in) :: dummy, ent_type, ierr + integer(in) :: ngv, ent_type, ierr character(*),parameter :: subName = '(seq_io_read_moab_tags) ' @@ -2726,13 +2726,13 @@ subroutine seq_io_read_moab_tags(filename, mbxid, dname, tag_list, matrix, nx) endif ! find out the number of global cells, needed for defining the variables length - ierr = iMOAB_GetGlobalInfo( mbxid, dummy, ng) + ierr = iMOAB_GetGlobalInfo( mbxid, ngv, ng) lny = 1 ! do we need 2 var, or just 1 ierr = iMOAB_GetMeshInfo ( mbxid, nvert, nvise, nbl, nsurf, nvisBC ) if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then ent_type = 0 ns = nvert(1) ! local vertices - lnx = dummy ! number of global vertices + lnx = ngv ! number of global vertices else ent_type = 1 ns = nvise(1) ! local cells diff --git a/driver-moab/main/seq_rest_mod.F90 b/driver-moab/main/seq_rest_mod.F90 index 25fbad42bbf6..f20c4c0a185c 100644 --- a/driver-moab/main/seq_rest_mod.F90 +++ b/driver-moab/main/seq_rest_mod.F90 @@ -386,7 +386,7 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) integer (in), pointer :: l2racc_lm_cnt integer (in) :: nx_lnd ! will be used if land and atm are on same grid - integer (in) :: ierr, dummy + integer (in) :: ierr, ngv real(r8), dimension(:,:), pointer :: p_x2oacc_om real(r8), dimension(:,:), pointer :: p_o2racc_om @@ -453,12 +453,12 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) if (lnd_present) then if(samegrid_al) then ! nx for land will be from global nb atmosphere - ierr = iMOAB_GetGlobalInfo(mbaxid, dummy, nx_lnd) ! max id for land will come from atm + ierr = iMOAB_GetGlobalInfo(mbaxid, ngv, nx_lnd) ! max id for land will come from atm call seq_io_read(moab_rest_file, mblxid, 'fractions_lx', & 'afrac:lfrac:lfrin', nx=nx_lnd) else if(samegrid_lr) then ! nx for land will be from global nb rof - ierr = iMOAB_GetGlobalInfo(mbrxid, dummy, nx_lnd) ! max id for land will come from rof + ierr = iMOAB_GetGlobalInfo(mbrxid, ngv, nx_lnd) ! max id for land will come from rof call seq_io_read(moab_rest_file, mblxid, 'fractions_lx', & 'afrac:lfrac:lfrin', nx=nx_lnd) else ! is this ever true ? @@ -474,13 +474,13 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) p_l2racc_lm => prep_rof_get_l2racc_lm() if(samegrid_al) then ! nx for land will be from global nb atmosphere - ierr = iMOAB_GetGlobalInfo(mbaxid, dummy, nx_lnd) ! max id for land will come from atm + ierr = iMOAB_GetGlobalInfo(mbaxid, ngv, nx_lnd) ! max id for land will come from atm call seq_io_read(moab_rest_file, mblxid, 'l2racc_lx', & trim(tagname), & matrix = p_l2racc_lm, nx=nx_lnd) else if(samegrid_lr) then ! nx for land will be from global nb rof - ierr = iMOAB_GetGlobalInfo(mbrxid, dummy, nx_lnd) ! max id for land will come from rof + ierr = iMOAB_GetGlobalInfo(mbrxid, ngv, nx_lnd) ! max id for land will come from rof call seq_io_read(moab_rest_file, mblxid, 'l2racc_lx', & trim(tagname), & matrix = p_l2racc_lm, nx=nx_lnd) @@ -1013,7 +1013,7 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & integer (in), pointer :: l2racc_lm_cnt integer (in) :: nx_lnd ! will be used if land and atm are on same grid - integer (in) :: ierr, dummy + integer (in) :: ierr, ngv real(r8), dimension(:,:), pointer :: p_x2oacc_om real(r8), dimension(:,:), pointer :: p_o2racc_om @@ -1192,13 +1192,13 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & if (lnd_present) then if(samegrid_al) then ! nx for land will be from global nb atmosphere - ierr = iMOAB_GetGlobalInfo(mbaxid, dummy, nx_lnd) ! max id for land will come from atm + ierr = iMOAB_GetGlobalInfo(mbaxid, ngv, nx_lnd) ! max id for land will come from atm call seq_io_write(rest_file, mblxid, 'fractions_lx', & 'afrac:lfrac:lfrin', & ! seq_frac_mod: character(*),parameter :: fraclist_l = 'afrac:lfrac:lfrin' whead=whead, wdata=wdata, nx=nx_lnd) else if(samegrid_lr) then ! nx for land will be from global nb atmosphere - ierr = iMOAB_GetGlobalInfo(mbrxid, dummy, nx_lnd) ! max id for land will come from rof + ierr = iMOAB_GetGlobalInfo(mbrxid, ngv, nx_lnd) ! max id for land will come from rof call seq_io_write(rest_file, mblxid, 'fractions_lx', & 'afrac:lfrac:lfrin', & ! seq_frac_mod: character(*),parameter :: fraclist_l = 'afrac:lfrac:lfrin' whead=whead, wdata=wdata, nx=nx_lnd) @@ -1220,13 +1220,13 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & p_l2racc_lm => prep_rof_get_l2racc_lm() if(samegrid_al) then ! nx for land will be from global nb atmosphere - ierr = iMOAB_GetGlobalInfo(mbaxid, dummy, nx_lnd) ! max id for land will come from atm + ierr = iMOAB_GetGlobalInfo(mbaxid, ngv, nx_lnd) ! max id for land will come from atm call seq_io_write(rest_file, mblxid, 'l2racc_lx', & trim(tagname), & whead=whead, wdata=wdata, matrix = p_l2racc_lm, nx=nx_lnd) else if(samegrid_lr) then ! nx for land will be from global nb atmosphere - ierr = iMOAB_GetGlobalInfo(mbrxid, dummy, nx_lnd) ! max id for land will come from rof + ierr = iMOAB_GetGlobalInfo(mbrxid, ngv, nx_lnd) ! max id for land will come from rof call seq_io_write(rest_file, mblxid, 'l2racc_lx', & trim(tagname), & whead=whead, wdata=wdata, matrix = p_l2racc_lm, nx=nx_lnd) From 393ba2536ba4b07344abd93957f5fc7181651fa0 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Wed, 8 Oct 2025 09:53:53 -0500 Subject: [PATCH 15/29] scm land case fixes introduce mb_scm_land logical it is false, except when land is identified as scm (lnd_nx and lnd_ny are both 1) it is turn on if needed in lnd moab cx init it will trigger land migrate of PC cloud, not a read from domain file also, Sa2l map will involve type vertex for land history and restart for land still have to be fixed --- driver-moab/main/cplcomp_exchange_mod.F90 | 104 +++++++++++++++------- driver-moab/main/prep_lnd_mod.F90 | 7 +- driver-moab/shr/seq_comm_mct.F90 | 1 + 3 files changed, 81 insertions(+), 31 deletions(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 95c34f3e39a5..d13630b9ce13 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -24,6 +24,7 @@ module cplcomp_exchange_mod use seq_comm_mct, only : mhpgid ! iMOAB app id for atm pgx grid, on atm pes use seq_comm_mct, only : atm_pg_active ! flag if PG mesh instanced use seq_comm_mct, only : mlnid , mblxid ! iMOAB app id for land , on land pes and coupler pes + use seq_comm_mct, only : mb_scm_land ! logical used to identify land scm case; moab will migrate land then use seq_comm_mct, only : mphaid ! iMOAB app id for phys atm; comp atm is 5, phys 5+200 use seq_comm_mct, only : MPSIID, mbixid ! sea-ice on comp pes and on coupler pes use seq_comm_mct, only : mrofid, mbrxid ! iMOAB id of moab rof app on comp pes and on coupler too @@ -1066,6 +1067,7 @@ subroutine cplcomp_moab_Init(infodata,comp) character(CXX) :: newlist integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) logical :: rof_present, lnd_prognostic + integer :: land_nx, land_ny ! try to identify if scm case; then land mesh will be migrated, not read real(r8), allocatable :: tagValues(:) ! used for setting aream tags for atm domain read case integer :: arrsize ! for the size of tagValues type(mct_list) :: temp_list @@ -1128,7 +1130,8 @@ subroutine cplcomp_moab_Init(infodata,comp) ! send mesh to coupler !!!! FULL ATM if ( trim(atm_mesh) == 'none' ) then ! full model - if (atm_pg_active) then ! change : send the pg2 mesh, not coarse mesh, when atm pg active + if (atm_pg_active) then ! change : send the point cloud phys grid mesh, not coarse mesh, + ! when atm pg active ierr = iMOAB_SendMesh(mhpgid, mpicom_join, mpigrp_cplid, id_join, partMethod) else ! still use the mhid, original coarse mesh @@ -1556,7 +1559,13 @@ subroutine cplcomp_moab_Init(infodata,comp) call seq_comm_getinfo(cplid ,mpigrp=mpigrp_cplid) ! receiver group call seq_comm_getinfo(id_old,mpigrp=mpigrp_old) ! component group pes call seq_infodata_GetData(infodata,rof_present=rof_present, lnd_prognostic=lnd_prognostic) + call seq_infodata_GetData(infodata,lnd_nx=land_nx, lnd_ny=land_ny) + if (land_nx .eq. 1 .and. land_ny .eq.1 ) then + ! turn on mb_scm_land + mb_scm_land = .true. ! we identified a scm case for land, we will migrate mesh, not read + ! the domain file anymore + endif ! use land full mesh if (MPI_COMM_NULL /= mpicom_new ) then ! we are on the coupler pes appname = "COUPLE_LAND"//C_NULL_CHAR @@ -1566,27 +1575,60 @@ subroutine cplcomp_moab_Init(infodata,comp) write(logunit,*) subname,' error in registering coupler land ' call shr_sys_abort(subname//' ERROR in registering coupler land') endif - ! do not receive the mesh anymore, read it from file, then pair it with mlnid, component land PC mesh - ! similar to rof mosart mesh - ! do not cull in case of data land, like all other data models - ! for regular land model, cull, because the lnd component culls too - if (lnd_prognostic) then - ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;REPARTITION'//C_NULL_CHAR - else - ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;REPARTITION;NO_CULLING'//C_NULL_CHAR - endif - call seq_infodata_GetData(infodata,lnd_domain=lnd_domain) - outfile = trim(lnd_domain)//C_NULL_CHAR - nghlay = 0 ! no ghost layers - if (seq_comm_iamroot(CPLID) ) then - write(logunit, *) "loading land domain file from file: ", trim(lnd_domain), & - " with options: ", trim(ropts) + endif + + if (mb_scm_land) then + ! change : send the point cloud land mesh (1 point usually) + ! when mb_scm_land + if (MPI_COMM_NULL /= mpicom_old ) then ! it means we are on the component pes (land) + ! send mesh to coupler then + ierr = iMOAB_SendMesh(mlnid, mpicom_join, mpigrp_cplid, id_join, partMethod) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in sending mesh from lnd comp ' + call shr_sys_abort(subname//' ERROR in sending mesh from lnd comp') + endif endif - ierr = iMOAB_LoadMesh(mblxid, outfile, ropts, nghlay) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in reading land coupler mesh from ', trim(lnd_domain) - call shr_sys_abort(subname//' ERROR in reading land coupler mesh') + if (MPI_COMM_NULL /= mpicom_new ) then ! we are on the coupler pes + ierr = iMOAB_ReceiveMesh(mblxid, mpicom_join, mpigrp_old, id_old) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in receiving mesh from lnd comp ' + call shr_sys_abort(subname//' ERROR in receiving mesh from lnd comp') + endif endif + if (MPI_COMM_NULL /= mpicom_old) then ! we are on component lnd pes again, release buffers + context_id = id_join + ierr = iMOAB_FreeSenderBuffers(mlnid, context_id) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in freeing buffers ' + call shr_sys_abort(subname//' ERROR in freeing buffers ') + endif + endif + else + if (MPI_COMM_NULL /= mpicom_new ) then ! we are on the coupler pes + ! do not receive the mesh anymore, read it from file, then pair it with mlnid, component land PC mesh + ! similar to rof mosart mesh + ! do not cull in case of data land, like all other data models + ! for regular land model, cull, because the lnd component culls too + if (lnd_prognostic) then + ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;REPARTITION'//C_NULL_CHAR + else + ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;REPARTITION;NO_CULLING'//C_NULL_CHAR + endif + call seq_infodata_GetData(infodata,lnd_domain=lnd_domain) + outfile = trim(lnd_domain)//C_NULL_CHAR + nghlay = 0 ! no ghost layers + if (seq_comm_iamroot(CPLID) ) then + write(logunit, *) "loading land domain file from file: ", trim(lnd_domain), & + " with options: ", trim(ropts) + endif + ierr = iMOAB_LoadMesh(mblxid, outfile, ropts, nghlay) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in reading land coupler mesh from ', trim(lnd_domain) + call shr_sys_abort(subname//' ERROR in reading land coupler mesh') + endif + endif + endif + if (MPI_COMM_NULL /= mpicom_new ) then ! we are on the coupler pes ! need to add global id tag to the app, it will be used in restart tagtype = 0 ! dense, integer numco = 1 @@ -1645,22 +1687,24 @@ subroutine cplcomp_moab_Init(infodata,comp) endif ! end of coupler pes - ! we are now on joint pes, compute comm graph between lnd and coupler model - typeA = 2 ! point cloud on component PEs, land - typeB = 3 ! full mesh on coupler pes, we just read it + if (mlnid >= 0) then ierr = iMOAB_GetMeshInfo ( mlnid, nvert, nvise, nbl, nsurf, nvisBC ) - comp%mbApCCid = mlnid ! phys atm + comp%mbApCCid = mlnid ! land comp%mbGridType = typeA - 2 ! 0 or 1, pc or cells comp%mblsize = nvert(1) ! vertices endif - ierr = iMOAB_ComputeCommGraph( mlnid, mblxid, mpicom_join, mpigrp_old, mpigrp_cplid, & - typeA, typeB, id_old, id_join) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for lnd model ' - call shr_sys_abort(subname//' ERROR in computing comm graph for lnd model ') + if ( .not. mb_scm_land ) then + ! we are now on joint pes, compute comm graph between lnd and coupler model + typeA = 2 ! point cloud on component PEs, land + typeB = 3 ! full mesh on coupler pes, we just read it + ierr = iMOAB_ComputeCommGraph( mlnid, mblxid, mpicom_join, mpigrp_old, mpigrp_cplid, & + typeA, typeB, id_old, id_join) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for lnd model ' + call shr_sys_abort(subname//' ERROR in computing comm graph for lnd model ') + endif endif - tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR call component_exch_moab(comp, mlnid, mblxid, 0, tagname, context_exch='doml') diff --git a/driver-moab/main/prep_lnd_mod.F90 b/driver-moab/main/prep_lnd_mod.F90 index bd737a1fccf5..4120b87b5bf5 100644 --- a/driver-moab/main/prep_lnd_mod.F90 +++ b/driver-moab/main/prep_lnd_mod.F90 @@ -15,6 +15,7 @@ module prep_lnd_mod use seq_comm_mct, only: mphaid ! iMOAB id for phys atm on atm pes use seq_comm_mct, only: mhpgid ! iMOAB id for atm pgx grid, on atm pes; created with se and gll grids use seq_comm_mct, only: mblxid ! iMOAB id for land migrated mesh to coupler pes + use seq_comm_mct, only: mb_scm_land ! flag that identifies PC for land; use seq_comm_mct, only: mbrxid ! iMOAB id of moab rof on coupler pes (FV now) use seq_comm_mct, only: mbintxal ! iMOAB id for intx mesh between atm and lnd use seq_comm_mct, only: mbintxrl ! iMOAB id for intx mesh between river and land @@ -594,7 +595,11 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln else type1 = 2 ! in this case, atm is just PC endif - type2 = 3; ! FV mesh on coupler land + if (mb_scm_land) then + type2 = 2 ! point cloud for land too, on coupler side; just one point, actually + else + type2 = 3; ! FV mesh on coupler land + endif ierr = iMOAB_ComputeCommGraph( mbaxid, mblxid, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & atm(1)%cplcompid, lnd(1)%cplcompid) if (ierr .ne. 0) then diff --git a/driver-moab/shr/seq_comm_mct.F90 b/driver-moab/shr/seq_comm_mct.F90 index d48aa3e20f76..86e9e4db48c3 100644 --- a/driver-moab/shr/seq_comm_mct.F90 +++ b/driver-moab/shr/seq_comm_mct.F90 @@ -227,6 +227,7 @@ module seq_comm_mct integer, public :: mbintxao ! iMOAB id for intersection mesh between ocean and atmosphere integer, public :: mbintxoa ! iMOAB id for intersection mesh between atmosphere and ocean integer, public :: mblxid ! iMOAB id for land mesh migrated to coupler pes + logical, public :: mb_scm_land = .false. ! land will be migrated if this is true, for scm case; usually one point only !!#ifdef MOABDEBUG integer, public :: mblx2id ! iMOAB id for land mesh instanced from MCT on coupler pes integer, public :: mbox2id ! iMOAB id for ocn mesh instanced from MCT on coupler pes From 25380703143c390589515cfbf601d1333a8a6f6c Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Wed, 8 Oct 2025 10:00:25 -0500 Subject: [PATCH 16/29] fix for conus case problem appeared for conus, debug mode case: SMS_D_Ln5.conusx4v1pg2_r05_IcoswISC30E3r5.F2010.chrysalis_gnu euler formula assumes there is only one connected region, when number of edges is computed HSFC does not guarantee that, especially for RRM meshes select arbitrarily 3 as more connected regions could be possible we could actually allocate much more, nedges_c to be number of quads * 4 Also, when skip match needs to be checked first, before mapid triggered in debug mode only --- components/homme/src/share/semoab_mod.F90 | 1 + driver-moab/main/seq_map_mod.F90 | 10 +++++----- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/components/homme/src/share/semoab_mod.F90 b/components/homme/src/share/semoab_mod.F90 index 4461c2fb39a7..0c5725edd2b2 100644 --- a/components/homme/src/share/semoab_mod.F90 +++ b/components/homme/src/share/semoab_mod.F90 @@ -530,6 +530,7 @@ subroutine create_moab_pg_mesh(par, elem, fv_nphys) ! first count the number of edges in the coarse mesh; ! use euler: v-m+f = 2 => m = v + f - 2 nedges_c = nverts_c + nelemd - 1 ! could be more, if unconnected regions ? + nedges_c = nedges_c + 3 ! assume more than one connected region if ( nedges_c < 0 ) nedges_c = 0 ! it cannot be negative internal_edges = 0 boundary_edges = 0 diff --git a/driver-moab/main/seq_map_mod.F90 b/driver-moab/main/seq_map_mod.F90 index df2489b446cd..845e5269af0b 100644 --- a/driver-moab/main/seq_map_mod.F90 +++ b/driver-moab/main/seq_map_mod.F90 @@ -112,7 +112,7 @@ subroutine seq_map_init_rcfile( mapper, comp_s, comp_d, & if (mct_gsmap_Identical(gsmap_s,gsmap_d)) then if(.not.skip_match) call seq_map_mapmatch(mapid,gsmap_s=gsmap_s,gsmap_d=gsmap_d,strategy="copy") - if (mapid > 0 .and. .not. skip_match) then + if (.not. skip_match .and. mapid > 0 ) then call seq_map_mappoint(mapid,mapper) else if(skip_match .and. seq_comm_iamroot(CPLID)) then @@ -128,7 +128,7 @@ subroutine seq_map_init_rcfile( mapper, comp_s, comp_d, & elseif (samegrid) then if(.not.skip_match) call seq_map_mapmatch(mapid,gsmap_s=gsmap_s,gsmap_d=gsmap_d,strategy="rearrange") - if (mapid > 0 .and. .not. skip_match) then + if (.not. skip_match .and. mapid > 0 ) then call seq_map_mappoint(mapid,mapper) else if(skip_match .and. seq_comm_iamroot(CPLID)) then @@ -151,7 +151,7 @@ subroutine seq_map_init_rcfile( mapper, comp_s, comp_d, & if(.not.skip_match) call seq_map_mapmatch(mapid,gsMap_s=gsMap_s,gsMap_d=gsMap_d,mapfile=mapfile,strategy=maptype) - if (mapid > 0 .and. .not. skip_match) then + if (.not. skip_match .and. mapid > 0) then call seq_map_mappoint(mapid,mapper) else if(skip_match .and. seq_comm_iamroot(CPLID)) then @@ -302,7 +302,7 @@ subroutine seq_map_init_rearrolap(mapper, comp_s, comp_d, string, no_match) if (mct_gsmap_Identical(gsmap_s,gsmap_d)) then if(.not.skip_match) call seq_map_mapmatch(mapid,gsmap_s=gsmap_s,gsmap_d=gsmap_d,strategy="copy") - if (mapid > 0 .and. .not.skip_match) then + if (.not.skip_match .and. mapid > 0) then call seq_map_mappoint(mapid,mapper) else if(skip_match .and. seq_comm_iamroot(CPLID)) then @@ -318,7 +318,7 @@ subroutine seq_map_init_rearrolap(mapper, comp_s, comp_d, string, no_match) else if(.not.skip_match) call seq_map_mapmatch(mapid,gsmap_s=gsmap_s,gsmap_d=gsmap_d,strategy="rearrange") - if (mapid > 0 .and. .not.skip_match) then + if (.not.skip_match .and. mapid > 0) then call seq_map_mappoint(mapid,mapper) else if(skip_match .and. seq_comm_iamroot(CPLID)) then From 85ed4e51ccfbe59211760ad9675355fd9436a995 Mon Sep 17 00:00:00 2001 From: Iulian Grindeanu Date: Wed, 8 Oct 2025 15:47:36 -0500 Subject: [PATCH 17/29] retrieve moab matrices only when it makes sense matrices o2x_am, used for atm merge, should not be retrieved if ocean is not active similar for ice, and xao fluxes --- driver-moab/main/prep_atm_mod.F90 | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/driver-moab/main/prep_atm_mod.F90 b/driver-moab/main/prep_atm_mod.F90 index 295efb8aeb87..a5492c58f1dd 100644 --- a/driver-moab/main/prep_atm_mod.F90 +++ b/driver-moab/main/prep_atm_mod.F90 @@ -1355,33 +1355,41 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) call shr_sys_abort(subname//' error in getting fractions_am from atm instance ') endif + if (mboxid > 0) then ! retrieve projection only when ox is active ? tagname = trim(seq_flds_o2x_fields)//C_NULL_CHAR arrsize = noflds * lsize ! allocate (o2x_am (lsize, noflds)) ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, o2x_am) if (ierr .ne. 0) then call shr_sys_abort(subname//' error in getting o2x_am array ') endif - + endif + + if (mbixid > 0) then tagname = trim(seq_flds_i2x_fields)//C_NULL_CHAR arrsize = niflds * lsize ! allocate (i2x_am (lsize, niflds)) ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, i2x_am) if (ierr .ne. 0) then call shr_sys_abort(subname//' error in getting i2x_am array ') endif + endif + if (mblxid > 0) then ! tagname = trim(seq_flds_l2x_fields)//C_NULL_CHAR arrsize = nlflds * lsize ! allocate (l2x_am (lsize, nlflds)) ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, l2x_am) if (ierr .ne. 0) then call shr_sys_abort(subname//' error in getting l2x_am array ') endif + endif + if (mboxid > 0) then tagname = trim(seq_flds_xao_fields)//C_NULL_CHAR arrsize = nxflds * lsize ! allocate (xao_am (lsize, nxflds)) ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, xao_am) if (ierr .ne. 0) then call shr_sys_abort(subname//' error in getting xao_om array ') endif + endif do n = 1,lsize From 1b130da6f4387aa3bf4307729bd2695947016a10 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Tue, 14 Oct 2025 14:59:19 -0500 Subject: [PATCH 18/29] Isolate atm ops in prep_rof_mrg_moab with rof_heat Ops involving a2x_rm in prep_rof_mrg_moab should only happen when rof_heat is true. --- driver-moab/main/prep_rof_mod.F90 | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/driver-moab/main/prep_rof_mod.F90 b/driver-moab/main/prep_rof_mod.F90 index 29d2be67bd82..ee68e44f99f1 100644 --- a/driver-moab/main/prep_rof_mod.F90 +++ b/driver-moab/main/prep_rof_mod.F90 @@ -1703,15 +1703,19 @@ subroutine prep_rof_mrg_moab (infodata, cime_model) ! (name of the fields) ! need these always, not only the first time l2x_r => l2r_rx(1) - a2x_r => a2r_rx(1) x2r_r => component_get_x2c_cx(rof(1)) nflds = mct_aVect_nRattr(x2r_r) ! these are saved after first time - naflds = mct_aVect_nRattr(a2x_r) + + if (rof_heat) then + a2x_r => a2r_rx(1) + naflds = mct_aVect_nRattr(a2x_r) + allocate(a2x_rm (lsize, naflds)) + endif + nlflds = mct_aVect_nRattr(l2x_r) allocate(x2r_rm (lsize, nflds)) - allocate(a2x_rm (lsize, naflds)) allocate(l2x_rm (lsize, nlflds)) ! allocate fractions too ! use the fraclist fraclist_r = 'lfrac:lfrin:rfrac' @@ -1867,11 +1871,13 @@ subroutine prep_rof_mrg_moab (infodata, cime_model) endif ! a2x_rm (lsize, naflds)) - tagname = trim(seq_flds_a2x_fields_to_rof)//C_NULL_CHAR - arrsize = naflds * lsize ! allocate (a2x_rm (lsize, naflds)) - ierr = iMOAB_GetDoubleTagStorage ( mbrxid, tagname, arrsize , ent_type, a2x_rm) - if (ierr .ne. 0) then - call shr_sys_abort(subname//' error in getting a2x_rm array ') + if (rof_heat) then + tagname = trim(seq_flds_a2x_fields_to_rof)//C_NULL_CHAR + arrsize = naflds * lsize ! allocate (a2x_rm (lsize, naflds)) + ierr = iMOAB_GetDoubleTagStorage ( mbrxid, tagname, arrsize , ent_type, a2x_rm) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' error in getting a2x_rm array ') + endif endif ! l2x_rm tagname = trim(seq_flds_l2x_fluxes_to_rof)//C_NULL_CHAR From 172c10eb6984b5005a328a1ceafcabb83bfae228 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Wed, 29 Oct 2025 21:46:05 -0500 Subject: [PATCH 19/29] add areafact_samegrid to component_init_areacor_moab Add areafact_samegrid to component_init_areacor_moab and when its true, set the area correction factors to 1. As done in driver-mct. --- driver-moab/main/cime_comp_mod.F90 | 10 +++--- driver-moab/main/component_mod.F90 | 56 ++++++++++++++++-------------- 2 files changed, 35 insertions(+), 31 deletions(-) diff --git a/driver-moab/main/cime_comp_mod.F90 b/driver-moab/main/cime_comp_mod.F90 index 6e14dfac203a..c23659131297 100644 --- a/driver-moab/main/cime_comp_mod.F90 +++ b/driver-moab/main/cime_comp_mod.F90 @@ -2178,28 +2178,28 @@ subroutine cime_init() call mpi_barrier(mpicom_GLOID,ierr) if (atm_present) call component_init_areacor(atm, areafact_samegrid, seq_flds_a2x_fluxes) ! send initial data to coupler - if (atm_present) call component_init_areacor_moab(atm, mphaid, mbaxid, seq_flds_a2x_fluxes, seq_flds_a2x_fields) + if (atm_present) call component_init_areacor_moab(atm, areafact_samegrid, mphaid, mbaxid, seq_flds_a2x_fluxes, seq_flds_a2x_fields) ! component_exch_moab(atm(1), mphaid, mbaxid, 0, seq_flds_a2x_fields) call mpi_barrier(mpicom_GLOID,ierr) if (lnd_present) call component_init_areacor(lnd, areafact_samegrid, seq_flds_l2x_fluxes) ! MOABTODO : lnd is vertex or cell ? - if (lnd_present) call component_init_areacor_moab(lnd, mlnid, mblxid, seq_flds_l2x_fluxes, seq_flds_l2x_fields) + if (lnd_present) call component_init_areacor_moab(lnd, areafact_samegrid, mlnid, mblxid, seq_flds_l2x_fluxes, seq_flds_l2x_fields) !component_exch_moab(lnd(1), mlnid, mblxid, 0, seq_flds_l2x_fields) call mpi_barrier(mpicom_GLOID,ierr) if (rof_present) call component_init_areacor(rof, areafact_samegrid, seq_flds_r2x_fluxes) - if (rof_present) call component_init_areacor_moab(rof, mrofid, mbrxid, seq_flds_r2x_fluxes, seq_flds_r2x_fields) + if (rof_present) call component_init_areacor_moab(rof, areafact_samegrid, mrofid, mbrxid, seq_flds_r2x_fluxes, seq_flds_r2x_fields) !component_exch_moab(rof(1), mrofid, mbrxid, 0, seq_flds_r2x_fields) call mpi_barrier(mpicom_GLOID,ierr) if (ocn_present) call component_init_areacor(ocn, areafact_samegrid, seq_flds_o2x_fluxes) - if (ocn_present) call component_init_areacor_moab(ocn, mpoid, mboxid, seq_flds_o2x_fluxes, seq_flds_o2x_fields) + if (ocn_present) call component_init_areacor_moab(ocn, areafact_samegrid, mpoid, mboxid, seq_flds_o2x_fluxes, seq_flds_o2x_fields) ! component_exch_moab(ocn(1), mpoid, mboxid, 0, seq_flds_o2x_fields) call mpi_barrier(mpicom_GLOID,ierr) if (ice_present) call component_init_areacor(ice, areafact_samegrid, seq_flds_i2x_fluxes) - if (ice_present) call component_init_areacor_moab(ice, mpsiid, mbixid, seq_flds_i2x_fluxes, seq_flds_i2x_fields) + if (ice_present) call component_init_areacor_moab(ice, areafact_samegrid, mpsiid, mbixid, seq_flds_i2x_fluxes, seq_flds_i2x_fields) !component_exch_moab(ice(1), mpsiid, mbixid, 0, seq_flds_i2x_fields) call mpi_barrier(mpicom_GLOID,ierr) diff --git a/driver-moab/main/component_mod.F90 b/driver-moab/main/component_mod.F90 index 68c378e92bbb..0e81834d9691 100644 --- a/driver-moab/main/component_mod.F90 +++ b/driver-moab/main/component_mod.F90 @@ -740,7 +740,7 @@ subroutine component_init_areacor(comp, samegrid, seq_flds_c2x_fluxes) end subroutine component_init_areacor -subroutine component_init_areacor_moab (comp, mbccid, mbcxid, seq_flds_c2x_fluxes, seq_flds_c2x_fields) +subroutine component_init_areacor_moab (comp, samegrid, mbccid, mbcxid, seq_flds_c2x_fluxes, seq_flds_c2x_fields) !--------------------------------------------------------------- ! COMPONENT PES and CPL/COMPONENT (for exchange only) ! @@ -754,6 +754,7 @@ subroutine component_init_areacor_moab (comp, mbccid, mbcxid, seq_flds_c2x_fluxe ! ! Arguments type(component_type) , intent(inout) :: comp(:) + logical , intent(in) :: samegrid integer , intent(in) :: mbccid ! comp side integer , intent(in) :: mbcxid ! coupler side ! point cloud or FV type, to use vertices or cells for setting/getting the area tags and corrections @@ -790,32 +791,35 @@ subroutine component_init_areacor_moab (comp, mbccid, mbcxid, seq_flds_c2x_fluxe lsize = comp(1)%mblsize allocate(areas (lsize, 3)) ! lsize is along grid; read mask too allocate(factors (lsize, 2)) - factors = 1.0 ! initialize with 1.0 all factors; then maybe correct them - ! get areas - tagname='area:aream:mask'//C_NULL_CHAR - arrsize = 3 * lsize - ierr = iMOAB_GetDoubleTagStorage ( mbccid, tagname, arrsize , comp(1)%mbGridType, areas(1,1) ) - if (ierr .ne. 0) then - call shr_sys_abort(subname//' cannot get areas ') - endif - ! now compute the factors - do i=1,lsize - rmask = areas(i,3) - rarea = areas(i, 1) - raream = areas(i, 2) - if ( abs(rmask) >= 1.0e-06) then - if (rarea * raream /= 0.0_R8) then - factors(i,1) = rarea/raream - factors(i,2)= 1.0_R8/factors(i,1) - else - write(logunit,*) trim(subname),' ERROR area,aream= ', & - rarea,raream,' in ',i,lsize - call shr_sys_flush(logunit) - call shr_sys_abort() - endif - endif - enddo + factors = 1.0_r8 ! initialize with 1.0 all factors; + if (.not.samegrid) then ! only correct if not monogrid (SCM case) + ! get areas + tagname='area:aream:mask'//C_NULL_CHAR + arrsize = 3 * lsize + ierr = iMOAB_GetDoubleTagStorage ( mbccid, tagname, arrsize , comp(1)%mbGridType, areas(1,1) ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' cannot get areas ') + endif + ! now compute the factors + do i=1,lsize + rmask = areas(i,3) + + rarea = areas(i, 1) + raream = areas(i, 2) + if ( abs(rmask) >= 1.0e-06) then + if (rarea * raream /= 0.0_R8) then + factors(i,1) = rarea/raream + factors(i,2)= 1.0_R8/factors(i,1) + else + write(logunit,*) trim(subname),' ERROR area,aream= ', & + rarea,raream,' in ',i,lsize + call shr_sys_flush(logunit) + call shr_sys_abort() + endif + endif + enddo + endif ! set factors as tags ! define the tags mdl2drv and drv2mdl on component sides, and compute them based on area and aream tagname = 'mdl2drv:drv2mdl'//C_NULL_CHAR From 62024f30c83e5564fdc26a6278b663de1e5b5f8c Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Sun, 2 Nov 2025 21:55:20 -0600 Subject: [PATCH 20/29] Add more handling of single column mode for land Add more handling of single column mode for land. --- driver-moab/main/cplcomp_exchange_mod.F90 | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index d13630b9ce13..ba8ee2be71e7 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -998,8 +998,8 @@ subroutine copy_aream_from_area(mbappid) ! copy aream from area if (mbappid >= 0) then ! coupler procs ierr = iMOAB_GetMeshInfo ( mbappid, nvert, nvise, nbl, nsurf, nvisBC ) - if ( (mbappid .eq. mbaxid) .and. (.not. atm_pg_active)) then - ! this is the spectral case, for atm only + if (.not.atm_pg_active) then + ! this is the spectral monogrid case arrSize = nvert(1) ! cells ent_type = 0 ! vertices else @@ -1662,11 +1662,16 @@ subroutine cplcomp_moab_Init(infodata,comp) nfields=mct_list_nitem (temp_list) if (nfields > 0) then ierr = iMOAB_GetMeshInfo ( mblxid, nvert, nvise, nbl, nsurf, nvisBC ) - arrsize = nvise(1)*nfields + if (mb_scm_land) then + arrsize = nvert(1)*nfields + ent_type = 0 ! cell + else + arrsize = nvise(1)*nfields + ent_type = 1 ! cell + endif allocate(tagValues(arrsize)) tagname = trim(newlist)//C_NULL_CHAR tagValues = 0.0_r8 - ent_type = 1 ! cells ierr = iMOAB_SetDoubleTagStorage ( mblxid, tagname, arrsize , ent_type, tagValues) if (ierr .ne. 0) then write(logunit,*) subname,' error in zeroing Flrr tags on land', ierr @@ -1691,6 +1696,8 @@ subroutine cplcomp_moab_Init(infodata,comp) if (mlnid >= 0) then ierr = iMOAB_GetMeshInfo ( mlnid, nvert, nvise, nbl, nsurf, nvisBC ) comp%mbApCCid = mlnid ! land + ! MOAB TODO: check this logic. Seems to assume typeA might be + ! 3 or 2 but only typeB has that possibility. comp%mbGridType = typeA - 2 ! 0 or 1, pc or cells comp%mblsize = nvert(1) ! vertices endif @@ -1707,6 +1714,9 @@ subroutine cplcomp_moab_Init(infodata,comp) endif tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR call component_exch_moab(comp, mlnid, mblxid, 0, tagname, context_exch='doml') + if (mblxid > 0) then ! on coupler pes only + call copy_aream_from_area(mblxid) + endif #ifdef MOABDEBUG outfile = 'recMeshLand.h5m'//C_NULL_CHAR From 4055dfed124fd222c6745f057df1bf358f2ca372 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Sun, 2 Nov 2025 21:58:48 -0600 Subject: [PATCH 21/29] Cover scm for land in mb routines Also check for single column model when land mb id is passed in and set ent_type appropriately. --- driver-moab/shr/shr_moab_mod.F90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/driver-moab/shr/shr_moab_mod.F90 b/driver-moab/shr/shr_moab_mod.F90 index b51a0def95bd..9fde5f4d5c71 100644 --- a/driver-moab/shr/shr_moab_mod.F90 +++ b/driver-moab/shr/shr_moab_mod.F90 @@ -7,7 +7,7 @@ module shr_moab_mod use shr_sys_mod, only: shr_sys_abort use seq_comm_mct, only: logunit - use seq_comm_mct, only: atm_pg_active, mbaxid + use seq_comm_mct, only: atm_pg_active, mbaxid,mb_scm_land,mblxid use shr_kind_mod, only: CXX => shr_kind_CXX use shr_kind_mod , only: R8 => SHR_KIND_R8 use iso_c_binding @@ -55,7 +55,8 @@ integer function mbGetnCells(moabid) endif ! only case on coupler side that we actually want number of vertices is when we use spectral atm - if ( .not. atm_pg_active .and. mbaxid .eq. moabid) then + if ((.not. atm_pg_active .and. (mbaxid .eq. moabid)) .or. & + (mb_scm_land .and. (mblxid .eq. moabid))) then mbGetnCells = nvert(1) else mbGetnCells = nvise(1) @@ -90,8 +91,8 @@ subroutine mbGetCellTagVals(mbid, intag,inarray,nMax) character(*), parameter :: subname = '(mbGetCellTagVals) ' !----------------------------------------------------------------------- ! - - if ( .not. atm_pg_active .and. mbaxid .eq. mbid) then + if ((.not. atm_pg_active .and. (mbaxid .eq. mbid)) .or. & + (mb_scm_land .and. (mblxid .eq. mbid))) then ent_type = 0 else ent_type = 1 @@ -134,7 +135,8 @@ subroutine mbSetCellTagVals(mbid, intag,inarray,nMax) !----------------------------------------------------------------------- ! - if ( .not. atm_pg_active .and. mbaxid .eq. mbid) then + if ((.not. atm_pg_active .and. (mbaxid .eq. mbid)) .or. & + (mb_scm_land .and. (mblxid .eq. mbid))) then ent_type = 0 else ent_type = 1 From 666ce4510f87cf3fd9d322934f744a86b8bdef62 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Sun, 2 Nov 2025 22:00:54 -0600 Subject: [PATCH 22/29] Zero Si_snowh for single_column In single_column case, zero out Si_snowh to match eam history with driver-mct --- driver-moab/main/prep_atm_mod.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/driver-moab/main/prep_atm_mod.F90 b/driver-moab/main/prep_atm_mod.F90 index a5492c58f1dd..52bd01820347 100644 --- a/driver-moab/main/prep_atm_mod.F90 +++ b/driver-moab/main/prep_atm_mod.F90 @@ -1046,6 +1046,7 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) integer :: n,ka,ki,kl,ko,kx,kof,kif,klf,klf_st,i,i1,o1 integer, save :: lsize integer, save :: index_x2a_Sf_lfrac, index_x2a_Sf_ifrac, index_x2a_Sf_ofrac + integer, save :: index_x2a_Si_snowh character(CL) :: atm_gnam, lnd_gnam character(CL) :: fracstr, fracstr_st @@ -1076,6 +1077,7 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) ! for moab info character(CXX) ::tagname, mct_field integer :: ent_type, ierr, arrsize + logical :: single_column #ifdef MOABDEBUG character*32 :: outfile, wopts, lnum #endif @@ -1090,6 +1092,7 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) !----------------------------------------------------------------------- ! call seq_comm_getdata(CPLID, iamroot=iamroot) + call seq_infodata_getdata(infodata,single_column=single_column) if (first_time) then @@ -1122,6 +1125,7 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) index_x2a_Sf_lfrac = mct_aVect_indexRA(x2a_a,'Sf_lfrac') index_x2a_Sf_ifrac = mct_aVect_indexRA(x2a_a,'Sf_ifrac') index_x2a_Sf_ofrac = mct_aVect_indexRA(x2a_a,'Sf_ofrac') + index_x2a_Si_snowh = mct_aVect_indexRA(x2a_a,'Si_snowh') !ngflds = mct_aVect_nRattr(g2x_o) allocate(fractions_am(lsize,5)) ! there are 5 fractions 'afrac:ifrac:ofrac:lfrac:lfrin' @@ -1396,6 +1400,8 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) x2a_am(n, index_x2a_Sf_lfrac) = fractions_am(n, klf) ! x2a_a%rAttr(index_x2a_Sf_lfrac,n) = fractions_a%Rattr(klf,n) x2a_am(n, index_x2a_Sf_ifrac) = fractions_am(n, kif) ! x2a_a%rAttr(index_x2a_Sf_ifrac,n) = fractions_a%Rattr(kif,n) x2a_am(n, index_x2a_Sf_ofrac) = fractions_am(n, kof) ! x2a_a%rAttr(index_x2a_Sf_ofrac,n) = fractions_a%Rattr(kof,n) + ! zero out Si_snowh for single column to match MCT + if(single_column) x2a_am(n, index_x2a_Si_snowh) = 0.0_r8 end do !--- document fraction operations --- From 4630bc6d8bbdbe541b7e967528fa25a837efd5ff Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Mon, 3 Nov 2025 10:56:19 -0600 Subject: [PATCH 23/29] Set ent_type correctly for scm land Set ent_type correctly for scm land --- driver-moab/main/seq_io_mod.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/driver-moab/main/seq_io_mod.F90 b/driver-moab/main/seq_io_mod.F90 index a50726a02eb3..c03d525ba18e 100644 --- a/driver-moab/main/seq_io_mod.F90 +++ b/driver-moab/main/seq_io_mod.F90 @@ -33,7 +33,7 @@ module seq_io_mod use shr_sys_mod, only: shr_sys_abort use seq_comm_mct, only: logunit, CPLID, seq_comm_setptrs use seq_comm_mct, only: seq_comm_namelen, seq_comm_name - use seq_comm_mct, only: mbaxid, atm_pg_active + use seq_comm_mct, only: mbaxid, atm_pg_active,mblxid,mb_scm_land use seq_flds_mod, only : seq_flds_lookup use mct_mod ! mct wrappers use pio @@ -1719,7 +1719,8 @@ subroutine seq_io_write_moab_tags(filename, mbxid, dname, tag_list, whead,wdata, ! get the local size ns ierr = iMOAB_GetMeshInfo ( mbxid, nvert, nvise, nbl, nsurf, nvisBC ) - if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then + if ((.not. atm_pg_active .and. (mbaxid .eq. mbxid)) .or. & + (mb_scm_land .and. (mblxid .eq. mbxid))) then ent_type = 0 ns = nvert(1) ! local vertices lnx = ngv ! number of global vertices @@ -2729,7 +2730,8 @@ subroutine seq_io_read_moab_tags(filename, mbxid, dname, tag_list, matrix, nx) ierr = iMOAB_GetGlobalInfo( mbxid, ngv, ng) lny = 1 ! do we need 2 var, or just 1 ierr = iMOAB_GetMeshInfo ( mbxid, nvert, nvise, nbl, nsurf, nvisBC ) - if (mbxid .eq. mbaxid .and. .not. atm_pg_active) then + if ((.not. atm_pg_active .and. (mbaxid .eq. mbxid)) .or. & + (mb_scm_land .and. (mblxid .eq. mbxid))) then ent_type = 0 ns = nvert(1) ! local vertices lnx = ngv ! number of global vertices From c07c3ae521be4ead8689282088946822f445a9a3 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Mon, 3 Nov 2025 10:57:31 -0600 Subject: [PATCH 24/29] Set ent_type correctly with land is single column Set ent_type correctly with land is in single column mode for section that zeros river variables. --- driver-moab/main/prep_lnd_mod.F90 | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/driver-moab/main/prep_lnd_mod.F90 b/driver-moab/main/prep_lnd_mod.F90 index 4120b87b5bf5..b73080a10f3a 100644 --- a/driver-moab/main/prep_lnd_mod.F90 +++ b/driver-moab/main/prep_lnd_mod.F90 @@ -22,6 +22,7 @@ module prep_lnd_mod use seq_comm_mct, only: mbaxid ! iMOAB id for atm migrated mesh to coupler pes use seq_comm_mct, only: atm_pg_active ! whether the atm uses FV mesh or not ; made true if fv_nphys > 0 + use seq_comm_mct, only: mb_scm_land ! use dimensions_mod, only: np ! for atmosphere use seq_comm_mct, only: seq_comm_getinfo => seq_comm_setptrs use seq_map_type_mod @@ -400,9 +401,15 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln write(logunit,*) subname,' cant get size of land mesh' call shr_sys_abort(subname//' ERROR in getting size of land mesh') endif - ! land is now cell mesh on coupler side - mlsize = nvise(1) - ent_type = 1 ! cell + ! land is usually cell on coupler but could be point + if(mb_scm_land) then + mlsize = nvert(1) + ent_type = 0 ! point cloud + else + mlsize = nvise(1) + ent_type = 1 ! cell + endif + ! set to 0 all fields that are projected from river nrflds = mct_aVect_nRattr(r2x_lx(1)) ! these are the numbers of fields in seq_flds_r2x_fields arrsize = nrflds*mlsize From 2aa691ef1f65c47646fbaefc2315849e3b89d73a Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Wed, 5 Nov 2025 20:10:16 -0600 Subject: [PATCH 25/29] Change how mb_scm_land is set Use single_column and scm_multcol to determine if mb_scm_land should be set. DP case will not have land_nx=land_ny=1 --- driver-moab/main/cplcomp_exchange_mod.F90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index ba8ee2be71e7..ee6e6dda2f7a 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -1066,8 +1066,7 @@ subroutine cplcomp_moab_Init(infodata,comp) character(CXX) :: tagname character(CXX) :: newlist integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) - logical :: rof_present, lnd_prognostic - integer :: land_nx, land_ny ! try to identify if scm case; then land mesh will be migrated, not read + logical :: rof_present, lnd_prognostic, single_column, scm_multcols real(r8), allocatable :: tagValues(:) ! used for setting aream tags for atm domain read case integer :: arrsize ! for the size of tagValues type(mct_list) :: temp_list @@ -1559,9 +1558,9 @@ subroutine cplcomp_moab_Init(infodata,comp) call seq_comm_getinfo(cplid ,mpigrp=mpigrp_cplid) ! receiver group call seq_comm_getinfo(id_old,mpigrp=mpigrp_old) ! component group pes call seq_infodata_GetData(infodata,rof_present=rof_present, lnd_prognostic=lnd_prognostic) - call seq_infodata_GetData(infodata,lnd_nx=land_nx, lnd_ny=land_ny) - - if (land_nx .eq. 1 .and. land_ny .eq.1 ) then + call seq_infodata_GetData(infodata,single_column=single_column, & + scm_multcols=scm_multcols) + if (single_column .or. scm_multcols) then ! turn on mb_scm_land mb_scm_land = .true. ! we identified a scm case for land, we will migrate mesh, not read ! the domain file anymore From 88f92c7533a0684b84503cbecce59e29daa2f8bf Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Fri, 7 Nov 2025 14:36:52 -0600 Subject: [PATCH 26/29] Set mbGridType to 0 for mlnid Set mbGridType to 0 for mlnid. Some code was moved for the scm case which broke the calculation that was here. --- driver-moab/main/cplcomp_exchange_mod.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index ee6e6dda2f7a..e14aadcb2dbd 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -1695,9 +1695,7 @@ subroutine cplcomp_moab_Init(infodata,comp) if (mlnid >= 0) then ierr = iMOAB_GetMeshInfo ( mlnid, nvert, nvise, nbl, nsurf, nvisBC ) comp%mbApCCid = mlnid ! land - ! MOAB TODO: check this logic. Seems to assume typeA might be - ! 3 or 2 but only typeB has that possibility. - comp%mbGridType = typeA - 2 ! 0 or 1, pc or cells + comp%mbGridType = 0 ! 0 or 1, pc or cells comp%mblsize = nvert(1) ! vertices endif if ( .not. mb_scm_land ) then From 0771e5e5ccb8f1cc2f99cbe4c65a0e54be68f4ba Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Sun, 9 Nov 2025 23:19:54 -0600 Subject: [PATCH 27/29] Make sure areas are initialized in areacor Make sure areas are initialized in areacor by moving the iMOAB_GetDouble call that does it outside of the if(.not.samegrid) block. aream is needed later. --- driver-moab/main/component_mod.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/driver-moab/main/component_mod.F90 b/driver-moab/main/component_mod.F90 index 0e81834d9691..18a80bfcc1bb 100644 --- a/driver-moab/main/component_mod.F90 +++ b/driver-moab/main/component_mod.F90 @@ -793,14 +793,14 @@ subroutine component_init_areacor_moab (comp, samegrid, mbccid, mbcxid, seq_flds allocate(factors (lsize, 2)) factors = 1.0_r8 ! initialize with 1.0 all factors; + ! get areas + tagname='area:aream:mask'//C_NULL_CHAR + arrsize = 3 * lsize + ierr = iMOAB_GetDoubleTagStorage ( mbccid, tagname, arrsize , comp(1)%mbGridType, areas(1,1) ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' cannot get areas ') + endif if (.not.samegrid) then ! only correct if not monogrid (SCM case) - ! get areas - tagname='area:aream:mask'//C_NULL_CHAR - arrsize = 3 * lsize - ierr = iMOAB_GetDoubleTagStorage ( mbccid, tagname, arrsize , comp(1)%mbGridType, areas(1,1) ) - if (ierr .ne. 0) then - call shr_sys_abort(subname//' cannot get areas ') - endif ! now compute the factors do i=1,lsize rmask = areas(i,3) @@ -819,7 +819,7 @@ subroutine component_init_areacor_moab (comp, samegrid, mbccid, mbcxid, seq_flds endif endif enddo - endif + endif ! if .not.samegrid ! set factors as tags ! define the tags mdl2drv and drv2mdl on component sides, and compute them based on area and aream tagname = 'mdl2drv:drv2mdl'//C_NULL_CHAR From 17d4f6e85ab1757f1164d919a8dd97fa60d4eb04 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Sun, 9 Nov 2025 23:40:33 -0600 Subject: [PATCH 28/29] In atm merge, make sure arrays for other components are initialized xao_am, l2x_am, i2x_am, and o2x_am are only set using GetDouble if the corresponding model is present. For at least o2x_am, it was still used in a merge even if not present. So 0 those arrays if the corresponding moab app id is not set. This zero's the Faoo_h2otemp field for the moab coupler. --- driver-moab/main/prep_atm_mod.F90 | 62 ++++++++++++++++--------------- 1 file changed, 32 insertions(+), 30 deletions(-) diff --git a/driver-moab/main/prep_atm_mod.F90 b/driver-moab/main/prep_atm_mod.F90 index 52bd01820347..17b6a1105d4e 100644 --- a/driver-moab/main/prep_atm_mod.F90 +++ b/driver-moab/main/prep_atm_mod.F90 @@ -1360,40 +1360,48 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) endif if (mboxid > 0) then ! retrieve projection only when ox is active ? - tagname = trim(seq_flds_o2x_fields)//C_NULL_CHAR - arrsize = noflds * lsize ! allocate (o2x_am (lsize, noflds)) - ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, o2x_am) - if (ierr .ne. 0) then - call shr_sys_abort(subname//' error in getting o2x_am array ') - endif + tagname = trim(seq_flds_o2x_fields)//C_NULL_CHAR + arrsize = noflds * lsize ! allocate (o2x_am (lsize, noflds)) + ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, o2x_am) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' error in getting o2x_am array ') + endif + else + o2x_am(:,:)=0.0_r8 endif if (mbixid > 0) then - tagname = trim(seq_flds_i2x_fields)//C_NULL_CHAR - arrsize = niflds * lsize ! allocate (i2x_am (lsize, niflds)) - ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, i2x_am) - if (ierr .ne. 0) then - call shr_sys_abort(subname//' error in getting i2x_am array ') + tagname = trim(seq_flds_i2x_fields)//C_NULL_CHAR + arrsize = niflds * lsize ! allocate (i2x_am (lsize, niflds)) + ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, i2x_am) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' error in getting i2x_am array ') + endif + else + i2x_am(:,:)=0.0_r8 endif - endif if (mblxid > 0) then ! - tagname = trim(seq_flds_l2x_fields)//C_NULL_CHAR - arrsize = nlflds * lsize ! allocate (l2x_am (lsize, nlflds)) - ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, l2x_am) - if (ierr .ne. 0) then - call shr_sys_abort(subname//' error in getting l2x_am array ') + tagname = trim(seq_flds_l2x_fields)//C_NULL_CHAR + arrsize = nlflds * lsize ! allocate (l2x_am (lsize, nlflds)) + ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, l2x_am) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' error in getting l2x_am array ') + endif + else + l2x_am(:,:)=0.0_r8 endif - endif if (mboxid > 0) then - tagname = trim(seq_flds_xao_fields)//C_NULL_CHAR - arrsize = nxflds * lsize ! allocate (xao_am (lsize, nxflds)) - ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, xao_am) - if (ierr .ne. 0) then - call shr_sys_abort(subname//' error in getting xao_om array ') + tagname = trim(seq_flds_xao_fields)//C_NULL_CHAR + arrsize = nxflds * lsize ! allocate (xao_am (lsize, nxflds)) + ierr = iMOAB_GetDoubleTagStorage ( mbaxid, tagname, arrsize , ent_type, xao_am) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' error in getting xao_om array ') + endif + else + xao_am(:,:)=0.0_r8 endif - endif do n = 1,lsize @@ -1456,12 +1464,6 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) enddo endif ! first time - ! we need to do something equivalent, to copy in a2x_am the tags from those shared indices - ! call mct_aVect_copy(aVin=l2x_a, aVout=x2a_a, vector=mct_usevector, sharedIndices=l2x_SharedIndices) - !call mct_aVect_copy(aVin=o2x_a, aVout=x2a_a, vector=mct_usevector, sharedIndices=o2x_SharedIndices) - !call mct_aVect_copy(aVin=i2x_a, aVout=x2a_a, vector=mct_usevector, sharedIndices=i2x_SharedIndices) - !call mct_aVect_copy(aVin=xao_a, aVout=x2a_a, vector=mct_usevector, sharedIndices=xao_SharedIndices) - ! If flux to atm is coming only from the ocean (based on field being in o2x_a) - ! -- then scale by both ocean and ice fraction ! If flux to atm is coming only from the land or ice or coupler From 6e346b10982583b42ddd11406e1fc66496ffd721 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Tue, 11 Nov 2025 22:58:35 -0600 Subject: [PATCH 29/29] Make sure land fluxes are zero when no land model 'Fall_' fluxes that are not merged have to be set to 0 in MOAB when there's no land model. Otherwise they will keep the default value. Found while testing the aquaplanet case and dust fluxes were -1e+10 --- driver-moab/main/prep_atm_mod.F90 | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/driver-moab/main/prep_atm_mod.F90 b/driver-moab/main/prep_atm_mod.F90 index 17b6a1105d4e..7937568f2b12 100644 --- a/driver-moab/main/prep_atm_mod.F90 +++ b/driver-moab/main/prep_atm_mod.F90 @@ -1077,7 +1077,7 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) ! for moab info character(CXX) ::tagname, mct_field integer :: ent_type, ierr, arrsize - logical :: single_column + logical :: single_column,lnd_present #ifdef MOABDEBUG character*32 :: outfile, wopts, lnum #endif @@ -1092,7 +1092,7 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) !----------------------------------------------------------------------- ! call seq_comm_getdata(CPLID, iamroot=iamroot) - call seq_infodata_getdata(infodata,single_column=single_column) + call seq_infodata_getdata(infodata,single_column=single_column,lnd_present=lnd_present) if (first_time) then @@ -1531,6 +1531,13 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) end if end if endif + ! if no land, make sure any 'Fall' fluxes that are not merged are 0. + ! e.g. dust fluxes + if(.not.lnd_present) then + if((lindx(ka)>0) .and. .not.lstate(ka) .and. .not.lmerge(ka)) then + x2a_am(n,ka) = 0.0_r8 + endif + endif if (iindx(ka) > 0 .and. fraci > 0._r8) then if (imerge(ka)) then x2a_am(n, ka) = x2a_am(n, ka) + i2x_am(n, iindx(ka)) * fraci