Merge branch 'mpaiao-pr-zsoifl' (PR #7319)

Add an alternate way of defining soil layers and depth to bedrock in ELM.
These changes were implemented with the intention of being bit-for-bit with previously generated surface files.
1. The code no longer expects nlevsoi in the surface file to be 10 (the default in ELM). Instead, it will use nlevsoi provided in the surface file, as long as nlevsoi <= nlevgrnd.
2. If ZSOI(nlevsoi) is provided in the surface file, ELM will use these layers instead of the default, and append layers
between nlevsoi+1 and nlevgrnd if needed using the default parameters, plus an offset to ensure that zsoi always increases.
4. The scaling parameters for the soil layers are now declared in elm_varpar.F90, to avoid magic numbers.
5. Slightly revised the code that reads the grid-cell dependent depth to bedrock, so it works with either aveDTB (the original name in ELM) and zbedrock (the name in CLM), so surface files borrowed from CLM would still work with ELM and generate the same soil layers.
6. Deleted some unused variables and (I think) unecessary calls to check_dim for nlevsoi so the surface file can have different values.

[BFB]

Author: bishtgautam

components/elm/src/biogeophys/SoilStateType.F90

@@ -11,6 +11,7 @@ module SoilStateType
   use ncdio_pio       , only : ncd_pio_openfile, ncd_inqfdims, ncd_pio_closefile, ncd_inqdid, ncd_inqdlen
   use elm_varpar      , only : more_vertlayers, numpft, numrad
   use elm_varpar      , only : nlevsoi, nlevgrnd, nlevlak, nlevsoifl, nlayer, nlayert, nlevurb, nlevsno
+  use elm_varpar      , only : scalez, zecoeff
   use landunit_varcon , only : istice, istdlak, istwet, istsoil, istcrop, istice_mec
   use column_varcon   , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv 
   use elm_varcon      , only : zsoi, dzsoi, zisoi, spval, namet, grlnd
@@ -440,6 +441,11 @@ subroutine InitCold(this, bounds)
     call getfil (fsurdat, locfn, 0)
     call ncd_pio_openfile (ncid, locfn, 0)
 
+    ! --------------------------------------------------------------------
+    !    Make sure nlevsoifl and nlevsoi match.  At this point, we keep this test, but
+    ! initVertical should have taken care of nlevsoi when the  value in the input file
+    ! differs from the default (10 layers).
+    ! --------------------------------------------------------------------
     call ncd_inqdlen(ncid,dimid,nlevsoifl,name='nlevsoi')
     if ( .not. more_vertlayers )then
        if ( nlevsoifl /= nlevsoi )then
@@ -450,6 +456,38 @@ subroutine InitCold(this, bounds)
        ! read in layers, interpolate to high resolution grid later
     end if
 
+
+
+    ! --------------------------------------------------------------------
+    !    Define the soil layers from the input file.  We first check if ZSOI is available
+    ! in the file, in which case we read the information directly from the file. Otherwise,
+    ! we assume the original ELM parameters. The input soil depths will be used for 
+    ! interpolating soil properties (e.g., sand and clay).
+    ! --------------------------------------------------------------------
+    allocate(zsoifl(1:nlevsoifl), zisoifl(0:nlevsoifl), dzsoifl(1:nlevsoifl))
+    ! Try to read soil information from the file.
+    call ncd_io(ncid=ncid, varname='ZSOI', flag='read', data=zsoifl, dim1name=grlnd, readvar=readvar)
+    if (.not. readvar ) then
+       !    Variable ZSOI not found, use the ELM parameters.
+       do j = 1, nlevsoifl
+          zsoifl(j) = scalez*(exp(zecoeff*(j-0.5_r8))-1._r8)    !node depths
+       end do
+    end if
+
+    dzsoifl(1) = 0.5_r8*(zsoifl(1)+zsoifl(2))             !thickness b/n two interfaces
+    do j = 2,nlevsoifl-1
+       dzsoifl(j)= 0.5_r8*(zsoifl(j+1)-zsoifl(j-1))
+    enddo
+    dzsoifl(nlevsoifl) = zsoifl(nlevsoifl)-zsoifl(nlevsoifl-1)
+
+    zisoifl(0) = 0._r8
+    do j = 1, nlevsoifl-1
+       zisoifl(j) = 0.5_r8*(zsoifl(j)+zsoifl(j+1))         !interface depths
+    enddo
+    zisoifl(nlevsoifl) = zsoifl(nlevsoifl) + 0.5_r8*dzsoifl(nlevsoifl)
+
+
+
     ! Read in organic matter dataset
 
     organic_max = ParamsShareInst%organic_max
@@ -548,27 +586,6 @@ subroutine InitCold(this, bounds)
 
     call ncd_pio_closefile(ncid)
 
-    ! --------------------------------------------------------------------
-    ! get original soil depths to be used in interpolation of sand and clay
-    ! --------------------------------------------------------------------
-
-    allocate(zsoifl(1:nlevsoifl), zisoifl(0:nlevsoifl), dzsoifl(1:nlevsoifl))
-    do j = 1, nlevsoifl
-       zsoifl(j) = 0.025*(exp(0.5_r8*(j-0.5_r8))-1._r8)    !node depths
-    enddo
-
-    dzsoifl(1) = 0.5_r8*(zsoifl(1)+zsoifl(2))             !thickness b/n two interfaces
-    do j = 2,nlevsoifl-1
-       dzsoifl(j)= 0.5_r8*(zsoifl(j+1)-zsoifl(j-1))
-    enddo
-    dzsoifl(nlevsoifl) = zsoifl(nlevsoifl)-zsoifl(nlevsoifl-1)
-
-    zisoifl(0) = 0._r8
-    do j = 1, nlevsoifl-1
-       zisoifl(j) = 0.5_r8*(zsoifl(j)+zsoifl(j+1))         !interface depths
-    enddo
-    zisoifl(nlevsoifl) = zsoifl(nlevsoifl) + 0.5_r8*dzsoifl(nlevsoifl)
-
     ! --------------------------------------------------------------------
     ! Set soil hydraulic and thermal properties: non-lake
     ! --------------------------------------------------------------------

components/elm/src/main/elm_varpar.F90

@@ -72,6 +72,17 @@ module elm_varpar
 
   integer, parameter :: nlevslp = 11          ! number of slope percentile levels
 
+  !   Parameters that define the thickness of the soil layers. By default, we define the 
+  ! mid-point of the soil layers as:
+  !
+  !    zsoi(j) = scalez * ( exp( zecoeff * (j - 0.5_r8) ) - 1.0_r8 )
+  !
+  !   The default values scalez = 0.025_r8 and zecoeff = 0.50_r8 trace back to the original
+  ! ELM configuration, but depending on the specific needs by the user, these can be
+  ! modified.
+  real(r8), parameter :: scalez  = 0.025_r8
+  real(r8), parameter :: zecoeff = 0.50_r8
+
   ! constants for decomposition cascade
 
   integer :: i_met_lit 

components/elm/src/main/initVerticalMod.F90

@@ -14,7 +14,8 @@ module initVerticalMod
   use spmdMod        , only : masterproc
   use elm_varpar     , only : more_vertlayers, nlevsno, nlevgrnd, nlevlak
   use elm_varpar     , only : toplev_equalspace, nlev_equalspace
-  use elm_varpar     , only : nlevsoi, nlevsoifl, nlevurb, nlevslp 
+  use elm_varpar     , only : nlevsoi, nlevsoifl, nlevurb, nlevslp
+  use elm_varpar     , only : nlevdecomp, scalez, zecoeff
   use elm_varctl     , only : fsurdat, iulog, use_var_soil_thick
   use elm_varctl     , only : use_vancouver, use_mexicocity, use_vertsoilc, use_extralakelayers, use_extrasnowlayers
   use elm_varctl     , only : use_erosion, use_polygonal_tundra
@@ -26,7 +27,7 @@ module initVerticalMod
   use ColumnType     , only : col_pp
   use ColumnDataType , only : col_ws
   use SnowHydrologyMod, only : InitSnowLayers
-  use ncdio_pio
+  use ncdio_pio       , only : file_desc_t, ncd_io, ncd_pio_openfile, ncd_pio_closefile , ncd_inqdlen
   use topounit_varcon  , only : max_topounits
   use GridcellType     , only : grc_pp
   !
@@ -50,12 +51,16 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
     real(r8)            , intent(in)    :: thick_wall(bounds%begl:)
     real(r8)            , intent(in)    :: thick_roof(bounds%begl:)
     !
-    ! LOCAL VARAIBLES:
-    integer               :: c,l,t,ti,topi,g,i,j,lev     ! indices 
+    ! LOCAL PARAMETERS
+    integer, parameter :: ndtbname = 2 ! Number of names to try for depth to bedrock
+    !
+    ! LOCAL VARIABLES:
+    integer               :: c,l,t,ti,topi,g,i,j,lev,n     ! indices 
     type(file_desc_t)     :: ncid              ! netcdf id
-    logical               :: readvar 
+    logical               :: readvar           ! Flag: variable was successfully read
     integer               :: dimid             ! dimension id
     character(len=256)    :: locfn             ! local filename
+    real(r8) ,pointer     :: zsoi_in(:)        ! read in - soil information
     real(r8) ,pointer     :: std (:)           ! read in - topo_std
     real(r8) ,pointer     :: tslope (:)        ! read in - topo_slope
     real(r8) ,pointer     :: gradz(:)          ! read in - gradz (polygonal tundra only)
@@ -68,7 +73,6 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
     real(r8)              :: slopebeta         ! temporary
     real(r8)              :: slopemax          ! temporary
     integer               :: ier               ! error status
-    real(r8)              :: scalez = 0.025_r8 ! Soil layer thickness discretization (m)
     real(r8)              :: thick_equal = 0.2
     real(r8) ,pointer     :: lakedepth_in(:,:)   ! read in - lakedepth 
     real(r8), allocatable :: zurb_wall(:,:)    ! wall (layer node depth)
@@ -80,6 +84,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
     real(r8)              :: depthratio        ! ratio of lake depth to standard deep lake depth 
     integer               :: begc, endc
     integer               :: begl, endl
+    character(len=256), dimension(ndtbname) :: dtbname ! List of possible names for depth to bedrock
     !------------------------------------------------------------------------
 
     begc = bounds%begc; endc= bounds%endc
@@ -95,47 +100,81 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
     call ncd_pio_openfile (ncid, locfn, 0)
 
     call ncd_inqdlen(ncid, dimid, nlevsoifl, name='nlevsoi')
-    if ( .not. more_vertlayers )then
-       if ( nlevsoifl /= nlevsoi )then
-          call shr_sys_abort(' ERROR: Number of soil layers on file does NOT match the number being used'//&
+
+
+    if ( .not. more_vertlayers ) then
+       !    Removed the requirement for nlevsoifl to match nlevsoi, but we make sure the
+       ! number of input layers do not exceed the maximum allocated, to avoid segmentation
+       ! violation errors.
+       if ( nlevsoifl > nlevgrnd ) then
+          call shr_sys_abort(' ERROR: Number of soil layers on file exceeds the maximum number of layers allowed (nlevgrnd)'//&
                errMsg(__FILE__, __LINE__))
+       elseif ( nlevsoifl /= nlevsoi ) then
+          !   Surface file has a different number of soil levels, update the simulation to
+          ! match the surface file.
+          nlevsoi    = nlevsoifl
+          if (use_vertsoilc) nlevdecomp = nlevsoifl
        end if
     else
        ! read in layers, interpolate to high resolution grid later
     end if
 
     ! --------------------------------------------------------------------
-    ! Define layer structure for soil, lakes, urban walls and roof 
-    ! Vertical profile of snow is not initialized here - but below
-    ! --------------------------------------------------------------------
-    
-    ! Soil layers and interfaces (assumed same for all non-lake patches)
+    !    Define layer structure for soil, lakes, urban walls and roof. We first check
+    ! whether or not soil layers exist in the surfacefile. If so, we use the layers from
+    ! the file. Otherwise, we define the layers using the default parameters, but
+    ! further checking if the run should include intermediate layers with equal 
+    ! thicknesses.
+    !     In soil layers and interfaces (assumed same for all non-lake patches),
     ! "0" refers to soil surface and "nlevsoi" refers to the bottom of model soil
-    
-    if ( more_vertlayers )then
+    ! 
+    ! Note: vertical profile of snow is not initialized here - but below
+    ! --------------------------------------------------------------------
+    ! Try to read soil information from the file.
+    allocate (zsoi_in(nlevsoi))
+    call ncd_io(ncid=ncid, varname='ZSOI', flag='read', data=zsoi_in, dim1name=grlnd, readvar=readvar)
+    if ( readvar ) then
+       ! -----------------------------------------------------------------
+       !    File contains soil information. We must complete the soil depth information for
+       ! layers beneath nlevsoi, using the original scaling parameters for increasing the
+       ! depth of the layers, but acknowledging that the layer must be beneath the deepest
+       ! input soil layer.
+       ! -----------------------------------------------------------------
+       zsoi(1:nlevsoi) = zsoi_in(1:nlevsoi)
+       do j = nlevsoi+1, nlevgrnd
+          zsoi(j) = zsoi(nlevsoi) + &
+             scalez*(exp(zecoeff*(j -0.5_r8))-exp(zecoeff*(nlevsoi-0.5_r8)))
+       end do
+
+
+    elseif ( more_vertlayers )then
        ! replace standard exponential grid with a grid that starts out exponential, 
        ! then has several evenly spaced layers, then finishes off exponential. 
        ! this allows the upper soil to behave as standard, but then continues 
        ! with higher resolution to a deeper depth, so that, for example, permafrost
        ! dynamics are not lost due to an inability to resolve temperature, moisture, 
        ! and biogeochemical dynamics at the base of the active layer
        do j = 1, toplev_equalspace
-          zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8)    !node depths
+          zsoi(j) = scalez*(exp(zecoeff*(j-0.5_r8))-1._r8)    !node depths
        enddo
 
        do j = toplev_equalspace+1,toplev_equalspace + nlev_equalspace
           zsoi(j) = zsoi(j-1) + thick_equal
        enddo
 
        do j = toplev_equalspace + nlev_equalspace +1, nlevgrnd
-          zsoi(j) = scalez*(exp(0.5_r8*((j - nlev_equalspace)-0.5_r8))-1._r8) + nlev_equalspace * thick_equal
+          zsoi(j) = scalez*(exp(zecoeff*((j - nlev_equalspace)-0.5_r8))-1._r8) + nlev_equalspace * thick_equal
        enddo
     else
-
+       ! -----------------------------------------------------------------
+       !    Soil layers not available from the input, and no additional layers needed. Use the
+       ! default soil thickness settings.
+       ! -----------------------------------------------------------------
        do j = 1, nlevgrnd
-          zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8)    !node depths
+          zsoi(j) = scalez*(exp(zecoeff*(j-0.5_r8))-1._r8)    !node depths
        enddo
     end if
+    deallocate(zsoi_in)
 
     dzsoi(1) = 0.5_r8*(zsoi(1)+zsoi(2))             !thickness b/n two interfaces
     do j = 2,nlevgrnd-1
@@ -629,9 +668,22 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
 
       if (use_var_soil_thick) then
          allocate(dtb(bounds%begg:bounds%endg,1:max_topounits))
-         call ncd_io(ncid=ncid, varname='aveDTB', flag='read', data=dtb, dim1name=grlnd, readvar=readvar)
+
+         !    ELM and CLM use different names for depth to bedrock in the surface file
+         ! ('aveDTB' and 'zbedrock', respectively).  To keep cross-model compatibility, we 
+         ! test both names before falling back to the default number of layers. 
+         dtbname(1) = 'aveDTB'
+         dtbname(2) = 'zbedrock'
+         readvar = .false.
+         do n=1,ndtbname
+            if (.not. readvar) then
+               call ncd_io(ncid=ncid, varname=dtbname(n), flag='read', data=dtb, dim1name=grlnd, readvar=readvar)
+            end if
+         end do
+
+
          if (.not. readvar) then
-            write(iulog,*) 'aveDTB not in surfdata: reverting to default 10 layers.'
+            write(iulog,fmt='(a,i5,a)') 'aveDTB not in surfdata: reverting to default ',nlevsoi,' layers.'
             do c = begc,endc
                col_pp%nlevbed(c) = nlevsoi
 	       col_pp%zibed(c) = zisoi(nlevsoi)

components/elm/src/main/surfrdMod.F90

@@ -9,7 +9,7 @@ module surfrdMod
   use shr_kind_mod    , only : r8 => shr_kind_r8
   use shr_log_mod     , only : errMsg => shr_log_errMsg
   use abortutils      , only : endrun
-  use elm_varpar      , only : nlevsoifl, numpft, numcft
+  use elm_varpar      , only : numpft, numcft
   use landunit_varcon , only : numurbl
   use elm_varcon      , only : grlnd
   use elm_varctl      , only : iulog, scmlat, scmlon, single_column, firrig_data
@@ -712,7 +712,6 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
     ! !LOCAL VARIABLES:
     integer  :: n,nl,nurb,g, t,tm,ti                ! indices
     integer  :: dimid,varid                ! netCDF id's
-    real(r8) :: nlevsoidata(nlevsoifl)
     logical  :: found                      ! temporary for error check
     integer  :: nindx                      ! temporary for error check
     integer  :: ier                        ! error status
@@ -750,7 +749,6 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
     allocate(pctglc_mec_tot(begg:endg,1:max_topounits))
     allocate(pctspec(begg:endg,1:max_topounits))
 
-    call check_dim(ncid, 'nlevsoi', nlevsoifl)
 
        ! Obtain non-grid surface properties of surface dataset other than percent pft
 
@@ -1533,7 +1531,6 @@ subroutine surfrd_topounit_data(begg, endg, lfsurdat)
     call getfil( lfsurdat, locfn, 0 )
     call ncd_pio_openfile (ncid, trim(locfn), 0)
 
-    !call check_dim(ncid, 'nlevsoi', nlevsoifl)
     call check_var(ncid=ncid, varname='MaxTopounitElv', vardesc=vardesc, readvar=readvar)
     if (readvar) then
        call ncd_io(ncid=ncid, varname='MaxTopounitElv', flag='read', data=maxTopoElv, &