From 17fc106322f2abc3e2bca76a8121edc3663a8db9 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Fri, 21 Aug 2026 00:13:11 -0500 Subject: [PATCH 1/7] Fix e3sm.exe link order so PnetCDF comes after SCORPIO's piof/pioc Linking the netcdf interface target (which bundles libpnetcdf.a) directly onto e3sm.exe pinned it ahead of piof/pioc on the static link line, since CMake collapses repeated references to the same target to one position. GNU ld resolves static archives in a single pass, so piof/pioc's calls into PnetCDF (ncmpi_*) went unresolved, failing the link. Only pull netcdf's include dirs onto the cpl target directly, and let its link libraries reach the executable exclusively through csm_share -> spio, where PnetCDF is already ordered after piof/pioc. Co-Authored-By: Claude Sonnet 5 --- components/cmake/build_model.cmake | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/components/cmake/build_model.cmake b/components/cmake/build_model.cmake index 18d5ab728bdd..d9c489640ccb 100644 --- a/components/cmake/build_model.cmake +++ b/components/cmake/build_model.cmake @@ -263,19 +263,31 @@ macro(build_model COMP_CLASS COMP_NAME) add_executable(${TARGET_NAME}) target_sources(${TARGET_NAME} PRIVATE ${REAL_SOURCES}) - # driver-mct/main sources (e.g. cime_comp_mod.F90) use netcdf directly, but - # the component libraries only link netcdf PRIVATEly (via csm_share), so - # its usage requirements (e.g. include dirs for netcdf.mod) do not - # propagate up to this exe target. Find/link it explicitly here too. - find_package(NETCDF REQUIRED) - target_link_libraries(${TARGET_NAME} netcdf) - foreach(ITEM IN LISTS COMP_CLASSES) if (NOT ITEM STREQUAL "cpl") target_link_libraries(${TARGET_NAME} ${ITEM}) endif() endforeach() + # driver-mct/main sources (e.g. cime_comp_mod.F90) use netcdf directly, but + # the component libraries only link netcdf PRIVATEly (via csm_share), so + # its usage requirements (e.g. include dirs for netcdf.mod) do not + # propagate up to this exe target. Pull in just the include dirs here. + # + # Deliberately NOT using target_link_libraries(${TARGET_NAME} netcdf): + # CMake treats "netcdf" as a single graph node, so linking it directly to + # this executable pins its (and PnetCDF's) position to wherever it's + # first encountered, ahead of the component libraries below that also + # transitively link it via csm_share -> spio. But piof/pioc (also pulled + # in via csm_share -> spio) reference PnetCDF symbols, so libpnetcdf.a + # must come after them on the link line, or the linker drops its symbols + # before piof/pioc need them. Only requesting the include dirs avoids + # adding that extra, mis-ordering link edge, while the link libraries + # still reach this target correctly ordered via csm_share -> spio. + find_package(NETCDF REQUIRED) + get_target_property(NETCDF_INTERFACE_INCLUDE_DIRS netcdf INTERFACE_INCLUDE_DIRECTORIES) + target_include_directories(${TARGET_NAME} PRIVATE ${NETCDF_INTERFACE_INCLUDE_DIRS}) + if (USE_MOAB) target_link_libraries(${TARGET_NAME} ${MOAB_LIBRARIES}) target_include_directories(${TARGET_NAME} PRIVATE ${MOAB_INCLUDE_DIRS}) From f058f9eb14ba8c3dbb997a2c13478f2d09f8168f Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Fri, 21 Aug 2026 14:00:32 -0500 Subject: [PATCH 2/7] Support planar MPAS meshes in the MOAB mesh instance Handle MPAS meshes that are not on a sphere (on_a_sphere false or sphere_radius == 0), such as the MALI planar Greenland mesh. Build vertex coordinates from latVertex/lonVertex projected onto the unit sphere, and normalize areaCell by the Earth radius squared so cell areas are in steradians, consistent with every other coupler mesh. Spherical meshes are unaffected. Co-authored-by: Iulian Grindeanu Co-Authored-By: Claude Fable 5 --- .../src/framework/mpas_moabmesh.F | 41 +++++++++++++++---- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/components/mpas-framework/src/framework/mpas_moabmesh.F b/components/mpas-framework/src/framework/mpas_moabmesh.F index 7079f0351e3a..3987869e717c 100644 --- a/components/mpas-framework/src/framework/mpas_moabmesh.F +++ b/components/mpas-framework/src/framework/mpas_moabmesh.F @@ -30,6 +30,7 @@ subroutine init_moab_mpas(domain, ext_comp_id, pidmoab) iMOAB_ResolveSharedEntities, iMOAB_DetermineGhostEntities, & iMOAB_DefineTagStorage, iMOAB_SetIntTagStorage , & iMOAB_UpdateMeshInfo, iMOAB_SetDoubleTagStorage + use shr_const_mod, only : SHR_CONST_REARTH type (domain_type), intent(inout) :: domain integer , intent(in) :: ext_comp_id @@ -43,6 +44,7 @@ subroutine init_moab_mpas(domain, ext_comp_id, pidmoab) integer, dimension(:), pointer :: nEdgesOnCell integer, dimension(:), pointer :: indexToVertexID, indexToCellID real(kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex + real(kind=RKIND), dimension(:), pointer :: latVertex, lonVertex real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, areaCell logical, pointer :: on_a_sphere, is_periodic real(kind=RKIND), pointer :: x_period, y_period @@ -58,6 +60,8 @@ subroutine init_moab_mpas(domain, ext_comp_id, pidmoab) character*100 tagname, lnum integer tagtype, numco, tag_sto_len, ent_type, tagindex, currentVertex real (kind=RKIND), pointer :: sphere_radius + real (kind=RKIND) :: moab_radius, latvc, lonvc + logical :: planar_mesh c_comm = domain % dminfo % comm write(lnum,"(I0.2)")ext_comp_id @@ -98,6 +102,14 @@ subroutine init_moab_mpas(domain, ext_comp_id, pidmoab) ! call mpas_pool_get_array(meshPool, 'yCell', yCell) ! call mpas_pool_get_array(meshPool, 'zCell', zCell) + ! planar meshes (e.g. MALI land ice) have sphere_radius 0 and x/y/zVertex in + ! projected meters; build unit-sphere coordinates from latVertex/lonVertex instead + planar_mesh = (.not. on_a_sphere) .or. (sphere_radius <= 0.0_RKIND) + if (planar_mesh) then + call mpas_pool_get_array(meshPool, 'latVertex', latVertex) + call mpas_pool_get_array(meshPool, 'lonVertex', lonVertex) + endif + call mpas_log_write(' MOAB instance: number of vertices:: $i number of cells:: $i solve: v:$i c:$i', intArgs=(/nVertices, nCells, nVerticesSolve, nCellsSolve/) ) !! allocate(indexUsed(nVertices), invMap(nVertices) ) ! conservative, invMap should be smaller @@ -136,12 +148,27 @@ subroutine init_moab_mpas(domain, ext_comp_id, pidmoab) all_connects(i1) = indexUsed( all_connects(i1) ) enddo allocate(moab_vert_coords(3*currentVertex)) - do i1 =1, currentVertex - moab_vert_coords(3*i1-2) = xVertex(invMap(i1))/sphere_radius - moab_vert_coords(3*i1-1) = yVertex(invMap(i1))/sphere_radius - moab_vert_coords(3*i1 ) = zVertex(invMap(i1))/sphere_radius - ! call mpas_log_write('i:: $i coords:: $r $r $r $r', intArgs=(/i1/), realArgs=(/moab_vert_coords(3*i1-2),moab_vert_coords(3*i1-1), moab_vert_coords(3*i1)/) ) - enddo + if (planar_mesh) then + ! unit-sphere coordinates from vertex latitude/longitude; areas stay + ! consistent with the coupler (steradians) by dividing by Earth radius^2, + ! matching the areaCell*0.24635127E-13 conversion in glc_domain_mct + moab_radius = SHR_CONST_REARTH + do i1 =1, currentVertex + latvc = latVertex(invMap(i1)) + lonvc = lonVertex(invMap(i1)) + moab_vert_coords(3*i1-2) = cos( latvc ) * cos( lonvc ) ! x coordinate + moab_vert_coords(3*i1-1) = cos( latvc ) * sin( lonvc ) ! y coordinate + moab_vert_coords(3*i1 ) = sin( latvc ) ! z coordinate + enddo + else + moab_radius = sphere_radius + do i1 =1, currentVertex + moab_vert_coords(3*i1-2) = xVertex(invMap(i1))/sphere_radius + moab_vert_coords(3*i1-1) = yVertex(invMap(i1))/sphere_radius + moab_vert_coords(3*i1 ) = zVertex(invMap(i1))/sphere_radius + ! call mpas_log_write('i:: $i coords:: $r $r $r $r', intArgs=(/i1/), realArgs=(/moab_vert_coords(3*i1-2),moab_vert_coords(3*i1-1), moab_vert_coords(3*i1)/) ) + enddo + endif dimcoord = 3*currentVertex dimen = 3 ierr = iMOAB_CreateVertices(pid, dimcoord, dimen, moab_vert_coords) @@ -183,7 +210,7 @@ subroutine init_moab_mpas(domain, ext_comp_id, pidmoab) n=0 do ic=1, nCellsSolve n= n+1 - data(n)=areaCell(ic) / (sphere_radius * sphere_radius) + data(n)=areaCell(ic) / (moab_radius * moab_radius) enddo tagname='area'//C_NULL_CHAR From 1dbf142345ea1bc463ba8ded0b87ba8a39f282cf Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Fri, 21 Aug 2026 14:00:44 -0500 Subject: [PATCH 3/7] Create a MOAB mesh instance for MALI and refactor its import/export Register the MALI mesh with iMOAB on the component PEs (app id mbglid) via init_moab_mpas, define and zero the g2x, x2g, and domain tags, and add glc_domain_moab to set lon/lat/area/mask/frac on the mesh. MALI planar-mesh cell areas are converted to steradians. Refactor glc_import_mct and glc_export_mct to take raw real(:,:) arrays instead of MCT attribute vectors, following the mpas-ocean pattern, so a single copy of the field logic serves both the MCT and MOAB drivers. The MOAB path exchanges data with the tag storage through transposed twin arrays; the MCT path passes x2g%rAttr and g2x%rAttr directly. Declare the GLC iMOAB app ids (mbglid, mbgxid, mbintxlg, mbintxgl) in seq_comm_mct and initialize them to -1. Co-authored-by: Iulian Grindeanu Co-Authored-By: Claude Fable 5 --- .../mpas-albany-landice/driver/glc_comp_mct.F | 324 ++++++++++++++++-- driver-moab/shr/seq_comm_mct.F90 | 9 + 2 files changed, 305 insertions(+), 28 deletions(-) diff --git a/components/mpas-albany-landice/driver/glc_comp_mct.F b/components/mpas-albany-landice/driver/glc_comp_mct.F index 5454606b71b0..f3ded9f74274 100644 --- a/components/mpas-albany-landice/driver/glc_comp_mct.F +++ b/components/mpas-albany-landice/driver/glc_comp_mct.F @@ -48,6 +48,18 @@ module glc_comp_mct use iso_c_binding, only : c_char, c_loc, c_ptr, c_int use mpas_c_interfacing, only : mpas_f_to_c_string, mpas_c_to_f_string +#ifdef HAVE_MOAB + use shr_kind_mod , only: cxx => SHR_KIND_CXX + use mpas_moabmesh + use seq_comm_mct, only: MBGLID + use seq_comm_mct, only: num_moab_exports + use iMOAB, only: iMOAB_DefineTagStorage, iMOAB_SetDoubleTagStorage, & + iMOAB_GetMeshInfo, iMOAB_GetDoubleTagStorage +#ifdef MOABDEBUG + use iMOAB, only: iMOAB_WriteMesh +#endif +#endif + ! MALI modules use li_core use li_core_interface @@ -66,6 +78,14 @@ module glc_comp_mct private :: glc_domain_mct ! !PRIVATE MODULE VARIABLES +#ifdef HAVE_MOAB + private :: glc_domain_moab + integer , private :: mblsize, totalmbls, totalmbls_r + ! moab arrays for exchange with the coupler; tag storage is (cell, field) ordered, + ! while import/export routines expect (field, cell), hence the transpose twins + real (kind=RKIND) , allocatable, private :: g2x_gm(:,:), g2x_gm2(:,:) + real (kind=RKIND) , allocatable, private :: x2g_gm(:,:), x2g_gm2(:,:) +#endif integer, private :: my_task @@ -185,6 +205,14 @@ subroutine glc_init_mct( EClock, cdata_g, x2g_g, g2x_g, NLFilename )!{{{ logical, pointer :: tempLogicalConfig, config_create_all_logs_in_e3sm character(len=StrKIND), pointer :: tempCharConfig +#ifdef HAVE_MOAB + integer :: ierrmb, numco, tagtype, tagindex, ent_type + character(CXX) :: tagname +#ifdef MOABDEBUG + character(CXX) :: outfile, wopts +#endif +#endif + interface subroutine xml_stream_parser(xmlname, mgr_p, comm, ierr) bind(c) @@ -548,6 +576,18 @@ end subroutine xml_stream_get_attributes !----------------------------------------------------------------------- call t_stopf ('mali_init') +#ifdef HAVE_MOAB + call init_moab_mpas(domain, glcID, MBGLID) + call mpas_log_write('initialized MOAB MPAS land ice instance... ') +#ifdef MOABDEBUG + if (MBGLID .ge. 0 ) then ! we are on glc comp pes + outfile = 'GLC_mpas.h5m'//C_NULL_CHAR + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR + ierrmb = iMOAB_WriteMesh(MBGLID, trim(outfile), trim(wopts)) + endif +#endif +#endif + !----------------------------------------------------------------------- ! ! check for consistency of mali and sync clock initial time @@ -602,6 +642,54 @@ end subroutine xml_stream_get_attributes nsend = mct_avect_nRattr(g2x_g) nrecv = mct_avect_nRattr(x2g_g) +#ifdef HAVE_MOAB + ! initialize moab tag storage for the coupling fields, on cells; + ! also allocate the arrays used to move data between MALI and the tags + mblsize = lsize + totalmbls = mblsize * nsend ! size of the double array for g2x fields + allocate ( g2x_gm (lsize, nsend) ) + allocate ( g2x_gm2(nsend, lsize) ) + g2x_gm = 0.0_RKIND + g2x_gm2 = 0.0_RKIND + + tagtype = 1 ! dense, double + numco = 1 ! one value per cell + ent_type = 1 ! cells + tagname = trim(seq_flds_g2x_fields)//C_NULL_CHAR + ierrmb = iMOAB_DefineTagStorage(MBGLID, tagname, tagtype, numco, tagindex) + if ( ierrmb /= 0 ) then + call mpas_log_write('cannot define tags for MOAB g2x fields', MPAS_LOG_ERR) + endif + ierrmb = iMOAB_SetDoubleTagStorage ( MBGLID, tagname, totalmbls, ent_type, g2x_gm ) + if ( ierrmb /= 0 ) then + call mpas_log_write('fail to zero MOAB g2x fields', MPAS_LOG_ERR) + endif + + totalmbls_r = mblsize * nrecv ! size of the double array for x2g fields + allocate ( x2g_gm (lsize, nrecv) ) + allocate ( x2g_gm2(nrecv, lsize) ) + x2g_gm = 0.0_RKIND + x2g_gm2 = 0.0_RKIND + + tagname = trim(seq_flds_x2g_fields)//C_NULL_CHAR + ierrmb = iMOAB_DefineTagStorage(MBGLID, tagname, tagtype, numco, tagindex) + if ( ierrmb /= 0 ) then + call mpas_log_write('cannot define tags for MOAB x2g fields', MPAS_LOG_ERR) + endif + ierrmb = iMOAB_SetDoubleTagStorage ( MBGLID, tagname, totalmbls_r, ent_type, x2g_gm ) + if ( ierrmb /= 0 ) then + call mpas_log_write('fail to zero MOAB x2g fields', MPAS_LOG_ERR) + endif + + ! add domain tags and fill them with the same values used for the mct domain + tagname = trim(seq_flds_dom_fields)//C_NULL_CHAR + ierrmb = iMOAB_DefineTagStorage(MBGLID, tagname, tagtype, numco, tagindex) + if ( ierrmb /= 0 ) then + call mpas_log_write('cannot define tags for MOAB dom fields', MPAS_LOG_ERR) + endif + call glc_domain_moab(MBGLID) +#endif + !----------------------------------------------------------------------- ! ! initialize necessary coupling info @@ -736,10 +824,26 @@ end subroutine xml_stream_get_attributes ! !----------------------------------------------------------------------- - call glc_export_mct(g2x_g, errorCode) +#ifdef HAVE_MOAB + g2x_gm2 = 0.0_RKIND + call glc_export_mct(g2x_gm2, errorCode) if (errorCode /= 0) then call mpas_log_write('Error in glc_export_mct', MPAS_LOG_CRIT) endif + ! set the initial g2x fields into the MOAB tags on cells + g2x_gm = transpose(g2x_gm2) + tagname = trim(seq_flds_g2x_fields)//C_NULL_CHAR + ent_type = 1 ! cells + ierrmb = iMOAB_SetDoubleTagStorage ( MBGLID, tagname, totalmbls, ent_type, g2x_gm ) + if ( ierrmb /= 0 ) then + call mpas_log_write('fail to set MOAB g2x fields', MPAS_LOG_ERR) + endif +#else + call glc_export_mct(g2x_g%rAttr, errorCode) + if (errorCode /= 0) then + call mpas_log_write('Error in glc_export_mct', MPAS_LOG_CRIT) + endif +#endif ! Setup clock for initial runs ! Some MPASO ocean pools handled here, but we don't have them. @@ -829,6 +933,14 @@ subroutine glc_run_mct( EClock, cdata_g, x2g_g, g2x_g)!{{{ character(len=StrKIND) :: timeStamp, streamName integer :: err, err_tmp, globalErr, streamDirection, iam logical :: solveVelo, streamActive +#ifdef HAVE_MOAB + integer :: ierrmb, ent_type + character(CXX) :: tagname +#ifdef MOABDEBUG + integer :: cur_glc_stepno + character(CXX) :: outfile, wopts, lnum +#endif +#endif iam = domain % dminfo % my_proc_id @@ -860,7 +972,18 @@ subroutine glc_run_mct( EClock, cdata_g, x2g_g, g2x_g)!{{{ err = ior(err, err_tmp) ! Import state from coupler - call glc_import_mct(x2g_g, err_tmp) +#ifdef HAVE_MOAB + tagname = trim(seq_flds_x2g_fields)//C_NULL_CHAR + ent_type = 1 ! cells + ierrmb = iMOAB_GetDoubleTagStorage ( MBGLID, tagname, totalmbls_r, ent_type, x2g_gm ) + if ( ierrmb /= 0 ) then + call mpas_log_write('fail to get MOAB x2g fields', MPAS_LOG_ERR) + endif + x2g_gm2 = transpose(x2g_gm) + call glc_import_mct(x2g_gm2, err_tmp) +#else + call glc_import_mct(x2g_g%rAttr, err_tmp) +#endif err = ior(err,err_tmp) ! Initialize time average fields @@ -1026,7 +1149,26 @@ subroutine glc_run_mct( EClock, cdata_g, x2g_g, g2x_g)!{{{ ! end do ! Export state to coupler - call glc_export_mct(g2x_g, err_tmp) +#ifdef HAVE_MOAB + g2x_gm2 = 0.0_RKIND + call glc_export_mct(g2x_gm2, err_tmp) + g2x_gm = transpose(g2x_gm2) + tagname = trim(seq_flds_g2x_fields)//C_NULL_CHAR + ent_type = 1 ! cells + ierrmb = iMOAB_SetDoubleTagStorage ( MBGLID, tagname, totalmbls, ent_type, g2x_gm ) + if ( ierrmb /= 0 ) then + call mpas_log_write('fail to set MOAB g2x fields', MPAS_LOG_ERR) + endif +#ifdef MOABDEBUG + call seq_timemgr_EClockGetData( EClock, stepno=cur_glc_stepno ) + write(lnum,"(I0.2)") cur_glc_stepno + outfile = 'glc_export_'//trim(lnum)//'.h5m'//C_NULL_CHAR + wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR + ierrmb = iMOAB_WriteMesh(MBGLID, outfile, wopts) +#endif +#else + call glc_export_mct(g2x_g%rAttr, err_tmp) +#endif err = ior(err,err_tmp) ! Move time levels back the way MPAS will expect them on the next time step @@ -1115,6 +1257,13 @@ subroutine glc_final_mct( EClock, cdata_g, x2g_g, g2x_g)!{{{ call mpas_framework_finalize(domain % dminfo, domain, io_system) +#ifdef HAVE_MOAB + if (allocated(g2x_gm)) deallocate(g2x_gm) + if (allocated(g2x_gm2)) deallocate(g2x_gm2) + if (allocated(x2g_gm)) deallocate(x2g_gm) + if (allocated(x2g_gm2)) deallocate(x2g_gm2) +#endif + ! Reset I/O logs call shr_file_setLogUnit (shrlogunit) call shr_file_setLogLevel(shrloglev) @@ -1370,6 +1519,120 @@ subroutine glc_domain_mct( lsize, gsMap_g, dom_g )!{{{ !----------------------------------------------------------------------- end subroutine glc_domain_mct!}}} +#ifdef HAVE_MOAB +!*********************************************************************** +! +! !IROUTINE: glc_domain_moab +! +! !INTERFACE: + subroutine glc_domain_moab( mbid )!{{{ + +! !DESCRIPTION: +! Set the domain tags (lat, lon, area, mask, frac) on the MOAB mesh for MALI; +! same values as the mct domain set by glc_domain_mct + + use shr_const_mod, only: SHR_CONST_PI + use iMOAB, only: iMOAB_SetDoubleTagStorage, iMOAB_GetMeshInfo + + implicit none + integer, intent(in) :: mbid ! moab app id for the MALI mesh on component pes + + real(kind=RKIND), allocatable :: data(:) + real(kind=RKIND) :: r2d + integer :: i, n, ierr, arrsize, ent_type + character(CXX) :: tagname + + type (block_type), pointer :: block_ptr + type (mpas_pool_type), pointer :: meshPool + integer, pointer :: nCellsSolve + real (kind=RKIND), dimension(:), pointer :: lonCell, latCell, areaCell + + integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) + + r2d = 180.0_RKIND / SHR_CONST_PI + + ! find out the number of local cells in the moab mesh + ierr = iMOAB_GetMeshInfo ( mbid, nvert, nvise, nbl, nsurf, nvisBC ) + if ( ierr /= 0 ) then + call mpas_log_write('fail to get moab mesh info for MALI', MPAS_LOG_ERR) + endif + arrsize = nvise(1) + allocate(data(arrsize)) + ent_type = 1 ! cells + + n = 0 + block_ptr => domain % blocklist + do while(associated(block_ptr)) + call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(meshPool, 'lonCell', lonCell) + do i = 1, nCellsSolve + n = n + 1 + data(n) = lonCell(i) * r2d + end do + block_ptr => block_ptr % next + end do + tagname = 'lon'//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage ( mbid, tagname, arrsize, ent_type, data ) + if ( ierr /= 0 ) then + call mpas_log_write('fail to set moab lon tag for MALI', MPAS_LOG_ERR) + endif + + n = 0 + block_ptr => domain % blocklist + do while(associated(block_ptr)) + call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(meshPool, 'latCell', latCell) + do i = 1, nCellsSolve + n = n + 1 + data(n) = latCell(i) * r2d + end do + block_ptr => block_ptr % next + end do + tagname = 'lat'//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage ( mbid, tagname, arrsize, ent_type, data ) + if ( ierr /= 0 ) then + call mpas_log_write('fail to set moab lat tag for MALI', MPAS_LOG_ERR) + endif + + n = 0 + block_ptr => domain % blocklist + do while(associated(block_ptr)) + call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) + ! Note: only supporting planar meshes, same as glc_domain_mct + do i = 1, nCellsSolve + n = n + 1 + data(n) = areaCell(i) * 0.24635127E-13_RKIND ! This is the conversion factor to rad^2 + end do + block_ptr => block_ptr % next + end do + tagname = 'area'//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage ( mbid, tagname, arrsize, ent_type, data ) + if ( ierr /= 0 ) then + call mpas_log_write('fail to set moab area tag for MALI', MPAS_LOG_ERR) + endif + + ! For now, assume mask and frac are 1 everywhere, same as glc_domain_mct + data(:) = 1.0_RKIND + tagname = 'mask'//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage ( mbid, tagname, arrsize, ent_type, data ) + if ( ierr /= 0 ) then + call mpas_log_write('fail to set moab mask tag for MALI', MPAS_LOG_ERR) + endif + tagname = 'frac'//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage ( mbid, tagname, arrsize, ent_type, data ) + if ( ierr /= 0 ) then + call mpas_log_write('fail to set moab frac tag for MALI', MPAS_LOG_ERR) + endif + + deallocate(data) + +!----------------------------------------------------------------------- + end subroutine glc_domain_moab!}}} +#endif !================================================================================= @@ -1378,7 +1641,9 @@ subroutine glc_import_mct(x2g_g, errorCode) ! !INPUT/OUTPUT PARAMETERS: - type(mct_aVect) , intent(inout) :: x2g_g + ! raw field data, shaped (nfields, ncells); the mct caller passes the + ! attribute vector rAttr array, the moab caller passes the transposed tag data + real (kind=RKIND), dimension(:,:), intent(in) :: x2g_g ! !OUTPUT PARAMETERS: @@ -1452,22 +1717,22 @@ subroutine glc_import_mct(x2g_g, errorCode) do i = 1, nCellsSolve n = n + 1 if (trim(config_smb_source) == 'coupler') then - sfcMassBal(i) = x2g_g % rAttr(index_x2g_Flgl_qice, n) + sfcMassBal(i) = x2g_g(index_x2g_Flgl_qice, n) elseif (trim(config_smb_source) == 'file') then ! do nothing, use data read-in from file else call mpas_log_write("Unknown value for 'config_smb_source'", MPAS_LOG_ERR) errorCode = ior(errorCode, 1) endif - floatingBasalMassBal(i) = x2g_g % rAttr(index_x2g_Fogx_qiceli, n) + floatingBasalMassBal(i) = x2g_g(index_x2g_Fogx_qiceli, n) if (nISMIP6OceanLayers > 0) then do iLev = 1, nISMIP6OceanLayers - fractionalMaskVal = x2g_g % rAttr(index_x2g_So_tf3d_mask(iLev), n) + fractionalMaskVal = x2g_g(index_x2g_So_tf3d_mask(iLev), n) if (fractionalMaskVal > 0.5_RKIND) then ! Only use MALI grid cells that have at least 50% overlap with MPAS-Ocean cells ! if valid, mark our integer mask as 1 and scale the TF value by the fractional mask orig3dOceanMask(iLev,i) = 1 - ismip6shelfMelt_3dThermalForcing(iLev, i) = x2g_g % rAttr(index_x2g_So_tf3d(iLev), n) / fractionalMaskVal + ismip6shelfMelt_3dThermalForcing(iLev, i) = x2g_g(index_x2g_So_tf3d(iLev), n) / fractionalMaskVal else ! if not using, insert 0 in mask and bad val in TF field orig3dOceanMask(iLev,i) = 0 @@ -1475,10 +1740,10 @@ subroutine glc_import_mct(x2g_g, errorCode) endif enddo endif -! surfaceTemperature(i) = x2g_g % rAttr(index_x2g_Sl_tsrf, n) -!JW basalOceanHeatflx(i) = x2g_g % rAttr(index_x2g_Fogo_qiceh, n) -! basalOceanHeatflx(i) = x2g_g % rAttr(index_x2g_Fogx_qicehi, n) - !OceanDensity (i) = x2g_g % rAttr(, n) +! surfaceTemperature(i) = x2g_g(index_x2g_Sl_tsrf, n) +!JW basalOceanHeatflx(i) = x2g_g(index_x2g_Fogo_qiceh, n) +! basalOceanHeatflx(i) = x2g_g(index_x2g_Fogx_qicehi, n) + !OceanDensity (i) = x2g_g(, n) end do block => block % next @@ -1496,7 +1761,10 @@ subroutine glc_export_mct(g2x_g, errorCode) !------------------------------------------------------------------- - type(mct_aVect) , intent(inout) :: g2x_g + ! raw field data, shaped (nfields, ncells); the mct caller passes the + ! attribute vector rAttr array, the moab caller passes an array that is + ! transposed afterwards and set into the moab tags + real (kind=RKIND), dimension(:,:), intent(inout) :: g2x_g integer, intent(out) :: errorCode !------------------------------------------------------------------- @@ -1560,45 +1828,45 @@ subroutine glc_export_mct(g2x_g, errorCode) n = n + 1 ! Fogg_rofl - g2x_g % rAttr(index_g2x_Fogg_rofl,n) = avgBareIceAblationApplied(i) + g2x_g(index_g2x_Fogg_rofl,n) = avgBareIceAblationApplied(i) ! Figg_rofi - g2x_g % rAttr(index_g2x_Figg_rofi,n) = 0.0 ! placeholder + g2x_g(index_g2x_Figg_rofi,n) = 0.0 ! placeholder ! Fogg_rofi - g2x_g % rAttr(index_g2x_Fogg_rofi,n) = avgCalvingFlux(i) - g2x_g % rAttr(index_g2x_Fogg_rofi,n) = g2x_g % rAttr(index_g2x_Fogg_rofi,n) + avgFaceMeltFlux(i) + g2x_g(index_g2x_Fogg_rofi,n) = avgCalvingFlux(i) + g2x_g(index_g2x_Fogg_rofi,n) = g2x_g(index_g2x_Fogg_rofi,n) + avgFaceMeltFlux(i) if (trim(config_basal_mass_bal_float) == 'ismip6') then ! if MALI is calculating ISMF, add that to rofi ! In some configurations, ISMF will be calculated in coupler or MPAS-Ocean ! Note avgFloatingBMBFlux is + for mass gain to ice sheet, so sign should be flipped for rofi - g2x_g % rAttr(index_g2x_Fogg_rofi,n) = g2x_g % rAttr(index_g2x_Fogg_rofi,n) - avgFloatingBMBFlux(i) + g2x_g(index_g2x_Fogg_rofi,n) = g2x_g(index_g2x_Fogg_rofi,n) - avgFloatingBMBFlux(i) endif - g2x_g % rAttr(index_g2x_Sg_topo, n) = max(0.0, upperSurface(i)) !updated to avoid warning for values below sea level - g2x_g % rAttr(index_g2x_Sg_tbot, n) = temperature(nVertlevels,i) - SHR_CONST_TKTRIP + g2x_g(index_g2x_Sg_topo, n) = max(0.0, upperSurface(i)) !updated to avoid warning for values below sea level + g2x_g(index_g2x_Sg_tbot, n) = temperature(nVertlevels,i) - SHR_CONST_TKTRIP ! layerThickness is not populated on init, so calculating like this instead: - g2x_g % rAttr(index_g2x_Sg_dztbot, n) = layerThicknessFractions(nVertLevels) * thickness(i) / 2.0 - g2x_g % rAttr(index_g2x_Sg_lithop, n) = Thickness(i)*SHR_CONST_RHOICE*SHR_CONST_G + g2x_g(index_g2x_Sg_dztbot, n) = layerThicknessFractions(nVertLevels) * thickness(i) / 2.0 + g2x_g(index_g2x_Sg_lithop, n) = Thickness(i)*SHR_CONST_RHOICE*SHR_CONST_G !!!Mask configurations - g2x_g % rAttr(index_g2x_Sg_ice_covered, n) = li_mask_is_ice_int(cellMask(i)) + g2x_g(index_g2x_Sg_ice_covered, n) = li_mask_is_ice_int(cellMask(i)) ! Sg_icemask should be 1 in locations where GLC will use a SMB value ! This is cells with ice or bare land (no ocean cells) ! This can easily be checked by any cells where upperSurface is ! above sea level! if (upperSurface(i) > config_sea_level) then - g2x_g % rAttr(index_g2x_Sg_icemask, n) = 1 + g2x_g(index_g2x_Sg_icemask, n) = 1 else - g2x_g % rAttr(index_g2x_Sg_icemask, n) = 0 + g2x_g(index_g2x_Sg_icemask, n) = 0 endif ! icemask_coupled_fluxes for now should be the same as icemask ! If we add in the ability to prevent the ice sheet from evolving ! (i.e., to ignore these fluxes) then this should be set to 0 everywhere. - g2x_g % rAttr(index_g2x_Sg_icemask_coupled_fluxes, n) = g2x_g % rAttr(index_g2x_Sg_icemask, n) + g2x_g(index_g2x_Sg_icemask_coupled_fluxes, n) = g2x_g(index_g2x_Sg_icemask, n) ! Add masks used for ice shelf basal melt calculation - g2x_g % rAttr(index_g2x_Sg_icemask_floating, n) = li_mask_is_floating_ice_int(cellMask(i)) - g2x_g % rAttr(index_g2x_Sg_icemask_grounded, n) = li_mask_is_grounded_ice_int(cellMask(i)) + g2x_g(index_g2x_Sg_icemask_floating, n) = li_mask_is_floating_ice_int(cellMask(i)) + g2x_g(index_g2x_Sg_icemask_grounded, n) = li_mask_is_grounded_ice_int(cellMask(i)) end do diff --git a/driver-moab/shr/seq_comm_mct.F90 b/driver-moab/shr/seq_comm_mct.F90 index 7937e69cf969..126ee0316e65 100644 --- a/driver-moab/shr/seq_comm_mct.F90 +++ b/driver-moab/shr/seq_comm_mct.F90 @@ -243,6 +243,10 @@ module seq_comm_mct integer, public :: mbintxlr ! iMOAB id for intersection mesh between land and river integer, public :: mbintxrl ! iMOAB id for intersection mesh between river and land integer, public :: mbintxri ! iMOAB id for intersection mesh between river and ice + integer, public :: mbglid ! iMOAB id for glc (MALI land ice), on glc component pes + integer, public :: mbgxid ! iMOAB id for glc migrated mesh to coupler pes + integer, public :: mbintxlg ! iMOAB id for read map between land and glc + integer, public :: mbintxgl ! iMOAB id for read map between glc and land integer, public :: num_moab_exports ! iMOAB id for atm phys grid, on atm pes @@ -695,6 +699,11 @@ subroutine seq_comm_init(global_comm_in, driver_comm_in, nmlfile, drv_comm_id) mbintxar = -1 ! iMOAB id for intx mesh between atm and river mbintxlr = -1 ! iMOAB id for intx mesh between land and river mbintxrl = -1 ! iMOAB id for intx mesh between river and land + mbintxri = -1 ! iMOAB id for intx mesh between river and ice + mbglid = -1 ! iMOAB id for glc (MALI land ice) on component pes + mbgxid = -1 ! iMOAB id for glc migrated to coupler pes + mbintxlg = -1 ! iMOAB id of moab instance of map read from lnd2glc map file + mbintxgl = -1 ! iMOAB id of moab instance of map read from glc2lnd map file num_moab_exports = 0 ! mostly used in debugging deallocate(comps,comms) From 89fd6b601a09fd6a0e23a608eea77c9992ed547f Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Fri, 21 Aug 2026 14:00:59 -0500 Subject: [PATCH 4/7] Migrate the GLC mesh to the coupler in driver-moab Add cplcomp_moab_init_glc to send the MALI mesh from the component PEs and receive it on the coupler PEs as app COUPLE_MPASLNDICE (mbgxid), define the g2x, x2g, and domain tags there, and exchange the domain tags so aream defaults to area until a map file provides it. Dispatch to it from cplcomp_moab_init for the 'g' component. Initialize the glc fractions on the coupler mesh: define gfrac and lfrac tags on mbgxid and set gfrac from the migrated domain frac. Enable the glc domain, fraction, g2x, and x2g writes in the coupler history file, and fix seq_diag_glc_moab to read from mbgxid. After this change the GLC mesh reaches the coupler and appears in cpl history files; no fields flow yet. Co-Authored-By: Claude Fable 5 --- driver-moab/main/cplcomp_exchange_mod.F90 | 91 ++++++++++++++++++++++- driver-moab/main/seq_diag_moab.F90 | 4 +- driver-moab/main/seq_frac_mct.F90 | 26 ++++++- driver-moab/main/seq_hist_mod.F90 | 24 +++--- 4 files changed, 128 insertions(+), 17 deletions(-) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 0e0a9d9d8db2..a0a14ab6ff82 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -16,6 +16,7 @@ module cplcomp_exchange_mod 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 use seq_flds_mod, only: seq_flds_l2x_fields, seq_flds_x2l_fields ! use seq_flds_mod, only: seq_flds_r2x_fields, seq_flds_x2r_fields, seq_flds_r2x_fluxes + use seq_flds_mod, only: seq_flds_g2x_fields, seq_flds_x2g_fields ! needed for MOAB init of glc fields on coupler side use seq_comm_mct, only: cplid, logunit use seq_comm_mct, only: seq_comm_getinfo => seq_comm_setptrs, seq_comm_iamin use cplcomp_moab_helpers_mod, only: moab_register_app, moab_send_mesh, moab_receive_mesh, & @@ -30,6 +31,7 @@ module cplcomp_exchange_mod 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 + use seq_comm_mct, only : mbglid, mbgxid ! iMOAB id of moab glc app on comp pes and on coupler too use shr_mpi_mod, only: shr_mpi_max ! use dimensions_mod, only : np ! for atmosphere use iso_c_binding @@ -60,6 +62,7 @@ module cplcomp_exchange_mod private :: cplcomp_moab_init_lnd private :: cplcomp_moab_init_ice private :: cplcomp_moab_init_rof + private :: cplcomp_moab_init_glc private :: cplcomp_moab_resolve_comm_types private :: cplcomp_moab_compute_comm_graph private :: cplcomp_moab_atm_phys_cid @@ -807,6 +810,86 @@ subroutine cplcomp_moab_init_ice(infodata, comp, id_old, id_join, mpicom_old, mp end subroutine cplcomp_moab_init_ice + subroutine cplcomp_moab_init_glc(infodata, comp, id_old, id_join, mpicom_old, mpicom_new, mpicom_join, dead_comps, partMethod, subname) + + ! Migrate the MALI (land ice) MOAB mesh from glc component pes to the + ! coupler pes, and define the coupling field tags on the coupler side. + ! There is no data-glc variant; the glc mesh is always a real MPAS mesh + ! created on the component pes by init_moab_mpas. + + use iMOAB, only: iMOAB_WriteMesh, iMOAB_GetMeshInfo + use seq_infodata_mod, only: seq_infodata_type + + type(seq_infodata_type), intent(in) :: infodata + type(component_type), intent(inout) :: comp + integer, intent(in) :: id_old, id_join + integer, intent(in) :: mpicom_old, mpicom_new, mpicom_join + logical, intent(in) :: dead_comps + integer, intent(in) :: partMethod + character(len=*), intent(in) :: subname + + integer :: mpigrp_cplid, mpigrp_old + integer :: ierr + character*200 :: appname, outfile, wopts + integer :: nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) + + call seq_comm_getinfo(cplid ,mpigrp=mpigrp_cplid) ! receiver group + call seq_comm_getinfo(id_old,mpigrp=mpigrp_old) ! component group pes + if (MPI_COMM_NULL /= mpicom_old ) then ! it means we are on the component pes +#ifdef MOABDEBUG + outfile = 'wholeGlc.h5m'//C_NULL_CHAR + wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbglid, outfile, wopts) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing glc mesh' + call shr_sys_abort(subname//' ERROR in writing glc mesh') + endif +#endif + if (mbglid >= 0) then + ierr = iMOAB_GetMeshInfo ( mbglid, nvert, nvise, nbl, nsurf, nvisBC ) + comp%mbApCCid = mbglid ! glc imoab app id + comp%mbGridType = 1 ! 0 or 1, pc or cells + comp%mblsize = nvise(1) ! cells + endif + ! send glc mesh to coupler + call moab_send_mesh(mbglid, mpicom_join, mpigrp_cplid, id_join, partMethod, subname) + endif + if (MPI_COMM_NULL /= mpicom_new ) then ! we are on the coupler pes + appname = "COUPLE_MPASLNDICE" + ! migrated mesh gets another app id, moab glc to coupler (mbgxid) + call moab_register_app(appname, mpicom_new, id_join, mbgxid, subname) + call moab_receive_mesh(mbgxid, mpicom_join, mpigrp_old, id_old, subname) + + call moab_define_double_tag(mbgxid, trim(seq_flds_g2x_fields), subname) + call moab_define_double_tag(mbgxid, trim(seq_flds_x2g_fields), subname) + call moab_define_double_tag(mbgxid, trim(seq_flds_dom_fields)//":norm8wt", subname) + endif + + if (mbglid .ge. 0) then ! we are on component glc pes + call moab_free_sender_buffers(mbglid, id_join, subname) + endif + + ! transport the domain tags to the coupler side, and set aream = area there + ! as a fallback; when land is present, glc aream is overwritten later from the + ! lnd2glc map file area_b (prep_glc_init), matching the mct driver behavior + call moab_exchange_domain_tags(comp, mbglid, mbgxid, 'lat:lon:area:frac:mask', 'domg') + +#ifdef MOABDEBUG + if (mbgxid >= 0) then ! coupler pes only +! debug test + outfile = 'recLndIce.h5m'//C_NULL_CHAR + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! +! write out the mesh file to disk + ierr = iMOAB_WriteMesh(mbgxid, trim(outfile), trim(wopts)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing glc mesh on coupler ' + call shr_sys_abort(subname//' ERROR in writing glc mesh on coupler ') + endif + endif +#endif + + end subroutine cplcomp_moab_init_glc + subroutine cplcomp_moab_pick_load_options(filename, ropts, ierr) ! Inspect filename's netCDF contents and choose MOAB load options that match ! the file's grid format. MOAB's parallel readers couple format to partition @@ -1023,7 +1106,7 @@ subroutine cplcomp_moab_Init(infodata,comp) character(len=*),parameter :: subname = "(cplcomp_moab_Init) " - integer :: maxMH, maxMPO, maxMLID, maxMSID, maxMRID + integer :: maxMH, maxMPO, maxMLID, maxMSID, maxMRID, maxMGID integer :: partMethod logical :: dead_comps @@ -1044,6 +1127,7 @@ subroutine cplcomp_moab_Init(infodata,comp) call shr_mpi_max(mlnid, maxMLID, mpicom_join, all=.true.) call shr_mpi_max(MPSIID, maxMSID, mpicom_join, all=.true.) call shr_mpi_max(mrofid, maxMRID, mpicom_join, all=.true.) + call shr_mpi_max(mbglid, maxMGID, mpicom_join, all=.true.) if (seq_comm_iamroot(CPLID) ) then write(logunit, *) "MOAB coupling for ", comp%oneletterid,' ', comp%ntype @@ -1078,6 +1162,11 @@ subroutine cplcomp_moab_Init(infodata,comp) call cplcomp_moab_init_rof(infodata, comp, id_old, id_join, mpicom_old, mpicom_new, mpicom_join, dead_comps, partMethod, subname) endif ! end for rof coupler set up +!!!!!!!!!!!!!!!! LAND ICE (GLC) + if (comp%oneletterid == 'g' .and. maxMGID /= -1) then + call cplcomp_moab_init_glc(infodata, comp, id_old, id_join, mpicom_old, mpicom_new, mpicom_join, dead_comps, partMethod, subname) + endif ! end for glc coupler set up + end subroutine cplcomp_moab_Init diff --git a/driver-moab/main/seq_diag_moab.F90 b/driver-moab/main/seq_diag_moab.F90 index bcc0b8d48d5d..8dba6b47646f 100644 --- a/driver-moab/main/seq_diag_moab.F90 +++ b/driver-moab/main/seq_diag_moab.F90 @@ -1234,14 +1234,14 @@ end subroutine seq_diag_rof_moab subroutine seq_diag_glc_moab( glc, infodata) + use seq_comm_mct, only : mbgxid + type(component_type) , intent(in) :: glc ! component type for instance1 type(seq_infodata_type) , intent(in) :: infodata !EOP !----- local ----- - !TODO change this to seq_comm_mct when GLC is ported - integer(in) :: mbgxid integer(in) :: n,ic,nf,ip ! generic index integer(in) :: lSize ! size of mesh real(r8) :: ca_g ! area of a grid cell diff --git a/driver-moab/main/seq_frac_mct.F90 b/driver-moab/main/seq_frac_mct.F90 index 3332afe45af8..249e030d0cae 100644 --- a/driver-moab/main/seq_frac_mct.F90 +++ b/driver-moab/main/seq_frac_mct.F90 @@ -172,6 +172,7 @@ module seq_frac_mct ! for tri grid, sameg_al would be false use seq_comm_mct, only : mbrxid ! iMOAB id of moab rof migrated to coupler pes + use seq_comm_mct, only : mbgxid ! iMOAB id of moab glc migrated to coupler pes use iMOAB, only : iMOAB_DefineTagStorage, iMOAB_SetDoubleTagStorage, & iMOAB_WriteMesh, iMOAB_SendElementTag, iMOAB_ReceiveElementTag, & @@ -408,11 +409,32 @@ subroutine seq_frac_init( infodata, & ! Initialize fractions on glc grid decomp, just an initial "guess", updated later - ! MOABTODO: add capabiility for MOAB if (glc_present) then call mct_aVect_init(fractions_g,rList=fraclist_g,lsize=0) call mct_aVect_zero(fractions_g) + if (mbgxid .ge. 0 ) then + arrSize = mbGetnCells(mbgxid) + tagname = trim(fraclist_g)//C_NULL_CHAR ! 'gfrac:lfrac' + tagtype = 1 ! dense, double + numco = 1 ! + ierr = iMOAB_DefineTagStorage(mbgxid, tagname, tagtype, numco, tagindex ) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in defining fraction tags on glc mesh on cpl ' + call shr_sys_abort(subname//' ERROR in defining fraction tags on glc mesh on cpl') + endif + + ! zero the fractions on glc + allocate(tagValues(arrSize*2) ) + tagValues = 0.0_r8 + call mbSetCellTagVals(mbgxid, fraclist_g,tagValues,arrSize*2) + deallocate(tagValues) + ! set gfrac to the domain frac + allocate(tagValues(arrSize)) + call mbGetCellTagVals(mbgxid, 'frac',tagValues,arrSize) + call mbSetCellTagVals(mbgxid, 'gfrac',tagValues,arrSize) + deallocate(tagValues) + endif end if ! Initialize fractions on land grid decomp, just an initial "guess", updated later @@ -765,7 +787,7 @@ subroutine seq_frac_init( infodata, & if (lnd_present) call seq_frac_check(fractions_l,mblxid,'lnd init') - if (glc_present) call seq_frac_check(fractions_g,-1,'glc init') + if (glc_present) call seq_frac_check(fractions_g,mbgxid,'glc init') if (rof_present) call seq_frac_check(fractions_r,mbrxid,'rof init') if (wav_present) call seq_frac_check(fractions_w,-1,'wav init') if (iac_present) call seq_frac_check(fractions_z,-1,'iac init') diff --git a/driver-moab/main/seq_hist_mod.F90 b/driver-moab/main/seq_hist_mod.F90 index 8e8ce4959afd..25c3d97cc482 100644 --- a/driver-moab/main/seq_hist_mod.F90 +++ b/driver-moab/main/seq_hist_mod.F90 @@ -42,11 +42,12 @@ module seq_hist_mod use prep_aoflux_mod, only: prep_aoflux_get_xao_ox use prep_aoflux_mod, only: prep_aoflux_get_xao_ax - use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid + use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid, mbgxid use seq_flds_mod, only: seq_flds_a2x_fields, seq_flds_xao_fields, seq_flds_o2x_fields, seq_flds_x2o_fields 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 seq_flds_mod, only: seq_flds_g2x_fields, seq_flds_x2g_fields use shr_moab_mod, only: mbGetnCells,mbGetCellTagVals,mbSetCellTagVals use component_type_mod @@ -390,17 +391,16 @@ subroutine seq_hist_write(infodata, EClock_d, & deallocate(mask) endif -! ! Ready to uncomment once mbgxid exists -! if (glc_present) then -! call seq_io_write(hist_file, mbgxid, 'domg', & -! trim(seq_flds_dom_fields), whead=whead, wdata=wdata, dims2do=latlonid) -! call seq_io_write(hist_file, mbgxid, 'fracg', & -! 'gfrac:lfrac', whead=whead, wdata=wdata, dims2din=latlonid) -! call seq_io_write(hist_file, mbgxid, 'g2x', & -! trim(seq_flds_g2x_fields), nt=1, whead=whead, wdata=wdata, dims2din=latlonid) -! call seq_io_write(hist_file, mbgxid, 'x2g', & -! trim(seq_flds_x2g_fields), nt=1, whead=whead, wdata=wdata, dims2din=latlonid) -! endif + if (glc_present) then + call seq_io_write(hist_file, mbgxid, 'domg', & + trim(seq_flds_dom_fields), whead=whead, wdata=wdata, nx=glc_nx, ny=glc_ny, nt=1, dims2do=latlonid) + call seq_io_write(hist_file, mbgxid, 'fracg', & + 'gfrac:lfrac', whead=whead, wdata=wdata, nx=glc_nx, ny=glc_ny, nt=1, dims2din=latlonid) + call seq_io_write(hist_file, mbgxid, 'g2x', & + trim(seq_flds_g2x_fields), whead=whead, wdata=wdata, nx=glc_nx, ny=glc_ny, nt=1, dims2din=latlonid) + call seq_io_write(hist_file, mbgxid, 'x2g', & + trim(seq_flds_x2g_fields), whead=whead, wdata=wdata, nx=glc_nx, ny=glc_ny, nt=1, dims2din=latlonid) + endif !MOAB TODO: convert this ! if (wav_present) then From 42394724d6ace3a434385e9906eeabb46c46a935 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Fri, 21 Aug 2026 14:01:26 -0500 Subject: [PATCH 5/7] Add two-way lnd-glc coupling with elevation classes to driver-moab Refresh prep_glc_mod from driver-mct (prep_glc_mrg_ocn, the thermal forcing/ice shelf split, l2gacc averaging counters, the lfrin fix) and add the MOAB implementation of the lnd->glc path: accumulate the per-elevation-class land fields in an array accumulator following the prep_rof pattern, average and set them back into the land coupler tags (which also keeps the l2x1yrg auxiliary history working), map all fields in one batched lfrac-normalized seq_map_map, downscale with array-based vertical interpolation between elevation classes, pre-multiply qice by area/aream to preserve conservation through the component-side area correction, and renormalize SMB so the global integral over the ice sheet matches the land integral. Results are written directly into the x2g tag names so the merge is a no-op. Add map_glc2lnd_ec_moab for the glc->lnd direction: because each elevation class needs a different weight field, build pre-multiplied numerator tags (frac_n*icemask, topo*w_n, icemask) on the glc mesh, map them all raw in one call, and divide on the land mesh, filling virtual columns with the class mean elevation where a class has no ice. Wire it into prep_lnd, which registers the GLC_LND_COU map and reads glc2lnd_fmapname through iMOAB. Read the lnd2glc and glc2lnd map files with moab_map_init_rcfile (weight ids scalar_l2g, flux_l2g, flux_g2l), taking aream on the glc coupler mesh from the lnd2glc map area_b as driver-mct does. Replace the inherited seq_map_init_rcfile calls for the glc mappers with seq_map_mapinit, since driver-moab never populates the coupler-side MCT gsmaps those need (they segfaulted in MCT sparse-matrix init). Restart the lnd->glc accumulator, the glc fractions, and g2x through the coupler restart file, and wire the accumulate/average/map/merge, component exchanges, area correction, zeroing, and glc budget call into the driver run loop. Validated on anlgce-ub22 (gnu): ERS.ne30pg2_r05_IcoswISC30E3r5_gis4to40.MALISIA and ERS_Ld5.ne30pg2_r05_IcoswISC30E3r5_gis4to40.IGELM_MLI (with MALI_USE_ALBANY=FALSE) both pass with bit-for-bit restarts; IGELM SMB renormalization factors are 1.006-1.018 and MALI globalStats shows the coupler SMB arriving (avgNetAccumulation ~0.33 m/yr). The ocn->glc thermal forcing and ice shelf paths, glc->ocn/ice runoff mapping, and the x2gacc restart are deferred to a follow-up. Co-Authored-By: Claude Fable 5 --- driver-moab/main/cime_comp_mod.F90 | 55 +- driver-moab/main/map_glc2lnd_mod.F90 | 256 +++++++ driver-moab/main/map_lnd2glc_mod.F90 | 93 ++- driver-moab/main/prep_glc_mod.F90 | 1056 ++++++++++++++++++++++++-- driver-moab/main/prep_lnd_mod.F90 | 110 ++- driver-moab/main/seq_rest_mod.F90 | 120 ++- 6 files changed, 1547 insertions(+), 143 deletions(-) diff --git a/driver-moab/main/cime_comp_mod.F90 b/driver-moab/main/cime_comp_mod.F90 index 594a2f5456da..4274634489bb 100644 --- a/driver-moab/main/cime_comp_mod.F90 +++ b/driver-moab/main/cime_comp_mod.F90 @@ -156,6 +156,7 @@ module cime_comp_mod use seq_diag_moab, only : seq_diag_zero_moab, seq_diag_lnd_moab use seq_diag_moab, only : seq_diag_rof_moab , seq_diag_ocn_moab, seq_diag_atm_moab use seq_diag_moab, only : seq_diag_ice_moab , seq_diag_accum_moab, seq_diag_print_moab + use seq_diag_moab, only : seq_diag_glc_moab use seq_diagBGC_moab, only : seq_diagBGC_zero_moab, seq_diagBGC_lnd_moab use seq_diagBGC_moab, only : seq_diagBGC_rof_moab , seq_diagBGC_ocn_moab, seq_diagBGC_atm_moab use seq_diagBGC_moab, only : seq_diagBGC_ice_moab , seq_diagBGC_accum_moab @@ -1500,7 +1501,9 @@ end subroutine cime_pre_init2 subroutine cime_init() use seq_flds_mod , only : seq_flds_x2a_fields, seq_flds_a2x_fields, seq_flds_l2x_fields, & seq_flds_o2x_fields, seq_flds_r2x_fields, seq_flds_i2x_fields + use seq_flds_mod , only : seq_flds_g2x_fluxes, seq_flds_g2x_fields use seq_comm_mct , only : mphaid, mbaxid, mlnid, mblxid, mrofid, mbrxid, mpoid, mboxid, mpsiid, mbixid + use seq_comm_mct , only : mbglid, mbgxid use seq_comm_mct, only: num_moab_exports ! used to count the steps for moab files integer :: nfields, numpts @@ -2144,7 +2147,7 @@ subroutine cime_init() ! MOABTODO: l2r intx, a2r intx call prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) - call prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) + call prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glctf, ocn_c2_glcshelf) call prep_wav_init(infodata, atm_c2_wav, ocn_c2_wav, ice_c2_wav) @@ -2253,9 +2256,8 @@ subroutine cime_init() call mpi_barrier(mpicom_GLOID,ierr) if (ice_present) call component_init_areacor_moab(ice, areafact_samegrid, mpsiid, mbixid, seq_flds_i2x_fluxes, seq_flds_i2x_fields) - ! MOAB TODO: convert these to moab call mpi_barrier(mpicom_GLOID,ierr) - !if (glc_present) call component_init_areacor(glc, areafact_samegrid, seq_flds_g2x_fluxes) + if (glc_present) call component_init_areacor_moab(glc, areafact_samegrid, mbglid, mbgxid, seq_flds_g2x_fluxes, seq_flds_g2x_fields) call mpi_barrier(mpicom_GLOID,ierr) !if (wav_present) call component_init_areacor(wav, areafact_samegrid, seq_flds_w2x_fluxes) @@ -3310,8 +3312,12 @@ subroutine cime_run() !---------------------------------------------------------- !| GLC SETUP-SEND !---------------------------------------------------------- - if (glc_present .and. glcrun_alarm) then - call cime_run_glc_setup_send(lnd2glc_averaged_now, prep_glc_accum_avg_called) + if (glc_present) then + if (glcrun_alarm) then + call cime_run_glc_setup_send(lnd2glc_averaged_now, prep_glc_accum_avg_called) + else + if (iamin_CPLID) call prep_glc_zero_fields_moab() + endif endif ! ------------------------------------------------------------------------ @@ -4467,12 +4473,16 @@ subroutine cime_run_ocnglc_coupling() if (glc_present) then + ! create o2x_gx for either ocn-glc coupling or ocn-glc shelf coupling + if (ocn_c2_glctf .or. (ocn_c2_glcshelf .and. glcshelf_c2_ocn)) then + call prep_glc_calc_o2x_gx(ocn_c2_glctf, ocn_c2_glcshelf, timer='CPL:glcprep_ocn2glc') !remap ocean fields to o2x_g at ocean couping interval + endif + + ! if ice-shelf coupling is on, now proceed to handle those calculations here in the coupler if (ocn_c2_glcshelf .and. glcshelf_c2_ocn) then ! the boundary flux calculations done in the coupler require inputs from both GLC and OCN, ! so they will only be valid if both OCN->GLC and GLC->OCN - call prep_glc_calc_o2x_gx(timer='CPL:glcprep_ocn2glc') !remap ocean fields to o2x_g at ocean couping interval - call prep_glc_calculate_subshelf_boundary_fluxes ! this is actual boundary layer flux calculation !this outputs !x2g_g/g2x_g, where latter is going @@ -4576,7 +4586,7 @@ subroutine cime_run_lnd_recv_post() ! Accumulate rof and glc inputs (module variables in prep_rof_mod and prep_glc_mod) if (lnd_c2_rof) call prep_rof_accum_lnd_moab(timer='CPL:lndpost_accl2r') - if (lnd_c2_glc .or. do_hist_l2x1yrg) call prep_glc_accum_lnd(timer='CPL:lndpost_accl2g' ) + if (lnd_c2_glc .or. do_hist_l2x1yrg) call prep_glc_accum_lnd_moab(timer='CPL:lndpost_accl2g' ) if (lnd_c2_iac) call prep_iac_accum(timer='CPL:lndpost_accl2z') if (drv_threading) call seq_comm_setnthreads(nthreads_GLOID) @@ -4589,6 +4599,9 @@ end subroutine cime_run_lnd_recv_post subroutine cime_run_glc_setup_send(lnd2glc_averaged_now, prep_glc_accum_avg_called) + use seq_flds_mod , only : seq_flds_x2g_fields + use seq_comm_mct , only : mbglid, mbgxid + logical, intent(inout) :: lnd2glc_averaged_now ! Set to .true. if lnd2glc averages are taken this timestep (otherwise left unchanged) logical, intent(inout) :: prep_glc_accum_avg_called ! Set to .true. if prep_glc_accum_avg is called here (otherwise left unchanged) @@ -4603,22 +4616,23 @@ subroutine cime_run_glc_setup_send(lnd2glc_averaged_now, prep_glc_accum_avg_call ! NOTE - only create appropriate input to glc if the avg_alarm is on if (lnd_c2_glc .or. ocn_c2_glcshelf) then if (glcrun_avg_alarm) then - call prep_glc_accum_avg(timer='CPL:glcprep_avg', & + call prep_glc_accum_avg_moab(timer='CPL:glcprep_avg', & lnd2glc_averaged_now=lnd2glc_averaged_now) prep_glc_accum_avg_called = .true. if (lnd_c2_glc) then - ! Note that l2x_gx is obtained from mapping the module variable l2gacc_lx - call prep_glc_calc_l2x_gx(fractions_lx, timer='CPL:glcprep_lnd2glc') + ! Note that the mapped fields are obtained from the accumulated lnd + ! fields (set back into the land mesh tags by prep_glc_accum_avg_moab) + call prep_glc_calc_l2x_gx_moab(fractions_lx, timer='CPL:glcprep_lnd2glc') - call prep_glc_mrg_lnd(infodata, fractions_gx, timer_mrg='CPL:glcprep_mrgx2g') + call prep_glc_mrg_lnd_moab(infodata, timer_mrg='CPL:glcprep_mrgx2g') endif call component_diag(infodata, glc, flow='x2c', comment='send glc', & info_debug=info_debug, timer_diag='CPL:glcprep_diagav') else - call prep_glc_zero_fields() + call prep_glc_zero_fields_moab() endif ! glcrun_avg_alarm end if ! lnd_c2_glc or ocn_c2_glcshelf @@ -4640,7 +4654,7 @@ subroutine cime_run_glc_setup_send(lnd2glc_averaged_now, prep_glc_accum_avg_call !| cpl -> glc !---------------------------------------------------- if (iamin_CPLALLGLCID .and. glc_prognostic) then - call component_exch(glc, flow='x2c', & + call component_exch_moab(glc(1), mbgxid, mbglid, 'x2c', seq_flds_x2g_fields, & infodata=infodata, infodata_string='cpl2glc_run', & mpicom_barrier=mpicom_CPLALLGLCID, run_barriers=run_barriers, & timer_barrier='CPL:C2G_BARRIER', timer_comp_exch='CPL:C2G', & @@ -4661,7 +4675,7 @@ subroutine cime_run_glc_accum_avg(lnd2glc_averaged_now, prep_glc_accum_avg_calle call t_drvstartf ('CPL:AVG_L2X1YRG',cplrun=.true.,barrier=mpicom_CPLID) if (drv_threading) call seq_comm_setnthreads(nthreads_CPLID) - call prep_glc_accum_avg(timer='CPL:glcprep_avg', & + call prep_glc_accum_avg_moab(timer='CPL:glcprep_avg', & lnd2glc_averaged_now=lnd2glc_averaged_now) prep_glc_accum_avg_called = .true. @@ -4673,11 +4687,15 @@ end subroutine cime_run_glc_accum_avg subroutine cime_run_glc_recv_post() + use seq_flds_mod , only : seq_flds_g2x_fields + use seq_comm_mct , only : mbglid, mbgxid + !---------------------------------------------------------- ! glc -> cpl !---------------------------------------------------------- if (iamin_CPLALLGLCID) then - call component_exch(glc, flow='c2x', infodata=infodata, infodata_string='glc2cpl_run', & + call component_exch_moab(glc(1), mbglid, mbgxid, 'c2x', seq_flds_g2x_fields, & + infodata=infodata, infodata_string='glc2cpl_run', & mpicom_barrier=mpicom_CPLALLGLCID, run_barriers=run_barriers, & timer_barrier='CPL:G2C_BARRIER', timer_comp_exch='CPL:G2C', & timer_map_exch='CPL:g2c_glcg2glcx', timer_infodata_exch='CPL:g2c_infoexch') @@ -5044,6 +5062,11 @@ subroutine cime_run_calc_budgets2(in_cplrun) if (ice_present) then call seq_diag_ice_moab(ice(ens1), infodata, do_i2x=.true.) endif + if (glc_present) then + ! covers the g2x (runoff to ocn/ice) budget terms; the x2g budget terms of + ! the mct driver are deferred with the rest of the ocn<->glc port + call seq_diag_glc_moab(glc(ens1), infodata) + endif if (do_bgc_budgets) then if (atm_present) then call seq_diagBGC_atm_moab(atm(ens1), infodata, do_a2x=.true., do_x2a=.true.) diff --git a/driver-moab/main/map_glc2lnd_mod.F90 b/driver-moab/main/map_glc2lnd_mod.F90 index 4ae638ed71cf..7072ea1335b8 100644 --- a/driver-moab/main/map_glc2lnd_mod.F90 +++ b/driver-moab/main/map_glc2lnd_mod.F90 @@ -13,6 +13,7 @@ module map_glc2lnd_mod #include "shr_assert.h" use seq_comm_mct, only : logunit use shr_kind_mod, only : r8 => shr_kind_r8 + use shr_kind_mod, only : cxx => SHR_KIND_CXX use glc_elevclass_mod, only : glc_get_num_elevation_classes, glc_get_elevation_class, & glc_mean_elevation_virtual, glc_elevclass_as_string, & GLC_ELEVCLASS_ERR_NONE, GLC_ELEVCLASS_ERR_TOO_LOW, & @@ -33,6 +34,7 @@ module map_glc2lnd_mod !-------------------------------------------------------------------------- public :: map_glc2lnd_ec ! map all fields from GLC -> LND grid that need to be separated by elevation class + public :: map_glc2lnd_ec_moab ! moab (tag-based) version of map_glc2lnd_ec !-------------------------------------------------------------------------- ! Private interfaces @@ -220,6 +222,260 @@ subroutine map_glc2lnd_ec(g2x_g, & end subroutine map_glc2lnd_ec + !----------------------------------------------------------------------- + subroutine map_glc2lnd_ec_moab(mapper, & + frac_field, topo_field, icemask_field, extra_fields, frac_l_out) + ! + ! !DESCRIPTION: + ! moab version of map_glc2lnd_ec: maps fields from the GLC mesh (mapper%src_mbid) + ! to the LND mesh (mapper%tgt_mbid), separated by elevation class, operating on + ! moab tags. The per-EC results are written into the tags named + ! on the land mesh (these are members of seq_flds_x2l_fields, already + ! defined there). + ! + ! The mct version does, for each elevation class n, two normalized maps: + ! frac_n_l = M(icemask*frac_n) / M(icemask) + ! field_n_l = M(icemask*frac_n*field) / M(icemask*frac_n) + ! Because the same conservative map M is used throughout, this is reproduced here + ! with ONE raw (norm=.false.) map of all pre-multiplied numerators, followed by + ! explicit divisions on the land side (the raw row sums cancel in the ratios). + ! The numerators are staged in the per-EC destination tag names on the glc mesh + ! (defined at init), plus one scratch tag 'Sg_icemsk_num' for M(icemask). + ! + ! Assumes the following tags exist: + ! - on the glc mesh: frac_field, topo_field, icemask_field, each extra field, + ! the per-EC names for frac/topo/extras, and Sg_icemsk_num + ! - on the lnd mesh: the same per-EC names and Sg_icemsk_num + ! + ! !USES: + use iMOAB, only : iMOAB_GetMeshInfo, iMOAB_GetDoubleTagStorage, iMOAB_SetDoubleTagStorage + use iso_c_binding, only : C_NULL_CHAR + use glc_elevclass_mod, only : glc_all_elevclass_strings, GLC_ELEVCLASS_STRLEN + ! + ! !ARGUMENTS: + type(seq_map), intent(inout) :: mapper ! conservative glc->lnd mapper with moab context + character(len=*), intent(in) :: frac_field ! name of glc field containing glc ice fraction + character(len=*), intent(in) :: topo_field ! name of glc field containing glc topo + character(len=*), intent(in) :: icemask_field ! name of glc field containing ice mask + character(len=*), intent(in) :: extra_fields ! colon-delimited additional fields ('' or ' ' for none) + real(r8), optional, intent(out) :: frac_l_out(:,0:) ! normalized per-EC fractions on the land mesh + ! + ! !LOCAL VARIABLES: + integer :: lsize_g, lsize_l + integer :: nEC, nextra, nnum + integer :: n, i, k, kfrac0, ktopo0, kextra0, kmask + integer :: ierr, ent_type, arrsize + integer :: nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) + + real(r8), allocatable :: glc_frac(:) ! total ice fraction in each glc cell + real(r8), allocatable :: glc_topo(:) ! topographic height of each glc cell + real(r8), allocatable :: glc_icemask(:) ! icemask of each glc cell + real(r8), allocatable :: glc_extra(:,:) ! extra fields on the glc mesh + real(r8), allocatable :: frac_this_ec(:) ! ice fraction in one elevation class + integer , allocatable :: glc_elevclass(:) + real(r8), allocatable :: num_g(:,:) ! numerators on the glc mesh + real(r8), allocatable :: num_l(:,:) ! mapped numerators on the land mesh + real(r8), allocatable :: out_l(:,:) ! normalized per-EC fields on the land mesh + real(r8) :: denom, topo_virtual + + character(len=GLC_ELEVCLASS_STRLEN), allocatable :: ec_strings(:) + character(CXX) :: numlist ! colon-separated list of all numerator tag names + character(CXX) :: tagname + type(mct_list) :: temp_list + type(mct_string) :: mctOStr + character(CXX) :: extra_names(20) ! names of the extra fields (assumed few) + + ! dummy zero-size attribute vectors to satisfy the seq_map_map interface + type(mct_aVect) :: av_dum_s, av_dum_d + + character(len=*), parameter :: scratch_icemask_num = 'Sg_icemsk_num' + character(len=*), parameter :: subname = 'map_glc2lnd_ec_moab' + !----------------------------------------------------------------------- + + if (mapper%src_mbid < 0 .or. mapper%tgt_mbid < 0) return + + ent_type = 1 ! cells on both meshes + + ierr = iMOAB_GetMeshInfo ( mapper%src_mbid, nvert, nvise, nbl, nsurf, nvisBC ) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting glc mesh info') + lsize_g = nvise(1) + ierr = iMOAB_GetMeshInfo ( mapper%tgt_mbid, nvert, nvise, nbl, nsurf, nvisBC ) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting lnd mesh info') + lsize_l = nvise(1) + + nEC = glc_get_num_elevation_classes() + + ! parse the extra field names + nextra = 0 + if (len_trim(extra_fields) > 0) then + call mct_list_init(temp_list, extra_fields) + nextra = mct_list_nitem(temp_list) + if (nextra > size(extra_names)) call shr_sys_abort(subname//' ERROR too many extra fields') + do i = 1, nextra + call mct_list_get(mctOStr, i, temp_list) + extra_names(i) = mct_string_toChar(mctOStr) + call mct_string_clean(mctOStr) + end do + call mct_list_clean(temp_list) + end if + + ! ------------------------------------------------------------------------ + ! Extract needed fields from the glc mesh tags + ! ------------------------------------------------------------------------ + + allocate(glc_frac(lsize_g), glc_topo(lsize_g), glc_icemask(lsize_g)) + allocate(frac_this_ec(lsize_g), glc_elevclass(lsize_g)) + tagname = trim(frac_field)//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mapper%src_mbid, tagname, lsize_g, ent_type, glc_frac) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting '//trim(frac_field)) + tagname = trim(topo_field)//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mapper%src_mbid, tagname, lsize_g, ent_type, glc_topo) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting '//trim(topo_field)) + tagname = trim(icemask_field)//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mapper%src_mbid, tagname, lsize_g, ent_type, glc_icemask) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting '//trim(icemask_field)) + if (nextra > 0) then + allocate(glc_extra(lsize_g, nextra)) + do i = 1, nextra + tagname = trim(extra_names(i))//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mapper%src_mbid, tagname, lsize_g, ent_type, glc_extra(:,i)) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting '//trim(extra_names(i))) + end do + end if + + call get_glc_elevation_classes(glc_topo, glc_elevclass) + + ! ------------------------------------------------------------------------ + ! Build the numerator fields and the combined tag list + ! layout per elevation class n (0..nEC): + ! NN = frac_n * icemask (this is M's input w_n) + ! NN = topo * w_n + ! NN = extra * w_n + ! plus one extra column: Sg_icemsk_num = icemask + ! ------------------------------------------------------------------------ + + nnum = (nEC+1)*(2+nextra) + 1 + allocate(num_g(lsize_g, nnum)) + allocate(ec_strings(0:nEC)) + ec_strings = glc_all_elevclass_strings(include_zero = .true.) + + numlist = '' + k = 0 + kfrac0 = 1 ! columns kfrac0 + n*(2+nextra) hold w_n, etc. + do n = 0, nEC + call get_frac_this_ec(glc_frac, glc_elevclass, n, frac_this_ec) + ! frac numerator: w_n = frac_n * icemask + k = k + 1 + num_g(:,k) = frac_this_ec(:) * glc_icemask(:) + call add_to_list(numlist, trim(frac_field)//trim(ec_strings(n))) + ! topo numerator: topo * w_n + k = k + 1 + num_g(:,k) = glc_topo(:) * num_g(:,k-1) + call add_to_list(numlist, trim(topo_field)//trim(ec_strings(n))) + ! extra numerators + do i = 1, nextra + k = k + 1 + num_g(:,k) = glc_extra(:,i) * num_g(:,k-1-i) + call add_to_list(numlist, trim(extra_names(i))//trim(ec_strings(n))) + end do + end do + ! icemask numerator (for the frac normalization) + k = k + 1 + kmask = k + num_g(:,k) = glc_icemask(:) + call add_to_list(numlist, scratch_icemask_num) + + ! set the numerators on the glc mesh tags + arrsize = nnum * lsize_g + tagname = trim(numlist)//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage(mapper%src_mbid, tagname, arrsize, ent_type, num_g) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR setting numerator tags on glc mesh') + + ! ------------------------------------------------------------------------ + ! One raw map of all numerators glc -> lnd + ! ------------------------------------------------------------------------ + + call mct_aVect_init(av_dum_s, rList = frac_field, lsize = 0) + call mct_aVect_init(av_dum_d, rList = frac_field, lsize = 0) + call seq_map_map(mapper, av_dum_s, av_dum_d, fldlist=trim(numlist), norm=.false.) + call mct_aVect_clean(av_dum_s) + call mct_aVect_clean(av_dum_d) + + ! ------------------------------------------------------------------------ + ! Get mapped numerators on the land mesh and normalize + ! ------------------------------------------------------------------------ + + allocate(num_l(lsize_l, nnum)) + allocate(out_l(lsize_l, nnum-1)) + num_l = 0.0_r8 + arrsize = nnum * lsize_l + ierr = iMOAB_GetDoubleTagStorage(mapper%tgt_mbid, tagname, arrsize, ent_type, num_l) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting numerator tags on lnd mesh') + + do n = 0, nEC + kfrac0 = 1 + n*(2+nextra) ! column of w_n + ktopo0 = kfrac0 + 1 + topo_virtual = glc_mean_elevation_virtual(n) + do i = 1, lsize_l + ! frac_n_l = M(w_n) / M(icemask) + denom = num_l(i, kmask) + if (denom /= 0.0_r8) then + out_l(i, kfrac0) = num_l(i, kfrac0) / denom + else + out_l(i, kfrac0) = 0.0_r8 + end if + ! field_n_l = M(field*w_n) / M(w_n) + denom = num_l(i, kfrac0) + if (denom /= 0.0_r8) then + out_l(i, ktopo0) = num_l(i, ktopo0) / denom + do k = 1, nextra + out_l(i, ktopo0+k) = num_l(i, ktopo0+k) / denom + end do + else + out_l(i, ktopo0) = 0.0_r8 + do k = 1, nextra + out_l(i, ktopo0+k) = 0.0_r8 + end do + end if + ! set the topo field for virtual columns (no contributing glc cells) + if (out_l(i, kfrac0) <= 0.0_r8) then + out_l(i, ktopo0) = topo_virtual + end if + end do + if (present(frac_l_out)) then + frac_l_out(:, n) = out_l(:, kfrac0) + end if + end do + + ! ------------------------------------------------------------------------ + ! Store the normalized per-EC fields back into the land mesh tags + ! (all list entries except the trailing Sg_icemsk_num scratch) + ! ------------------------------------------------------------------------ + + i = len_trim(numlist) - len(scratch_icemask_num) - 1 ! strip ':Sg_icemsk_num' + tagname = numlist(1:i)//C_NULL_CHAR + arrsize = (nnum-1) * lsize_l + ierr = iMOAB_SetDoubleTagStorage(mapper%tgt_mbid, tagname, arrsize, ent_type, out_l) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR setting per-EC tags on lnd mesh') + + deallocate(glc_frac, glc_topo, glc_icemask, frac_this_ec, glc_elevclass) + if (allocated(glc_extra)) deallocate(glc_extra) + deallocate(num_g, num_l, out_l, ec_strings) + + contains + + subroutine add_to_list(list, name) + ! append name to a colon-separated list + character(len=*), intent(inout) :: list + character(len=*), intent(in) :: name + if (len_trim(list) == 0) then + list = trim(name) + else + list = trim(list)//':'//trim(name) + end if + end subroutine add_to_list + + end subroutine map_glc2lnd_ec_moab !----------------------------------------------------------------------- subroutine get_glc_elevation_classes(glc_topo, glc_elevclass) diff --git a/driver-moab/main/map_lnd2glc_mod.F90 b/driver-moab/main/map_lnd2glc_mod.F90 index 1b418515e912..60bf7b9efe71 100644 --- a/driver-moab/main/map_lnd2glc_mod.F90 +++ b/driver-moab/main/map_lnd2glc_mod.F90 @@ -32,12 +32,17 @@ module map_lnd2glc_mod !-------------------------------------------------------------------------- public :: map_lnd2glc ! map one field from LND -> GLC grid + ! array-based pieces of the algorithm, used by the moab path in prep_glc_mod + ! (the horizontal map is done there with one batched seq_map_map call on tags; + ! these routines provide the elevation-class assignment and the vertical + ! interpolation on plain arrays fetched from the glc mesh tags) + public :: get_glc_elevation_classes ! get the elevation class of each point on the glc grid + public :: map_lnd2glc_vertical_interp ! vertically interpolate per-EC data to the ice sheet topography !-------------------------------------------------------------------------- ! Private interfaces !-------------------------------------------------------------------------- - private :: get_glc_elevation_classes ! get the elevation class of each point on the glc grid private :: map_bare_land ! remap the field of interest for the bare land "elevation class" private :: map_ice_covered ! remap the field of interest for all elevation classes (excluding bare land) @@ -192,6 +197,92 @@ subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, & end subroutine map_lnd2glc + !----------------------------------------------------------------------- + subroutine map_lnd2glc_vertical_interp(topo_g, topo_g_EC, data_g_EC, data_g_bareland, & + glc_elevclass, data_g) + ! + ! !DESCRIPTION: + ! Vertically interpolate a field, already horizontally mapped to the glc grid in + ! each elevation class, onto the actual ice sheet topography. This is the + ! array-based equivalent of the combination of map_bare_land / map_ice_covered + ! (minus the horizontal maps, which the caller has already done): the output is + ! the bare-land (EC 0) value where the glc cell is ice-free, and the linear + ! vertical interpolation between bounding elevation classes where ice-covered. + ! + ! All arrays are on the glc grid decomposition. data_g_EC and topo_g_EC hold the + ! horizontally-mapped field and Sl_topo for elevation classes 1..nEC. + ! + ! Note: the mct path (map_ice_covered) stores the per-EC arrays in default real + ! precision; here everything stays in r8, a roundoff-level difference. + ! + ! !ARGUMENTS: + real(r8), intent(in) :: topo_g(:) ! ice topographic height on the glc grid + real(r8), intent(in) :: topo_g_EC(:,:) ! mapped per-EC topo (lsize_g, nEC) + real(r8), intent(in) :: data_g_EC(:,:) ! mapped per-EC field (lsize_g, nEC) + real(r8), intent(in) :: data_g_bareland(:) ! mapped bare-land (EC 0) field + integer , intent(in) :: glc_elevclass(:) ! elevation class of each glc point (0 = bare) + real(r8), intent(out) :: data_g(:) ! result on the glc grid + ! + ! !LOCAL VARIABLES: + integer :: lsize_g ! number of cells on glc grid + integer :: nEC ! number of elevation classes + integer :: n, ec + real(r8) :: elev_l, elev_u ! lower and upper elevations in interpolation range + real(r8) :: d_elev ! elev_u - elev_l + + character(len=*), parameter :: subname = 'map_lnd2glc_vertical_interp' + !----------------------------------------------------------------------- + + lsize_g = size(data_g) + nEC = size(data_g_EC, 2) + SHR_ASSERT_FL((size(topo_g) == lsize_g), __FILE__, __LINE__) + SHR_ASSERT_FL((size(glc_elevclass) == lsize_g), __FILE__, __LINE__) + + do n = 1, lsize_g + + if (glc_elevclass(n) == 0) then + ! bare land (ice-free) point: use the bare-land value + data_g(n) = data_g_bareland(n) + + ! For each ice sheet point, find bounding EC values... + else if (topo_g(n) < topo_g_EC(n,1)) then + ! lower than lowest mean EC elevation value + data_g(n) = data_g_EC(n,1) + + else if (topo_g(n) >= topo_g_EC(n,nEC)) then + ! higher than highest mean EC elevation value + data_g(n) = data_g_EC(n,nEC) + + else + ! do linear interpolation of data in the vertical + do ec = 2, nEC + if (topo_g(n) < topo_g_EC(n, ec)) then + elev_l = topo_g_EC(n, ec-1) + elev_u = topo_g_EC(n, ec) + d_elev = elev_u - elev_l + if (d_elev <= 0) then + ! This shouldn't happen, but handle it in case it does. In this case, + ! let's arbitrarily use the mean of the two elevation classes, rather + ! than the weighted mean. + write(logunit,*) subname//' WARNING: topo diff between elevation classes <= 0' + write(logunit,*) 'n, ec, elev_l, elev_u = ', n, ec, elev_l, elev_u + write(logunit,*) 'Simply using mean of the two elevation classes,' + write(logunit,*) 'rather than the weighted mean.' + data_g(n) = data_g_EC(n,ec-1) * 0.5_r8 & + + data_g_EC(n,ec) * 0.5_r8 + else + data_g(n) = data_g_EC(n,ec-1) * (elev_u - topo_g(n)) / d_elev & + + data_g_EC(n,ec) * (topo_g(n) - elev_l) / d_elev + end if + + exit + end if + end do + end if ! topo_g(n) + end do ! lsize_g + + end subroutine map_lnd2glc_vertical_interp + !----------------------------------------------------------------------- subroutine get_glc_elevation_classes(glc_ice_covered, glc_topo, glc_elevclass) ! diff --git a/driver-moab/main/prep_glc_mod.F90 b/driver-moab/main/prep_glc_mod.F90 index 339acaaaac04..8828b060e849 100644 --- a/driver-moab/main/prep_glc_mod.F90 +++ b/driver-moab/main/prep_glc_mod.F90 @@ -3,11 +3,13 @@ module prep_glc_mod #include "shr_assert.h" use shr_kind_mod , only: r8 => SHR_KIND_R8 use shr_kind_mod , only: cl => SHR_KIND_CL + use shr_kind_mod , only: CXX => SHR_KIND_CXX use shr_sys_mod , only: shr_sys_abort, shr_sys_flush use seq_comm_mct , only: num_inst_glc, num_inst_lnd, num_inst_frc, & num_inst_ocn use seq_comm_mct , only: CPLID, GLCID, logunit use seq_comm_mct , only: seq_comm_getData=>seq_comm_setptrs + use seq_comm_mct , only: mblxid, mbgxid, mbintxlg, mbintxgl ! iMOAB app ids: lnd and glc on coupler, l2g and g2l map holders use seq_infodata_mod, only: seq_infodata_type, seq_infodata_getdata use seq_map_type_mod use seq_map_mod @@ -31,20 +33,31 @@ module prep_glc_mod public :: prep_glc_init public :: prep_glc_mrg_lnd + public :: prep_glc_mrg_ocn + public :: prep_glc_mrg_lnd_moab public :: prep_glc_accum_lnd public :: prep_glc_accum_ocn public :: prep_glc_accum_avg + public :: prep_glc_accum_lnd_moab + public :: prep_glc_accum_avg_moab public :: prep_glc_calc_l2x_gx public :: prep_glc_calc_o2x_gx + public :: prep_glc_calc_l2x_gx_moab public :: prep_glc_zero_fields + public :: prep_glc_zero_fields_moab + + public :: prep_glc_get_l2gacc_lm + public :: prep_glc_get_l2gacc_lm_cnt + public :: prep_glc_get_sharedFieldsLndGlc public :: prep_glc_get_l2x_gx public :: prep_glc_get_l2gacc_lx public :: prep_glc_get_l2gacc_lx_one_instance public :: prep_glc_get_l2gacc_lx_cnt + public :: prep_glc_get_l2gacc_lx_cnt_avg public :: prep_glc_get_o2x_gx public :: prep_glc_get_x2gacc_gx @@ -53,8 +66,8 @@ module prep_glc_mod public :: prep_glc_get_mapper_Sl2g public :: prep_glc_get_mapper_Fl2g - public :: prep_glc_get_mapper_So2g - public :: prep_glc_get_mapper_Fo2g + public :: prep_glc_get_mapper_So2g_shelf + public :: prep_glc_get_mapper_Fo2g_shelf public :: prep_glc_calculate_subshelf_boundary_fluxes @@ -68,6 +81,7 @@ module prep_glc_mod private :: prep_glc_map_one_state_field_lnd2glc private :: prep_glc_map_qice_conservative_lnd2glc private :: prep_glc_renormalize_smb + private :: prep_glc_renormalize_smb_moab !-------------------------------------------------------------------------- ! Private data @@ -76,8 +90,9 @@ module prep_glc_mod ! mappers type(seq_map), pointer :: mapper_Sl2g type(seq_map), pointer :: mapper_Fl2g - type(seq_map), pointer :: mapper_So2g - type(seq_map), pointer :: mapper_Fo2g + type(seq_map), pointer :: mapper_So2g_shelf + type(seq_map), pointer :: mapper_Fo2g_shelf + type(seq_map), pointer :: mapper_So2g_tf type(seq_map), pointer :: mapper_Fg2l ! attribute vectors @@ -91,6 +106,7 @@ module prep_glc_mod type(mct_aVect), pointer :: l2gacc_lx(:) ! Lnd export, lnd grid, cpl pes - allocated in driver integer , target :: l2gacc_lx_cnt ! l2gacc_lx: number of time samples accumulated + integer , target :: l2gacc_lx_cnt_avg ! l2gacc_lx: number of time samples averaged ! other module variables integer :: mpicom_CPLID ! MPI cpl communicator @@ -129,22 +145,39 @@ module prep_glc_mod real(r8), allocatable :: outOceanHeatFlux(:) real(r8), allocatable :: outIceHeatFlux(:) + ! moab support: accumulation of the per-elevation-class lnd fields on the coupler + ! land mesh, and bookkeeping for the lnd->glc downscaling done on moab tags + character(CXX) :: sharedFieldsLndGlc ! = seq_flds_l2x_fields_to_glc, the tag list + real(r8), allocatable, target :: l2gacc_lm(:,:) ! accumulated lnd fields, (lsize_lm, nflds_lg) + real(r8), allocatable :: l2x_lm2(:,:) ! scratch for reading the instantaneous lnd tags + integer , target :: l2gacc_lm_cnt ! l2gacc_lm: number of time samples accumulated + integer , target :: l2gacc_lm_cnt_avg ! l2gacc_lm: number of samples in last average + integer :: lsize_lm = 0 ! number of local cells, land coupler mesh + integer :: lsize_gm = 0 ! number of local cells, glc coupler mesh + integer :: nflds_lg = 0 ! number of fields in sharedFieldsLndGlc + !================================================================================================ contains !================================================================================================ - subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) + subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glctf, ocn_c2_glcshelf) !--------------------------------------------------------------- ! Description ! Initialize module attribute vectors and mapping variables ! + use iMOAB, only : iMOAB_RegisterApplication, iMOAB_MigrateMapMesh, & + iMOAB_ComputeCommGraph, iMOAB_DefineTagStorage, iMOAB_GetMeshInfo + use iso_c_binding, only : C_NULL_CHAR + use shr_string_mod, only : shr_string_listGetNum + ! ! Arguments type (seq_infodata_type) , intent(inout) :: infodata logical , intent(in) :: lnd_c2_glc ! .true. => lnd to glc coupling on - logical , intent(in) :: ocn_c2_glcshelf ! .true. => ocn to glc coupling on + logical , intent(in) :: ocn_c2_glctf ! .true. => ocn to glc thermal forcing coupling on + logical , intent(in) :: ocn_c2_glcshelf ! .true. => ocn to glc shelf coupling on ! ! Local Variables integer :: eli, egi, eoi @@ -164,6 +197,15 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) type(mct_avect), pointer :: x2g_gx type(mct_avect), pointer :: o2x_ox + ! moab locals + integer :: ierr, idintx, type_grid, arearead + integer :: tagtype, numco, tagindex + integer :: mpigrp_CPLID ! coupler pes group + integer :: nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) + character(CL) :: appname + character(CL) :: wgtIdSl2g, wgtIdFl2g, wgtIdFg2l + character(CXX) :: tagname + character(*), parameter :: subname = '(prep_glc_init)' character(*), parameter :: F00 = "('"//subname//" : ', 4A )" !--------------------------------------------------------------- @@ -178,8 +220,9 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) allocate(mapper_Sl2g) allocate(mapper_Fl2g) - allocate(mapper_So2g) - allocate(mapper_Fo2g) + allocate(mapper_So2g_shelf) + allocate(mapper_So2g_tf) + allocate(mapper_Fo2g_shelf) allocate(mapper_Fg2l) smb_renormalize = prep_glc_do_renormalize_smb(infodata) @@ -195,6 +238,25 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) call mct_aVect_zero(l2gacc_lx(eli)) end do l2gacc_lx_cnt = 0 + l2gacc_lx_cnt_avg = 0 + + ! moab accumulator: the per-elevation-class lnd fields are accumulated in a + ! plain array read from the coupler land mesh tags (prep_rof pattern) + if (mblxid >= 0) then + sharedFieldsLndGlc = trim(seq_flds_l2x_fields_to_glc) + nflds_lg = shr_string_listGetNum(sharedFieldsLndGlc) + ierr = iMOAB_GetMeshInfo ( mblxid, nvert, nvise, nbl, nsurf, nvisBC ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR getting land coupler mesh info') + endif + lsize_lm = nvise(1) + allocate(l2gacc_lm(lsize_lm, nflds_lg)) + allocate(l2x_lm2(lsize_lm, nflds_lg)) + l2gacc_lm = 0._r8 + l2x_lm2 = 0._r8 + l2gacc_lm_cnt = 0 + l2gacc_lm_cnt_avg = 0 + endif end if if (glc_present .and. lnd_c2_glc) then @@ -216,41 +278,161 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) samegrid_lg = .true. if (trim(lnd_gnam) /= trim(glc_gnam)) samegrid_lg = .false. - if (iamroot_CPLID) then - write(logunit,*) ' ' - write(logunit,F00) 'Initializing mapper_Sl2g' - end if - call seq_map_init_rcfile(mapper_Sl2g, lnd(1), glc(1), & - 'seq_maps.rc', 'lnd2glc_smapname:', 'lnd2glc_smaptype:', samegrid_lg, & - 'mapper_Sl2g initialization', esmf_map_flag) - - if (iamroot_CPLID) then - write(logunit,*) ' ' - write(logunit,F00) 'Initializing mapper_Fl2g' - end if - call seq_map_init_rcfile(mapper_Fl2g, lnd(1), glc(1), & - 'seq_maps.rc', 'lnd2glc_fmapname:', 'lnd2glc_fmaptype:', samegrid_lg, & - 'mapper_Fl2g initialization', esmf_map_flag) - + ! the mct sparse-matrix init (seq_map_init_rcfile) is not usable in + ! driver-moab (coupler-side gsmaps are not populated); the mappers get + ! their real (moab) context below + call seq_map_mapinit(mapper_Sl2g, mpicom_CPLID) + call seq_map_mapinit(mapper_Fl2g, mpicom_CPLID) ! We need to initialize our own Fg2l mapper because in some cases (particularly ! TG compsets - dlnd forcing CISM) the system doesn't otherwise create a Fg2l ! mapper. - if (iamroot_CPLID) then - write(logunit,*) ' ' - write(logunit,F00) 'Initializing mapper_Fg2l' - end if - call seq_map_init_rcfile(mapper_Fg2l, glc(1), lnd(1), & - 'seq_maps.rc', 'glc2lnd_fmapname:', 'glc2lnd_fmaptype:', samegrid_lg, & - 'mapper_Fg2l initialization', esmf_map_flag) + call seq_map_mapinit(mapper_Fg2l, mpicom_CPLID) call prep_glc_set_g2x_lx_fields() + + ! now give the mappers moab context and load the mapping weights with iMOAB + if ((mblxid >= 0) .and. (mbgxid >= 0)) then + + if (samegrid_lg) then + call shr_sys_abort(subname// & + ' ERROR: moab lnd->glc coupling requires distinct lnd and glc grids with map files') + endif + + call seq_comm_getData(CPLID, mpigrp=mpigrp_CPLID) + type_grid = 3 ! FV-FV for both lnd and glc coupler meshes + + if (iamroot_CPLID) then + write(logunit,*) ' ' + write(logunit,F00) 'Initializing MOAB mapper_Sl2g and mapper_Fl2g' + end if + appname = "LND_GLC_COU"//C_NULL_CHAR + ! unique external id for the moab app holding the lnd->glc read map + idintx = 100*lnd(1)%cplcompid + glc(1)%cplcompid + ierr = iMOAB_RegisterApplication(trim(appname), mpicom_CPLID, idintx, mbintxlg) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in registering lnd glc map app' + call shr_sys_abort(subname//' ERROR in registering lnd glc map app') + endif + + ! scalar (bilinear) map, used for all the lnd->glc downscaling maps + mapper_Sl2g%src_mbid = mblxid + mapper_Sl2g%tgt_mbid = mbgxid + mapper_Sl2g%intx_mbid = mbintxlg + mapper_Sl2g%src_context = lnd(1)%cplcompid + mapper_Sl2g%intx_context = idintx + wgtIdSl2g = 'scalar_l2g' + mapper_Sl2g%weight_identifier = wgtIdSl2g + mapper_Sl2g%mbname = 'mapper_Sl2g' + arearead = 0 ! no need for areas + call moab_map_init_rcfile( mapper_Sl2g, type_grid, & + 'seq_maps.rc', 'lnd2glc_smapname:', 'lnd2glc_smaptype:', samegrid_lg, & + arearead, wgtIdSl2g, 'mapper_Sl2g MOAB initialization', esmf_map_flag) + + ! flux (conservative) map; arearead=2 loads the map file area_b into the glc + ! mesh aream, reproducing the mct driver's seq_map_readdata of area_b -> aream + mapper_Fl2g%src_mbid = mblxid + mapper_Fl2g%tgt_mbid = mbgxid + mapper_Fl2g%intx_mbid = mbintxlg + mapper_Fl2g%src_context = lnd(1)%cplcompid + mapper_Fl2g%intx_context = idintx + wgtIdFl2g = 'flux_l2g' + mapper_Fl2g%weight_identifier = wgtIdFl2g + mapper_Fl2g%mbname = 'mapper_Fl2g' + arearead = 2 ! area_b for glc aream + call moab_map_init_rcfile( mapper_Fl2g, type_grid, & + 'seq_maps.rc', 'lnd2glc_fmapname:', 'lnd2glc_fmaptype:', samegrid_lg, & + arearead, wgtIdFl2g, 'mapper_Fl2g MOAB initialization', esmf_map_flag) + + ! one mesh migration and one comm graph cover both maps (same coverage) + ierr = iMOAB_MigrateMapMesh (mblxid, mbintxlg, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type_grid, lnd(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating lnd mesh for map lnd 2 glc' + call shr_sys_abort(subname//' ERROR in migrating lnd mesh for map lnd 2 glc') + endif + ierr = iMOAB_ComputeCommGraph( mblxid, mbintxlg, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, & + type_grid, type_grid, lnd(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, LND-GLC' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, LND-GLC') + endif + + ! define the per-elevation-class projection target tags on the glc mesh + tagtype = 1 ! dense, double + numco = 1 + tagname = trim(seq_flds_l2x_fields_to_glc)//C_NULL_CHAR + ierr = iMOAB_DefineTagStorage(mbgxid, tagname, tagtype, numco, tagindex ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR defining lnd per-EC tags on the glc mesh') + endif + + ! glc coupler mesh size, used by the downscaling + ierr = iMOAB_GetMeshInfo ( mbgxid, nvert, nvise, nbl, nsurf, nvisBC ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR getting glc coupler mesh info') + endif + lsize_gm = nvise(1) + + ! moab context for the glc->lnd conservative map (used here by the smb + ! renormalization). If prep_lnd has set it up already (glc_c2_lnd), just + ! point our mapper at the same map app and weights; otherwise (e.g. TG + ! compsets) register and load it ourselves. + wgtIdFg2l = 'flux_g2l' + if (mbintxgl < 0) then + if (iamroot_CPLID) then + write(logunit,*) ' ' + write(logunit,F00) 'Initializing MOAB mapper_Fg2l (glc->lnd map app)' + end if + appname = "GLC_LND_COU"//C_NULL_CHAR + idintx = 100*glc(1)%cplcompid + lnd(1)%cplcompid + ierr = iMOAB_RegisterApplication(trim(appname), mpicom_CPLID, idintx, mbintxgl) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in registering glc lnd map app' + call shr_sys_abort(subname//' ERROR in registering glc lnd map app') + endif + mapper_Fg2l%src_mbid = mbgxid + mapper_Fg2l%tgt_mbid = mblxid + mapper_Fg2l%intx_mbid = mbintxgl + mapper_Fg2l%src_context = glc(1)%cplcompid + mapper_Fg2l%intx_context = idintx + mapper_Fg2l%weight_identifier = wgtIdFg2l + mapper_Fg2l%mbname = 'mapper_Fg2l' + arearead = 0 + call moab_map_init_rcfile( mapper_Fg2l, type_grid, & + 'seq_maps.rc', 'glc2lnd_fmapname:', 'glc2lnd_fmaptype:', samegrid_lg, & + arearead, wgtIdFg2l, 'mapper_Fg2l (prep_glc) MOAB initialization', esmf_map_flag) + ierr = iMOAB_MigrateMapMesh (mbgxid, mbintxgl, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type_grid, glc(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating glc mesh for map glc 2 lnd' + call shr_sys_abort(subname//' ERROR in migrating glc mesh for map glc 2 lnd') + endif + ierr = iMOAB_ComputeCommGraph( mbgxid, mbintxgl, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, & + type_grid, type_grid, glc(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, GLC-LND' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, GLC-LND') + endif + else + ! reuse the map app and weights loaded by prep_lnd_init + mapper_Fg2l%src_mbid = mbgxid + mapper_Fg2l%tgt_mbid = mblxid + mapper_Fg2l%intx_mbid = mbintxgl + mapper_Fg2l%src_context = glc(1)%cplcompid + mapper_Fg2l%intx_context = 100*glc(1)%cplcompid + lnd(1)%cplcompid + mapper_Fg2l%weight_identifier = wgtIdFg2l + mapper_Fg2l%mbname = 'mapper_Fg2l' + endif + + endif ! mblxid and mbgxid + end if call shr_sys_flush(logunit) end if - if (glc_present .and. ocn_c2_glcshelf) then - + ! setup needed for either kind of ocn2glc coupling + if (glc_present .and. (ocn_c2_glctf .or. ocn_c2_glcshelf)) then call seq_comm_getData(CPLID, & mpicom=mpicom_CPLID, iamroot=iamroot_CPLID) @@ -275,21 +457,18 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) x2gacc_gx_cnt = 0 samegrid_go = .true. if (trim(ocn_gnam) /= trim(glc_gnam)) samegrid_go = .false. - if (iamroot_CPLID) then - write(logunit,*) ' ' - write(logunit,F00) 'Initializing mapper_So2g' - end if - call seq_map_init_rcfile(mapper_So2g, ocn(1), glc(1), & - 'seq_maps.rc','ocn2glc_smapname:','ocn2glc_smaptype:',samegrid_go, & - 'mapper_So2g initialization',esmf_map_flag) - if (iamroot_CPLID) then - write(logunit,*) ' ' - write(logunit,F00) 'Initializing mapper_Fo2g' - end if - call seq_map_init_rcfile(mapper_Fo2g, ocn(1), glc(1), & - 'seq_maps.rc','ocn2glc_fmapname:','ocn2glc_fmaptype:',samegrid_go, & - 'mapper_Fo2g initialization',esmf_map_flag) + end if + + ! setup needed for ocn2glc TF coupling + ! MOABTODO: give these mappers moab context when the ocn<->glc coupling is ported + if (glc_present .and. ocn_c2_glctf) then + call seq_map_mapinit(mapper_So2g_tf, mpicom_CPLID) + end if + ! setup needed for ocn2glcshelf coupling + if (glc_present .and. ocn_c2_glcshelf) then + call seq_map_mapinit(mapper_So2g_shelf, mpicom_CPLID) + call seq_map_mapinit(mapper_Fo2g_shelf, mpicom_CPLID) !Initialize module-level arrays associated with compute_melt_fluxes allocate(oceanTemperature(lsize_g)) allocate(oceanSalinity(lsize_g)) @@ -307,10 +486,9 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) ! TODO: Can we allocate these only while used or are we worried about performance hit? ! TODO: add deallocates! - call shr_sys_flush(logunit) - end if + call shr_sys_flush(logunit) end subroutine prep_glc_init @@ -502,6 +680,7 @@ subroutine prep_glc_accum_avg(timer, lnd2glc_averaged_now) call mct_avect_avg(l2gacc_lx(eli), l2gacc_lx_cnt) end do end if + l2gacc_lx_cnt_avg = l2gacc_lx_cnt l2gacc_lx_cnt = 0 ! Accumulation for OCN @@ -521,6 +700,247 @@ subroutine prep_glc_accum_avg(timer, lnd2glc_averaged_now) end subroutine prep_glc_accum_avg + !================================================================================================ + + subroutine prep_glc_accum_lnd_moab(timer) + + !--------------------------------------------------------------- + ! Description + ! Accumulate the per-elevation-class land forcing for glc, by reading the + ! coupler land mesh tags into a plain array (prep_rof accumulation pattern) + ! + use iMOAB, only : iMOAB_GetDoubleTagStorage + use iso_c_binding, only : C_NULL_CHAR + ! + ! Arguments + character(len=*), intent(in) :: timer + ! + ! Local Variables + integer :: ierr, ent_type, arrsize + character(CXX) :: tagname + character(*), parameter :: subname = '(prep_glc_accum_lnd_moab)' + !--------------------------------------------------------------- + + if (.not. allocated(l2gacc_lm)) return ! moab accumulator not set up + + call t_drvstartf (trim(timer),barrier=mpicom_CPLID) + tagname = trim(sharedFieldsLndGlc)//C_NULL_CHAR + arrsize = nflds_lg * lsize_lm + ent_type = 1 ! cells + ierr = iMOAB_GetDoubleTagStorage ( mblxid, tagname, arrsize, ent_type, l2x_lm2 ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' error in getting per-EC lnd fields for glc accumulation') + endif + if (l2gacc_lm_cnt == 0) then + l2gacc_lm = l2x_lm2 + else + l2gacc_lm = l2gacc_lm + l2x_lm2 + endif + l2gacc_lm_cnt = l2gacc_lm_cnt + 1 + call t_drvstopf (trim(timer)) + + end subroutine prep_glc_accum_lnd_moab + + !================================================================================================ + + subroutine prep_glc_accum_avg_moab(timer, lnd2glc_averaged_now) + + !--------------------------------------------------------------- + ! Description + ! Finalize the accumulation of the land forcing for glc, and set the averaged + ! values back into the coupler land mesh tags, from where the lnd->glc maps + ! read them (prep_rof set-back pattern). + ! + ! Note: the mct version also averages the ocn accumulation (x2gacc); that part + ! of the moab port is deferred with the rest of the ocn<->glc coupling. + ! + use iMOAB, only : iMOAB_SetDoubleTagStorage + use iso_c_binding, only : C_NULL_CHAR + ! + ! Arguments + character(len=*), intent(in) :: timer + logical, intent(inout) :: lnd2glc_averaged_now ! Set to .true. if lnd2glc averages were taken this timestep (otherwise left unchanged) + ! + ! Local Variables + integer :: ierr, ent_type, arrsize + real(r8) :: ravg ! averaging factor + character(CXX) :: tagname + character(*), parameter :: subname = '(prep_glc_accum_avg_moab)' + !--------------------------------------------------------------- + + if (.not. allocated(l2gacc_lm)) return ! moab accumulator not set up + + call t_drvstartf (trim(timer),barrier=mpicom_CPLID) + if (l2gacc_lm_cnt > 0) then + lnd2glc_averaged_now = .true. + end if + if (l2gacc_lm_cnt > 1) then + ravg = 1.0_r8/real(l2gacc_lm_cnt, r8) + l2gacc_lm = l2gacc_lm * ravg + end if + l2gacc_lm_cnt_avg = l2gacc_lm_cnt + l2gacc_lm_cnt = 0 + + ! set the averaged values back into the land mesh tags + tagname = trim(sharedFieldsLndGlc)//C_NULL_CHAR + arrsize = nflds_lg * lsize_lm + ent_type = 1 ! cells + ierr = iMOAB_SetDoubleTagStorage ( mblxid, tagname, arrsize, ent_type, l2gacc_lm ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' error in setting accumulated per-EC lnd fields on land mesh') + endif + call t_drvstopf (trim(timer)) + + end subroutine prep_glc_accum_avg_moab + + !================================================================================================ + + subroutine prep_glc_mrg_ocn(infodata, fractions_gx, timer_mrg) + + !--------------------------------------------------------------- + ! Description + ! Merge glc inputs + ! + ! Arguments + type(seq_infodata_type) , intent(in) :: infodata + type(mct_aVect) , intent(in) :: fractions_gx(:) + character(len=*) , intent(in) :: timer_mrg + ! + ! Local Variables + integer :: egi, eoi, efi + type(mct_avect), pointer :: x2g_gx + character(*), parameter :: subname = '(prep_glc_mrg_ocn)' + !--------------------------------------------------------------- + + call t_drvstartf (trim(timer_mrg),barrier=mpicom_CPLID) + do egi = 1,num_inst_glc + ! Use fortran mod to address ensembles in merge + eoi = mod((egi-1),num_inst_ocn) + 1 + efi = mod((egi-1),num_inst_frc) + 1 + + x2g_gx => component_get_x2c_cx(glc(egi)) + call prep_glc_merge_ocn_forcing(o2x_gx(eoi), fractions_gx(efi), x2g_gx) + enddo + call t_drvstopf (trim(timer_mrg)) + + end subroutine prep_glc_mrg_ocn + + !================================================================================================ + + subroutine prep_glc_merge_ocn_forcing( o2x_g, fractions_g, x2g_g ) + + !----------------------------------------------------------------------- + ! Description + ! "Merge" ocean forcing for glc input. + ! + ! State fields are copied directly, meaning that averages are taken just over the + ! ocean-covered portion of the glc domain. + ! + ! Flux fields are downweighted by landfrac, which effectively sends a 0 flux from the + ! non-ocean-covered portion of the glc domain. + ! + ! Arguments + type(mct_aVect), intent(inout) :: o2x_g ! input + type(mct_aVect), intent(in) :: fractions_g + type(mct_aVect), intent(inout) :: x2g_g ! output + !----------------------------------------------------------------------- + + integer :: num_flux_fields + integer :: num_state_fields + integer :: nflds + integer :: i,n + integer :: mrgstr_index + integer :: index_o2x + integer :: index_x2g + integer :: index_ofrac + integer :: lsize + logical :: iamroot + logical, save :: first_time = .true. + character(CL),allocatable :: mrgstr(:) ! temporary string + character(CL) :: field ! string converted to char + character(*), parameter :: subname = '(prep_glc_merge_ocn_forcing) ' + + !----------------------------------------------------------------------- + + call seq_comm_getdata(CPLID, iamroot=iamroot) + lsize = mct_aVect_lsize(x2g_g) + + !num_flux_fields = shr_string_listGetNum(trim(seq_flds_x2g_fluxes_from_ocn)) + num_flux_fields = 0 + num_state_fields = shr_string_listGetNum(trim(seq_flds_x2g_tf_states_from_ocn)) + + if (first_time) then + nflds = num_flux_fields + num_state_fields + allocate(mrgstr(nflds)) + end if + + mrgstr_index = 1 + + do i = 1, num_state_fields + call seq_flds_getField(field, i, seq_flds_x2g_tf_states_from_ocn) + index_o2x = mct_aVect_indexRA(o2x_g, trim(field)) + index_x2g = mct_aVect_indexRA(x2g_g, trim(field)) + + if (first_time) then + mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // & + ' = o2x%'//trim(field) + end if + + do n = 1, lsize + x2g_g%rAttr(index_x2g,n) = o2x_g%rAttr(index_o2x,n) + end do + + mrgstr_index = mrgstr_index + 1 + enddo + + !index_lfrac = mct_aVect_indexRA(fractions_g,"lfrac") + !do i = 1, num_flux_fields + + ! call seq_flds_getField(field, i, seq_flds_x2g_fluxes_from_lnd) + ! index_l2x = mct_aVect_indexRA(l2x_g, trim(field)) + ! index_x2g = mct_aVect_indexRA(x2g_g, trim(field)) + + ! if (trim(field) == qice_fieldname) then + + ! if (first_time) then + ! mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // & + ! ' = l2x%'//trim(field) + ! end if + + ! ! treat qice as if it were a state variable, with a simple copy. + ! do n = 1, lsize + ! x2g_g%rAttr(index_x2g,n) = l2x_g%rAttr(index_l2x,n) + ! end do + + ! else + ! write(logunit,*) subname,' ERROR: Flux fields other than ', & + ! qice_fieldname, ' currently are not handled in lnd2glc remapping.' + ! write(logunit,*) '(Attempt to handle flux field <', trim(field), '>.)' + ! write(logunit,*) 'Substantial thought is needed to determine how to remap other fluxes' + ! write(logunit,*) 'in a smooth, conservative manner.' + ! call shr_sys_abort(subname//& + ! ' ERROR: Flux fields other than qice currently are not handled in lnd2glc remapping.') + ! endif ! qice_fieldname + + ! mrgstr_index = mrgstr_index + 1 + + !end do + + if (first_time) then + if (iamroot) then + write(logunit,'(A)') subname//' Summary:' + do i = 1,nflds + write(logunit,'(A)') trim(mrgstr(i)) + enddo + endif + deallocate(mrgstr) + endif + + first_time = .false. + + end subroutine prep_glc_merge_ocn_forcing + + !================================================================================================ subroutine prep_glc_mrg_lnd(infodata, fractions_gx, timer_mrg) @@ -604,7 +1024,7 @@ subroutine prep_glc_merge_lnd_forcing( l2x_g, fractions_g, x2g_g ) mrgstr_index = 1 do i = 1, num_state_fields - call seq_flds_getField(field, i, seq_flds_x2g_states) + call seq_flds_getField(field, i, seq_flds_x2g_states_from_lnd) index_l2x = mct_aVect_indexRA(l2x_g, trim(field)) index_x2g = mct_aVect_indexRA(x2g_g, trim(field)) @@ -668,13 +1088,15 @@ subroutine prep_glc_merge_lnd_forcing( l2x_g, fractions_g, x2g_g ) end subroutine prep_glc_merge_lnd_forcing - subroutine prep_glc_calc_o2x_gx(timer) + subroutine prep_glc_calc_o2x_gx(ocn_c2_glctf, ocn_c2_glcshelf, timer) !--------------------------------------------------------------- ! Description ! Create o2x_gx ! Arguments character(len=*), intent(in) :: timer + logical, intent(in) :: ocn_c2_glctf + logical, intent(in) :: ocn_c2_glcshelf character(*), parameter :: subname = '(prep_glc_calc_o2x_gx)' ! Local Variables @@ -684,15 +1106,14 @@ subroutine prep_glc_calc_o2x_gx(timer) call t_drvstartf (trim(timer),barrier=mpicom_CPLID) do eoi = 1,num_inst_ocn o2x_ox => component_get_c2x_cx(ocn(eoi)) -!MOABTODO: uncomment when porting glc -! if (ocn_c2_glctf) then -! call seq_map_map(mapper_So2g_tf, o2x_ox, o2x_gx(eoi), & -! fldlist=seq_flds_x2g_tf_states_from_ocn,norm=.true.) -! end if -! if (ocn_c2_glcshelf) then -! call seq_map_map(mapper_So2g_shelf, o2x_ox, o2x_gx(eoi), & -! fldlist=seq_flds_x2g_shelf_states_from_ocn,norm=.true.) -! end if + if (ocn_c2_glctf) then + call seq_map_map(mapper_So2g_tf, o2x_ox, o2x_gx(eoi), & + fldlist=seq_flds_x2g_tf_states_from_ocn,norm=.true.) + end if + if (ocn_c2_glcshelf) then + call seq_map_map(mapper_So2g_shelf, o2x_ox, o2x_gx(eoi), & + fldlist=seq_flds_x2g_shelf_states_from_ocn,norm=.true.) + end if enddo call t_drvstopf (trim(timer)) @@ -857,7 +1278,7 @@ subroutine prep_glc_calculate_subshelf_boundary_fluxes !Done here instead of in glc-frequency mapping so it happens within ocean coupling interval. ! Also could map o2x_ox->o2x_gx(1) but using x2g_gx as destination allows us to see ! these fields on the GLC grid of the coupler history file, which helps with debugging. - call seq_map_map(mapper_So2g, o2x_ox, x2g_gx, & + call seq_map_map(mapper_So2g_shelf, o2x_ox, x2g_gx, & fldlist=seq_flds_x2g_shelf_states_from_ocn,norm=.true.) ! inputs to melt flux calculation @@ -959,12 +1380,67 @@ subroutine prep_glc_zero_fields() do egi = 1,num_inst_glc x2g_gx => component_get_x2c_cx(glc(egi)) - call mct_aVect_zero(x2g_gx) + if (associated(x2g_gx)) then + call mct_aVect_zero(x2g_gx) + else + write(logunit,*) ' ' + write(logunit,*) 'Warning: x2g_gx not associated for glc(', egi, ')' + end if end do + end subroutine prep_glc_zero_fields !================================================================================================ + subroutine prep_glc_zero_fields_moab() + + !--------------------------------------------------------------- + ! Description + ! Set glc input tags (x2g fields on the coupler-side glc mesh) to zero + ! + ! This is the moab counterpart of prep_glc_zero_fields; see the note there + ! about why zeroing is needed for exact restart tests. + + use iMOAB, only : iMOAB_GetMeshInfo, iMOAB_SetDoubleTagStorage + use seq_comm_mct, only : mbgxid + use ISO_C_BINDING, only : C_NULL_CHAR + use shr_kind_mod, only : CXX => shr_kind_CXX + use shr_string_mod, only : shr_string_listGetNum + + ! Local Variables + integer :: ierr, ent_type, arrsize, nxflds, lsize_gm + integer :: nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) + real(r8), allocatable :: tmparray(:) + character(CXX) :: tagname + character(*), parameter :: subname = '(prep_glc_zero_fields_moab)' + !--------------------------------------------------------------- + + if (mbgxid < 0) return ! nothing to do if the glc coupler mesh does not exist + + ierr = iMOAB_GetMeshInfo ( mbgxid, nvert, nvise, nbl, nsurf, nvisBC ) + if (ierr .ne. 0) then + write(logunit,*) subname,' cant get size of glc mesh on coupler' + call shr_sys_abort(subname//' ERROR in getting size of glc mesh on coupler') + endif + lsize_gm = nvise(1) + ent_type = 1 ! cells + + nxflds = shr_string_listGetNum(seq_flds_x2g_fields) + arrsize = nxflds * lsize_gm + allocate (tmparray(arrsize)) + tmparray = 0._r8 + tagname = trim(seq_flds_x2g_fields)//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage(mbgxid, tagname, arrsize, ent_type, tmparray) + if (ierr .ne. 0) then + write(logunit,*) subname,' cant zero out x2g tags on glc coupler mesh' + call shr_sys_abort(subname//' cant zero out x2g tags on glc coupler mesh') + endif + deallocate (tmparray) + + end subroutine prep_glc_zero_fields_moab + + !================================================================================================ + subroutine prep_glc_map_qice_conservative_lnd2glc(egi, eli, fractions_lx, & mapper_Sl2g, mapper_Fg2l) @@ -1202,8 +1678,9 @@ subroutine prep_glc_renormalize_smb(eli, fractions_lx, g2x_gx, mapper_Fg2l, area aream_l(:) = dom_l%data%rAttr(km,:) ! Export land fractions from fractions_lx to a local array + ! Note that for E3SM we are using lfrin instead of lfrac allocate(lfrac(lsize_l)) - call mct_aVect_exportRattr(fractions_lx, "lfrac", lfrac) + call mct_aVect_exportRattr(fractions_lx, "lfrin", lfrac) ! Map Sg_icemask from the glc grid to the land grid. ! This may not be necessary, if Sg_icemask_l has already been mapped from Sg_icemask_g. @@ -1386,6 +1863,8 @@ subroutine prep_glc_renormalize_smb(eli, fractions_lx, g2x_gx, mapper_Fg2l, area endif if (iamroot) then + write(logunit,*) 'global_accum_on_land_grid = ', global_accum_on_land_grid + write(logunit,*) 'global_accum_on_glc_grid = ', global_accum_on_glc_grid write(logunit,*) 'accum_renorm_factor = ', accum_renorm_factor write(logunit,*) 'ablat_renorm_factor = ', ablat_renorm_factor endif @@ -1410,6 +1889,436 @@ end subroutine prep_glc_renormalize_smb !================================================================================================ + subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) + !--------------------------------------------------------------- + ! Description + ! moab version of prep_glc_calc_l2x_gx (+ the merge preparation): + ! - one batched horizontal map of all the per-elevation-class land fields to the + ! glc mesh, weighted by lfrac with normalization (the same semantics as the + ! per-field maps done through map_lnd2glc in the mct version) + ! - elevation-class vertical interpolation on arrays fetched from the glc mesh tags + ! - conservative correction of the qice flux (area/aream pre-adjustment and, if + ! enabled, global smb renormalization) + ! - results are written directly into the x2g tag names on the glc mesh, which + ! makes the merge a plain no-op, matching the mct merge (plain copies) + ! + use iMOAB, only : iMOAB_GetDoubleTagStorage, iMOAB_SetDoubleTagStorage + use iso_c_binding, only : C_NULL_CHAR + use map_lnd2glc_mod, only : get_glc_elevation_classes, map_lnd2glc_vertical_interp + use shr_string_mod, only : shr_string_listGetNum + ! + ! Arguments + type(mct_aVect) , intent(in) :: fractions_lx(:) + character(len=*), intent(in) :: timer + ! + ! Local Variables + integer :: num_flux_fields + integer :: num_state_fields + integer :: field_num + integer :: ierr, ent_type, arrsize, n, ec, kf, nEC + character(len=cl) :: fieldname + real(r8), allocatable :: data_lg(:,:) ! all mapped per-EC fields on the glc mesh + real(r8), allocatable :: glc_ice_covered(:) + real(r8), allocatable :: glc_topo(:) + real(r8), allocatable :: topo_g_EC(:,:) + real(r8), allocatable :: data_g_EC(:,:) + real(r8), allocatable :: data_g_bare(:) + real(r8), allocatable :: data_g(:) ! downscaled field on the glc mesh + real(r8), allocatable :: area_g(:) + real(r8), allocatable :: aream_g(:) + integer , allocatable :: glc_elevclass(:) + character(CXX) :: tagname + character(*), parameter :: subname = '(prep_glc_calc_l2x_gx_moab)' + !--------------------------------------------------------------- + + if ((mblxid < 0) .or. (mbgxid < 0)) return + + call t_drvstartf (trim(timer),barrier=mpicom_CPLID) + + ent_type = 1 ! cells + + ! Horizontal map of all per-EC fields at once, weighted by lfrac and normalized. + ! The source tags hold the averaged accumulation, set back on the land mesh by + ! prep_glc_accum_avg_moab. The attribute vector arguments are metadata only. + call seq_map_map(mapper_Sl2g, l2gacc_lx(1), l2x_gx(1), & + fldlist=trim(sharedFieldsLndGlc), norm=.true., & + avwts_s=fractions_lx(1), avwtsfld_s='lfrac') + + ! fetch the mapped per-EC fields and the glc state needed for the vertical interpolation + allocate(data_lg(lsize_gm, nflds_lg)) + data_lg = 0._r8 + tagname = trim(sharedFieldsLndGlc)//C_NULL_CHAR + arrsize = nflds_lg * lsize_gm + ierr = iMOAB_GetDoubleTagStorage(mbgxid, tagname, arrsize, ent_type, data_lg) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR getting mapped per-EC fields on the glc mesh') + endif + + allocate(glc_ice_covered(lsize_gm), glc_topo(lsize_gm), glc_elevclass(lsize_gm)) + tagname = trim(Sg_frac_field)//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mbgxid, tagname, lsize_gm, ent_type, glc_ice_covered) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting '//Sg_frac_field) + tagname = trim(Sg_topo_field)//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mbgxid, tagname, lsize_gm, ent_type, glc_topo) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting '//Sg_topo_field) + + call get_glc_elevation_classes(glc_ice_covered, glc_topo, glc_elevclass) + + nEC = glc_get_num_elevation_classes() + allocate(topo_g_EC(lsize_gm, nEC)) + allocate(data_g_EC(lsize_gm, nEC)) + allocate(data_g_bare(lsize_gm)) + allocate(data_g(lsize_gm)) + do ec = 1, nEC + kf = mct_aVect_indexRA(l2gacc_lx(1), 'Sl_topo'//glc_elevclass_as_string(ec)) + topo_g_EC(:,ec) = data_lg(:,kf) + end do + + ! area arrays needed for the qice conservation correction + allocate(area_g(lsize_gm), aream_g(lsize_gm)) + tagname = 'area'//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mbgxid, tagname, lsize_gm, ent_type, area_g) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting area on the glc mesh') + tagname = 'aream'//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mbgxid, tagname, lsize_gm, ent_type, aream_g) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting aream on the glc mesh') + + num_flux_fields = shr_string_listGetNum(trim(seq_flds_x2g_fluxes_from_lnd)) + num_state_fields = shr_string_listGetNum(trim(seq_flds_x2g_states_from_lnd)) + + do field_num = 1, num_flux_fields + call seq_flds_getField(fieldname, field_num, seq_flds_x2g_fluxes_from_lnd) + + if (trim(fieldname) == qice_fieldname) then + + do ec = 1, nEC + kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(ec)) + data_g_EC(:,ec) = data_lg(:,kf) + end do + kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(0)) + data_g_bare(:) = data_lg(:,kf) + + call map_lnd2glc_vertical_interp(glc_topo, topo_g_EC, data_g_EC, data_g_bare, & + glc_elevclass, data_g) + + ! Preemptive adjustment for area differences between the glc model and the + ! coupler: the drv2mdl flux correction on the component side multiplies by + ! aream/area, so multiply here by area/aream (see the discussion in + ! prep_glc_map_qice_conservative_lnd2glc) + do n = 1, lsize_gm + if (aream_g(n) > 0.0_r8) then + data_g(n) = data_g(n) * area_g(n)/aream_g(n) + else + data_g(n) = 0.0_r8 + endif + enddo + + if (smb_renormalize) then + call prep_glc_renormalize_smb_moab(fractions_lx(1), aream_g, data_g) + end if + + tagname = trim(qice_fieldname)//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage(mbgxid, tagname, lsize_gm, ent_type, data_g) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR setting '//qice_fieldname) + + else + write(logunit,*) subname,' ERROR: Flux fields other than ', & + qice_fieldname, ' currently are not handled in lnd2glc remapping.' + call shr_sys_abort(subname//& + ' ERROR: Flux fields other than qice currently are not handled in lnd2glc remapping.') + endif ! qice_fieldname + end do + + do field_num = 1, num_state_fields + call seq_flds_getField(fieldname, field_num, seq_flds_x2g_states_from_lnd) + + do ec = 1, nEC + kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(ec)) + data_g_EC(:,ec) = data_lg(:,kf) + end do + kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(0)) + data_g_bare(:) = data_lg(:,kf) + + call map_lnd2glc_vertical_interp(glc_topo, topo_g_EC, data_g_EC, data_g_bare, & + glc_elevclass, data_g) + + tagname = trim(fieldname)//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage(mbgxid, tagname, lsize_gm, ent_type, data_g) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR setting '//trim(fieldname)) + end do + + deallocate(data_lg, glc_ice_covered, glc_topo, glc_elevclass) + deallocate(topo_g_EC, data_g_EC, data_g_bare, data_g) + deallocate(area_g, aream_g) + + call t_drvstopf (trim(timer)) + + end subroutine prep_glc_calc_l2x_gx_moab + + !================================================================================================ + + subroutine prep_glc_renormalize_smb_moab(fractions_lx, aream_g, qice_g) + + ! moab version of prep_glc_renormalize_smb: renormalizes the surface mass balance + ! so that the global integral on the glc grid equals the one on the land grid. + ! Same algorithm as the mct version, with the land- and glc-side fields fetched + ! from the coupler mesh tags and the per-EC land fractions recomputed through + ! map_glc2lnd_ec_moab. + + use iMOAB, only : iMOAB_GetDoubleTagStorage + use iso_c_binding, only : C_NULL_CHAR + use map_glc2lnd_mod, only : map_glc2lnd_ec_moab + use shr_mpi_mod, only : shr_mpi_sum, shr_mpi_bcast + + ! Arguments + type(mct_aVect) , intent(in) :: fractions_lx ! fractions on the land grid (metadata only) + real(r8) , intent(in) :: aream_g(:) ! cell areas on glc grid, for mapping + real(r8) , intent(inout) :: qice_g(:) ! qice data on glc grid + ! + ! Local Variables + logical :: iamroot + integer :: ierr, ent_type, n, ec, kf, nEC + + real(r8), allocatable :: aream_l(:) ! cell areas on land grid, for mapping + real(r8), allocatable :: lfrac(:) ! land fraction (lfrin) on land grid + real(r8), allocatable :: Sg_icemask_l(:) ! icemask on land grid + real(r8), allocatable :: Sg_icemask_g(:) ! icemask on glc grid + real(r8), allocatable :: qice_l(:,:) ! SMB (Flgl_qice) per EC on land grid + real(r8), allocatable :: frac_l(:,:) ! EC fractions (Sg_ice_covered) on land grid + + type(mct_aVect) :: av_dum ! zero-size av for the seq_map_map interface + + real(r8) :: local_accum_on_land_grid + real(r8) :: global_accum_on_land_grid + real(r8) :: local_accum_on_glc_grid + real(r8) :: global_accum_on_glc_grid + + real(r8) :: local_ablat_on_land_grid + real(r8) :: global_ablat_on_land_grid + real(r8) :: local_ablat_on_glc_grid + real(r8) :: global_ablat_on_glc_grid + + real(r8) :: accum_renorm_factor ! ratio between global accumulation on the two grids + real(r8) :: ablat_renorm_factor ! ratio between global ablation on the two grids + + real(r8) :: effective_area ! grid cell area multiplied by min(lfrac,Sg_icemask_l) + + character(CXX) :: tagname + character(*), parameter :: subname = '(prep_glc_renormalize_smb_moab)' + !--------------------------------------------------------------- + + call seq_comm_getdata(CPLID, iamroot=iamroot) + ent_type = 1 ! cells + nEC = glc_get_num_elevation_classes() + + ! land areas for the mapping and the land fraction basis; note that for E3SM we + ! are using lfrin instead of lfrac (see prep_glc_renormalize_smb) + allocate(aream_l(lsize_lm), lfrac(lsize_lm)) + tagname = 'aream'//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mblxid, tagname, lsize_lm, ent_type, aream_l) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting aream on the land mesh') + tagname = 'lfrin'//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mblxid, tagname, lsize_lm, ent_type, lfrac) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting lfrin on the land mesh') + + ! Map Sg_icemask from the glc mesh to the land mesh tag of the same name (see + ! prep_glc_renormalize_smb for why this mapping is redone here); the result is + ! identical to what prep_lnd_calc_g2x_lx produces from the same static g2x data. + call mct_aVect_init(av_dum, rList=Sg_icemask_field, lsize=0) + call seq_map_map(mapper_Fg2l, av_dum, av_dum, fldlist=Sg_icemask_field, norm=.true.) + call mct_aVect_clean(av_dum) + allocate(Sg_icemask_l(lsize_lm)) + tagname = trim(Sg_icemask_field)//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mblxid, tagname, lsize_lm, ent_type, Sg_icemask_l) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting Sg_icemask on the land mesh') + + ! Map Sg_ice_covered (per elevation class) from glc to land, and return the + ! normalized per-EC fractions (see prep_glc_renormalize_smb for the rationale) + allocate(frac_l(lsize_lm, 0:nEC)) + call map_glc2lnd_ec_moab(mapper_Fg2l, & + frac_field = Sg_frac_field, & + topo_field = Sg_topo_field, & + icemask_field = Sg_icemask_field, & + extra_fields = ' ', & ! no extra fields + frac_l_out = frac_l) + + ! qice per elevation class on the land grid, from the averaged accumulator + allocate(qice_l(lsize_lm, 0:nEC)) + do ec = 0, nEC + kf = mct_aVect_indexRA(l2gacc_lx(1), qice_fieldname//glc_elevclass_as_string(ec)) + qice_l(:,ec) = l2gacc_lm(:,kf) + end do + + ! Sum qice over local land grid cells + + local_accum_on_land_grid = 0.0_r8 + local_ablat_on_land_grid = 0.0_r8 + + do n = 1, lsize_lm + + effective_area = min(lfrac(n),Sg_icemask_l(n)) * aream_l(n) + + do ec = 0, nEC + + if (qice_l(n,ec) >= 0.0_r8) then + local_accum_on_land_grid = local_accum_on_land_grid & + + effective_area * frac_l(n,ec) * qice_l(n,ec) + else + local_ablat_on_land_grid = local_ablat_on_land_grid & + + effective_area * frac_l(n,ec) * qice_l(n,ec) + endif + + enddo ! ec + + enddo ! n + + call shr_mpi_sum(local_accum_on_land_grid, & + global_accum_on_land_grid, & + mpicom_CPLID, 'accum_l') + + call shr_mpi_sum(local_ablat_on_land_grid, & + global_ablat_on_land_grid, & + mpicom_CPLID, 'ablat_l') + + call shr_mpi_bcast(global_accum_on_land_grid, mpicom_CPLID) + call shr_mpi_bcast(global_ablat_on_land_grid, mpicom_CPLID) + + ! Sum qice_g over local glc grid cells (see prep_glc_renormalize_smb for the + ! discussion of the areas used here) + + allocate(Sg_icemask_g(size(qice_g))) + tagname = trim(Sg_icemask_field)//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mbgxid, tagname, size(qice_g), ent_type, Sg_icemask_g) + if (ierr .ne. 0) call shr_sys_abort(subname//' ERROR getting Sg_icemask on the glc mesh') + + local_accum_on_glc_grid = 0.0_r8 + local_ablat_on_glc_grid = 0.0_r8 + + do n = 1, size(qice_g) + + if (qice_g(n) >= 0.0_r8) then + local_accum_on_glc_grid = local_accum_on_glc_grid & + + Sg_icemask_g(n) * aream_g(n) * qice_g(n) + else + local_ablat_on_glc_grid = local_ablat_on_glc_grid & + + Sg_icemask_g(n) * aream_g(n) * qice_g(n) + endif + + enddo ! n + + call shr_mpi_sum(local_accum_on_glc_grid, & + global_accum_on_glc_grid, & + mpicom_CPLID, 'accum_g') + + call shr_mpi_sum(local_ablat_on_glc_grid, & + global_ablat_on_glc_grid, & + mpicom_CPLID, 'ablat_g') + + call shr_mpi_bcast(global_accum_on_glc_grid, mpicom_CPLID) + call shr_mpi_bcast(global_ablat_on_glc_grid, mpicom_CPLID) + + ! Renormalize + + if (global_accum_on_glc_grid > 0.0_r8) then + accum_renorm_factor = global_accum_on_land_grid / global_accum_on_glc_grid + else + accum_renorm_factor = 0.0_r8 + endif + + if (global_ablat_on_glc_grid < 0.0_r8) then ! negative by definition + ablat_renorm_factor = global_ablat_on_land_grid / global_ablat_on_glc_grid + else + ablat_renorm_factor = 0.0_r8 + endif + + if (iamroot) then + write(logunit,*) 'moab global_accum_on_land_grid = ', global_accum_on_land_grid + write(logunit,*) 'moab global_accum_on_glc_grid = ', global_accum_on_glc_grid + write(logunit,*) 'moab accum_renorm_factor = ', accum_renorm_factor + write(logunit,*) 'moab ablat_renorm_factor = ', ablat_renorm_factor + endif + + do n = 1, size(qice_g) + if (qice_g(n) >= 0.0_r8) then + qice_g(n) = qice_g(n) * accum_renorm_factor + else + qice_g(n) = qice_g(n) * ablat_renorm_factor + endif + enddo + + deallocate(aream_l) + deallocate(lfrac) + deallocate(Sg_icemask_l) + deallocate(Sg_icemask_g) + deallocate(qice_l) + deallocate(frac_l) + + end subroutine prep_glc_renormalize_smb_moab + + !================================================================================================ + + subroutine prep_glc_mrg_lnd_moab(infodata, timer_mrg) + + !--------------------------------------------------------------- + ! Description + ! moab counterpart of prep_glc_mrg_lnd. The mct merge is a plain copy of the + ! mapped fields into x2g; in the moab path prep_glc_calc_l2x_gx_moab already + ! writes the downscaled results directly into the x2g tag names on the glc + ! mesh, so there is nothing left to do here besides optional debug output. + ! +#ifdef MOABDEBUG + use iMOAB, only : iMOAB_WriteMesh + use seq_comm_mct, only : num_moab_exports + use iso_c_binding, only : C_NULL_CHAR +#endif + ! + ! Arguments + type(seq_infodata_type) , intent(in) :: infodata + character(len=*) , intent(in) :: timer_mrg + ! + ! Local Variables +#ifdef MOABDEBUG + integer :: ierr + character*32 :: outfile, wopts, lnum +#endif + character(*), parameter :: subname = '(prep_glc_mrg_lnd_moab)' + !--------------------------------------------------------------- + + call t_drvstartf (trim(timer_mrg), barrier=mpicom_CPLID) +#ifdef MOABDEBUG + if (mbgxid .ge. 0 ) then ! we are on coupler pes, for sure + write(lnum,"(I0.2)") num_moab_exports + outfile = 'GlcCplAftMrg'//trim(lnum)//'.h5m'//C_NULL_CHAR + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbgxid, trim(outfile), trim(wopts)) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' error in writing glc mesh after merge') + endif + endif +#endif + call t_drvstopf (trim(timer_mrg)) + + end subroutine prep_glc_mrg_lnd_moab + + !================================================================================================ + + function prep_glc_get_l2gacc_lm() + real(r8), pointer :: prep_glc_get_l2gacc_lm(:,:) + prep_glc_get_l2gacc_lm => l2gacc_lm + end function prep_glc_get_l2gacc_lm + + function prep_glc_get_l2gacc_lm_cnt() + integer, pointer :: prep_glc_get_l2gacc_lm_cnt + prep_glc_get_l2gacc_lm_cnt => l2gacc_lm_cnt + end function prep_glc_get_l2gacc_lm_cnt + + function prep_glc_get_sharedFieldsLndGlc() + character(CXX) :: prep_glc_get_sharedFieldsLndGlc + prep_glc_get_sharedFieldsLndGlc = sharedFieldsLndGlc + end function prep_glc_get_sharedFieldsLndGlc + + !================================================================================================ + function prep_glc_get_l2x_gx() type(mct_aVect), pointer :: prep_glc_get_l2x_gx(:) prep_glc_get_l2x_gx => l2x_gx(:) @@ -1431,6 +2340,11 @@ function prep_glc_get_l2gacc_lx_cnt() prep_glc_get_l2gacc_lx_cnt => l2gacc_lx_cnt end function prep_glc_get_l2gacc_lx_cnt + function prep_glc_get_l2gacc_lx_cnt_avg() + integer, pointer :: prep_glc_get_l2gacc_lx_cnt_avg + prep_glc_get_l2gacc_lx_cnt_avg => l2gacc_lx_cnt_avg + end function prep_glc_get_l2gacc_lx_cnt_avg + function prep_glc_get_o2x_gx() type(mct_aVect), pointer :: prep_glc_get_o2x_gx(:) prep_glc_get_o2x_gx => o2x_gx(:) @@ -1456,15 +2370,15 @@ function prep_glc_get_mapper_Fl2g() prep_glc_get_mapper_Fl2g => mapper_Fl2g end function prep_glc_get_mapper_Fl2g - function prep_glc_get_mapper_So2g() - type(seq_map), pointer :: prep_glc_get_mapper_So2g - prep_glc_get_mapper_So2g=> mapper_So2g - end function prep_glc_get_mapper_So2g + function prep_glc_get_mapper_So2g_shelf() + type(seq_map), pointer :: prep_glc_get_mapper_So2g_shelf + prep_glc_get_mapper_So2g_shelf=> mapper_So2g_shelf + end function prep_glc_get_mapper_So2g_shelf - function prep_glc_get_mapper_Fo2g() - type(seq_map), pointer :: prep_glc_get_mapper_Fo2g - prep_glc_get_mapper_Fo2g=> mapper_Fo2g - end function prep_glc_get_mapper_Fo2g + function prep_glc_get_mapper_Fo2g_shelf() + type(seq_map), pointer :: prep_glc_get_mapper_Fo2g_shelf + prep_glc_get_mapper_Fo2g_shelf=> mapper_Fo2g_shelf + end function prep_glc_get_mapper_Fo2g_shelf !*********************************************************************** ! diff --git a/driver-moab/main/prep_lnd_mod.F90 b/driver-moab/main/prep_lnd_mod.F90 index c81951b686d2..efbbea65339f 100644 --- a/driver-moab/main/prep_lnd_mod.F90 +++ b/driver-moab/main/prep_lnd_mod.F90 @@ -19,6 +19,8 @@ module prep_lnd_mod 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: mbgxid ! iMOAB id for glc migrated mesh to coupler pes + use seq_comm_mct, only: mbintxgl ! iMOAB id for read map between glc 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 @@ -34,6 +36,7 @@ module prep_lnd_mod use component_type_mod, only: component_get_x2c_cx, component_get_c2x_cx use component_type_mod, only: lnd, atm, rof, glc use map_glc2lnd_mod , only: map_glc2lnd_ec + use map_glc2lnd_mod , only: map_glc2lnd_ec_moab use iso_c_binding use iMOAB , only: iMOAB_ComputeCommGraph, iMOAB_ComputeMeshIntersectionOnSphere, & iMOAB_ComputeScalarProjectionWeights, iMOAB_DefineTagStorage, iMOAB_RegisterApplication, & @@ -153,7 +156,7 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln ! MOAB stuff integer :: ierr, idintx, rank character*32 :: appname - character*32 :: dm1, dm2, dofnameS, dofnameT, wgtIdFr2l, wgtIdFa2l, wgtIdSa2l + character*32 :: dm1, dm2, dofnameS, dofnameT, wgtIdFr2l, wgtIdFa2l, wgtIdSa2l, wgtIdFg2l integer :: orderS, orderT, volumetric, noConserve, validate, fInverseDistanceMap integer :: fNoBubble, monotonicity ! will do comm graph over coupler PES, in 2-hop strategy @@ -199,6 +202,8 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln wgtIdFr2l = 'flux_r2l' wgtIdFa2l = 'flux_a2l' wgtIdSa2l = 'scalar_a2l' + wgtIdFg2l = 'flux_g2l' ! must match the identifier used by prep_glc_init + compute_maps_online_r2l = cpl_compute_maps_online ! read from disk or compute online compute_maps_online_a2l = cpl_compute_maps_online ! read from disk or compute online ! compute_maps_online_a2l = .false. ! Explicitly force read from disk @@ -666,23 +671,85 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln call shr_sys_flush(logunit) if (glc_c2_lnd) then - if (iamroot_CPLID) then - write(logunit,*) ' ' - write(logunit,F00) 'Initializing mapper_Sg2l' - end if - call seq_map_init_rcfile(mapper_Sg2l, glc(1), lnd(1), & - 'seq_maps.rc','glc2lnd_smapname:','glc2lnd_smaptype:',samegrid_lg, & - 'mapper_Sg2l initialization',esmf_map_flag) - - if (iamroot_CPLID) then - write(logunit,*) ' ' - write(logunit,F00) 'Initializing mapper_Fg2l' - end if - call seq_map_init_rcfile(mapper_Fg2l, glc(1), lnd(1), & - 'seq_maps.rc','glc2lnd_fmapname:','glc2lnd_fmaptype:',samegrid_lg, & - 'mapper_Fg2l initialization',esmf_map_flag) + ! the mct sparse-matrix init (seq_map_init_rcfile) is not usable in + ! driver-moab (coupler-side gsmaps are not populated); the mappers get + ! their real (moab) context below + call seq_map_mapinit(mapper_Sg2l, mpicom_CPLID) + call seq_map_mapinit(mapper_Fg2l, mpicom_CPLID) call prep_lnd_set_glc2lnd_fields() + + ! moab context for the glc->lnd conservative map; read from the map file + if ((mbgxid .ge. 0) .and. (mblxid .ge. 0)) then + + if (samegrid_lg) then + call shr_sys_abort(subname// & + ' ERROR: moab glc->lnd coupling requires distinct lnd and glc grids with map files') + endif + + if (iamroot_CPLID) then + write(logunit,*) ' ' + write(logunit,F00) 'Initializing MOAB mapper_Fg2l' + end if + appname = "GLC_LND_COU" + ! unique external id for the moab app holding the glc->lnd read map + idintx = 100*glc(1)%cplcompid + lnd(1)%cplcompid + ierr = iMOAB_RegisterApplication(trim(appname)//C_NULL_CHAR, mpicom_CPLID, idintx, mbintxgl) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in registering glc lnd map app' + call shr_sys_abort(subname//' ERROR in registering glc lnd map app') + endif + + mapper_Fg2l%src_mbid = mbgxid + mapper_Fg2l%tgt_mbid = mblxid + mapper_Fg2l%intx_mbid = mbintxgl + mapper_Fg2l%src_context = glc(1)%cplcompid + mapper_Fg2l%intx_context = idintx + mapper_Fg2l%weight_identifier = wgtIdFg2l + mapper_Fg2l%mbname = 'mapper_Fg2l' + type1 = 3 ! FV-FV + arearead = 0 ! no need for areas + call moab_map_init_rcfile( mapper_Fg2l, type1, & + 'seq_maps.rc', 'glc2lnd_fmapname:', 'glc2lnd_fmaptype:', samegrid_lg, & + arearead, wgtIdFg2l, 'mapper_Fg2l MOAB initialization', esmf_map_flag) + + call seq_comm_getinfo(CPLID ,mpigrp=mpigrp_CPLID) + ierr = iMOAB_MigrateMapMesh (mbgxid, mbintxgl, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, glc(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating glc mesh for map glc 2 lnd' + call shr_sys_abort(subname//' ERROR in migrating glc mesh for map glc 2 lnd') + endif + ierr = iMOAB_ComputeCommGraph( mbgxid, mbintxgl, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, & + type1, type1, glc(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, GLC-LND' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, GLC-LND') + endif + + ! The per-EC destination tags exist on the land mesh (part of x2l fields). + ! Define the same names on the glc mesh, where they stage the pre-multiplied + ! numerators for the elevation-class mapping (map_glc2lnd_ec_moab), and + ! define the scratch tag used for the icemask normalization on both meshes. + tagtype = 1 ! dense, double + numco = 1 + tagname = trim(seq_flds_x2l_fields_from_glc)//C_NULL_CHAR + ierr = iMOAB_DefineTagStorage(mbgxid, tagname, tagtype, numco, tagindex ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR defining glc per-EC staging tags on the glc mesh') + endif + tagname = 'Sg_icemsk_num'//C_NULL_CHAR + ierr = iMOAB_DefineTagStorage(mbgxid, tagname, tagtype, numco, tagindex ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR defining Sg_icemsk_num on the glc mesh') + endif + ierr = iMOAB_DefineTagStorage(mblxid, tagname, tagtype, numco, tagindex ) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR defining Sg_icemsk_num on the land mesh') + endif + + endif ! mbgxid and mblxid + endif call shr_sys_flush(logunit) @@ -950,15 +1017,14 @@ subroutine prep_lnd_calc_g2x_lx(timer) call seq_map_map(mapper_Fg2l, g2x_gx, g2x_lx(egi), & fldlist = glc2lnd_non_ec_fields, norm=.true.) - ! Map fields that are separated by elevation class on the land grid - call map_glc2lnd_ec( & - g2x_g = g2x_gx, & + ! Map fields that are separated by elevation class on the land grid; the moab + ! version operates on the coupler mesh tags (result lands directly in the + ! per-EC x2l tag names on the land mesh) and returns early without moab context + call map_glc2lnd_ec_moab(mapper_Fg2l, & frac_field = glc_frac_field, & topo_field = glc_topo_field, & icemask_field = glc_icemask_field, & - extra_fields = glc2lnd_ec_extra_fields, & - mapper = mapper_Fg2l, & - g2x_l = g2x_lx(egi)) + extra_fields = glc2lnd_ec_extra_fields) enddo call t_drvstopf (trim(timer)) diff --git a/driver-moab/main/seq_rest_mod.F90 b/driver-moab/main/seq_rest_mod.F90 index 1c8ae005f391..d65abfecef3d 100644 --- a/driver-moab/main/seq_rest_mod.F90 +++ b/driver-moab/main/seq_rest_mod.F90 @@ -75,12 +75,16 @@ module seq_rest_mod use seq_flds_mod, only: seq_flds_a2x_fields, seq_flds_xao_fields, seq_flds_o2x_fields, seq_flds_x2o_fields use seq_flds_mod, only: seq_flds_i2x_fields, seq_flds_r2x_fields + use seq_flds_mod, only: seq_flds_g2x_fields use prep_rof_mod, only: prep_rof_get_o2racc_om ! return a pointer to a moab matrix use prep_rof_mod, only: prep_rof_get_l2racc_lm_cnt use prep_rof_mod, only: prep_rof_get_l2racc_lm use prep_rof_mod, only: prep_rof_get_sharedFieldsLndRof + use prep_glc_mod, only: prep_glc_get_l2gacc_lm_cnt + use prep_glc_mod, only: prep_glc_get_l2gacc_lm + use prep_glc_mod, only: prep_glc_get_sharedFieldsLndGlc implicit none private @@ -163,7 +167,7 @@ module seq_rest_mod subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) - use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid ! coupler side instances + use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid, mbgxid ! coupler side instances use seq_comm_mct , only: num_moab_exports ! it is used only as a counter for moab h5m files use seq_comm_mct, only: atm_pg_active ! whether the atm/lnd mesh is cells (pg2/FV) or a point cloud (np4) use iMOAB, only: iMOAB_GetGlobalInfo @@ -185,6 +189,7 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) integer (in), pointer :: x2oacc_om_cnt ! replacement, moab version for x2oacc_ox_cnt integer (in), pointer :: l2racc_lm_cnt + integer (in), pointer :: l2gacc_lm_cnt integer (in) :: nx_lnd ! will be used if land and atm are on same grid integer (in) :: ngv, nge ! global num vertices / elements from iMOAB_GetGlobalInfo integer (in) :: ierr @@ -192,6 +197,7 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) real(r8), dimension(:,:), pointer :: p_x2oacc_om real(r8), dimension(:,:), pointer :: p_o2racc_om real(r8), dimension(:,:), pointer :: p_l2racc_lm + real(r8), dimension(:,:), pointer :: p_l2gacc_lm character(len=*), parameter :: subname = "(seq_rest_mb_read) " @@ -311,15 +317,37 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) matrix = p_o2racc_om ) call seq_io_read(moab_rest_file, o2racc_om_cnt, 'o2racc_ox_cnt') end if -! MOABTODO -! if (lnd_present .and. glc_prognostic) then -! -! l2gacc_lx => prep_glc_get_l2gacc_lx() -! l2gacc_lx_cnt => prep_glc_get_l2gacc_lx_cnt() -! call seq_io_read(rest_file, gsmap, l2gacc_lx, 'l2gacc_lx') -! call seq_io_read(rest_file, l2gacc_lx_cnt ,'l2gacc_lx_cnt') -! end if + if (lnd_present .and. glc_prognostic) then + tagname = prep_glc_get_sharedFieldsLndGlc() + l2gacc_lm_cnt => prep_glc_get_l2gacc_lm_cnt() + p_l2gacc_lm => prep_glc_get_l2gacc_lm() + if(samegrid_al) then + ! land is on the atm mesh; max global land id comes from atm. + ierr = iMOAB_GetGlobalInfo(mbaxid, ngv, nge) + if (atm_pg_active) then + nx_lnd = nge ! atm mesh is cells + else + nx_lnd = ngv ! atm mesh is a point cloud + endif + call seq_io_read(moab_rest_file, mblxid, 'l2gacc_lx', & + trim(tagname), & + matrix = p_l2gacc_lm, nx=nx_lnd) + else if(samegrid_lr) then + ! land is on the rof mesh; max global land id comes from rof. + ierr = iMOAB_GetGlobalInfo(mbrxid, ngv, nge) + nx_lnd = nge + call seq_io_read(moab_rest_file, mblxid, 'l2gacc_lx', & + trim(tagname), & + matrix = p_l2gacc_lm, nx=nx_lnd) + else + call seq_io_read(moab_rest_file, mblxid, 'l2gacc_lx', & + trim(tagname), & + matrix = p_l2gacc_lm ) + endif + call seq_io_read(rest_file, l2gacc_lm_cnt ,'l2gacc_lx_cnt') + end if +! MOABTODO (deferred with the rest of the ocn<->glc coupling) ! if (ocn_c2_glcshelf) then ! gsmap => component_get_gsmap_cx(glc(1)) ! x2gacc_gx => prep_glc_get_x2gacc_gx() @@ -359,12 +387,13 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) call seq_io_read(moab_rest_file, mbrxid, 'r2x_rx', & trim(seq_flds_r2x_fields) ) endif + if (glc_present) then + call seq_io_read(moab_rest_file, mbgxid, 'fractions_gx', & + 'gfrac:lfrac') ! fraclist_g = 'gfrac:lfrac' + call seq_io_read(moab_rest_file, mbgxid, 'g2x_gx', & + trim(seq_flds_g2x_fields) ) + endif !MOABTODO -! if (glc_present) then -! gsmap => component_get_gsmap_cx(glc(1)) -! call seq_io_read(rest_file, gsmap, fractions_gx, 'fractions_gx') -! call seq_io_read(rest_file, glc, 'c2x', 'g2x_gx') -! endif ! if (wav_present) then ! gsmap => component_get_gsmap_cx(wav(1)) ! call seq_io_read(rest_file, gsmap, fractions_wx, 'fractions_wx') @@ -425,7 +454,7 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & atm, lnd, ice, ocn, rof, glc, wav, esp, iac, & tag, samegrid_al, samegrid_lr, rest_file) - use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid ! coupler side instances + use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid, mbgxid ! coupler side instances use seq_comm_mct , only: num_moab_exports ! it is used only as a counter for moab h5m files use seq_comm_mct, only: atm_pg_active ! whether the atm/lnd mesh is cells (pg2/FV) or a point cloud (np4) use iMOAB, only: iMOAB_GetGlobalInfo @@ -475,6 +504,7 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & integer (in), pointer :: x2oacc_om_cnt ! replacement, moab version for x2oacc_ox_cnt integer (in), pointer :: l2racc_lm_cnt + integer (in), pointer :: l2gacc_lm_cnt integer (in) :: nx_lnd ! will be used if land and atm are on same grid integer (in) :: lnd_nx, lnd_ny ! global land grid dims (from infodata), for non-samegrid land restart I/O integer (in) :: ngv, nge ! global num vertices / elements from iMOAB_GetGlobalInfo @@ -483,6 +513,7 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & real(r8), dimension(:,:), pointer :: p_x2oacc_om real(r8), dimension(:,:), pointer :: p_o2racc_om real(r8), dimension(:,:), pointer :: p_l2racc_lm + real(r8), dimension(:,:), pointer :: p_l2gacc_lm character(len=*),parameter :: subname = "(seq_rest_mb_write) " !------------------------------------------------------------------------------- @@ -723,16 +754,38 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & call seq_io_write(rest_file, o2racc_om_cnt, 'o2racc_ox_cnt', & whead=whead, wdata=wdata) end if -! MOABTODO -! if (lnd_present .and. glc_prognostic) then -! gsmap => component_get_gsmap_cx(lnd(1)) -! l2gacc_lx => prep_glc_get_l2gacc_lx() -! l2gacc_lx_cnt => prep_glc_get_l2gacc_lx_cnt() -! call seq_io_write(rest_file, gsmap, l2gacc_lx, 'l2gacc_lx', & -! whead=whead, wdata=wdata) -! call seq_io_write(rest_file, l2gacc_lx_cnt, 'l2gacc_lx_cnt', & -! whead=whead, wdata=wdata) -! end if + if (lnd_present .and. glc_prognostic) then + tagname = prep_glc_get_sharedFieldsLndGlc() + l2gacc_lm_cnt => prep_glc_get_l2gacc_lm_cnt() + p_l2gacc_lm => prep_glc_get_l2gacc_lm() + if(samegrid_al) then + ! land is on the atm mesh; max global land id comes from atm. + ierr = iMOAB_GetGlobalInfo(mbaxid, ngv, nge) + if (atm_pg_active) then + nx_lnd = nge ! atm mesh is cells + else + nx_lnd = ngv ! atm mesh is a point cloud + endif + call seq_io_write(rest_file, mblxid, 'l2gacc_lx', & + trim(tagname), & + whead=whead, wdata=wdata, matrix = p_l2gacc_lm, nx=nx_lnd) + else if(samegrid_lr) then + ! land is on the rof mesh; max global land id comes from rof. + ierr = iMOAB_GetGlobalInfo(mbrxid, ngv, nge) + nx_lnd = nge + call seq_io_write(rest_file, mblxid, 'l2gacc_lx', & + trim(tagname), & + whead=whead, wdata=wdata, matrix = p_l2gacc_lm, nx=nx_lnd) + else + ! land on its own (masked) grid: size PIO global dim from infodata land dims + call seq_io_write(rest_file, mblxid, 'l2gacc_lx', & + trim(tagname), & + whead=whead, wdata=wdata, matrix = p_l2gacc_lm, nx=lnd_nx, ny=lnd_ny ) + endif + call seq_io_write(rest_file, l2gacc_lm_cnt, 'l2gacc_lx_cnt', & + whead=whead, wdata=wdata) + end if +! MOABTODO (deferred with the rest of the ocn<->glc coupling) ! if (ocn_c2_glcshelf) then ! gsmap => component_get_gsmap_cx(glc(1)) ! x2gacc_gx => prep_glc_get_x2gacc_gx() @@ -781,14 +834,15 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & trim(seq_flds_r2x_fields), & whead=whead, wdata=wdata) endif + if (glc_present) then + call seq_io_write(rest_file, mbgxid, 'fractions_gx', & + 'gfrac:lfrac', & ! fraclist_g = 'gfrac:lfrac' + whead=whead, wdata=wdata) + call seq_io_write(rest_file, mbgxid, 'g2x_gx', & + trim(seq_flds_g2x_fields), & + whead=whead, wdata=wdata) + endif ! MOABTODO -! if (glc_present) then -! gsmap => component_get_gsmap_cx(glc(1)) -! call seq_io_write(rest_file, gsmap, fractions_gx, 'fractions_gx', & -! whead=whead, wdata=wdata) -! call seq_io_write(rest_file, glc, 'c2x', 'g2x_gx', & -! whead=whead, wdata=wdata) -! endif ! if (wav_present) then ! gsmap => component_get_gsmap_cx(wav(1)) ! call seq_io_write(rest_file, gsmap, fractions_wx, 'fractions_wx', & @@ -818,7 +872,7 @@ end subroutine seq_rest_mb_write #ifdef MOABDEBUG subroutine write_moab_state ( before_reading ) ! debug, write files - use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid ! coupler side instances + use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid, mbgxid ! coupler side instances use seq_comm_mct, only: num_moab_exports use iso_c_binding use iMOAB, only: iMOAB_WriteMesh From cf429914744952082bbcb2871e935a040707449e Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Sat, 22 Aug 2026 00:24:55 -0500 Subject: [PATCH 6/7] Match driver-mct bit for bit in the moab lnd-glc history output Store the per-elevation-class arrays used by the lnd->glc vertical interpolation in default real precision, as the mct map_ice_covered does, instead of keeping them in r8. The demotion is what mct does today, so reproducing it is what makes the downscaled fields (notably x2g_Sl_tsrf) agree; with the same weights, the iMOAB ApplyScalarProjectionWeights result then matches mct_sMat_avMult exactly. Stop setting the averaged lnd->glc accumulation back into the coupler land mesh tags. Instead, prep_glc_calc_l2x_gx_moab saves the instantaneous l2x values, stages the averages into the tags for the duration of the batched map, and restores them right after, so the coupler history keeps showing instantaneous l2x fields the way the mct driver does. Add an optional matrix argument to seq_hist_writeaux so the l2x1yrg auxiliary file, whose data is the accumulator rather than the tags, can still be written from the array. Do the glc->lnd elevation-class normalization as a reciprocal multiply with a zero denominator passing through as a multiply by zero, matching the mct seq_map_avNormArr idiom, which removes the ~1e-15 differences in the per-class frac and topo fields. Validated on anlgce-ub22 (gnu) by comparing 2-day IGELM_MLI twins with cprnc: of 286 coupler history fields, all are bit for bit except x2g_Flgl_qice (~4e-15 normalized) and the two MALI g2x fields that echo it. That residual is a uniform 1-ulp factor across all 17352 glc cells, i.e. shr_mpi_sum ordering in the SMB renormalization over the different moab and mct decompositions, not a difference in the coupling itself; glc_renormalize_smb='off' gives strict bit-for-bit. Co-Authored-By: Claude Opus 5 (1M context) --- driver-moab/main/cime_comp_mod.F90 | 7 +++- driver-moab/main/map_glc2lnd_mod.F90 | 22 +++++------ driver-moab/main/map_lnd2glc_mod.F90 | 9 +++-- driver-moab/main/prep_glc_mod.F90 | 58 ++++++++++++++++------------ driver-moab/main/seq_hist_mod.F90 | 17 +++++++- 5 files changed, 70 insertions(+), 43 deletions(-) diff --git a/driver-moab/main/cime_comp_mod.F90 b/driver-moab/main/cime_comp_mod.F90 index 4274634489bb..6871f628d073 100644 --- a/driver-moab/main/cime_comp_mod.F90 +++ b/driver-moab/main/cime_comp_mod.F90 @@ -5146,6 +5146,7 @@ subroutine cime_run_write_history(lnd2glc_averaged_now) type(ESMF_Time) :: etime_curr ! Current model time real(r8) :: tbnds1_offset ! Time offset for call to seq_hist_writeaux + real(r8), pointer :: p_l2gacc_lm(:,:) ! averaged lnd->glc accumulator (moab) if (iamin_CPLID) then @@ -5336,6 +5337,9 @@ subroutine cime_run_write_history(lnd2glc_averaged_now) years_offset = -1) call t_startf('CPL:seq_hist_writeaux-l2x1yrg') + ! the averaged data lives in the prep_glc accumulator array, not in the + ! (instantaneous) l2x tags; the attribute vector supplies the field list + p_l2gacc_lm => prep_glc_get_l2gacc_lm() do eli = 1,num_inst_lnd inst_suffix = component_get_suffix(lnd(eli)) ! Use yr_offset=-1 so the file with fields from year 1 has time stamp @@ -5344,7 +5348,8 @@ subroutine cime_run_write_history(lnd2glc_averaged_now) aname='l2x1yr_glc',dname='doml',inst_suffix=trim(inst_suffix), & nx=lnd_nx, ny=lnd_ny, nt=1, write_now=.true., & tbnds1_offset = tbnds1_offset, yr_offset=-1, & - av_to_write=prep_glc_get_l2gacc_lx_one_instance(eli)) + av_to_write=prep_glc_get_l2gacc_lx_one_instance(eli), & + matrix=p_l2gacc_lm) enddo call t_stopf('CPL:seq_hist_writeaux-l2x1yrg') diff --git a/driver-moab/main/map_glc2lnd_mod.F90 b/driver-moab/main/map_glc2lnd_mod.F90 index 7072ea1335b8..08dc226cc1ad 100644 --- a/driver-moab/main/map_glc2lnd_mod.F90 +++ b/driver-moab/main/map_glc2lnd_mod.F90 @@ -417,26 +417,24 @@ subroutine map_glc2lnd_ec_moab(mapper, & ktopo0 = kfrac0 + 1 topo_virtual = glc_mean_elevation_virtual(n) do i = 1, lsize_l + ! The divisions are done as reciprocal multiplies (with a zero denominator + ! passing through as a multiply by zero), exactly like the normalization in + ! the mct seq_map_avNormArr, so the results match the mct driver bit for bit. ! frac_n_l = M(w_n) / M(icemask) denom = num_l(i, kmask) if (denom /= 0.0_r8) then - out_l(i, kfrac0) = num_l(i, kfrac0) / denom - else - out_l(i, kfrac0) = 0.0_r8 + denom = 1.0_r8/denom end if + out_l(i, kfrac0) = num_l(i, kfrac0) * denom ! field_n_l = M(field*w_n) / M(w_n) denom = num_l(i, kfrac0) if (denom /= 0.0_r8) then - out_l(i, ktopo0) = num_l(i, ktopo0) / denom - do k = 1, nextra - out_l(i, ktopo0+k) = num_l(i, ktopo0+k) / denom - end do - else - out_l(i, ktopo0) = 0.0_r8 - do k = 1, nextra - out_l(i, ktopo0+k) = 0.0_r8 - end do + denom = 1.0_r8/denom end if + out_l(i, ktopo0) = num_l(i, ktopo0) * denom + do k = 1, nextra + out_l(i, ktopo0+k) = num_l(i, ktopo0+k) * denom + end do ! set the topo field for virtual columns (no contributing glc cells) if (out_l(i, kfrac0) <= 0.0_r8) then out_l(i, ktopo0) = topo_virtual diff --git a/driver-moab/main/map_lnd2glc_mod.F90 b/driver-moab/main/map_lnd2glc_mod.F90 index 60bf7b9efe71..a0442e6b5fe8 100644 --- a/driver-moab/main/map_lnd2glc_mod.F90 +++ b/driver-moab/main/map_lnd2glc_mod.F90 @@ -212,13 +212,14 @@ subroutine map_lnd2glc_vertical_interp(topo_g, topo_g_EC, data_g_EC, data_g_bare ! All arrays are on the glc grid decomposition. data_g_EC and topo_g_EC hold the ! horizontally-mapped field and Sl_topo for elevation classes 1..nEC. ! - ! Note: the mct path (map_ice_covered) stores the per-EC arrays in default real - ! precision; here everything stays in r8, a roundoff-level difference. + ! Note: the per-EC arrays are stored in default real precision, exactly like the + ! mct path (map_ice_covered), so that the interpolation reproduces the mct + ! answers bit for bit (see E3SM issue #8657 about the demotion itself). ! ! !ARGUMENTS: real(r8), intent(in) :: topo_g(:) ! ice topographic height on the glc grid - real(r8), intent(in) :: topo_g_EC(:,:) ! mapped per-EC topo (lsize_g, nEC) - real(r8), intent(in) :: data_g_EC(:,:) ! mapped per-EC field (lsize_g, nEC) + real , intent(in) :: topo_g_EC(:,:) ! mapped per-EC topo (lsize_g, nEC), default real like mct + real , intent(in) :: data_g_EC(:,:) ! mapped per-EC field (lsize_g, nEC), default real like mct real(r8), intent(in) :: data_g_bareland(:) ! mapped bare-land (EC 0) field integer , intent(in) :: glc_elevclass(:) ! elevation class of each glc point (0 = bare) real(r8), intent(out) :: data_g(:) ! result on the glc grid diff --git a/driver-moab/main/prep_glc_mod.F90 b/driver-moab/main/prep_glc_mod.F90 index 8828b060e849..876580c4bd65 100644 --- a/driver-moab/main/prep_glc_mod.F90 +++ b/driver-moab/main/prep_glc_mod.F90 @@ -747,24 +747,21 @@ subroutine prep_glc_accum_avg_moab(timer, lnd2glc_averaged_now) !--------------------------------------------------------------- ! Description - ! Finalize the accumulation of the land forcing for glc, and set the averaged - ! values back into the coupler land mesh tags, from where the lnd->glc maps - ! read them (prep_rof set-back pattern). + ! Finalize the accumulation of the land forcing for glc. The averaged values + ! stay in the l2gacc_lm array (matching the mct l2gacc_lx attribute vector); + ! prep_glc_calc_l2x_gx_moab stages them into the land mesh tags only for the + ! duration of the horizontal map, so the l2x tags keep their instantaneous + ! values for the coupler history. ! ! Note: the mct version also averages the ocn accumulation (x2gacc); that part ! of the moab port is deferred with the rest of the ocn<->glc coupling. ! - use iMOAB, only : iMOAB_SetDoubleTagStorage - use iso_c_binding, only : C_NULL_CHAR - ! ! Arguments character(len=*), intent(in) :: timer logical, intent(inout) :: lnd2glc_averaged_now ! Set to .true. if lnd2glc averages were taken this timestep (otherwise left unchanged) ! ! Local Variables - integer :: ierr, ent_type, arrsize real(r8) :: ravg ! averaging factor - character(CXX) :: tagname character(*), parameter :: subname = '(prep_glc_accum_avg_moab)' !--------------------------------------------------------------- @@ -780,15 +777,6 @@ subroutine prep_glc_accum_avg_moab(timer, lnd2glc_averaged_now) end if l2gacc_lm_cnt_avg = l2gacc_lm_cnt l2gacc_lm_cnt = 0 - - ! set the averaged values back into the land mesh tags - tagname = trim(sharedFieldsLndGlc)//C_NULL_CHAR - arrsize = nflds_lg * lsize_lm - ent_type = 1 ! cells - ierr = iMOAB_SetDoubleTagStorage ( mblxid, tagname, arrsize, ent_type, l2gacc_lm ) - if (ierr .ne. 0) then - call shr_sys_abort(subname//' error in setting accumulated per-EC lnd fields on land mesh') - endif call t_drvstopf (trim(timer)) end subroutine prep_glc_accum_avg_moab @@ -1920,8 +1908,10 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) real(r8), allocatable :: data_lg(:,:) ! all mapped per-EC fields on the glc mesh real(r8), allocatable :: glc_ice_covered(:) real(r8), allocatable :: glc_topo(:) - real(r8), allocatable :: topo_g_EC(:,:) - real(r8), allocatable :: data_g_EC(:,:) + ! the per-EC arrays are default real, matching the mct map_ice_covered storage, + ! so the vertical interpolation reproduces the mct answers bit for bit + real , allocatable :: topo_g_EC(:,:) + real , allocatable :: data_g_EC(:,:) real(r8), allocatable :: data_g_bare(:) real(r8), allocatable :: data_g(:) ! downscaled field on the glc mesh real(r8), allocatable :: area_g(:) @@ -1938,12 +1928,32 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) ent_type = 1 ! cells ! Horizontal map of all per-EC fields at once, weighted by lfrac and normalized. - ! The source tags hold the averaged accumulation, set back on the land mesh by - ! prep_glc_accum_avg_moab. The attribute vector arguments are metadata only. + ! The averaged accumulation (l2gacc_lm) is staged into the land mesh tags only + ! for the duration of the map; the instantaneous l2x values are saved first and + ! restored right after, so the coupler history still shows instantaneous l2x + ! fields exactly like the mct driver. The attribute vector arguments are + ! metadata only. + arrsize = nflds_lg * lsize_lm + tagname = trim(sharedFieldsLndGlc)//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage(mblxid, tagname, arrsize, ent_type, l2x_lm2) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR saving instantaneous per-EC lnd fields') + endif + ierr = iMOAB_SetDoubleTagStorage(mblxid, tagname, arrsize, ent_type, l2gacc_lm) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR staging averaged per-EC lnd fields') + endif + call seq_map_map(mapper_Sl2g, l2gacc_lx(1), l2x_gx(1), & fldlist=trim(sharedFieldsLndGlc), norm=.true., & avwts_s=fractions_lx(1), avwtsfld_s='lfrac') + ! restore the instantaneous l2x values on the land mesh + ierr = iMOAB_SetDoubleTagStorage(mblxid, tagname, arrsize, ent_type, l2x_lm2) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' ERROR restoring instantaneous per-EC lnd fields') + endif + ! fetch the mapped per-EC fields and the glc state needed for the vertical interpolation allocate(data_lg(lsize_gm, nflds_lg)) data_lg = 0._r8 @@ -1971,7 +1981,7 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) allocate(data_g(lsize_gm)) do ec = 1, nEC kf = mct_aVect_indexRA(l2gacc_lx(1), 'Sl_topo'//glc_elevclass_as_string(ec)) - topo_g_EC(:,ec) = data_lg(:,kf) + topo_g_EC(:,ec) = real(data_lg(:,kf)) end do ! area arrays needed for the qice conservation correction @@ -1993,7 +2003,7 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) do ec = 1, nEC kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(ec)) - data_g_EC(:,ec) = data_lg(:,kf) + data_g_EC(:,ec) = real(data_lg(:,kf)) end do kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(0)) data_g_bare(:) = data_lg(:,kf) @@ -2034,7 +2044,7 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) do ec = 1, nEC kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(ec)) - data_g_EC(:,ec) = data_lg(:,kf) + data_g_EC(:,ec) = real(data_lg(:,kf)) end do kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(0)) data_g_bare(:) = data_lg(:,kf) diff --git a/driver-moab/main/seq_hist_mod.F90 b/driver-moab/main/seq_hist_mod.F90 index 25c3d97cc482..cf5ef437743c 100644 --- a/driver-moab/main/seq_hist_mod.F90 +++ b/driver-moab/main/seq_hist_mod.F90 @@ -1085,7 +1085,7 @@ end subroutine seq_hist_writeavg !=============================================================================== subroutine seq_hist_writeaux(infodata, EClock_d, comp, flow, aname, dname, inst_suffix, & - nx, ny, nt, write_now, flds, tbnds1_offset, yr_offset, av_to_write) + nx, ny, nt, write_now, flds, tbnds1_offset, yr_offset, av_to_write, matrix) implicit none @@ -1124,10 +1124,16 @@ subroutine seq_hist_writeaux(infodata, EClock_d, comp, flow, aname, dname, inst_ integer , optional, intent(in) :: yr_offset ! If av_to_write is provided, use it to get the list of tags to write. - ! The data is always written from MOAB tags, not from the attribute vector. + ! The data is written from MOAB tags unless 'matrix' is provided. ! Otherwise, get the tag list from 'comp', based on 'flow'. type(mct_avect), target , optional, intent(in) :: av_to_write + ! If matrix is provided (non-averaged writes only), write the data directly + ! from this (numpts, nflds) array instead of reading the MOAB tags; the columns + ! must follow the tag-list order. Used for the l2x1yrg aux file, whose data + ! comes from the prep_glc accumulator rather than the instantaneous l2x tags. + real(r8), dimension(:,:), pointer, optional :: matrix + !--- local --- type(mct_avect), pointer :: av character(CL) :: case_name ! case name @@ -1422,6 +1428,13 @@ subroutine seq_hist_writeaux(infodata, EClock_d, comp, flow, aname, dname, inst_ nx=nx, ny=ny, nt=ncnt(found), & file_ind=found, & use_float=(.not. use_double)) + else if (present(matrix)) then + ! Non-averaged, all fields, data supplied by the caller + call seq_io_write(hist_file(found), mbxid, trim(aname), & + trim(tag_list), whead=whead, wdata=wdata, & + nx=nx, ny=ny, nt=ncnt(found), & + matrix=matrix, file_ind=found, & + use_float=(.not. use_double)) else ! Non-averaged, all fields: write from MOAB tags call seq_io_write(hist_file(found), mbxid, trim(aname), & From b604502d3e715e24283f0aabc2569c69a1fcec75 Mon Sep 17 00:00:00 2001 From: Robert Jacob Date: Mon, 31 Aug 2026 00:32:34 -0500 Subject: [PATCH 7/7] Keep the lnd-glc elevation-class downscaling in double precision The per-elevation-class staging arrays used by the lnd->glc vertical downscaling were declared with a bare real and the horizontally mapped r8 values were explicitly demoted into them, so every per-EC field value and per-EC topography was rounded to single precision before the vertical interpolation ran. Everything upstream and downstream of that interpolation is r8, and the driver is not built with -fdefault-real-8, so this cost the downscaled fields sent to GLC most of their precision. It looks like an oversight inherited from CESM rather than a memory optimization: the arrays are only lsize_g * nEC. Declare data_g_EC and topo_g_EC as real(r8) and drop the real() demotion, in map_ice_covered in both drivers. driver-moab needs it in two more places, because the moab path does not go through map_ice_covered: the topo_g_EC and data_g_EC dummy arguments of map_lnd2glc_vertical_interp, and the matching allocatables and demotions in prep_glc_calc_l2x_gx_moab that feed it. Answer-changing for compsets with prognostic GLC coupling. Measured on chrysalis (gnu, IGELM_MLI, Albany FO solver) by comparing driver-mct coupler history before and after: 27 of 286 fields move, the surface mass balance x2g_Flgl_qice by 1.0e-7 normalized against a single-precision epsilon of 1.2e-7, and x2g_Sl_tsrf by 8.8e-9. Both drivers are changed together so they stay matched to each other; the moab-vs-mct comparison improves from 9 differing fields to 7, with l2x_Sl_topo01 and x2l_Sg_topo01 becoming bit for bit. Fixes #8657 Co-Authored-By: Claude Opus 5 (1M context) --- driver-mct/main/map_lnd2glc_mod.F90 | 8 ++++---- driver-moab/main/map_lnd2glc_mod.F90 | 16 ++++++---------- driver-moab/main/prep_glc_mod.F90 | 12 +++++------- 3 files changed, 15 insertions(+), 21 deletions(-) diff --git a/driver-mct/main/map_lnd2glc_mod.F90 b/driver-mct/main/map_lnd2glc_mod.F90 index 1b418515e912..d3114fe251cc 100644 --- a/driver-mct/main/map_lnd2glc_mod.F90 +++ b/driver-mct/main/map_lnd2glc_mod.F90 @@ -348,8 +348,8 @@ subroutine map_ice_covered(l2x_l, landfrac_l, fieldname, & type(mct_aVect) :: l2x_g_temp ! temporary attribute vector holding the remapped fields for this elevation class real(r8), pointer :: tmp_field_g(:) ! must be a pointer to satisfy the MCT interface - real, pointer :: data_g_EC(:,:) ! remapped field in each glc cell, in each EC - real, pointer :: topo_g_EC(:,:) ! remapped topo in each glc cell, in each EC + real(r8), pointer :: data_g_EC(:,:) ! remapped field in each glc cell, in each EC + real(r8), pointer :: topo_g_EC(:,:) ! remapped topo in each glc cell, in each EC ! 1 is probably enough, but use 10 to be safe, in case the length of the delimiter ! changes @@ -416,9 +416,9 @@ subroutine map_ice_covered(l2x_l, landfrac_l, fieldname, & fieldname_ec = fieldname // elevclass_as_string toponame_ec = toponame // elevclass_as_string call mct_aVect_exportRattr(l2x_g_temp, fieldname_ec, tmp_field_g) - data_g_EC(:,ec) = real(tmp_field_g) + data_g_EC(:,ec) = tmp_field_g call mct_aVect_exportRattr(l2x_g_temp, toponame_ec, tmp_field_g) - topo_g_EC(:,ec) = real(tmp_field_g) + topo_g_EC(:,ec) = tmp_field_g enddo ! ------------------------------------------------------------------------ diff --git a/driver-moab/main/map_lnd2glc_mod.F90 b/driver-moab/main/map_lnd2glc_mod.F90 index a0442e6b5fe8..23b6e1b0da08 100644 --- a/driver-moab/main/map_lnd2glc_mod.F90 +++ b/driver-moab/main/map_lnd2glc_mod.F90 @@ -212,14 +212,10 @@ subroutine map_lnd2glc_vertical_interp(topo_g, topo_g_EC, data_g_EC, data_g_bare ! All arrays are on the glc grid decomposition. data_g_EC and topo_g_EC hold the ! horizontally-mapped field and Sl_topo for elevation classes 1..nEC. ! - ! Note: the per-EC arrays are stored in default real precision, exactly like the - ! mct path (map_ice_covered), so that the interpolation reproduces the mct - ! answers bit for bit (see E3SM issue #8657 about the demotion itself). - ! ! !ARGUMENTS: real(r8), intent(in) :: topo_g(:) ! ice topographic height on the glc grid - real , intent(in) :: topo_g_EC(:,:) ! mapped per-EC topo (lsize_g, nEC), default real like mct - real , intent(in) :: data_g_EC(:,:) ! mapped per-EC field (lsize_g, nEC), default real like mct + real(r8), intent(in) :: topo_g_EC(:,:) ! mapped per-EC topo (lsize_g, nEC) + real(r8), intent(in) :: data_g_EC(:,:) ! mapped per-EC field (lsize_g, nEC) real(r8), intent(in) :: data_g_bareland(:) ! mapped bare-land (EC 0) field integer , intent(in) :: glc_elevclass(:) ! elevation class of each glc point (0 = bare) real(r8), intent(out) :: data_g(:) ! result on the glc grid @@ -440,8 +436,8 @@ subroutine map_ice_covered(l2x_l, landfrac_l, fieldname, & type(mct_aVect) :: l2x_g_temp ! temporary attribute vector holding the remapped fields for this elevation class real(r8), pointer :: tmp_field_g(:) ! must be a pointer to satisfy the MCT interface - real, pointer :: data_g_EC(:,:) ! remapped field in each glc cell, in each EC - real, pointer :: topo_g_EC(:,:) ! remapped topo in each glc cell, in each EC + real(r8), pointer :: data_g_EC(:,:) ! remapped field in each glc cell, in each EC + real(r8), pointer :: topo_g_EC(:,:) ! remapped topo in each glc cell, in each EC ! 1 is probably enough, but use 10 to be safe, in case the length of the delimiter ! changes @@ -508,9 +504,9 @@ subroutine map_ice_covered(l2x_l, landfrac_l, fieldname, & fieldname_ec = fieldname // elevclass_as_string toponame_ec = toponame // elevclass_as_string call mct_aVect_exportRattr(l2x_g_temp, fieldname_ec, tmp_field_g) - data_g_EC(:,ec) = real(tmp_field_g) + data_g_EC(:,ec) = tmp_field_g call mct_aVect_exportRattr(l2x_g_temp, toponame_ec, tmp_field_g) - topo_g_EC(:,ec) = real(tmp_field_g) + topo_g_EC(:,ec) = tmp_field_g enddo ! ------------------------------------------------------------------------ diff --git a/driver-moab/main/prep_glc_mod.F90 b/driver-moab/main/prep_glc_mod.F90 index 876580c4bd65..24c849c62efd 100644 --- a/driver-moab/main/prep_glc_mod.F90 +++ b/driver-moab/main/prep_glc_mod.F90 @@ -1908,10 +1908,8 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) real(r8), allocatable :: data_lg(:,:) ! all mapped per-EC fields on the glc mesh real(r8), allocatable :: glc_ice_covered(:) real(r8), allocatable :: glc_topo(:) - ! the per-EC arrays are default real, matching the mct map_ice_covered storage, - ! so the vertical interpolation reproduces the mct answers bit for bit - real , allocatable :: topo_g_EC(:,:) - real , allocatable :: data_g_EC(:,:) + real(r8), allocatable :: topo_g_EC(:,:) + real(r8), allocatable :: data_g_EC(:,:) real(r8), allocatable :: data_g_bare(:) real(r8), allocatable :: data_g(:) ! downscaled field on the glc mesh real(r8), allocatable :: area_g(:) @@ -1981,7 +1979,7 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) allocate(data_g(lsize_gm)) do ec = 1, nEC kf = mct_aVect_indexRA(l2gacc_lx(1), 'Sl_topo'//glc_elevclass_as_string(ec)) - topo_g_EC(:,ec) = real(data_lg(:,kf)) + topo_g_EC(:,ec) = data_lg(:,kf) end do ! area arrays needed for the qice conservation correction @@ -2003,7 +2001,7 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) do ec = 1, nEC kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(ec)) - data_g_EC(:,ec) = real(data_lg(:,kf)) + data_g_EC(:,ec) = data_lg(:,kf) end do kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(0)) data_g_bare(:) = data_lg(:,kf) @@ -2044,7 +2042,7 @@ subroutine prep_glc_calc_l2x_gx_moab(fractions_lx, timer) do ec = 1, nEC kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(ec)) - data_g_EC(:,ec) = real(data_lg(:,kf)) + data_g_EC(:,ec) = data_lg(:,kf) end do kf = mct_aVect_indexRA(l2gacc_lx(1), trim(fieldname)//glc_elevclass_as_string(0)) data_g_bare(:) = data_lg(:,kf)