diff --git a/cmake/FindOASIS3MCT.cmake b/cmake/FindOASIS3MCT.cmake new file mode 100644 index 0000000000..e09e08f02c --- /dev/null +++ b/cmake/FindOASIS3MCT.cmake @@ -0,0 +1,47 @@ +# Finds the OASIS3-MCT library. +# This will define the following variables: +# +# OASIS3MCT_FOUND - True if the system has the MCT library. +# OASIS3MCT_INCLUDE_DIRS - Include directories needed to use MCT. +# OASIS3MCT_LIBRARIES - Libraries needed to link to MCT. +# OASIS3MCT::OASIS3MCT - A target to use with `target_link_libraries`. + + +find_path(PSMILE_INCLUDE_DIR NAMES mod_oasis.mod PATH_SUFFIXES psmile.MPI1 lib/psmile.MPI1 build/lib/psmile.MPI1) +find_path(MCT_INCLUDE_DIR NAMES oas_mct_mod.mod PATH_SUFFIXES mct lib/mct build/lib/mct) +find_path(SCRIP_INCLUDE_DIR NAMES remap_conservative.mod PATH_SUFFIXES scrip lib/scrip build/lib/scrip) + +find_library(PSMILE_LIB NAMES psmile.MPI1 PATH_SUFFIXES lib) +find_library(MCT_LIB NAMES mct PATH_SUFFIXES lib) +find_library(MPEU_LIB NAMES mpeu PATH_SUFFIXES lib) +find_library(SCRIP_LIB NAMES scrip PATH_SUFFIXES lib) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(OASIS3MCT DEFAULT_MSG + PSMILE_INCLUDE_DIR + MCT_INCLUDE_DIR + SCRIP_INCLUDE_DIR + PSMILE_LIB + MCT_LIB + MPEU_LIB + SCRIP_LIB +) + +if(OASIS3MCT_FOUND) + set(OASIS3MCT_LIBRARIES ${PSMILE_LIB} ${MCT_LIB} ${MPEU_LIB} ${SCRIP_LIB}) + set(OASIS3MCT_INCLUDE_DIRS ${PSMILE_INCLUDE_DIR} ${MCT_INCLUDE_DIR} ${SCRIP_INCLUDE_DIR}) + if(NOT TARGET OASIS3MCT::OASIS3MCT) + add_library(OASIS3MCT::OASIS3MCT INTERFACE IMPORTED) + target_include_directories(OASIS3MCT::OASIS3MCT INTERFACE ${OASIS3MCT_INCLUDE_DIRS}) + target_link_libraries(OASIS3MCT::OASIS3MCT INTERFACE ${OASIS3MCT_LIBRARIES}) + target_link_options(OASIS3MCT::OASIS3MCT INTERFACE "-fopenmp") #TODO: see if this could be moved to SetBuildOptions.cmake + endif() +endif() + +unset(PSMILE_INCLUDE_DIR) +unset(MCT_INCLUDE_DIR) +unset(SCRIP_INCLUDE_DIR) +unset(PSMILE_LIB) +unset(MCT_LIB) +unset(MPEU_LIB) +unset(SCRIP_LIB) \ No newline at end of file diff --git a/src/clm5/CMakeLists.txt b/src/clm5/CMakeLists.txt index 719f97a131..0e32f30af5 100644 --- a/src/clm5/CMakeLists.txt +++ b/src/clm5/CMakeLists.txt @@ -30,6 +30,14 @@ foreach(f90in_file IN LISTS F90IN_FILES) endforeach() target_sources(${PROJECT_NAME} PRIVATE ${PREPROCESSED_F90_FILES}) +if (USE_OASIS) + target_sources(${PROJECT_NAME} + PRIVATE + oasis3/oas_defineMod.F90 + oasis3/oas_vardefMod.F90 + oasis3/oas_sendReceiveMod.F90 + ) +endif() target_sources(${PROJECT_NAME} PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/ncdio_pio.F90 diff --git a/src/clm5/biogeophys/HydrologyDrainageMod.F90 b/src/clm5/biogeophys/HydrologyDrainageMod.F90 index 3dd434661b..5a15cd59a2 100644 --- a/src/clm5/biogeophys/HydrologyDrainageMod.F90 +++ b/src/clm5/biogeophys/HydrologyDrainageMod.F90 @@ -50,7 +50,7 @@ subroutine HydrologyDrainage(bounds, & use column_varcon , only : icol_roof, icol_road_imperv, icol_road_perv, icol_sunwall, icol_shadewall use clm_varcon , only : denh2o, denice use clm_varctl , only : use_vichydro - use clm_varpar , only : nlevgrnd, nlevurb + use clm_varpar , only : nlevgrnd, nlevurb, nlevsoi use clm_time_manager , only : get_step_size, get_nstep use SoilHydrologyMod , only : CLMVICMap, Drainage, PerchedLateralFlow, LateralFlowPowerLaw use SoilWaterMovementMod , only : use_aquifer_layer @@ -106,7 +106,9 @@ subroutine HydrologyDrainage(bounds, & qflx_rsub_sat => waterflux_inst%qflx_rsub_sat_col , & ! soil saturation excess [mm h2o/s] qflx_drain => waterflux_inst%qflx_drain_col , & ! sub-surface runoff (mm H2O /s) qflx_surf => waterflux_inst%qflx_surf_col , & ! surface runoff (mm H2O /s) - qflx_infl => waterflux_inst%qflx_infl_col , & ! infiltration (mm H2O /s) + qflx_infl => waterflux_inst%qflx_infl_col , & ! infiltration (mm H2O /s) + qflx_rootsoi => waterflux_inst%qflx_rootsoi_col , & ! root and soil water exchange [mm H2O/s] [+ into root] + qflx_parflow => waterflux_inst%qflx_parflow_col , & ! source/sink flux passed to ParFlow for each soil layer qflx_qrgwl => waterflux_inst%qflx_qrgwl_col , & ! qflx_surf at glaciers, wetlands, lakes qflx_runoff => waterflux_inst%qflx_runoff_col , & ! total runoff ! (qflx_drain+qflx_surf+qflx_qrgwl) (mm H2O /s) @@ -114,7 +116,7 @@ subroutine HydrologyDrainage(bounds, & qflx_runoff_r => waterflux_inst%qflx_runoff_r_col , & ! Rural total runoff ! (qflx_drain+qflx_surf+qflx_qrgwl) (mm H2O /s) qflx_ice_runoff_snwcp => waterflux_inst%qflx_ice_runoff_snwcp_col, & ! solid runoff from snow capping (mm H2O /s) - qflx_irrig => irrigation_inst%qflx_irrig_col & ! irrigation flux (mm H2O /s) + qflx_irrig => irrigation_inst%qflx_irrig_col & ! irrigation flux (mm H2O /s) ) ! Determine time step and step size @@ -218,6 +220,17 @@ subroutine HydrologyDrainage(bounds, & end do + ! Calculate here the source/sink term for ParFlow + do j = 1, nlevsoi + do fc = 1, num_hydrologyc + c = filter_hydrologyc(fc) + if (j == 1) then + qflx_parflow(c,j) = qflx_infl(c) - qflx_rootsoi(c,j) + else + qflx_parflow(c,j) = -qflx_rootsoi(c,j) + end if + end do + end do end associate end subroutine HydrologyDrainage diff --git a/src/clm5/biogeophys/WaterfluxType.F90 b/src/clm5/biogeophys/WaterfluxType.F90 index 16b7217873..1f7806385a 100644 --- a/src/clm5/biogeophys/WaterfluxType.F90 +++ b/src/clm5/biogeophys/WaterfluxType.F90 @@ -72,6 +72,7 @@ module WaterfluxType real(r8), pointer :: qflx_adv_col (:,:) ! col advective flux across different soil layer interfaces [mm H2O/s] [+ downward] real(r8), pointer :: qflx_rootsoi_col (:,:) ! col root and soil water exchange [mm H2O/s] [+ into root] + real(r8), pointer :: qflx_parflow_col (:,:) ! col source/sink flux per soil layer sent to ParFlow [mm H2O/s] [- out from root] real(r8), pointer :: qflx_infl_col (:) ! col infiltration (mm H2O /s) real(r8), pointer :: qflx_surf_col (:) ! col surface runoff (mm H2O /s) real(r8), pointer :: qflx_drain_col (:) ! col sub-surface runoff (mm H2O /s) @@ -214,6 +215,7 @@ subroutine InitAllocate(this, bounds) allocate(this%qflx_drain_vr_col (begc:endc,1:nlevsoi)) ; this%qflx_drain_vr_col (:,:) = nan allocate(this%qflx_adv_col (begc:endc,0:nlevsoi)) ; this%qflx_adv_col (:,:) = nan allocate(this%qflx_rootsoi_col (begc:endc,1:nlevsoi)) ; this%qflx_rootsoi_col (:,:) = nan + allocate(this%qflx_parflow_col (begc:endc,1:nlevsoi)) ; this%qflx_parflow_col (:,:) = nan allocate(this%qflx_infl_col (begc:endc)) ; this%qflx_infl_col (:) = nan allocate(this%qflx_surf_col (begc:endc)) ; this%qflx_surf_col (:) = nan allocate(this%qflx_drain_col (begc:endc)) ; this%qflx_drain_col (:) = nan @@ -419,6 +421,11 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='water flux from soil to root in each soil-layer', & ptr_col=this%qflx_rootsoi_col, set_spec=spval, l2g_scale_type='veg', default='inactive') + this%qflx_parflow_col(begc:endc,:) = spval + call hist_addfld2d (fname='QPARFLOW', units='mm/s', type2d='levsoi', & + avgflag='A', long_name='Water fluxes per soil layer sent to Parflow', & + ptr_col=this%qflx_parflow_col, set_spec=spval, l2g_scale_type='veg', default='inactive') + this%qflx_evap_can_patch(begp:endp) = spval call hist_addfld1d (fname='QVEGE', units='mm/s', & avgflag='A', long_name='canopy evaporation', & diff --git a/src/clm5/cpl/lnd_comp_mct.F90 b/src/clm5/cpl/lnd_comp_mct.F90 index 923d6f6c34..37e979f144 100644 --- a/src/clm5/cpl/lnd_comp_mct.F90 +++ b/src/clm5/cpl/lnd_comp_mct.F90 @@ -64,6 +64,9 @@ subroutine lnd_init_mct( EClock, cdata_l, x2l_l, l2x_l, NLFilename ) use clm_varctl , only : nsrStartup, nsrContinue, nsrBranch use clm_cpl_indices , only : clm_cpl_indices_set use mct_mod , only : mct_aVect_init, mct_aVect_zero, mct_gsMap_lsize +#if defined(USE_OASIS) + use oas_defineMod , only : oas_definitions_init +#endif use ESMF ! ! !ARGUMENTS: @@ -228,6 +231,10 @@ subroutine lnd_init_mct( EClock, cdata_l, x2l_l, l2x_l, NLFilename ) call mct_aVect_init(l2x_l, rList=seq_flds_l2x_fields, lsize=lsize) call mct_aVect_zero(l2x_l) +#if defined(USE_OASIS) + ! Initialize OASIS3-MCT + call oas_definitions_init(bounds) +#endif ! Finish initializing clm call initialize2() diff --git a/src/clm5/cpl/lnd_import_export.F90 b/src/clm5/cpl/lnd_import_export.F90 index 4923b721af..b758484b76 100644 --- a/src/clm5/cpl/lnd_import_export.F90 +++ b/src/clm5/cpl/lnd_import_export.F90 @@ -8,6 +8,9 @@ module lnd_import_export use atm2lndType , only: atm2lnd_type use glc2lndMod , only: glc2lnd_type use clm_cpl_indices +#if defined(USE_OASIS) + use oas_sendReceiveMod +#endif ! implicit none !=============================================================================== @@ -276,7 +279,9 @@ subroutine lnd_import( bounds, x2l, glc_present, atm2lnd_inst, glc2lnd_inst) index_x2l_Flgg_hflx = index_x2l_Flgg_hflx, & index_x2l_Sg_icemask = index_x2l_Sg_icemask, & index_x2l_Sg_icemask_coupled_fluxes = index_x2l_Sg_icemask_coupled_fluxes) - +#if defined(USE_OASIS) + call oas_receive(bounds, atm2lnd_inst) +#endif end subroutine lnd_import !=============================================================================== @@ -423,7 +428,9 @@ subroutine lnd_export( bounds, lnd2atm_inst, lnd2glc_inst, l2x) end if end do - +#if defined(USE_OASIS) + call oas_send(bounds, lnd2atm_inst) +#endif end subroutine lnd_export end module lnd_import_export diff --git a/src/clm5/main/atm2lndType.F90 b/src/clm5/main/atm2lndType.F90 index 33452d5b2a..40548a4811 100644 --- a/src/clm5/main/atm2lndType.F90 +++ b/src/clm5/main/atm2lndType.F90 @@ -8,7 +8,7 @@ module atm2lndType use shr_kind_mod , only : r8 => shr_kind_r8 use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use shr_log_mod , only : errMsg => shr_log_errMsg - use clm_varpar , only : numrad, ndst, nlevgrnd !ndst = number of dust bins. + use clm_varpar , only : numrad, ndst, nlevgrnd, nlevsoi !ndst = number of dust bins. use clm_varcon , only : rair, grav, cpair, hfus, tfrz, spval use clm_varctl , only : iulog, use_c13, use_cn, use_lch4, use_cndv, use_fates, use_luna use decompMod , only : bounds_type @@ -116,6 +116,7 @@ module atm2lndType real(r8), pointer :: forc_flood_grc (:) => null() ! rof flood (mm/s) real(r8), pointer :: volr_grc (:) => null() ! rof volr total volume (m3) real(r8), pointer :: volrmch_grc (:) => null() ! rof volr main channel (m3) + real(r8), pointer :: parflow_psi_grc (:,:) => null() ! Parflow pressure head ! anomaly forcing real(r8), pointer :: af_precip_grc (:) => null() ! anomaly forcing @@ -540,9 +541,10 @@ subroutine InitAllocate(this, bounds) allocate(this%forc_snow_downscaled_col (begc:endc)) ; this%forc_snow_downscaled_col (:) = ival ! rof->lnd - allocate(this%forc_flood_grc (begg:endg)) ; this%forc_flood_grc (:) = ival - allocate(this%volr_grc (begg:endg)) ; this%volr_grc (:) = ival - allocate(this%volrmch_grc (begg:endg)) ; this%volrmch_grc (:) = ival + allocate(this%forc_flood_grc (begg:endg)) ; this%forc_flood_grc (:) = ival + allocate(this%volr_grc (begg:endg)) ; this%volr_grc (:) = ival + allocate(this%volrmch_grc (begg:endg)) ; this%volrmch_grc (:) = ival + allocate(this%parflow_psi_grc (begg:endg,1:nlevsoi)); this%parflow_psi_grc (:,:) = ival ! anomaly forcing allocate(this%bc_precip_grc (begg:endg)) ; this%bc_precip_grc (:) = ival @@ -577,7 +579,7 @@ end subroutine InitAllocate subroutine InitHistory(this, bounds) ! ! !USES: - use histFileMod, only : hist_addfld1d + use histFileMod, only : hist_addfld1d, hist_addfld2d ! ! !ARGUMENTS: class(atm2lnd_type) :: this @@ -608,6 +610,11 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='river channel main channel water storage', & ptr_lnd=this%volrmch_grc) + this%parflow_psi_grc(begg:endg, :) = 0._r8 + call hist_addfld2d (fname='PARFLOW_PSI', units='MPa', type2d='levsoi', & + avgflag='A', long_name='Parflow pressure head', & + ptr_lnd=this%parflow_psi_grc, default='inactive') + this%forc_wind_grc(begg:endg) = spval call hist_addfld1d (fname='WIND', units='m/s', & avgflag='A', long_name='atmospheric wind velocity magnitude', & diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index 5ef55973b1..7c975db627 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -11,7 +11,7 @@ module lnd2atmMod use shr_log_mod , only : errMsg => shr_log_errMsg use shr_megan_mod , only : shr_megan_mechcomps_n use shr_fire_emis_mod , only : shr_fire_emis_mechcomps_n - use clm_varpar , only : numrad, ndst, nlevgrnd !ndst = number of dust bins. + use clm_varpar , only : numrad, ndst, nlevgrnd, nlevsoi !ndst = number of dust bins. use clm_varcon , only : rair, grav, cpair, hfus, tfrz, spval use clm_varctl , only : iulog, use_lch4 use seq_drydep_mod , only : n_drydep, drydep_method, DD_XLND @@ -156,7 +156,7 @@ subroutine lnd2atm(bounds, & real(r8) , intent(in) :: net_carbon_exchange_grc( bounds%begg: ) ! net carbon exchange between land and atmosphere, positive for source (gC/m2/s) ! ! !LOCAL VARIABLES: - integer :: c, g ! indices + integer :: c, g, j ! indices real(r8) :: qflx_ice_runoff_col(bounds%begc:bounds%endc) ! total column-level ice runoff real(r8) :: eflx_sh_ice_to_liq_grc(bounds%begg:bounds%endg) ! sensible heat flux generated from the ice to liquid conversion, averaged to gridcell real(r8), parameter :: amC = 12.0_r8 ! Atomic mass number for Carbon @@ -406,6 +406,18 @@ subroutine lnd2atm(bounds, & waterstate_inst%tws_grc(g) = waterstate_inst%tws_grc(g) + atm2lnd_inst%volr_grc(g) / grc%area(g) * 1.e-3_r8 enddo + ! Calculate Parflow water fluxes + call c2g( bounds, nlevsoi, & + waterflux_inst%qflx_parflow_col (bounds%begc:bounds%endc, :), & + lnd2atm_inst%qflx_parflow_grc (bounds%begg:bounds%endg, :), & + c2l_scale_type= 'unity', l2g_scale_type='unity' ) +! TODO: Verify questionable factor 3.6/dz(g,j) ! +! do j = 1, nlevsoi +! do g = bounds%begg, bounds%endg +! lnd2atm_inst%qflx_parflow_grc(g) = lnd2atm_inst%qflx_parflow_grc(g)*3.6_r8/col%dz(g,j) +! enddo +! enddo + end subroutine lnd2atm !----------------------------------------------------------------------- diff --git a/src/clm5/main/lnd2atmType.F90 b/src/clm5/main/lnd2atmType.F90 index 023608a4a7..6159e24c73 100644 --- a/src/clm5/main/lnd2atmType.F90 +++ b/src/clm5/main/lnd2atmType.F90 @@ -10,7 +10,7 @@ module lnd2atmType use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun use decompMod , only : bounds_type - use clm_varpar , only : numrad, ndst, nlevgrnd !ndst = number of dust bins. + use clm_varpar , only : numrad, ndst, nlevgrnd, nlevsoi !ndst = number of dust bins. use clm_varcon , only : spval use clm_varctl , only : iulog, use_lch4 use shr_megan_mod , only : shr_megan_mechcomps_n @@ -68,9 +68,10 @@ module lnd2atmType real(r8), pointer :: qflx_rofliq_qgwl_grc (:) => null() ! rof liq -- glacier, wetland and lakes water balance residual component real(r8), pointer :: qflx_rofliq_h2osfc_grc (:) => null() ! rof liq -- surface water runoff component real(r8), pointer :: qflx_rofliq_drain_perched_grc (:) => null() ! rof liq -- perched water table runoff component - real(r8), pointer :: qflx_rofice_grc (:) => null() ! rof ice forcing - real(r8), pointer :: qflx_liq_from_ice_col(:) => null() ! liquid runoff from converted ice runoff - real(r8), pointer :: qirrig_grc (:) => null() ! irrigation flux + real(r8), pointer :: qflx_rofice_grc (:) => null() ! rof ice forcing + real(r8), pointer :: qflx_liq_from_ice_col (:) => null() ! liquid runoff from converted ice runoff + real(r8), pointer :: qflx_parflow_grc (:,:) => null() ! source/sink flux per soil layer sent to ParFlow [mm H2O/m^2/s] [- out from root] + real(r8), pointer :: qirrig_grc (:) => null() ! irrigation flux contains @@ -174,9 +175,10 @@ subroutine InitAllocate(this, bounds) allocate(this%qflx_rofliq_qgwl_grc (begg:endg)) ; this%qflx_rofliq_qgwl_grc (:) =ival allocate(this%qflx_rofliq_h2osfc_grc (begg:endg)) ; this%qflx_rofliq_h2osfc_grc (:) =ival allocate(this%qflx_rofliq_drain_perched_grc (begg:endg)) ; this%qflx_rofliq_drain_perched_grc (:) =ival - allocate(this%qflx_rofice_grc (begg:endg)) ; this%qflx_rofice_grc (:) =ival - allocate(this%qflx_liq_from_ice_col(begc:endc)) ; this%qflx_liq_from_ice_col(:) = ival - allocate(this%qirrig_grc (begg:endg)) ; this%qirrig_grc (:) =ival + allocate(this%qflx_rofice_grc (begg:endg)) ; this%qflx_rofice_grc (:) =ival + allocate(this%qflx_liq_from_ice_col(begc:endc)) ; this%qflx_liq_from_ice_col(:) =ival + allocate(this%qflx_parflow_grc (begg:endg,1:nlevsoi)); this%qflx_parflow_grc (:,:) =ival + allocate(this%qirrig_grc (begg:endg)) ; this%qirrig_grc (:) =ival if (shr_megan_mechcomps_n>0) then allocate(this%flxvoc_grc(begg:endg,1:shr_megan_mechcomps_n)); this%flxvoc_grc(:,:)=ival @@ -260,7 +262,7 @@ end subroutine ReadNamelist subroutine InitHistory(this, bounds) ! ! !USES: - use histFileMod, only : hist_addfld1d + use histFileMod, only : hist_addfld1d, hist_addfld2d ! ! !ARGUMENTS: class(lnd2atm_type) :: this @@ -324,6 +326,11 @@ subroutine InitHistory(this, bounds) ptr_lnd=this%nem_grc) end if + this%qflx_parflow_grc(begg:endg, :) = 0._r8 + call hist_addfld2d (fname='QPARFLOW_TO_OASIS', units='mm H2O/m2/s', type2d='levsoi', & + avgflag='A', long_name='source/sink flux per soil layer sent to ParFlow', & + ptr_lnd=this%qflx_parflow_grc, default='inactive') + end subroutine InitHistory end module lnd2atmType diff --git a/src/clm5/oasis3/oas_defineMod.F90 b/src/clm5/oasis3/oas_defineMod.F90 new file mode 100644 index 0000000000..a861fd9feb --- /dev/null +++ b/src/clm5/oasis3/oas_defineMod.F90 @@ -0,0 +1,169 @@ +module oas_defineMod + use mod_oasis + implicit none + save + private + + public :: oas_definitions_init + integer :: ierror +contains + + subroutine oas_definitions_init(bounds) + use spmdMod , only : masterproc + use clm_varpar , only : nlevsoi + use decompMod , only : bounds_type + use oas_vardefMod + + type(bounds_type) , intent(in) :: bounds + integer :: partition(3) + integer :: grid_id ! id returned by oasis_def_partition + integer :: var_nodims(2) + integer :: var_shape(1) ! not used by oasis_def_var + character(len=2) :: soil_layer + integer :: i, n_grid_points + + if (masterproc) then + call define_grid() + end if + + ! ----------------------------------------------------------------- + ! ... Define partition + ! ----------------------------------------------------------------- + n_grid_points = (bounds%endg - bounds%begg) + 1 + partition(1) = 1 ! Apple style partition + partition(2) = bounds%begg - 1 ! Global offset + partition(3) = n_grid_points ! # of grid cells allocated to this MPI task + call oasis_def_partition(grid_id, partition, ierror) + + ! ----------------------------------------------------------------- + ! ... Define coupling fields + ! ----------------------------------------------------------------- + var_nodims(1) = 1 ! var_nodims(1) is not used anymore in OASIS + var_nodims(2) = 1 ! number of fields in a bundle + + allocate(et_loss(nlevsoi)) + allocate(watsat(nlevsoi)) + allocate(psi(nlevsoi)) + + do i = 1, nlevsoi + write (soil_layer, '(I2.2)') i ! soil layer index (01-10) + + et_loss(i)%name = 'CLMFLX'//soil_layer ! Evapotranspiration fluxes sent to Parflow + watsat(i)%name = 'CLMSAT'//soil_layer ! Water saturation received from Parflow + psi(i)%name = 'CLMPSI'//soil_layer ! Pressure head received from Parflow + + call oasis_def_var(et_loss(i)%id, et_loss(i)%name, grid_id, var_nodims, OASIS_Out, var_shape, OASIS_Real, ierror) + call oasis_def_var(watsat(i)%id, watsat(i)%name, grid_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror) + call oasis_def_var(psi(i)%id, psi(i)%name, grid_id, var_nodims, OASIS_In, var_shape, OASIS_Real, ierror) + end do + + ! End definition phase + call oasis_enddef(ierror) + end subroutine oas_definitions_init + + subroutine define_grid() + use shr_kind_mod , only : r8 => shr_kind_r8 + use domainMod , only : ldomain + character(len=4), parameter :: grid_name='gclm' + integer, parameter :: SOUTH = 1 + integer, parameter :: NORTH = 2 + integer, parameter :: WEST = 1 + integer, parameter :: EAST = 2 + real(kind=r8), allocatable :: corner_lon(:,:), corner_lat(:,:) + real(kind=r8), allocatable :: oas_lon(:,:), oas_lat(:,:) + real(kind=r8), allocatable :: oas_corner_lon(:,:,:), oas_corner_lat(:,:,:) + integer, allocatable :: oas_mask(:,:) + real(kind=r8) :: center, offset, neighbor + integer :: write_grid_files + integer :: i, j + integer :: n_lons, n_lats, n_gridcells + + call oasis_start_grids_writing(write_grid_files) + + if (write_grid_files == 1) then + n_lons = ldomain%ni + n_lats = ldomain%nj + n_gridcells = ldomain%ns + + ! ----------------------------------------------------------------- + ! ... Define centers + ! ----------------------------------------------------------------- + allocate(oas_lon(n_gridcells, 1)) + allocate(oas_lat(n_gridcells, 1)) + do j = 1, n_lats + i = (j-1)*n_lons + 1 + oas_lon(i:n_lons, 1) = ldomain%lonc(:) + oas_lat(i:n_lons, 1) = ldomain%latc(j) + enddo + call oasis_write_grid(grid_name, n_gridcells, 1, oas_lon, oas_lat) + + ! ----------------------------------------------------------------- + ! ... Define corners + ! ----------------------------------------------------------------- + allocate(corner_lon(n_lons, WEST:EAST)) + allocate(corner_lat(n_lats, SOUTH:NORTH)) + allocate(oas_corner_lon(n_gridcells, 1, 4)) + allocate(oas_corner_lat(n_gridcells, 1, 4)) + + do j = 1, n_lats + if (j == 1) then + neighbor = ldomain%latc(j+1) + else + neighbor = ldomain%latc(j-1) + end if + center = ldomain%latc(j) + offset = abs(center - neighbor) / 2.0_r8 + corner_lat(j, SOUTH:NORTH) = [center-offset, center+offset] + enddo + + do i = 1, n_lons + if (i == 1) then + neighbor = ldomain%lonc(i+1) + else + neighbor = ldomain%lonc(i-1) + end if + center = ldomain%lonc(i) + offset = abs(center - neighbor) / 2.0_r8 + corner_lon(i, WEST:EAST) = [center-offset, center+offset] + enddo + + ! oas_corner indexing scheme + ! 4 +-----+ 3 NORTH + ! | | + ! | | + ! 1 +-----+ 2 SOUTH + ! WEST EAST + do j = 1, n_lats + i = (j-1)*n_lons + 1 + oas_corner_lon(i:n_lons,1,1:2) = corner_lon(:, WEST:EAST) ! bottom side + oas_corner_lat(i:n_lons,1,2) = corner_lat(j, SOUTH) ! right side + oas_corner_lat(i:n_lons,1,3) = corner_lat(j, NORTH) ! right side + enddo + + ! Fill missing longitudes and latitudes + oas_corner_lon(:,1,[4,3]) = oas_corner_lon(:,1,[1,2]) ! West-East for top side + oas_corner_lat(:,1,[1,4]) = oas_corner_lat(:,1,[2,3]) ! South-North for left side + call oasis_write_corner(grid_name, n_gridcells, 1, 4, oas_corner_lon, oas_corner_lat) + + ! ----------------------------------------------------------------- + ! ... Define mask + ! ----------------------------------------------------------------- + ! CLM5 landmask convention: 0 = ocean, 1 = land + allocate(oas_mask(n_gridcells, 1)) + do j = 1, n_lats + i = (j-1)*n_lons + 1 + oas_mask(i:n_lons,1) = ldomain%mask(:) + enddo + + ! Invert mask to conform to OASIS convention: 0 = not masked, 1 = masked + where (oas_mask(:,1) == 0) + oas_mask(:,1) = 1 + else where + oas_mask(:,1) = 0 + end where + call oasis_write_mask(grid_name, n_gridcells, 1, oas_mask) + call oasis_terminate_grids_writing() + deallocate(oas_lon, oas_lat, oas_corner_lon, oas_corner_lat, oas_mask, corner_lon, corner_lat) + end if + end subroutine define_grid +end module oas_defineMod \ No newline at end of file diff --git a/src/clm5/oasis3/oas_sendReceiveMod.F90 b/src/clm5/oasis3/oas_sendReceiveMod.F90 new file mode 100644 index 0000000000..2ccf57e1ba --- /dev/null +++ b/src/clm5/oasis3/oas_sendReceiveMod.F90 @@ -0,0 +1,50 @@ +module oas_sendReceiveMod + use shr_kind_mod , only: r8 => shr_kind_r8 + use clm_time_manager , only: get_curr_time, get_prev_time + use decompMod , only: bounds_type + use clm_varpar , only: nlevsoi + use oas_vardefMod + use mod_oasis + implicit none + save + private + + public :: oas_send + public :: oas_receive + integer :: days_elapsed, seconds_elapsed + integer :: n_grid_points + integer :: i, ierror + +contains + + subroutine oas_receive(bounds, atm2lnd_inst) + use atm2lndType, only: atm2lnd_type + type(bounds_type), intent(in) :: bounds + type(atm2lnd_type), intent(inout) :: atm2lnd_inst + real(kind=r8), allocatable :: buffer(:) + + n_grid_points = (bounds%endg - bounds%begg) + 1 + allocate(buffer(n_grid_points)) + call get_curr_time(days_elapsed, seconds_elapsed) + do i = 1, nlevsoi + call oasis_get(psi(i)%id, seconds_elapsed, buffer, ierror) + atm2lnd_inst%parflow_psi_grc(:,i) = buffer + end do + end subroutine oas_receive + + subroutine oas_send(bounds, lnd2atm_inst) + use lnd2atmType, only : lnd2atm_type + type(bounds_type), intent(in) :: bounds + type(lnd2atm_type), intent(inout) :: lnd2atm_inst + real(kind=r8), allocatable :: buffer(:) + + n_grid_points = (bounds%endg - bounds%begg) + 1 + allocate(buffer(n_grid_points)) + call get_prev_time(days_elapsed, seconds_elapsed) + do i = 1, nlevsoi + buffer = lnd2atm_inst%qflx_parflow_grc(:,i) + call oasis_put(et_loss(i)%id, seconds_elapsed, buffer, ierror) + end do + end subroutine oas_send + +end module oas_sendReceiveMod diff --git a/src/clm5/oasis3/oas_vardefMod.F90 b/src/clm5/oasis3/oas_vardefMod.F90 new file mode 100644 index 0000000000..12d5d9b87b --- /dev/null +++ b/src/clm5/oasis3/oas_vardefMod.F90 @@ -0,0 +1,18 @@ +module oas_vardefMod + implicit none + save + + ! Type for coupling field information + type :: oas_var + character(len=12) :: name ! Name of the coupling field + integer :: id ! Id of the field + end type oas_var + + ! Sent fields + type(oas_var), allocatable :: et_loss(:) ! soil ET loss + + ! Received fields from Parflow + type(oas_var), allocatable :: watsat(:) ! water saturation + type(oas_var), allocatable :: psi(:) ! pressure head + +end module oas_vardefMod \ No newline at end of file diff --git a/src/csm_share/CMakeLists.txt b/src/csm_share/CMakeLists.txt index b9f4d7bcec..6309ed82e4 100644 --- a/src/csm_share/CMakeLists.txt +++ b/src/csm_share/CMakeLists.txt @@ -112,4 +112,11 @@ target_sources(${PROJECT_NAME} PRIVATE ) target_link_libraries(${PROJECT_NAME} PUBLIC pio mct gptl) +if (USE_OASIS) + find_package(OASIS3MCT REQUIRED) + if (OASIS3MCT_FOUND) + target_link_libraries(${PROJECT_NAME} PUBLIC OASIS3MCT::OASIS3MCT) + target_compile_definitions(${PROJECT_NAME} PUBLIC USE_OASIS) + endif() +endif() install (TARGETS ${PROJECT_NAME} ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}) diff --git a/src/eclm/cime_comp_mod.F90 b/src/eclm/cime_comp_mod.F90 index 41da39c066..77edd82765 100644 --- a/src/eclm/cime_comp_mod.F90 +++ b/src/eclm/cime_comp_mod.F90 @@ -40,7 +40,9 @@ module cime_comp_mod use mct_mod ! mct_ wrappers for mct lib use perf_mod use ESMF - +#if defined(USE_OASIS) + use mod_oasis +#endif !---------------------------------------------------------------------------- ! component model interfaces (init, run, final methods) !---------------------------------------------------------------------------- @@ -583,6 +585,7 @@ module cime_comp_mod subroutine cime_pre_init1(esmf_log_option) use shr_pio_mod, only : shr_pio_init1, shr_pio_init2 use seq_comm_mct, only: num_inst_driver + !---------------------------------------------------------- !| Initialize MCT and MPI communicators and IO !---------------------------------------------------------- @@ -593,13 +596,25 @@ subroutine cime_pre_init1(esmf_log_option) logical :: comp_iamin(num_inst_total) character(len=seq_comm_namelen) :: comp_name(num_inst_total) integer :: it - integer :: driver_id + integer :: driver_id, oas_comp_id integer :: driver_comm call mpi_init(ierr) call shr_mpi_chkerr(ierr,subname//' mpi_init') + +#if defined(USE_OASIS) + call oasis_init_comp (oas_comp_id, 'eCLM', ierr) + if (ierr /= 0) then + call oasis_abort(oas_comp_id, 'cime_pre_init1', 'oasis_init_comp error') + end if + call oasis_get_localcomm(global_comm, ierr) + if (ierr /= 0) then + call oasis_abort(oas_comp_id, 'cime_pre_init1', 'oasis_get_localcomm error') + end if +#else call mpi_comm_dup(MPI_COMM_WORLD, global_comm, ierr) call shr_mpi_chkerr(ierr,subname//' mpi_comm_dup') +#endif comp_comm = MPI_COMM_NULL time_brun = mpi_wtime() @@ -4135,6 +4150,10 @@ subroutine cime_final() mpicom=mpicom_GLOID) endif +#if defined(USE_OASIS) + call oasis_terminate (ierr) + call shr_mpi_chkerr(ierr,subname//' oasis_terminate') +#endif call t_finalizef() end subroutine cime_final