OASIS3-MCT integration initial developments (ESCOMP#9)

Commit summary: * 6ede8a4 Avoid eCLM-Parflow deadlock by building internal MCT and depending on a patched version of OASIS3-MCT * fb58947 Fixed wrong loop index * 4160724 Added oas_sendReceiveMod and other changes * 3040a97 Added water flux QPARFLOW. This is equal to ET loss per soil layer * 3b18c34 Defined OASIS partition type and Parflow coupling fields * eafb1e1 Defined OASIS lnd grid at lnd_init_mct * b8f4a08 Prelim codes for oas_define_partition and oas_define_coupling_fields * dbcdf4e Prelim codes for oas_define_grid() * 3763f08 Scaffold for OASIS define routines * b7d87d0 Added OASIS init and terminate functions
HPSCTerrSys · Oct 13, 2021 · 2567de7 · 2567de7
1 parent 7251076
commit 2567de7
Show file tree

Hide file tree

Showing 14 changed files with 400 additions and 22 deletions.
diff --git a/cmake/FindOASIS3MCT.cmake 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)
diff --git 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
@@ -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,15 +106,17 @@ 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)
          qflx_runoff_u      => waterflux_inst%qflx_runoff_u_col      , & ! Urban total runoff (qflx_drain+qflx_surf) (mm H2O /s)
          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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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