Merge branch 'bishtgautam/lnd/mpp-v2-via-emi' into next (PR #1460)

Updates VSFM and adds EMI

This PR includes following changes:

1)Updates Variably Saturated Flow Model (VSFM) and converts to submodule.
a)Introduces Multi-Physics Problem (MPP) formulation.
VSFM is an example of an MPP for subsurface hydrology.
Another example of MPP is the PETSc-based Thermal Model (PTM).
for snow-soil-standing water thermal hydrology.
b)Upgrades compatibility with PETSc > 3.7.0.
c)Adds initial support for including lateral subsurface flow.

2)Adds External Model Interface (EMI)
EMI is a generic interface to couple ALM with multiple external models (EMs)
such as SBeTR, VSFM, FATES, PFLOTRAN, etc.

[BFB]

Conflicts:
	.gitmodules
	components/clm/bld/configure
	components/clm/src/biogeophys/BalanceCheckMod.F90
	components/clm/src/biogeophys/SoilWaterMovementMod.F90
	components/clm/src/main/clm_initializeMod.F90

.gitmodules

@@ -12,3 +12,7 @@
 [submodule "components/clm/src/external_models/fates"]
 	path = components/clm/src/external_models/fates
 	url = git@github.com:ACME-Climate/fates.git
+[submodule "alm-mpp"]
+	path = components/clm/src/external_models/mpp
+	url = git@github.com:ACME-Climate/mpp.git
+	branch = alm/develop

components/clm/bld/configure

@@ -651,13 +651,17 @@ sub write_filepath_cesmbld
 		     "betr/bgc_century",
 		     "betr/bgc_sminn",
 		     "biogeophys",
-                     "biogeophys/vsfm",
 		     "biogeochem", 
 		     "dyn_subgrid",
 		     "external_models/fates/main",
 		     "external_models/fates/biogeophys",
 		     "external_models/fates/biogeochem",
 		     "external_models/fates/fire",
+		     "external_models/mpp/src/mpp/dtypes",
+		     "external_models/mpp/src/mpp/thermal",
+		     "external_models/mpp/src/mpp/util",
+		     "external_models/mpp/src/mpp/vsfm",
+		     "external_models/emi/src/",
 		     "utils", 
 		     "cpl" );
 

components/clm/bld/namelist_files/namelist_defaults_clm4_5.xml

@@ -70,10 +70,19 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 
 <glc_snow_persistence_max_days>7300</glc_snow_persistence_max_days>
 
-<!-- VSFM default -->
+<!-- ========================================================================================  -->
+<!-- VSFM default                                                                              -->
+<!-- ========================================================================================  -->
 <use_vsfm>.false.</use_vsfm>
+<vsfm_include_seepage_bc>.false.</vsfm_include_seepage_bc>
 <vsfm_satfunc_type>smooth_brooks_corey_bz3</vsfm_satfunc_type>
 <vsfm_use_dynamic_linesearch>.false.</vsfm_use_dynamic_linesearch>
+<vsfm_lateral_model_type>none</vsfm_lateral_model_type>
+
+<!-- ========================================================================================  -->
+<!-- PETSc-based thermal model default                                                         -->
+<!-- ========================================================================================  -->
+<use_petsc_thermal_model>.false.</use_petsc_thermal_model>
 
 <!-- ================================================================== -->
 <!-- The default filenames are given relative to the root directory
@@ -1768,6 +1777,12 @@ this mask will have smb calculated over the entire global land surface
 <use_aereoxid_prog use_ed=".true.">.true.</use_aereoxid_prog>
 <use_aereoxid_prog use_cn=".true.">.true.</use_aereoxid_prog>
 
+<!-- =========================================  -->
+<!-- Defaults for lateral grid connectivity     -->
+<!-- =========================================  -->
+<lateral_connectivity>.false.</lateral_connectivity>
+<domain_decomposition_method>round_robin</domain_decomposition_method>
+
 <!-- =========================================  -->
 <!-- Defaults for clm_pflotran_inparm interface -->
 <!-- =========================================  -->

components/clm/bld/namelist_files/namelist_definition_clm4_5.xml

@@ -1553,6 +1553,19 @@ If TRUE (which is the default), check consistency between pct_nat_pft on the fla
 and pct_nat_pft read from the surface dataset.
 </entry>
 
+<!-- =========================================================================================  -->
+<!-- Namelist options related to lateral grid connectivity
+<!-- ========================================================================================  -->
+<entry id="lateral_connectivity" type="logical" category="clm_initialization"
+       group="clm_inparm" valid_values="" >
+If TRUE setup up lateral grid connectivity in CLM.
+</entry>
+
+<entry id="domain_decomp_type" type="char*256" category="clm_initialization"
+       group="clm_inparm" valid_values="round_robin,graph_partitioning" >
+Specifies the method for decomposing CLM grids across processors.
+</entry>
+
 <!-- =========================================================================================  -->
 <!-- Namelist options related to the bgc & pflotran interface
 <!-- ========================================================================================  -->
@@ -1604,6 +1617,29 @@ Type of saturation function used in VSFM.
 <entry id="vsfm_use_dynamic_linesearch"  type="logical" category="default_settings"
        group="clm_inparm" valid_values="">
 Runtime flag to dynamically modify PETSc SNES linesearch option when VSFM fails to converge before cutting the timestep.
+
+<entry id="vsfm_lateral_model_type" type="char*32" category="clm_physics"
+       group="clm_inparm"
+       valid_values="none,source_sink,three_dimensional" >
+Type of lateral flow formulation supported in VSFM.
+    none               = No lateral flow
+    source_sink        = 1D VSFM solution + Lateral flow modeled as source/sink term
+    three_dimensional  = 3D VSFM solution (only supported for serial run)
 </entry>
 
+<entry id=vsfm_include_seepage_bc" type="logical" category="default_settings"
+        group="clm_inparm" valid_values="">
+Runtime flag to add seepage boundary condition on the top soil layer to allow for exfiltration.
+</entry>
+
+<!-- ========================================================================================  -->
+<!-- Namelist options related to PETSc-based thermal model                                     -->
+<!-- ========================================================================================  -->
+
+<entry id="use_petsc_thermal_model" type="logical" category="default_settings"
+        group="clm_inparm" valid_values="">
+Runtime flag to turn on/off PETSc based thermal model.
+</entry>
+
+
 </namelist_definition>

components/clm/src/biogeophys/BalanceCheckMod.F90

@@ -239,6 +239,7 @@ subroutine BalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
           qflx_ice_dynbal            =>    waterflux_vars%qflx_ice_dynbal_grc         , & ! Input:  [real(r8) (:)   ]  ice runoff due to dynamic land cover change (mm H2O /s)
           snow_sources               =>    waterflux_vars%snow_sources_col            , & ! Output: [real(r8) (:)   ]  snow sources (mm H2O /s)  
           snow_sinks                 =>    waterflux_vars%snow_sinks_col              , & ! Output: [real(r8) (:)   ]  snow sinks (mm H2O /s)    
+          qflx_lateral               =>    waterflux_vars%qflx_lateral_col            , & ! Input:  [real(r8) (:)   ]  lateral flux of water to neighboring column (mm H2O /s)
 
           eflx_lwrad_out             =>    energyflux_vars%eflx_lwrad_out_patch       , & ! Input:  [real(r8) (:)   ]  emitted infrared (longwave) radiation (W/m**2)
           eflx_lwrad_net             =>    energyflux_vars%eflx_lwrad_net_patch       , & ! Input:  [real(r8) (:)   ]  net infrared (longwave) rad (W/m**2) [+ = to atm]
@@ -308,8 +309,10 @@ subroutine BalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
              errh2o(c) = endwb(c) - begwb(c) &
                   - (forc_rain_col(c) + forc_snow_col(c)  + qflx_floodc(c) + qflx_irrig(c) &
                   - qflx_evap_tot(c) - qflx_surf(c)  - qflx_h2osfc_surf(c) &
-                  - qflx_qrgwl(c) - qflx_drain(c) - qflx_drain_perched(c) - qflx_snwcp_ice(c)) * dtime
+                  - qflx_qrgwl(c) - qflx_drain(c) - qflx_drain_perched(c) - qflx_snwcp_ice(c) &
+                  - qflx_lateral(c) ) * dtime
              dwb(c) = (endwb(c)-begwb(c))/dtime
+
           else
 
              errh2o(c) = 0.0_r8
@@ -371,6 +374,7 @@ subroutine BalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
              write(iulog,*)'qflx_qrgwl                 = ',qflx_qrgwl(indexc)
              write(iulog,*)'qflx_drain                 = ',qflx_drain(indexc)
              write(iulog,*)'qflx_snwcp_ice             = ',qflx_snwcp_ice(indexc)
+             write(iulog,*)'qflx_lateral               = ',qflx_lateral(indexc)
              write(iulog,*)'total_plant_stored_h2o_col = ',total_plant_stored_h2o_col(indexc)
              write(iulog,*)'clm model is stopping'
              call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(__FILE__, __LINE__))
@@ -395,6 +399,7 @@ subroutine BalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
              write(iulog,*)'qflx_snwcp_ice             = ',qflx_snwcp_ice(indexc)
              write(iulog,*)'qflx_glcice_melt           = ',qflx_glcice_melt(indexc)
              write(iulog,*)'qflx_glcice_frz            = ',qflx_glcice_frz(indexc) 
+             write(iulog,*)'qflx_lateral               = ',qflx_lateral(indexc)
              write(iulog,*)'total_plant_stored_h2o_col = ',total_plant_stored_h2o_col(indexc)
              write(iulog,*)'clm model is stopping'
              call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(__FILE__, __LINE__))

components/clm/src/biogeophys/EnergyFluxType.F90

@@ -60,13 +60,18 @@ module EnergyFluxType
      real(r8), pointer :: eflx_traffic_lun        (:)   ! lun traffic sensible heat flux (W/m**2)
      real(r8), pointer :: eflx_wasteheat_lun      (:)   ! lun sensible heat flux from domestic heating/cooling sources of waste heat (W/m**2)
      real(r8), pointer :: eflx_heat_from_ac_lun   (:)   ! lun sensible heat flux to be put back into canyon due to removal by AC (W/m**2)
+     real(r8), pointer :: eflx_hs_h2osfc_col      (:)   ! heat flux on standing water [W/m2]
+     real(r8), pointer :: eflx_hs_top_snow_col    (:)   ! heat flux on top snow layer [W/m2]
+     real(r8), pointer :: eflx_hs_soil_col        (:)   ! heat flux on soil [W/m2
+     real(r8), pointer :: eflx_sabg_lyr_col       (:,:) ! absorbed solar radiation (col,lyr) [W/m2]
 
      ! Derivatives of energy fluxes
      real(r8), pointer :: dgnetdT_patch           (:)   ! patch derivative of net ground heat flux wrt soil temp  (W/m**2 K)
      real(r8), pointer :: netrad_patch            (:)   ! col net radiation (W/m**2) [+ = to sfc]
      real(r8), pointer :: cgrnd_patch             (:)   ! col deriv. of soil energy flux wrt to soil temp [W/m2/k]
      real(r8), pointer :: cgrndl_patch            (:)   ! col deriv. of soil latent heat flux wrt soil temp  [W/m**2/k]
      real(r8), pointer :: cgrnds_patch            (:)   ! col deriv. of soil sensible heat flux wrt soil temp [W/m2/k]
+     real(r8), pointer :: eflx_dhsdT_col          (:)   ! col deriv. of energy flux into surface layer wrt temp [W/m2/K]
 
      ! Canopy radiation
      real(r8), pointer :: dlrad_patch             (:)   ! col downward longwave radiation below the canopy [W/m2]
@@ -196,6 +201,12 @@ subroutine InitAllocate(this, bounds)
     allocate( this%eflx_wasteheat_lun      (begl:endl))             ; this%eflx_wasteheat_lun      (:)   = nan
     allocate( this%eflx_anthro_patch       (begp:endp))             ; this%eflx_anthro_patch       (:)   = nan
 
+    allocate( this%eflx_hs_top_snow_col    (begc:endc))             ; this%eflx_hs_top_snow_col    (:)   = nan
+    allocate( this%eflx_hs_h2osfc_col      (begc:endc))             ; this%eflx_hs_h2osfc_col      (:)   = nan
+    allocate( this%eflx_hs_soil_col        (begc:endc))             ; this%eflx_hs_soil_col        (:)   = nan
+    allocate( this%eflx_sabg_lyr_col       (begc:endc,-nlevsno+1:1)); this%eflx_sabg_lyr_col       (:,:) = nan
+    allocate( this%eflx_dhsdT_col          (begc:endc))             ; this%eflx_dhsdT_col          (:)   = nan
+
     allocate( this%dgnetdT_patch           (begp:endp))             ; this%dgnetdT_patch           (:)   = nan
     allocate( this%cgrnd_patch             (begp:endp))             ; this%cgrnd_patch             (:)   = nan
     allocate( this%cgrndl_patch            (begp:endp))             ; this%cgrndl_patch            (:)   = nan

components/clm/src/biogeophys/SoilStateType.F90

@@ -85,6 +85,7 @@ module SoilStateType
      procedure, private :: InitHistory  
      procedure, private :: InitCold    
      procedure, public  :: Restart
+     procedure, public  :: InitColdGhost
 
   end type soilstate_type
   !------------------------------------------------------------------------
@@ -117,11 +118,13 @@ subroutine InitAllocate(this, bounds)
     integer :: begp, endp
     integer :: begc, endc
     integer :: begg, endg
+    integer :: begc_all, endc_all
     !------------------------------------------------------------------------
 
-    begp = bounds%begp; endp= bounds%endp
-    begc = bounds%begc; endc= bounds%endc
-    begg = bounds%begg; endg= bounds%endg
+    begp     = bounds%begp    ; endp     = bounds%endp
+    begc     = bounds%begc    ; endc     = bounds%endc
+    begg     = bounds%begg    ; endg     = bounds%endg
+    begc_all = bounds%begc_all; endc_all = bounds%endc_all
 
     allocate(this%mss_frc_cly_vld_col  (begc:endc))                     ; this%mss_frc_cly_vld_col  (:)   = nan
     allocate(this%sandfrac_patch       (begp:endp))                     ; this%sandfrac_patch       (:)   = nan
@@ -131,14 +134,14 @@ subroutine InitAllocate(this, bounds)
     allocate(this%cellclay_col         (begc:endc,nlevsoi))             ; this%cellclay_col         (:,:) = nan 
     allocate(this%bd_col               (begc:endc,nlevgrnd))            ; this%bd_col               (:,:) = nan
 
-    allocate(this%hksat_col            (begc:endc,nlevgrnd))            ; this%hksat_col            (:,:) = spval
+    allocate(this%hksat_col            (begc_all:endc_all,nlevgrnd))    ; this%hksat_col            (:,:) = spval
     allocate(this%hksat_min_col        (begc:endc,nlevgrnd))            ; this%hksat_min_col        (:,:) = spval
     allocate(this%hk_l_col             (begc:endc,nlevgrnd))            ; this%hk_l_col             (:,:) = nan   
     allocate(this%smp_l_col            (begc:endc,nlevgrnd))            ; this%smp_l_col            (:,:) = nan   
     allocate(this%smpmin_col           (begc:endc))                     ; this%smpmin_col           (:)   = nan
 
-    allocate(this%bsw_col              (begc:endc,nlevgrnd))            ; this%bsw_col              (:,:) = nan
-    allocate(this%watsat_col           (begc:endc,nlevgrnd))            ; this%watsat_col           (:,:) = nan
+    allocate(this%bsw_col              (begc_all:endc_all,nlevgrnd))    ; this%bsw_col              (:,:) = nan
+    allocate(this%watsat_col           (begc_all:endc_all,nlevgrnd))    ; this%watsat_col           (:,:) = nan
     allocate(this%watdry_col           (begc:endc,nlevgrnd))            ; this%watdry_col           (:,:) = spval
     allocate(this%watopt_col           (begc:endc,nlevgrnd))            ; this%watopt_col           (:,:) = spval
     allocate(this%watfc_col            (begc:endc,nlevgrnd))            ; this%watfc_col            (:,:) = nan
@@ -855,4 +858,127 @@ subroutine Restart(this, bounds, ncid, flag)
 end if
   end subroutine Restart
 
+
+  !------------------------------------------------------------------------
+#ifdef USE_PETSC_LIB
+  subroutine InitColdGhost(this, bounds_proc)
+    !
+    ! !DESCRIPTION:
+    ! Assign soil properties for ghost/halo columns
+    !
+    ! !USES:
+    use domainLateralMod       , only : ExchangeColumnLevelGhostData
+    use shr_infnan_mod         , only : shr_infnan_isnan
+    use shr_infnan_mod         , only : isnan => shr_infnan_isnan
+    use landunit_varcon        , only : max_lunit
+    !
+    implicit none
+    !
+    ! !ARGUMENTS:
+    class(soilstate_type)            :: this
+    type(bounds_type), intent(in)    :: bounds_proc
+    !
+    integer             :: c,j                     ! indices
+    integer             :: nvals_col               ! number of values per subgrid category
+    integer             :: beg_idx, end_idx        ! begin/end index for accessing values in data_send/data_recv
+    real(r8) , parameter:: FILL_VALUE = -999999.d0 ! temporary
+    real(r8) , pointer  :: data_send_col(:)        ! data sent by local mpi rank
+    real(r8) , pointer  :: data_recv_col(:)        ! data received by local mpi rank
+
+    ! Number of values per soil column
+    nvals_col = 4*nlevgrnd ! (watsat + hksat + bsw + sucsat) * nlevgrnd
+
+    ! Allocate value
+    allocate(data_send_col((bounds_proc%endc     - bounds_proc%begc     + 1)*nvals_col))
+    allocate(data_recv_col((bounds_proc%endc_all - bounds_proc%begc_all + 1)*nvals_col))
+
+    ! Assemble the data to send
+    do c = bounds_proc%begc, bounds_proc%endc
+
+       beg_idx = (c - bounds_proc%begc)*nvals_col
+
+       do j = 1, nlevgrnd
+
+          beg_idx = beg_idx + 1
+          if (.not. isnan(this%watsat_col(c,j)) .and. this%watsat_col(c,j) /= spval) then
+             data_send_col(beg_idx) = this%watsat_col(c,j)
+          else
+             data_send_col(beg_idx) = FILL_VALUE
+          endif
+
+          beg_idx = beg_idx + 1
+          if (.not. isnan(this%hksat_col(c,j)) .and. this%hksat_col(c,j) /= spval) then
+             data_send_col(beg_idx) = this%hksat_col(c,j)
+          else
+             data_send_col(beg_idx) = FILL_VALUE
+          endif
+
+          beg_idx = beg_idx + 1
+          if (.not. isnan(this%bsw_col(c,j)) .and. this%bsw_col(c,j) /= spval) then
+             data_send_col(beg_idx) = this%bsw_col(c,j)
+          else
+             data_send_col(beg_idx) = FILL_VALUE
+          endif
+
+          beg_idx = beg_idx + 1
+          if (.not. isnan(this%sucsat_col(c,j)) .and. this%sucsat_col(c,j) /= spval) then
+             data_send_col(beg_idx) = this%sucsat_col(c,j)
+          else
+             data_send_col(beg_idx) = FILL_VALUE
+          endif
+       enddo
+    enddo
+
+    ! Send the data
+    call ExchangeColumnLevelGhostData(bounds_proc, nvals_col, data_send_col, data_recv_col)
+
+    ! Assign data corresponding to ghost/halo soil columns
+    do c = bounds_proc%endc + 1, bounds_proc%endc_all
+       beg_idx = (c - bounds_proc%begc)*nvals_col
+       do j = 1, nlevgrnd
+          beg_idx = beg_idx + 1
+          this%watsat_col(c,j) = data_recv_col(beg_idx)
+
+          beg_idx = beg_idx + 1
+          this%hksat_col(c,j) = data_recv_col(beg_idx)
+
+          beg_idx = beg_idx + 1
+          this%bsw_col(c,j) = data_recv_col(beg_idx)
+
+          beg_idx = beg_idx + 1
+          this%sucsat_col(c,j) = data_recv_col(beg_idx)
+       enddo
+    enddo
+
+    ! Free up memory
+    deallocate(data_send_col)
+    deallocate(data_recv_col)
+
+  end subroutine InitColdGhost
+
+#else
+
+  !------------------------------------------------------------------------
+  subroutine InitColdGhost(this, bounds_proc)
+    !
+    ! !DESCRIPTION:
+    ! Assign soil properties for ghost/halo columns
+    !
+    ! !USES:
+    implicit none
+    !
+    ! !ARGUMENTS:
+    class(soilstate_type)            :: this
+    type(bounds_type), intent(in)    :: bounds_proc
+
+    character(len=*), parameter :: subname = 'InitColdGhost'
+
+    call endrun(msg='ERROR ' // trim(subname) //': Requires '//&
+         'PETSc, but the code was compiled without -DUSE_PETSC_LIB')
+
+  end subroutine InitColdGhost
+
+#endif
+  !------------------------------------------------------------------------
+
 end module SoilStateType