move hydro solver switch from EDParamsMod to FatesInterfaceTypesMod

biogeophys/FatesPlantHydraulicsMod.F90

@@ -49,7 +49,6 @@ module FatesPlantHydraulicsMod
   use EDParamsMod       , only : hydr_psi0
   use EDParamsMod       , only : hydr_psicap
   use EDParamsMod       , only : hydr_htftype_node
-  use EDParamsMod       , only : hydr_solver
 
   use EDTypesMod        , only : ed_site_type
   use FatesPatchMod     , only : fates_patch_type
@@ -64,6 +63,7 @@ module FatesPlantHydraulicsMod
   use FatesInterfaceTypesMod  , only : hlm_ipedof
   use FatesInterfaceTypesMod  , only : numpft
   use FatesInterfaceTypesMod  , only : nlevsclass
+  use FatesInterfaceTypesMod  , only : hlm_hydr_solver
 
   use FatesAllometryMod, only    : bleaf
   use FatesAllometryMod, only    : bsap_allom
@@ -1419,7 +1419,7 @@ subroutine InitHydrSites(sites,bc_in)
      case(rhizlayer_aggmeth_none)
 
         csite_hydr%nlevrhiz = bc_in(s)%nlevsoil
-        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver,bc_in(s)%nlevsoil)
+        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hlm_hydr_solver,bc_in(s)%nlevsoil)
 
         do j=1,csite_hydr%nlevrhiz
            csite_hydr%map_r2s(j,1) = j
@@ -1431,7 +1431,7 @@ subroutine InitHydrSites(sites,bc_in)
      case(rhizlayer_aggmeth_combine12)
 
         csite_hydr%nlevrhiz     = max(1,bc_in(s)%nlevsoil-1)
-        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver,bc_in(s)%nlevsoil)
+        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hlm_hydr_solver,bc_in(s)%nlevsoil)
 
         csite_hydr%map_r2s(1,1) = 1
         j_bc                 = min(2,bc_in(s)%nlevsoil) ! this protects 1 soil layer
@@ -1449,7 +1449,7 @@ subroutine InitHydrSites(sites,bc_in)
      case(rhizlayer_aggmeth_balN)
 
         csite_hydr%nlevrhiz = min(aggN,bc_in(s)%nlevsoil)
-        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver,bc_in(s)%nlevsoil)
+        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hlm_hydr_solver,bc_in(s)%nlevsoil)
 
         ntoagg = int(ceiling(real(bc_in(s)%nlevsoil)/real(csite_hydr%nlevrhiz)-nearzero))
 
@@ -2583,21 +2583,21 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime)
               ! from leaf to the current soil layer.  This does NOT
               ! update cohort%th_*
 
-              if(hydr_solver == hydr_solver_2DNewton) then
+              if(hlm_hydr_solver == hydr_solver_2DNewton) then
 
                  call MatSolve2D(csite_hydr,ccohort,ccohort_hydr, &
                       dtime,qflx_tran_veg_indiv, &
                       sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, &
                       dth_layershell_col)
 
-              elseif(hydr_solver == hydr_solver_2DPicard) then
+              elseif(hlm_hydr_solver == hydr_solver_2DPicard) then
 
                  call PicardSolve2D(csite_hydr,ccohort,ccohort_hydr, &
                       dtime,qflx_tran_veg_indiv, &
                       sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, &
                       dth_layershell_col,csite_hydr%num_nodes)
 
-              elseif(hydr_solver == hydr_solver_1DTaylor ) then
+              elseif(hlm_hydr_solver == hydr_solver_1DTaylor ) then
 
                  ! ---------------------------------------------------------------------------------
                  ! Approach: do nlevsoi_hyd sequential solutions to Richards' equation,
@@ -4938,7 +4938,7 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, &
 
 
    ! This NaN's the scratch arrays
-   call csite_hydr%FlushSiteScratch(hydr_solver)
+   call csite_hydr%FlushSiteScratch(hlm_hydr_solver)
 
    ! This is the maximum number of iterations needed for this cohort
    ! (each soil layer has a different number, this saves the max)
@@ -5714,7 +5714,7 @@ subroutine PicardSolve2D(csite_hydr,cohort,cohort_hydr, &
        ft           => cohort%pft)
 
     ! This NaN's the scratch arrays
-    call csite_hydr%FlushSiteScratch(hydr_solver)
+    call csite_hydr%FlushSiteScratch(hlm_hydr_solver)
 
     ! This is the maximum number of iterations needed for this cohort
     ! (each soil layer has a different number, this saves the max)

main/EDParamsMod.F90

@@ -203,18 +203,6 @@ module EDParamsMod
    real(r8),protected,public :: hydr_psicap        !  sapwood water potential at which capillary reserves exhausted (MPa)
    character(len=param_string_length),parameter,public :: hydr_name_psicap = "fates_hydro_psicap"
 
-
-   ! Switch that defines which hydraulic solver to use
-   ! 1 = Taylor solution that solves plant fluxes with 1 layer
-   !     sequentially placing solution on top of previous layer solves
-   ! 2 = Picard solution that solves all fluxes in a plant and
-   !     the soil simultaneously, 2D: soil x (root + shell)
-   ! 3 = Newton-Raphson (Deprecated) solution that solves all fluxes in a plant and
-   !     the soil simultaneously, 2D: soil x (root + shell)
-
-   integer,protected,public :: hydr_solver        !  switch designating hydraulics numerical solver
-   character(len=param_string_length),parameter,public :: hydr_name_solver = "fates_hydro_solver"
-
    !Soil BGC parameters, mostly used for testing FATES when not coupled to the dynamics bgc hlm
    ! ----------------------------------------------------------------------------------------------
    real(r8),protected,public :: bgc_soil_salinity ! site-level soil salinity for FATES when not coupled to dynamic soil BGC of salinity
@@ -350,7 +338,6 @@ subroutine FatesParamsInit()
     hydr_kmax_rsurf2                      = nan
     hydr_psi0                             = nan
     hydr_psicap                           = nan
-    hydr_solver                           = -9
     bgc_soil_salinity                     = nan
     logging_dbhmin                        = nan
     logging_dbhmax                        = nan
@@ -506,9 +493,6 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=maxcohort_name, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
-    call fates_params%RegisterParameter(name=hydr_name_solver, dimension_shape=dimension_shape_scalar, &
-         dimension_names=dim_names_scalar)
-
     call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf1, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -745,10 +729,6 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetrieveParameter(name=hydr_name_psicap, &
           data=hydr_psicap)
 
-    call fates_params%RetrieveParameter(name=hydr_name_solver, &
-         data=tmpreal)
-    hydr_solver = nint(tmpreal)
-
     call fates_params%RetrieveParameter(name=bgc_name_soil_salinity, &
           data=bgc_soil_salinity)	  
 
@@ -905,7 +885,6 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),fmt0) 'hydro_kmax_rsurf2 = ',hydr_kmax_rsurf2  
         write(fates_log(),fmt0) 'hydro_psi0 = ',hydr_psi0
         write(fates_log(),fmt0) 'hydro_psicap = ',hydr_psicap
-        write(fates_log(),fmt0) 'hydro_solver = ',hydr_solver
         write(fates_log(),fmt0) 'bgc_soil_salinity = ', bgc_soil_salinity
         write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin
         write(fates_log(),fmt0) 'logging_dbhmax = ',logging_dbhmax

main/FatesInterfaceMod.F90

@@ -1475,6 +1475,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
          hlm_num_luh2_transitions  = unset_int
          hlm_use_cohort_age_tracking = unset_int
          hlm_daylength_factor_switch = unset_int
+         hlm_hydr_solver = unset_int
          hlm_use_logging   = unset_int
          hlm_use_ed_st3    = unset_int
          hlm_use_ed_prescribed_phys = unset_int
@@ -1764,6 +1765,11 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if
 
+         if(hlm_hydr_solver .eq. unset_int) then
+            write(fates_log(), *) 'FATES hydro solver is unset: hlm_hydr_solver, exiting'
+            call endrun(msg=errMsg(sourcefile, __LINE__))
+         end if
+
          if(hlm_use_sp.eq.itrue.and.hlm_use_nocomp.eq.ifalse)then
             write(fates_log(), *) 'SP cannot be on if nocomp mode is off. Exiting. '
             call endrun(msg=errMsg(sourcefile, __LINE__))
@@ -1973,6 +1979,12 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
                   write(fates_log(),*) 'Transfering hlm_daylength_factor_switch= ',ival,' to FATES'
                end if
 
+            case('hydr_solver')
+               hlm_hydr_solver = ival
+               if (fates_global_verbose()) then
+                  write(fates_log(),*) 'Transfering hlm_hydr_solver= ',ival,' to FATES'
+               end if
+
             case('use_logging')
                hlm_use_logging = ival
                if (fates_global_verbose()) then

main/FatesInterfaceTypesMod.F90

@@ -161,6 +161,14 @@ module FatesInterfaceTypesMod
    integer, public :: hlm_daylength_factor_switch ! This switch enables the use of the daylength factor from the HLM
                                                   ! 1 = TRUE, 0 = FALSE
 
+   integer, public :: hlm_hydr_solver             ! Switch that defines which hydraulic solver to use
+                                                  ! 1 = Taylor solution that solves plant fluxes with 1 layer
+                                                  !     sequentially placing solution on top of previous layer solves
+                                                  ! 2 = Picard solution that solves all fluxes in a plant and
+                                                  !     the soil simultaneously, 2D: soil x (root + shell)
+                                                  ! 3 = Newton-Raphson (Deprecated) solution that solves all fluxes in a plant and
+                                                  !     the soil simultaneously, 2D: soil x (root + shell)
+
    integer, public :: hlm_use_ed_st3              ! This flag signals whether or not to use
                                                   ! (ST)atic (ST)and (ST)ructure mode (ST3)
                                                   ! Essentially, this gives us the ability