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}) 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/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 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/cime_comp_mod.F90 b/driver-moab/main/cime_comp_mod.F90 index 594a2f5456da..6871f628d073 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.) @@ -5123,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 @@ -5313,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 @@ -5321,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/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/map_glc2lnd_mod.F90 b/driver-moab/main/map_glc2lnd_mod.F90 index 4ae638ed71cf..08dc226cc1ad 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,258 @@ 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 + ! 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 + 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 + 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 + 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..23b6e1b0da08 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,89 @@ 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. + ! + ! !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) ! @@ -348,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 @@ -416,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 339acaaaac04..24c849c62efd 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,235 @@ 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. 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. + ! + ! 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 + real(r8) :: ravg ! averaging factor + 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 + 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 +1012,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 +1076,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 +1094,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 +1266,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 +1368,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 +1666,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 +1851,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 +1877,456 @@ 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 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 + 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 +2348,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 +2378,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_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..cf5ef437743c 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 @@ -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), & 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 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)