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..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_meshes + use semoab_mod, only: create_moab_pg_mesh use iMOAB, only : iMOAB_RegisterApplication use iso_c_binding #endif @@ -625,45 +625,45 @@ subroutine phys_grid_init (pgN) call compute_global_area (pg) #ifdef HAVE_MOAB if (pgN > 0) then - appname="HM_COARSE"//C_NULL_CHAR - ATM_ID1 = 120 ! - ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, MHID) - if (ierr > 0 ) & - call abortmp('Error: cannot register moab app') - if(par%masterproc) then - write(iulog,*) " " - write(iulog,*) "register MOAB app:", trim(appname), " MHID=", MHID - write(iulog,*) " " - endif - appname="HM_FINE"//C_NULL_CHAR - ATM_ID1 = 119 ! this number should not conflict with other components IDs; how do we know? - ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, MHFID) - if (ierr > 0 ) & - call abortmp('Error: cannot register moab app for fine mesh') - if(par%masterproc) then - write(iulog,*) " " - write(iulog,*) "register MOAB app:", trim(appname), " MHFID=", MHFID - write(iulog,*) " " - endif - appname="HM_PGX"//C_NULL_CHAR - ATM_ID1 = ATMID(1) ! this number should not conflict with other components IDs; how do we know? - ! - ! in this case, we reuse the main atm id, mhid will not be used for intersection anymore - ! still, need to be careful - ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, mhpgid) - if (ierr > 0 ) & - call abortmp('Error: cannot register moab app for fine mesh') - if(par%masterproc) then - write(iulog,*) " " - write(iulog,*) "register MOAB app:", trim(appname), " MHPGID=", mhpgid - write(iulog,*) " " - endif + appname="HM_COARSE"//C_NULL_CHAR + ATM_ID1 = 120 ! + ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, MHID) + if (ierr > 0 ) & + call abortmp('Error: cannot register moab app') + if(par%masterproc) then + write(iulog,*) " " + write(iulog,*) "register MOAB app:", trim(appname), " MHID=", MHID + write(iulog,*) " " + endif + appname="HM_FINE"//C_NULL_CHAR + ATM_ID1 = 119 ! this number should not conflict with other components IDs; how do we know? + ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, MHFID) + if (ierr > 0 ) & + call abortmp('Error: cannot register moab app for fine mesh') + if(par%masterproc) then + write(iulog,*) " " + write(iulog,*) "register MOAB app:", trim(appname), " MHFID=", MHFID + write(iulog,*) " " + endif + appname="HM_PGX"//C_NULL_CHAR + ATM_ID1 = ATMID(1) ! this number should not conflict with other components IDs; how do we know? + ! + ! in this case, we reuse the main atm id, mhid will not be used for intersection anymore + ! still, need to be careful + ierr = iMOAB_RegisterApplication(appname, par%comm, ATM_ID1, mhpgid) + if (ierr > 0 ) & + call abortmp('Error: cannot register moab app for fine mesh') + if(par%masterproc) then + write(iulog,*) " " + write(iulog,*) "register MOAB app:", trim(appname), " MHPGID=", mhpgid + write(iulog,*) " " + endif ! instance distributed moab meshes from elem structures ! 1 ) spectral coarse mesh ! 2 ) GLL fine quad mesh (used mostly for visualization) ! 3 ) pgN FV type mesh, (most of the time pg2 mesh), used for coupling with other components; - call create_moab_meshes(par, elem, pgN) - endif + 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 4841cb84b543..0c5725edd2b2 100644 --- a/components/homme/src/share/semoab_mod.F90 +++ b/components/homme/src/share/semoab_mod.F90 @@ -76,6 +76,749 @@ integer function search_in(intarr, leng, value) end function search_in + 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 ? + 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 + 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 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 60f94f6992b3..18a80bfcc1bb 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 @@ -528,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 @@ -593,8 +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 - - call seq_map_map(mapper_Sa2l, av_s=dom_s%data, av_d=dom_d%data, fldlist='aream') + ! 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 @@ -734,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) ! @@ -748,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 @@ -784,7 +791,8 @@ 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 + + factors = 1.0_r8 ! initialize with 1.0 all factors; ! get areas tagname='area:aream:mask'//C_NULL_CHAR arrsize = 3 * lsize @@ -792,24 +800,26 @@ subroutine component_init_areacor_moab (comp, mbccid, mbcxid, seq_flds_c2x_fluxe 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 + if (.not.samegrid) then ! only correct if not monogrid (SCM case) + ! 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 ! 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 diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 6d5c4bdc638c..e14aadcb2dbd 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -10,7 +10,7 @@ module cplcomp_exchange_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_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 @@ -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 @@ -995,12 +996,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 (.not.atm_pg_active) then + ! this is the spectral monogrid case + 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 ) @@ -1059,7 +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 + 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 @@ -1122,11 +1129,12 @@ 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 - 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 +1153,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 +1192,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 +1213,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 +1221,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,21 +1232,18 @@ 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 + + tagname = trim(seq_flds_x2a_fields)//C_NULL_CHAR + numco = 1 ! usually 1 value per cell - 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 +1251,25 @@ subroutine cplcomp_moab_Init(infodata,comp) endif !add the normalization tag + 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 ) 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') + 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) + endif ! coupler pes #ifdef MOABDEBUG @@ -1541,7 +1558,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,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 + endif ! use land full mesh if (MPI_COMM_NULL /= mpicom_new ) then ! we are on the coupler pes appname = "COUPLE_LAND"//C_NULL_CHAR @@ -1551,27 +1574,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 @@ -1605,11 +1661,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 @@ -1630,24 +1691,29 @@ 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%mbGridType = typeA - 2 ! 0 or 1, pc or cells + comp%mbApCCid = mlnid ! land + comp%mbGridType = 0 ! 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') + if (mblxid > 0) then ! on coupler pes only + call copy_aream_from_area(mblxid) + endif #ifdef MOABDEBUG outfile = 'recMeshLand.h5m'//C_NULL_CHAR @@ -1898,7 +1964,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 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 d5a923722497..7937568f2b12 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 @@ -1033,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 @@ -1063,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,lnd_present #ifdef MOABDEBUG character*32 :: outfile, wopts, lnum #endif @@ -1077,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,lnd_present=lnd_present) if (first_time) then @@ -1088,7 +1104,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 @@ -1105,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' @@ -1296,7 +1317,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) @@ -1334,32 +1359,48 @@ subroutine prep_atm_mrg_moab(infodata, xao_ax) call shr_sys_abort(subname//' error in getting fractions_am from atm instance ') 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 ') + 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 + else + o2x_am(:,:)=0.0_r8 endif - - 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 ') + + 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 + else + i2x_am(:,:)=0.0_r8 endif - 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 ') + 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 + else + l2x_am(:,:)=0.0_r8 endif - 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 ') + 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 + else + xao_am(:,:)=0.0_r8 endif @@ -1367,6 +1408,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 --- @@ -1421,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 @@ -1494,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 @@ -1536,7 +1580,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) @@ -1549,7 +1597,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/main/prep_lnd_mod.F90 b/driver-moab/main/prep_lnd_mod.F90 index 08d208d1c60e..b73080a10f3a 100644 --- a/driver-moab/main/prep_lnd_mod.F90 +++ b/driver-moab/main/prep_lnd_mod.F90 @@ -15,12 +15,14 @@ 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 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 @@ -399,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 @@ -592,9 +600,13 @@ 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 + 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 - type2 = 3; ! FV mesh on coupler land 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/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 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 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 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() 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 6831e21dbce7..c03d525ba18e 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,mblxid,mb_scm_land 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 ) @@ -1670,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 @@ -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) @@ -1709,10 +1710,25 @@ 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 + ierr = iMOAB_GetGlobalInfo( mbxid, ngv, 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 ((.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 + 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 @@ -1724,11 +1740,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 ) - ns = nvise(1) ! local cells - +#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) @@ -2623,6 +2638,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 +2646,7 @@ end subroutine seq_io_read_char ! ! !REVISION HISTORY: ! 2025-07-20 - Cursor - initial documentation + ! 2025-09-24 spectral case atm ! ! !INTERFACE: ------------------------------------------------------------------ @@ -2677,7 +2694,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) ' @@ -2687,7 +2704,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) @@ -2711,8 +2727,19 @@ 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) - lnx = 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 ((.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 + 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 @@ -2721,9 +2748,6 @@ subroutine seq_io_read_moab_tags(filename, mbxid, dname, tag_list, matrix, nx) #endif lnx = 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 allocate(data1(ns)) allocate(data_reorder(ns)) allocate(dof(ns)) 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 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) 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 diff --git a/driver-moab/shr/seq_flds_mod.F90 b/driver-moab/shr/seq_flds_mod.F90 index f8abf9e4c396..9aa3be8d5b1a 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 @@ -260,7 +258,6 @@ module seq_flds_mod character(CXX) :: seq_flds_dom_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_x2a_fields character(CXX) :: seq_flds_i2x_fields @@ -390,7 +387,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 +420,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) @@ -4098,47 +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 - ! 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) - - - + 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 diff --git a/driver-moab/shr/shr_moab_mod.F90 b/driver-moab/shr/shr_moab_mod.F90 index f07d4da4ad6a..9fde5f4d5c71 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,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 @@ -53,7 +54,13 @@ 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)) .or. & + (mb_scm_land .and. (mblxid .eq. moabid))) then + mbGetnCells = nvert(1) + else + mbGetnCells = nvise(1) + endif end function mbGetnCells @@ -84,8 +91,13 @@ subroutine mbGetCellTagVals(mbid, intag,inarray,nMax) character(*), parameter :: subname = '(mbGetCellTagVals) ' !----------------------------------------------------------------------- ! - - ent_type=1 + 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 + endif + tagname = trim(intag)//C_NULL_CHAR ierr = iMOAB_GetDoubleTagStorage (mbid, tagname, nMax , ent_type, inarray) if (ierr .ne. 0) then @@ -123,7 +135,12 @@ subroutine mbSetCellTagVals(mbid, intag,inarray,nMax) !----------------------------------------------------------------------- ! - ent_type=1 + 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 + endif tagname = trim(intag)//C_NULL_CHAR ierr = iMOAB_SetDoubleTagStorage (mbid, tagname, nMax , ent_type, inarray) if (ierr .ne. 0) then