module initVerticalMod

#include "shr_assert.h"

  !-----------------------------------------------------------------------
  ! !DESCRIPTION:
  ! Initialize vertical components of column datatype
  !
  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 shr_sys_mod       , only : shr_sys_abort
  use decompMod         , only : bounds_type
  use spmdMod           , only : masterproc
  use clm_varpar        , only : nlevsno, nlevgrnd, nlevlak
  use clm_varpar        , only : toplev_equalspace, nlev_equalspace
  use clm_varpar        , only : nlevsoi, nlevsoifl, nlevurb 
  use clm_varctl        , only : fsurdat, iulog
  use clm_varctl        , only : use_vancouver, use_mexicocity, use_vertsoilc, use_extralakelayers
  use clm_varctl        , only : use_bedrock, soil_layerstruct
  use clm_varctl        , only : use_ed
  use clm_varcon        , only : zlak, dzlak, zsoi, dzsoi, zisoi, dzsoi_decomp, spval, ispval, grlnd 
  use column_varcon     , only : icol_roof, icol_sunwall, icol_shadewall, is_hydrologically_active
  use landunit_varcon   , only : istdlak, istice_mec
  use fileutils         , only : getfil
  use LandunitType      , only : lun                
  use GridcellType      , only : grc                
  use ColumnType        , only : col                
  use glcBehaviorMod    , only : glc_behavior_type
  use SnowHydrologyMod  , only : InitSnowLayers             
  use EDTypesMod        , only : ed_hist_scpfmaps
  use abortUtils        , only : endrun    
  use ncdio_pio
  !
  ! !PUBLIC TYPES:
  implicit none
  save
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  public :: initVertical
  ! !PRIVATE MEMBER FUNCTIONS:
  private :: ReadNL
  private :: hasBedrock  ! true if the given column type includes bedrock layers
  !

  character(len=*), parameter, private :: sourcefile = &
       __FILE__

  !
  !------------------------------------------------------------------------

contains

  !------------------------------------------------------------------------
  subroutine ReadNL( )
    !
    ! !DESCRIPTION:
    ! Read namelist for SoilStateType
    !
    ! !USES:
    use shr_mpi_mod    , only : shr_mpi_bcast
    use shr_log_mod    , only : errMsg => shr_log_errMsg
    use fileutils      , only : getavu, relavu, opnfil
    use clm_nlUtilsMod , only : find_nlgroup_name
    use clm_varctl     , only : iulog
    use spmdMod        , only : mpicom, masterproc
    use controlMod     , only : NLFilename
    !
    ! !ARGUMENTS:
    !
    ! !LOCAL VARIABLES:
    integer :: ierr                 ! error code
    integer :: unitn                ! unit for namelist file
    character(len=32) :: subname = 'InitVertical_readnl'  ! subroutine name
    !-----------------------------------------------------------------------

    character(len=*), parameter :: nl_name  = 'clm_inparm'  ! Namelist name
                                                                      
    ! MUST agree with name in namelist and read
    namelist /clm_inparm/ use_bedrock

    ! preset values

    use_bedrock = .false.

    if ( masterproc )then

       unitn = getavu()
       write(iulog,*) 'Read in '//nl_name//' namelist'
       call opnfil (NLFilename, unitn, 'F')
       call find_nlgroup_name(unitn, nl_name, status=ierr)
       if (ierr == 0) then
          read(unit=unitn, nml=clm_inparm, iostat=ierr)
          if (ierr /= 0) then
             call endrun(msg="ERROR reading '//nl_name//' namelist"//errmsg(sourcefile, __LINE__))
          end if
       else
          call endrun(msg="ERROR finding '//nl_name//' namelist"//errmsg(sourcefile, __LINE__))
       end if
       call relavu( unitn )

    end if

    call shr_mpi_bcast(use_bedrock, mpicom)

  end subroutine ReadNL

  !------------------------------------------------------------------------
  subroutine initVertical(bounds, glc_behavior, snow_depth, thick_wall, thick_roof)
    use clm_varcon, only : zmin_bedrock
    !
    ! !ARGUMENTS:
    type(bounds_type)   , intent(in)    :: bounds
    type(glc_behavior_type), intent(in) :: glc_behavior
    real(r8)            , intent(in)    :: snow_depth(bounds%begc:)
    real(r8)            , intent(in)    :: thick_wall(bounds%begl:)
    real(r8)            , intent(in)    :: thick_roof(bounds%begl:)
    !
    ! LOCAL VARAIBLES:
    integer               :: c,l,g,i,j,lev     ! indices 
    type(file_desc_t)     :: ncid              ! netcdf id
    logical               :: readvar 
    integer               :: dimid             ! dimension id
    character(len=256)    :: locfn             ! local filename
    real(r8) ,pointer     :: std (:)           ! read in - topo_std 
    real(r8) ,pointer     :: tslope (:)        ! read in - topo_slope 
    real(r8)              :: slope0            ! temporary
    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     :: zbedrock_in(:)   ! read in - z_bedrock
    real(r8) ,pointer     :: lakedepth_in(:)   ! read in - lakedepth 
    real(r8), allocatable :: zurb_wall(:,:)    ! wall (layer node depth)
    real(r8), allocatable :: zurb_roof(:,:)    ! roof (layer node depth)
    real(r8), allocatable :: dzurb_wall(:,:)   ! wall (layer thickness)
    real(r8), allocatable :: dzurb_roof(:,:)   ! roof (layer thickness)
    real(r8), allocatable :: ziurb_wall(:,:)   ! wall (layer interface)
    real(r8), allocatable :: ziurb_roof(:,:)   ! roof (layer interface)
    real(r8)              :: depthratio        ! ratio of lake depth to standard deep lake depth 
    integer               :: begc, endc
    integer               :: begl, endl
    integer               :: jmin_bedrock

    ! Possible values for levgrnd_class. The important thing is that, for a given column,
    ! layers that are fundamentally different (e.g., soil vs bedrock) have different
    ! values. This information is used in the vertical interpolation in init_interp.
    !
    ! IMPORTANT: These values should not be changed lightly. e.g., try to avoid changing
    ! the values assigned to LEVGRND_CLASS_STANDARD, LEVGRND_CLASS_DEEP_BEDROCK, etc.  The
    ! problem with changing these is that init_interp expects that layers with a value of
    ! (e.g.) 1 on the source file correspond to layers with a value of 1 on the
    ! destination file. So if you change the values of these constants, you either need to
    ! adequately inform users of this change, or build in some translation mechanism in
    ! init_interp (such as via adding more metadata to the restart file on the meaning of
    ! these different values).
    !
    ! The distinction between "shallow" and "deep" bedrock is not made explicitly
    ! elsewhere. But, since these classes have somewhat different behavior, they are
    ! distinguished explicitly here.
    integer, parameter :: LEVGRND_CLASS_STANDARD        = 1
    integer, parameter :: LEVGRND_CLASS_DEEP_BEDROCK    = 2
    integer, parameter :: LEVGRND_CLASS_SHALLOW_BEDROCK = 3
    !------------------------------------------------------------------------

    begc = bounds%begc; endc= bounds%endc
    begl = bounds%begl; endl= bounds%endl

    SHR_ASSERT_ALL((ubound(snow_depth)  == (/endc/)), errMsg(sourcefile, __LINE__))
    SHR_ASSERT_ALL((ubound(thick_wall)  == (/endl/)), errMsg(sourcefile, __LINE__))
    SHR_ASSERT_ALL((ubound(thick_roof)  == (/endl/)), errMsg(sourcefile, __LINE__))

    ! Open surface dataset to read in data below 

    call getfil (fsurdat, locfn, 0)
    call ncd_pio_openfile (ncid, locfn, 0)

    ! --------------------------------------------------------------------
    ! 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)
    ! "0" refers to soil surface and "nlevsoi" refers to the bottom of model soil

    if ( soil_layerstruct == '10SL_3.5m' ) then 
       do j = 1, nlevgrnd
          zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8)    !node depths
       enddo

       dzsoi(1) = 0.5_r8*(zsoi(1)+zsoi(2))             !thickness b/n two interfaces
       do j = 2,nlevgrnd-1
          dzsoi(j)= 0.5_r8*(zsoi(j+1)-zsoi(j-1))
       enddo
       dzsoi(nlevgrnd) = zsoi(nlevgrnd)-zsoi(nlevgrnd-1)

       zisoi(0) = 0._r8
       do j = 1, nlevgrnd-1
          zisoi(j) = 0.5_r8*(zsoi(j)+zsoi(j+1))         !interface depths
       enddo
       zisoi(nlevgrnd) = zsoi(nlevgrnd) + 0.5_r8*dzsoi(nlevgrnd)

    else if ( soil_layerstruct == '23SL_3.5m' )then
       ! Soil layer structure that starts with standard exponential
       ! and 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
       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
       enddo

       dzsoi(1) = 0.5_r8*(zsoi(1)+zsoi(2))             !thickness b/n two interfaces
       do j = 2,nlevgrnd-1
          dzsoi(j)= 0.5_r8*(zsoi(j+1)-zsoi(j-1))
       enddo
       dzsoi(nlevgrnd) = zsoi(nlevgrnd)-zsoi(nlevgrnd-1)

       zisoi(0) = 0._r8
       do j = 1, nlevgrnd-1
       zisoi(j) = 0.5_r8*(zsoi(j)+zsoi(j+1))         !interface depths
       enddo
       zisoi(nlevgrnd) = zsoi(nlevgrnd) + 0.5_r8*dzsoi(nlevgrnd)

    else if ( soil_layerstruct == '49SL_10m' ) then
       !scs: 10 meter soil column, nlevsoi set to 49 in clm_varpar
       do j = 1,10
          dzsoi(j)= 1.e-2_r8     !10mm layers
       enddo
       do j = 11,19
          dzsoi(j)= 1.e-1_r8     !100 mm layers
       enddo
       do j = 20,nlevsoi+1       !300 mm layers
          dzsoi(j)= 3.e-1_r8
       enddo
       do j = nlevsoi+2,nlevgrnd !10 meter bedrock layers
          dzsoi(j)= 10._r8
       enddo
       
       zisoi(0) = 0._r8
       do j = 1,nlevgrnd
          zisoi(j)= sum(dzsoi(1:j))
       enddo
       
       do j = 1, nlevgrnd
          zsoi(j) = 0.5*(zisoi(j-1) + zisoi(j))
       enddo

    else if ( soil_layerstruct == '20SL_8.5m' ) then
       do j = 1,4
          dzsoi(j)= j*0.02_r8          ! linear increase in layer thickness of 2cm each layer
       enddo
       do j = 5,13
          dzsoi(j)= dzsoi(4)+(j-4)*0.04_r8      ! linear increase in layer thickness of 2cm each layer
       enddo
       do j = 14,nlevsoi       
          dzsoi(j)= dzsoi(13)+(j-13)*0.10_r8     ! linear increase in layer thickness of 2cm each layer
       enddo
       do j = nlevsoi+1,nlevgrnd !bedrock layers
          dzsoi(j)= dzsoi(nlevsoi)+(((j-nlevsoi)*25._r8)**1.5_r8)/100._r8  ! bedrock layers
       enddo
       
       zisoi(0) = 0._r8
       do j = 1,nlevgrnd
          zisoi(j)= sum(dzsoi(1:j))
       enddo
       
       do j = 1, nlevgrnd
          zsoi(j) = 0.5*(zisoi(j-1) + zisoi(j))
       enddo
    end if

    ! define a vertical grid spacing such that it is the normal dzsoi if
    ! nlevdecomp =nlevgrnd, or else 1 meter
    if (use_vertsoilc) then
       dzsoi_decomp = dzsoi            !thickness b/n two interfaces
    else
       dzsoi_decomp(1) = 1.
    end if

    if (masterproc) then
       write(iulog, *) 'zsoi', zsoi(:) 
       write(iulog, *) 'zisoi: ', zisoi(:)
       write(iulog, *) 'dzsoi: ', dzsoi(:)
       write(iulog, *) 'dzsoi_decomp: ',dzsoi_decomp
    end if

    if (nlevurb > 0) then
       allocate(zurb_wall(bounds%begl:bounds%endl,nlevurb),    &
            zurb_roof(bounds%begl:bounds%endl,nlevurb),    &
            dzurb_wall(bounds%begl:bounds%endl,nlevurb),   &
            dzurb_roof(bounds%begl:bounds%endl,nlevurb),   &
            ziurb_wall(bounds%begl:bounds%endl,0:nlevurb), &
            ziurb_roof(bounds%begl:bounds%endl,0:nlevurb), &
            stat=ier)
       if (ier /= 0) then
          call shr_sys_abort(' ERROR allocation error for '//&
               'zurb_wall,zurb_roof,dzurb_wall,dzurb_roof,ziurb_wall,ziurb_roof'//&
               errMsg(sourcefile, __LINE__))
       end if
    end if

    if(use_ed)then
      call ed_hist_scpfmaps
    end if


    ! Column level initialization for urban wall and roof layers and interfaces
    do l = bounds%begl,bounds%endl

       ! "0" refers to urban wall/roof surface and "nlevsoi" refers to urban wall/roof bottom
       if (lun%urbpoi(l)) then
          if (use_vancouver) then       
             zurb_wall(l,1) = 0.010_r8/2._r8
             zurb_wall(l,2) = zurb_wall(l,1) + 0.010_r8/2._r8 + 0.020_r8/2._r8
             zurb_wall(l,3) = zurb_wall(l,2) + 0.020_r8/2._r8 + 0.070_r8/2._r8
             zurb_wall(l,4) = zurb_wall(l,3) + 0.070_r8/2._r8 + 0.070_r8/2._r8
             zurb_wall(l,5) = zurb_wall(l,4) + 0.070_r8/2._r8 + 0.030_r8/2._r8

             zurb_roof(l,1) = 0.010_r8/2._r8
             zurb_roof(l,2) = zurb_roof(l,1) + 0.010_r8/2._r8 + 0.010_r8/2._r8
             zurb_roof(l,3) = zurb_roof(l,2) + 0.010_r8/2._r8 + 0.010_r8/2._r8
             zurb_roof(l,4) = zurb_roof(l,3) + 0.010_r8/2._r8 + 0.010_r8/2._r8
             zurb_roof(l,5) = zurb_roof(l,4) + 0.010_r8/2._r8 + 0.030_r8/2._r8

             dzurb_wall(l,1) = 0.010_r8
             dzurb_wall(l,2) = 0.020_r8
             dzurb_wall(l,3) = 0.070_r8
             dzurb_wall(l,4) = 0.070_r8
             dzurb_wall(l,5) = 0.030_r8
             write(iulog,*)'Total thickness of wall: ',sum(dzurb_wall(l,:))
             write(iulog,*)'Wall layer thicknesses: ',dzurb_wall(l,:)

             dzurb_roof(l,1) = 0.010_r8
             dzurb_roof(l,2) = 0.010_r8
             dzurb_roof(l,3) = 0.010_r8
             dzurb_roof(l,4) = 0.010_r8
             dzurb_roof(l,5) = 0.030_r8
             write(iulog,*)'Total thickness of roof: ',sum(dzurb_roof(l,:))
             write(iulog,*)'Roof layer thicknesses: ',dzurb_roof(l,:)

             ziurb_wall(l,0) = 0.
             ziurb_wall(l,1) = dzurb_wall(l,1)
             do j = 2,nlevurb
                ziurb_wall(l,j) = sum(dzurb_wall(l,1:j))
             end do
             write(iulog,*)'Wall layer interface depths: ',ziurb_wall(l,:)

             ziurb_roof(l,0) = 0.
             ziurb_roof(l,1) = dzurb_roof(l,1)
             do j = 2,nlevurb
                ziurb_roof(l,j) = sum(dzurb_roof(l,1:j))
             end do
             write(iulog,*)'Roof layer interface depths: ',ziurb_roof(l,:)
          else if (use_mexicocity) then
             zurb_wall(l,1) = 0.015_r8/2._r8
             zurb_wall(l,2) = zurb_wall(l,1) + 0.015_r8/2._r8 + 0.120_r8/2._r8
             zurb_wall(l,3) = zurb_wall(l,2) + 0.120_r8/2._r8 + 0.150_r8/2._r8
             zurb_wall(l,4) = zurb_wall(l,3) + 0.150_r8/2._r8 + 0.150_r8/2._r8
             zurb_wall(l,5) = zurb_wall(l,4) + 0.150_r8/2._r8 + 0.015_r8/2._r8

             zurb_roof(l,1) = 0.010_r8/2._r8
             zurb_roof(l,2) = zurb_roof(l,1) + 0.010_r8/2._r8 + 0.050_r8/2._r8
             zurb_roof(l,3) = zurb_roof(l,2) + 0.050_r8/2._r8 + 0.050_r8/2._r8
             zurb_roof(l,4) = zurb_roof(l,3) + 0.050_r8/2._r8 + 0.050_r8/2._r8
             zurb_roof(l,5) = zurb_roof(l,4) + 0.050_r8/2._r8 + 0.025_r8/2._r8

             dzurb_wall(l,1) = 0.015_r8
             dzurb_wall(l,2) = 0.120_r8
             dzurb_wall(l,3) = 0.150_r8
             dzurb_wall(l,4) = 0.150_r8
             dzurb_wall(l,5) = 0.015_r8
             write(iulog,*)'Total thickness of wall: ',sum(dzurb_wall(l,:))
             write(iulog,*)'Wall layer thicknesses: ',dzurb_wall(l,:)

             dzurb_roof(l,1) = 0.010_r8
             dzurb_roof(l,2) = 0.050_r8
             dzurb_roof(l,3) = 0.050_r8
             dzurb_roof(l,4) = 0.050_r8
             dzurb_roof(l,5) = 0.025_r8
             write(iulog,*)'Total thickness of roof: ',sum(dzurb_roof(l,:))
             write(iulog,*)'Roof layer thicknesses: ',dzurb_roof(l,:)

             ziurb_wall(l,0) = 0.
             ziurb_wall(l,1) = dzurb_wall(l,1)
             do j = 2,nlevurb
                ziurb_wall(l,j) = sum(dzurb_wall(l,1:j))
             end do
             write(iulog,*)'Wall layer interface depths: ',ziurb_wall(l,:)

             ziurb_roof(l,0) = 0.
             ziurb_roof(l,1) = dzurb_roof(l,1)
             do j = 2,nlevurb
                ziurb_roof(l,j) = sum(dzurb_roof(l,1:j))
             end do
             write(iulog,*)'Roof layer interface depths: ',ziurb_roof(l,:)
          else
             do j = 1, nlevurb
                zurb_wall(l,j) = (j-0.5)*(thick_wall(l)/float(nlevurb))  !node depths
             end do
             do j = 1, nlevurb
                zurb_roof(l,j) = (j-0.5)*(thick_roof(l)/float(nlevurb))  !node depths
             end do

             dzurb_roof(l,1) = 0.5*(zurb_roof(l,1)+zurb_roof(l,2))    !thickness b/n two interfaces
             do j = 2,nlevurb-1
                dzurb_roof(l,j)= 0.5*(zurb_roof(l,j+1)-zurb_roof(l,j-1)) 
             enddo
             dzurb_roof(l,nlevurb) = zurb_roof(l,nlevurb)-zurb_roof(l,nlevurb-1)

             dzurb_wall(l,1) = 0.5*(zurb_wall(l,1)+zurb_wall(l,2))    !thickness b/n two interfaces
             do j = 2,nlevurb-1
                dzurb_wall(l,j)= 0.5*(zurb_wall(l,j+1)-zurb_wall(l,j-1)) 
             enddo
             dzurb_wall(l,nlevurb) = zurb_wall(l,nlevurb)-zurb_wall(l,nlevurb-1)

             ziurb_wall(l,0) = 0.
             do j = 1, nlevurb-1
                ziurb_wall(l,j) = 0.5*(zurb_wall(l,j)+zurb_wall(l,j+1))          !interface depths
             enddo
             ziurb_wall(l,nlevurb) = zurb_wall(l,nlevurb) + 0.5*dzurb_wall(l,nlevurb)

             ziurb_roof(l,0) = 0.
             do j = 1, nlevurb-1
                ziurb_roof(l,j) = 0.5*(zurb_roof(l,j)+zurb_roof(l,j+1))          !interface depths
             enddo
             ziurb_roof(l,nlevurb) = zurb_roof(l,nlevurb) + 0.5*dzurb_roof(l,nlevurb)
          end if
       end if
    end do

    do c = bounds%begc,bounds%endc
       l = col%landunit(c)

       if (lun%urbpoi(l)) then
          if (col%itype(c)==icol_sunwall .or. col%itype(c)==icol_shadewall) then
             col%z(c,1:nlevurb)  = zurb_wall(l,1:nlevurb)
             col%zi(c,0:nlevurb) = ziurb_wall(l,0:nlevurb)
             col%dz(c,1:nlevurb) = dzurb_wall(l,1:nlevurb)
             if (nlevurb < nlevgrnd) then
                col%z(c,nlevurb+1:nlevgrnd)  = spval
                col%zi(c,nlevurb+1:nlevgrnd) = spval
                col%dz(c,nlevurb+1:nlevgrnd) = spval
             end if
          else if (col%itype(c)==icol_roof) then
             col%z(c,1:nlevurb)  = zurb_roof(l,1:nlevurb)
             col%zi(c,0:nlevurb) = ziurb_roof(l,0:nlevurb)
             col%dz(c,1:nlevurb) = dzurb_roof(l,1:nlevurb)
             if (nlevurb < nlevgrnd) then
                col%z(c,nlevurb+1:nlevgrnd)  = spval
                col%zi(c,nlevurb+1:nlevgrnd) = spval
                col%dz(c,nlevurb+1:nlevgrnd) = spval
             end if
          else
             col%z(c,1:nlevgrnd)  = zsoi(1:nlevgrnd)
             col%zi(c,0:nlevgrnd) = zisoi(0:nlevgrnd)
             col%dz(c,1:nlevgrnd) = dzsoi(1:nlevgrnd)
          end if
       else if (lun%itype(l) /= istdlak) then
          col%z(c,1:nlevgrnd)  = zsoi(1:nlevgrnd)
          col%zi(c,0:nlevgrnd) = zisoi(0:nlevgrnd)
          col%dz(c,1:nlevgrnd) = dzsoi(1:nlevgrnd)
       end if
    end do

    if (nlevurb > 0) then
       deallocate(zurb_wall, zurb_roof, dzurb_wall, dzurb_roof, ziurb_wall, ziurb_roof)
    end if

    !-----------------------------------------------
    ! Set index defining depth to bedrock
    !-----------------------------------------------

    allocate(zbedrock_in(bounds%begg:bounds%endg))
    if (use_bedrock) then
       call ncd_io(ncid=ncid, varname='zbedrock', flag='read', data=zbedrock_in, dim1name=grlnd, readvar=readvar)
       if (.not. readvar) then
          if (masterproc) then
             call endrun( 'ERROR:: zbedrock not found on surface data set, and use_bedrock is true.'//errmsg(sourcefile, __LINE__) )
          end if
       end if

    !  if use_bedrock = false, set zbedrock to lowest layer bottom interface
    else
       if (masterproc) write(iulog,*) 'not using use_bedrock!!'
       zbedrock_in(:) = zisoi(nlevsoi)
    endif

    !  determine minimum index of minimum soil depth
    jmin_bedrock = 3
    do j = 3,nlevsoi 
       if (zisoi(j-1) < zmin_bedrock .and. zisoi(j) >= zmin_bedrock) then
          jmin_bedrock = j
       endif
    enddo

    if (masterproc) write(iulog,*) 'jmin_bedrock: ', jmin_bedrock

    !  Determine gridcell bedrock index
    do g = bounds%begg,bounds%endg
       grc%nbedrock(g) = nlevsoi
       do j = jmin_bedrock,nlevsoi 
          if (zisoi(j-1) < zbedrock_in(g) .and. zisoi(j) >= zbedrock_in(g)) then
             grc%nbedrock(g) = j
          end if
       end do
    end do

    !  Set column bedrock index
    do c = begc, endc
       g = col%gridcell(c)
       col%nbedrock(c) = grc%nbedrock(g) 
    end do

    deallocate(zbedrock_in)

    !-----------------------------------------------
    ! Set lake levels and layers (no interfaces)
    !-----------------------------------------------

    allocate(lakedepth_in(bounds%begg:bounds%endg))
    call ncd_io(ncid=ncid, varname='LAKEDEPTH', flag='read', data=lakedepth_in, dim1name=grlnd, readvar=readvar)
    if (.not. readvar) then
       if (masterproc) then
          write(iulog,*) 'WARNING:: LAKEDEPTH not found on surface data set. All lake columns will have lake depth', &
               ' set equal to default value.'
       end if
       lakedepth_in(:) = spval
    end if
    do c = begc, endc
       g = col%gridcell(c)
       col%lakedepth(c) = lakedepth_in(g)
    end do
    deallocate(lakedepth_in)

    ! Lake layers
    if (.not. use_extralakelayers) then
       dzlak(1) = 0.1_r8
       dzlak(2) = 1._r8
       dzlak(3) = 2._r8
       dzlak(4) = 3._r8
       dzlak(5) = 4._r8
       dzlak(6) = 5._r8
       dzlak(7) = 7._r8
       dzlak(8) = 7._r8
       dzlak(9) = 10.45_r8
       dzlak(10)= 10.45_r8

       zlak(1) =  0.05_r8
       zlak(2) =  0.6_r8
       zlak(3) =  2.1_r8
       zlak(4) =  4.6_r8
       zlak(5) =  8.1_r8
       zlak(6) = 12.6_r8
       zlak(7) = 18.6_r8
       zlak(8) = 25.6_r8
       zlak(9) = 34.325_r8
       zlak(10)= 44.775_r8
    else
       dzlak(1) =0.1_r8
       dzlak(2) =0.25_r8
       dzlak(3) =0.25_r8
       dzlak(4) =0.25_r8
       dzlak(5) =0.25_r8
       dzlak(6) =0.5_r8
       dzlak(7) =0.5_r8
       dzlak(8) =0.5_r8
       dzlak(9) =0.5_r8
       dzlak(10) =0.75_r8
       dzlak(11) =0.75_r8
       dzlak(12) =0.75_r8
       dzlak(13) =0.75_r8
       dzlak(14) =2_r8
       dzlak(15) =2_r8
       dzlak(16) =2.5_r8
       dzlak(17) =2.5_r8
       dzlak(18) =3.5_r8
       dzlak(19) =3.5_r8
       dzlak(20) =3.5_r8
       dzlak(21) =3.5_r8
       dzlak(22) =5.225_r8
       dzlak(23) =5.225_r8
       dzlak(24) =5.225_r8
       dzlak(25) =5.225_r8

       zlak(1) = dzlak(1)/2._r8
       do i=2,nlevlak
          zlak(i) = zlak(i-1) + (dzlak(i-1)+dzlak(i))/2._r8
       end do
    end if

    do c = bounds%begc,bounds%endc
       l = col%landunit(c)

       if (lun%itype(l) == istdlak) then

          if (col%lakedepth(c) == spval) then
             col%lakedepth(c)         = zlak(nlevlak) + 0.5_r8*dzlak(nlevlak)
             col%z_lake(c,1:nlevlak)  = zlak(1:nlevlak)
             col%dz_lake(c,1:nlevlak) = dzlak(1:nlevlak)

          else if (col%lakedepth(c) > 1._r8 .and. col%lakedepth(c) < 5000._r8) then

             depthratio                 = col%lakedepth(c) / (zlak(nlevlak) + 0.5_r8*dzlak(nlevlak)) 
             col%z_lake(c,1)            = zlak(1)
             col%dz_lake(c,1)           = dzlak(1)
             col%dz_lake(c,2:nlevlak-1) = dzlak(2:nlevlak-1)*depthratio
             col%dz_lake(c,nlevlak)     = dzlak(nlevlak)*depthratio - (col%dz_lake(c,1) - dzlak(1)*depthratio)
             do lev=2,nlevlak
                col%z_lake(c,lev) = col%z_lake(c,lev-1) + (col%dz_lake(c,lev-1)+col%dz_lake(c,lev))/2._r8
             end do

          else if (col%lakedepth(c) > 0._r8 .and. col%lakedepth(c) <= 1._r8) then

             col%dz_lake(c,:) = col%lakedepth(c) / nlevlak;
             col%z_lake(c,1)  = col%dz_lake(c,1) / 2._r8;
             do lev=2,nlevlak
                col%z_lake(c,lev) = col%z_lake(c,lev-1) + (col%dz_lake(c,lev-1)+col%dz_lake(c,lev))/2._r8
             end do

          else

             write(iulog,*)'Bad lake depth: lakedepth: ', col%lakedepth(c)
             call shr_sys_abort(errmsg(sourcefile, __LINE__))

          end if

          col%z(c,1:nlevgrnd)  = zsoi(1:nlevgrnd)
          col%zi(c,0:nlevgrnd) = zisoi(0:nlevgrnd)
          col%dz(c,1:nlevgrnd) = dzsoi(1:nlevgrnd)
       end if
    end do

    ! ------------------------------------------------------------------------
    ! Set classes of layers
    ! ------------------------------------------------------------------------

    do c = bounds%begc, bounds%endc
       l = col%landunit(c)
       if (hasBedrock(col_itype=col%itype(c), lun_itype=lun%itype(l))) then
          ! NOTE(wjs, 2015-10-17) We are assuming that points with bedrock have both
          ! "shallow" and "deep" bedrock. Currently, this is not true for lake columns:
          ! lakes do not distinguish between "shallow" bedrock and "normal" soil.
          ! However, that was just due to an oversight that is supposed to be corrected
          ! soon; so to keep things simple we assume that any point with bedrock
          ! potentially has both shallow and deep bedrock.
          col%levgrnd_class(c, 1:col%nbedrock(c)) = LEVGRND_CLASS_STANDARD
          if (col%nbedrock(c) < nlevsoi) then
             col%levgrnd_class(c, (col%nbedrock(c) + 1) : nlevsoi) = LEVGRND_CLASS_SHALLOW_BEDROCK
          end if
          col%levgrnd_class(c, (nlevsoi + 1) : nlevgrnd) = LEVGRND_CLASS_DEEP_BEDROCK
       else
          col%levgrnd_class(c, 1:nlevgrnd) = LEVGRND_CLASS_STANDARD
       end if
    end do

    do j = 1, nlevgrnd
       do c = bounds%begc, bounds%endc
          if (col%z(c,j) == spval) then
             col%levgrnd_class(c,j) = ispval
          end if
       end do
    end do

    !-----------------------------------------------
    ! Set cold-start values for snow levels, snow layers and snow interfaces 
    !-----------------------------------------------

    call InitSnowLayers(bounds, snow_depth(bounds%begc:bounds%endc))

    !-----------------------------------------------
    ! Read in topographic index and slope
    !-----------------------------------------------

    allocate(tslope(bounds%begg:bounds%endg))
    call ncd_io(ncid=ncid, varname='SLOPE', flag='read', data=tslope, dim1name=grlnd, readvar=readvar)
    if (.not. readvar) then
       call shr_sys_abort(' ERROR: TOPOGRAPHIC SLOPE NOT on surfdata file'//&
            errMsg(sourcefile, __LINE__)) 
    end if
    do c = begc,endc
       g = col%gridcell(c)
       ! check for near zero slopes, set minimum value
       col%topo_slope(c) = max(tslope(g), 0.2_r8)
    end do
    deallocate(tslope)

    allocate(std(bounds%begg:bounds%endg))
    call ncd_io(ncid=ncid, varname='STD_ELEV', flag='read', data=std, dim1name=grlnd, readvar=readvar)
    if (.not. readvar) then
       call shr_sys_abort(' ERROR: TOPOGRAPHIC STDdev (STD_ELEV) NOT on surfdata file'//&
            errMsg(sourcefile, __LINE__)) 
    end if
    do c = begc,endc
       g = col%gridcell(c)
       ! Topographic variables
       col%topo_std(c) = std(g)
    end do
    deallocate(std)

    !-----------------------------------------------
    ! SCA shape function defined
    !-----------------------------------------------

    do c = begc,endc
       l = col%landunit(c)
       g = col%gridcell(c)

       if (lun%itype(l)==istice_mec .and. glc_behavior%allow_multiple_columns_grc(g)) then
          ! ice_mec columns already account for subgrid topographic variability through
          ! their use of multiple elevation classes; thus, to avoid double-accounting for
          ! topographic variability in these columns, we ignore topo_std and use a value
          ! of n_melt that assumes little topographic variability within the column
          col%n_melt(c) = 10._r8
       else
          col%n_melt(c) = 200.0/max(10.0_r8, col%topo_std(c))
       end if

       ! microtopographic parameter, units are meters (try smooth function of slope)

       slopebeta = 3._r8
       slopemax = 0.4_r8
       slope0 = slopemax**(-1._r8/slopebeta)
       col%micro_sigma(c) = (col%topo_slope(c) + slope0)**(-slopebeta)
    end do

    call ncd_pio_closefile(ncid)

  end subroutine initVertical

  !-----------------------------------------------------------------------
  logical function hasBedrock(col_itype, lun_itype)
    !
    ! !DESCRIPTION:
    ! Returns true if the given column type has a representation of bedrock - i.e., a set
    ! of layers at the bottom of the column that are treated fundamentally differently
    ! from the upper layers.
    !
    ! !USES:
    use landunit_varcon, only : istice, istice_mec, isturb_MIN, isturb_MAX
    use column_varcon  , only : icol_road_perv, is_hydrologically_active
    !
    ! !ARGUMENTS:
    integer, intent(in) :: col_itype  ! col%itype value
    integer, intent(in) :: lun_itype  ! lun%itype value for the landunit on which this column sits
    ! If we had an easy way to figure out which landunit a column was on based on
    ! col_itype (which would be very helpful!), then we wouldn't need lun_itype.
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'hasBedrock'
    !-----------------------------------------------------------------------

    ! TODO(wjs, 2015-10-17) I don't like that the logic here implicitly duplicates logic
    ! elsewhere in the code. For example, if there were a change in the lake code so that
    ! it no longer treated the bottom layers as bedrock, then that change would need to be
    ! reflected here. One solution would be to set some has_bedrock flag in one central
    ! place, and then have the science code use that. But that could get messy in the
    ! science code. Another solution would be to decentralize the definition of
    ! hasBedrock, so that (for example) the lake code itself sets the value for lun_itype
    ! == istdlak - that way, hasBedrock(lake) would be more likely to get updated
    ! correctly if the lake logic changes.

    if (lun_itype == istice .or. lun_itype == istice_mec) then
       hasBedrock = .false.
    else if (lun_itype >= isturb_MIN .and. lun_itype <= isturb_MAX) then
       if (col_itype == icol_road_perv) then
          hasBedrock = .true.
       else
          hasBedrock = .false.
       end if
    else
       hasBedrock = .true.
    end if

    ! As an independent check of the above logic, assert that, at the very least, any
    ! hydrologically-active column is given hasBedrock = .true. This is to try to catch
    ! problems with new column types being added that aren't handled properly by the
    ! above logic, since (as noted in the todo note above) there is some implicit
    ! duplication of logic between this routine and other parts of the code, which is
    ! dangerous. For example, if a new "urban lawn" type is added, then it should have
    ! hasBedrock = .true. - and this omission will hopefully be caught by this assertion.
    if (is_hydrologically_active(col_itype=col_itype, lun_itype=lun_itype)) then
       SHR_ASSERT(hasBedrock, "hasBedrock should be true for all hydrologically-active columns")
    end if

  end function hasBedrock


end module initVerticalMod
