#include <misc.h>
#include <preproc.h>


module mkarbinitMod 1

contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: mkarbinit
!
! !INTERFACE:

  subroutine mkarbinit() 1,13
!
! !DESCRIPTION:
! Initializes the following time varying variables:
! water      : h2osno, h2ocan, h2osoi_liq, h2osoi_ice, h2osoi_vol
! snow       : snowdp, snl, dz, z, zi
! temperature: t_soisno, t_veg, t_grnd
!
! !USES:
    use shr_kind_mod , only : r8 => shr_kind_r8
    use shr_const_mod, only : SHR_CONST_TKFRZ
    use clmtype
    use clm_varpar   , only : nlevsoi, nlevgrnd, nlevsno, nlevlak, nlevurb
    use clm_varcon   , only : bdsno, istice, istwet, istsoil, isturb, &
                              denice, denh2o, spval, sb, icol_road_perv, &
                              icol_road_imperv, icol_roof, icol_sunwall, &
                              icol_shadewall
    use clm_varcon   , only : istice_mec, h2osno_max
    use clm_varctl   , only : iulog, pertlim
    use spmdMod      , only : masterproc
    use decompMod    , only : get_proc_bounds
    use shr_sys_mod  , only : shr_sys_flush
    use SNICARMod    , only : snw_rds_min

!
! !ARGUMENTS:
    implicit none
!
! !CALLED FROM:
! subroutine initialize in module initializeMod
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
! 3/07/08 Keith Oleson: initialize h2osoi_vol for all soil layers to 0.3
! 3/18/08 David Lawrence, initialize deep layers
! 03/28/08 Mark Flanner, initialize snow aerosols and grain size
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
    integer , pointer :: pcolumn(:)        ! column index associated with each pft
    integer , pointer :: ctype(:)          ! column type
    integer , pointer :: clandunit(:)      ! landunit index associated with each column
    integer , pointer :: ltype(:)          ! landunit type
    logical , pointer :: lakpoi(:)         ! true => landunit is a lake point
    integer , pointer :: plandunit(:)      ! landunit index associated with each pft
    logical , pointer :: urbpoi(:)         ! true => landunit is an urban point
    logical , pointer :: ifspecial(:)      ! true => landunit is not vegetated
    real(r8), pointer :: dz(:,:)           ! layer thickness depth (m)
    real(r8), pointer :: watsat(:,:)       ! volumetric soil water at saturation (porosity) (nlevgrnd)
    real(r8), pointer :: h2osoi_ice(:,:)   ! ice lens (kg/m2)
    real(r8), pointer :: h2osoi_liq(:,:)   ! liquid water (kg/m2)
    real(r8), pointer :: bsw2(:,:)         ! Clapp and Hornberger "b" for CN code
    real(r8), pointer :: psisat(:,:)       ! soil water potential at saturation for CN code (MPa)
    real(r8), pointer :: vwcsat(:,:)       ! volumetric water content at saturation for CN code (m3/m3)
    real(r8), pointer :: zi(:,:)           ! interface level below a "z" level (m)
    real(r8), pointer :: wa(:)             ! water in the unconfined aquifer (mm)
    real(r8), pointer :: wt(:)             ! total water storage (unsaturated soil water + groundwater) (mm)
    real(r8), pointer :: zwt(:)            ! water table depth (m)
!
! local pointers to implicit out arguments
!
    integer , pointer :: snl(:)             ! number of snow layers
    real(r8), pointer :: t_soisno(:,:)      ! soil temperature (Kelvin)  (-nlevsno+1:nlevgrnd)
    real(r8), pointer :: t_lake(:,:)        ! lake temperature (Kelvin)  (1:nlevlak)
    real(r8), pointer :: t_grnd(:)          ! ground temperature (Kelvin)
    real(r8), pointer :: t_veg(:)           ! vegetation temperature (Kelvin)
    real(r8), pointer :: t_ref2m(:)         ! 2 m height surface air temperature (Kelvin)
    real(r8), pointer :: t_ref2m_u(:)       ! Urban 2 m height surface air temperature (Kelvin)
    real(r8), pointer :: t_ref2m_r(:)       ! Rural 2 m height surface air temperature (Kelvin)
    real(r8), pointer :: h2osoi_vol(:,:)    ! volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3]
    real(r8), pointer :: h2ocan_col(:)      ! canopy water (mm H2O) (column-level)
    real(r8), pointer :: h2ocan_pft(:)      ! canopy water (mm H2O) (pft-level)
    real(r8), pointer :: h2osno(:)          ! snow water (mm H2O)
    real(r8), pointer :: snowdp(:)          ! snow height (m)
    real(r8), pointer :: eflx_lwrad_out(:)  ! emitted infrared (longwave) radiation (W/m**2)
    real(r8), pointer :: soilpsi(:,:)       ! soil water potential in each soil layer (MPa)
    real(r8), pointer :: snw_rds(:,:)       ! effective snow grain radius (col,lyr) [microns, m^-6]
    real(r8), pointer :: snw_rds_top(:)     ! snow grain size, top (col) [microns]
    real(r8), pointer :: sno_liq_top(:)     ! liquid water fraction (mass) in top snow layer (col) [frc]
    real(r8), pointer :: mss_bcpho(:,:)     ! mass of hydrophobic BC in snow (col,lyr) [kg]
    real(r8), pointer :: mss_bcphi(:,:)     ! mass of hydrophillic BC in snow (col,lyr) [kg]
    real(r8), pointer :: mss_bctot(:,:)     ! total mass of BC (pho+phi) (col,lyr) [kg]
    real(r8), pointer :: mss_bc_col(:)      ! total mass of BC in snow column (col) [kg]
    real(r8), pointer :: mss_bc_top(:)      ! total mass of BC in top snow layer (col) [kg]
    real(r8), pointer :: mss_cnc_bcphi(:,:) ! mass concentration of BC species 1 (col,lyr) [kg/kg]
    real(r8), pointer :: mss_cnc_bcpho(:,:) ! mass concentration of BC species 2 (col,lyr) [kg/kg]
    real(r8), pointer :: mss_ocpho(:,:)     ! mass of hydrophobic OC in snow (col,lyr) [kg]
    real(r8), pointer :: mss_ocphi(:,:)     ! mass of hydrophillic OC in snow (col,lyr) [kg]
    real(r8), pointer :: mss_octot(:,:)     ! total mass of OC (pho+phi) (col,lyr) [kg]
    real(r8), pointer :: mss_oc_col(:)      ! total mass of OC in snow column (col) [kg]
    real(r8), pointer :: mss_oc_top(:)      ! total mass of OC in top snow layer (col) [kg]
    real(r8), pointer :: mss_cnc_ocphi(:,:) ! mass concentration of OC species 1 (col,lyr) [kg/kg]
    real(r8), pointer :: mss_cnc_ocpho(:,:) ! mass concentration of OC species 2 (col,lyr) [kg/kg]
    real(r8), pointer :: mss_dst1(:,:)      ! mass of dust species 1 in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst2(:,:)      ! mass of dust species 2 in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst3(:,:)      ! mass of dust species 3 in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst4(:,:)      ! mass of dust species 4 in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dsttot(:,:)    ! total mass of dust in snow (col,lyr) [kg]
    real(r8), pointer :: mss_dst_col(:)     ! total mass of dust in snow column (col) [kg]
    real(r8), pointer :: mss_dst_top(:)     ! total mass of dust in top snow layer (col) [kg]
    real(r8), pointer :: mss_cnc_dst1(:,:)  ! mass concentration of dust species 1 (col,lyr) [kg/kg]
    real(r8), pointer :: mss_cnc_dst2(:,:)  ! mass concentration of dust species 2 (col,lyr) [kg/kg]
    real(r8), pointer :: mss_cnc_dst3(:,:)  ! mass concentration of dust species 3 (col,lyr) [kg/kg]
    real(r8), pointer :: mss_cnc_dst4(:,:)  ! mass concentration of dust species 4 (col,lyr) [kg/kg]
!
!
! !OTHER LOCAL VARIABLES:
!EOP
    integer :: j,l,c,p      ! indices
    integer :: nlevs        ! number of levels
    integer :: begp, endp   ! per-proc beginning and ending pft indices
    integer :: begc, endc   ! per-proc beginning and ending column indices
    integer :: begl, endl   ! per-proc beginning and ending landunit indices
    integer :: begg, endg   ! per-proc gridcell ending gridcell indices
    real(r8):: vwc,psi      ! for calculating soilpsi
    real(r8):: pertval      ! for calculating temperature perturbation
!-----------------------------------------------------------------------

    if ( masterproc )then
        write(iulog,*) 'Setting initial data to non-spun up values'
        if ( pertlim /= 0.0_r8 ) &
        write(iulog,*) 'Applying perturbation to initial surface temperature'
    end if

    ! Assign local pointers to derived subtypes components (landunit-level)

    ltype      => clm3%g%l%itype
    lakpoi     => clm3%g%l%lakpoi
    ifspecial  => clm3%g%l%ifspecial
    urbpoi     => clm3%g%l%urbpoi

    ! Assign local pointers to derived subtypes components (column-level)

    ctype            => clm3%g%l%c%itype
    clandunit        => clm3%g%l%c%landunit
    snl              => clm3%g%l%c%cps%snl
    dz               => clm3%g%l%c%cps%dz
    watsat           => clm3%g%l%c%cps%watsat
    bsw2             => clm3%g%l%c%cps%bsw2
    vwcsat           => clm3%g%l%c%cps%vwcsat
    psisat           => clm3%g%l%c%cps%psisat
    soilpsi          => clm3%g%l%c%cps%soilpsi
    h2osoi_ice       => clm3%g%l%c%cws%h2osoi_ice
    h2osoi_liq       => clm3%g%l%c%cws%h2osoi_liq
    h2osoi_vol       => clm3%g%l%c%cws%h2osoi_vol
    h2ocan_col       => clm3%g%l%c%cws%pws_a%h2ocan
    snowdp           => clm3%g%l%c%cps%snowdp
    h2osno           => clm3%g%l%c%cws%h2osno
    t_soisno         => clm3%g%l%c%ces%t_soisno
    t_lake           => clm3%g%l%c%ces%t_lake
    t_grnd           => clm3%g%l%c%ces%t_grnd
    zi               => clm3%g%l%c%cps%zi
    wa               => clm3%g%l%c%cws%wa
    wt               => clm3%g%l%c%cws%wt
    zwt              => clm3%g%l%c%cws%zwt
    snw_rds          => clm3%g%l%c%cps%snw_rds
    snw_rds_top      => clm3%g%l%c%cps%snw_rds_top
    sno_liq_top      => clm3%g%l%c%cps%sno_liq_top
    mss_bcpho        => clm3%g%l%c%cps%mss_bcpho
    mss_bcphi        => clm3%g%l%c%cps%mss_bcphi
    mss_bctot        => clm3%g%l%c%cps%mss_bctot
    mss_bc_col       => clm3%g%l%c%cps%mss_bc_col
    mss_bc_top       => clm3%g%l%c%cps%mss_bc_top
    mss_cnc_bcphi    => clm3%g%l%c%cps%mss_cnc_bcphi
    mss_cnc_bcpho    => clm3%g%l%c%cps%mss_cnc_bcpho
    mss_ocpho        => clm3%g%l%c%cps%mss_ocpho
    mss_ocphi        => clm3%g%l%c%cps%mss_ocphi
    mss_octot        => clm3%g%l%c%cps%mss_octot
    mss_oc_col       => clm3%g%l%c%cps%mss_oc_col
    mss_oc_top       => clm3%g%l%c%cps%mss_oc_top
    mss_cnc_ocphi    => clm3%g%l%c%cps%mss_cnc_ocphi
    mss_cnc_ocpho    => clm3%g%l%c%cps%mss_cnc_ocpho
    mss_dst1         => clm3%g%l%c%cps%mss_dst1
    mss_dst2         => clm3%g%l%c%cps%mss_dst2
    mss_dst3         => clm3%g%l%c%cps%mss_dst3
    mss_dst4         => clm3%g%l%c%cps%mss_dst4
    mss_dsttot       => clm3%g%l%c%cps%mss_dsttot
    mss_dst_col      => clm3%g%l%c%cps%mss_dst_col
    mss_dst_top      => clm3%g%l%c%cps%mss_dst_top
    mss_cnc_dst1     => clm3%g%l%c%cps%mss_cnc_dst1
    mss_cnc_dst2     => clm3%g%l%c%cps%mss_cnc_dst2
    mss_cnc_dst3     => clm3%g%l%c%cps%mss_cnc_dst3
    mss_cnc_dst4     => clm3%g%l%c%cps%mss_cnc_dst4

    ! Assign local pointers to derived subtypes components (pft-level)

    pcolumn        => clm3%g%l%c%p%column
    h2ocan_pft     => clm3%g%l%c%p%pws%h2ocan
    t_veg          => clm3%g%l%c%p%pes%t_veg
    t_ref2m        => clm3%g%l%c%p%pes%t_ref2m
    t_ref2m_u      => clm3%g%l%c%p%pes%t_ref2m_u
    t_ref2m_r      => clm3%g%l%c%p%pes%t_ref2m_r
    plandunit      => clm3%g%l%c%p%landunit
    eflx_lwrad_out => clm3%g%l%c%p%pef%eflx_lwrad_out  

    ! Determine subgrid bounds on this processor

    call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp)

    ! NOTE: h2ocan, h2osno, and snowdp has valid values everywhere
    ! canopy water (pft level)

    do p = begp, endp
       h2ocan_pft(p) = 0._r8
       
       ! added for canopy water mass balance under dynamic pft weights
       !clm3%g%l%c%p%pps%tlai(p) = 0._r8
       !clm3%g%l%c%p%pps%tsai(p) = 0._r8
       !clm3%g%l%c%p%pps%elai(p) = 0._r8
       !clm3%g%l%c%p%pps%esai(p) = 0._r8
       !clm3%g%l%c%p%pps%htop(p) = 0._r8
       !clm3%g%l%c%p%pps%hbot(p) = 0._r8
       !clm3%g%l%c%p%pps%frac_veg_nosno_alb(p) = 0._r8
    end do

    do c = begc,endc

       ! canopy water (column level)

       h2ocan_col(c) = 0._r8

       ! snow water

       l = clandunit(c)

       ! Note: Glacier_mec columns are initialized with half the maximum snow cover.
       ! This gives more realistic values of qflx_glcice sooner in the simulation
       !  for columns with net ablation, at the cost of delaying ice formation
       !  in columns with net accumulation.
       if (ltype(l)==istice) then
          h2osno(c) = h2osno_max
       elseif (ltype(l)==istice_mec) then
          h2osno(c) = 0.5_r8 * h2osno_max   ! 50 cm if h2osno_max = 1 m
       else
          h2osno(c) = 0._r8
       endif

       ! snow depth

       snowdp(c)  = h2osno(c) / bdsno

    end do

    ! Set snow layer number, depth and thickiness

    call snowdp2lev(begc, endc)

    ! Set snow/soil temperature, note:
    ! t_soisno only has valid values over non-lake
    ! t_lake   only has valid values over lake
    ! t_grnd has valid values over all land
    ! t_veg  has valid values over all land

    ! NOTE: THESE MEMORY COPIES ARE INEFFICIENT -- SINCE nlev LOOP IS NESTED FIRST!!!!
    do c = begc,endc

       t_soisno(c,-nlevsno+1:nlevgrnd) = spval
       t_lake(c,1:nlevlak) = spval

       l = clandunit(c)
       if (.not. lakpoi(l)) then  !not lake
          t_soisno(c,-nlevsno+1:0) = spval
          if (snl(c) < 0) then    !snow layer temperatures
             do j = snl(c)+1, 0
                t_soisno(c,j) = 250._r8
             enddo
          endif
          if (ltype(l)==istice .or. ltype(l)==istice_mec) then
             do j = 1, nlevgrnd
                t_soisno(c,j) = 250._r8

             end do
          else if (ltype(l) == istwet) then
             do j = 1, nlevgrnd
                t_soisno(c,j) = 277._r8
             end do
          else if (ltype(l) == isturb) then
#if (defined VANCOUVER)
             if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then 
               ! Set road top layer to initial air temperature and interpolate other
               ! layers down to 20C in bottom layer
               do j = 1, nlevurb
                  t_soisno(c,j) = 297.56 - (j-1) * ((297.56-293.16)/(nlevurb-1)) 
               end do
             ! Set wall and roof layers to initial air temperature
             else if (ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall .or. ctype(c) == icol_roof) then
               do j = 1, nlevurb
                  t_soisno(c,j) = 297.56
               end do
             else
               do j = 1, nlevurb
                  t_soisno(c,j) = 283._r8
               end do
             end if
#elif (defined MEXICOCITY)
             if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then 
               ! Set road top layer to initial air temperature and interpolate other
               ! layers down to 22C in bottom layer
               do j = 1, nlevurb
                  t_soisno(c,j) = 289.46 - (j-1) * ((289.46-295.16)/(nlevurb-1)) 
               end do
             ! Set wall and roof layers to initial air temperature
             else if (ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall .or. ctype(c) == icol_roof) then
               do j = 1, nlevurb
                  t_soisno(c,j) = 289.46
               end do
             else
               do j = 1, nlevurb
                  t_soisno(c,j) = 283._r8
               end do
             end if
#elif (defined GRANDVIEW)
             if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then 
             ! Set road layers to 18.5C
               do j = 1, nlevurb
                  t_soisno(c,j) = 291.66_r8
               end do
             else if (ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall) then
             ! Set wall layers to 18.35C
               do j = 1, nlevurb
                  t_soisno(c,j) = 291.51_r8
               end do
             else if (ctype(c) == icol_roof) then
             ! Set roof layers to 18.5C
               do j = 1, nlevurb
                  t_soisno(c,j) = 291.66_r8
               end do
             end if
#else
             if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then 
               do j = 1, nlevurb
                 t_soisno(c,j) = 274._r8
               end do
             ! Set sunwall, shadewall, roof to fairly high temperature to avoid initialization
             ! shock from large heating/air conditioning flux
             else if (ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall &
                      .or. ctype(c) == icol_roof) then
               do j = 1, nlevurb
                 t_soisno(c,j) = 292._r8
               end do
             end if
#endif
          else
             do j = 1, nlevgrnd
                t_soisno(c,j) = 274._r8
             end do
          endif
          t_grnd(c) = t_soisno(c,snl(c)+1)
       else                     !lake
          t_lake(c,1:nlevlak) = 277._r8
          t_grnd(c) = t_lake(c,1)
       endif
       if ( pertlim /= 0.0_r8 )then
          if (.not. lakpoi(l)) then  !not lake
             if (     ltype(l) == isturb) then
                nlevs = nlevurb
             else
                nlevs = nlevgrnd
             end if
             
             do j = 1, nlevs
                call random_number (pertval)
                pertval       = 2._r8*pertlim*(0.5_r8 - pertval)
                t_soisno(c,j) = t_soisno(c,j)*(1._r8 + pertval)
             end do
             t_grnd(c) = t_soisno(c,snl(c)+1)
          else                       !lake
             do j = 1, nlevlak
                call random_number (pertval)
                pertval     = 2._r8*pertlim*(0.5_r8 - pertval)
                t_lake(c,j) = t_lake(c,j)*(1._r8 + pertval)
             end do
             t_grnd(c) = t_lake(c,1)
          endif
       end if

    end do

    do p = begp, endp
       c = pcolumn(p)
       l = plandunit(p)

#if (defined VANCOUVER)
       t_veg(p) = 297.56
       t_ref2m(p) = 297.56
       if (urbpoi(l)) then
         t_ref2m_u(p) = 297.56
       else
         t_ref2m_u(p) = spval
       end if
       if (ifspecial(l)) then
         t_ref2m_r(p) = spval
       else
         t_ref2m_r(p) = 297.56
       end if 
#elif (defined MEXICOCITY)
       t_veg(p) = 289.46
       t_ref2m(p) = 289.46
       if (urbpoi(l)) then
         t_ref2m_u(p) = 289.46
       else
         t_ref2m_u(p) = spval
       end if
       if (ifspecial(l)) then
         t_ref2m_r(p) = spval
       else
         t_ref2m_r(p) = 289.46
       end if 
#elif (defined GRANDVIEW)
       ! Set to 19.0C
       t_veg(p) = 292.16
       t_ref2m(p) = 292.16
       if (urbpoi(l)) then
         t_ref2m_u(p) = 292.16
       else
         t_ref2m_u(p) = spval
       end if
       if (ifspecial(l)) then
         t_ref2m_r(p) = spval
       else
         t_ref2m_r(p) = 292.16
       end if 
#else
       t_veg(p) = 283._r8
       t_ref2m(p) = 283._r8
       if (urbpoi(l)) then
         t_ref2m_u(p) = 283._r8
       else
         t_ref2m_u(p) = spval
       end if
       if (ifspecial(l)) then
         t_ref2m_r(p) = spval
       else
         t_ref2m_r(p) = 283._r8
       end if 
#endif
       eflx_lwrad_out(p) = sb * (t_grnd(c))**4
    end do

    ! Set snow/soil ice and liquid mass

    ! volumetric water is set first and liquid content and ice lens are obtained
    ! NOTE: h2osoi_vol, h2osoi_liq and h2osoi_ice only have valid values over soil
    ! and urban pervious road (other urban columns have zero soil water)

    h2osoi_vol(begc:endc,         1:) = spval
    h2osoi_liq(begc:endc,-nlevsno+1:) = spval
    h2osoi_ice(begc:endc,-nlevsno+1:) = spval

    wa(begc:endc)  = 5000._r8
    wt(begc:endc)  = 5000._r8
    zwt(begc:endc) = 0._r8

    do c = begc,endc
       l = clandunit(c)
       if (.not. lakpoi(l)) then  !not lake
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_road_perv) then
                wa(c)  = 4800._r8
                wt(c)  = wa(c)
                zwt(c) = (25._r8 + zi(c,nlevsoi)) - wa(c)/0.2_r8 /1000._r8  ! One meter below soil column
             else
                wa(c)  = spval
                wt(c)  = spval
                zwt(c) = spval
             end if
          else
             wa(c)  = 4800._r8
             wt(c)  = wa(c)
             zwt(c) = (25._r8 + zi(c,nlevsoi)) - wa(c)/0.2_r8 /1000._r8  ! One meter below soil column
          end if
       end if
    end do

    do c = begc,endc
       l = clandunit(c)
       if (.not. lakpoi(l)) then  !not lake

          ! volumetric water
          if (ltype(l) == istsoil) then
             nlevs = nlevgrnd
             do j = 1, nlevs
                if (j > nlevsoi) then
                   h2osoi_vol(c,j) = 0.0_r8
                else
                   h2osoi_vol(c,j) = 0.3_r8
                endif
             end do
          else if (ltype(l) == isturb) then 
             nlevs = nlevurb
             do j = 1, nlevs
                if (ctype(c) == icol_road_perv .and. j <= nlevsoi) then
#if (defined GRANDVIEW)
                   h2osoi_vol(c,j) = 0.0_r8
#else
                   h2osoi_vol(c,j) = 0.3_r8
#endif
                else
                   h2osoi_vol(c,j) = 0.0_r8
                end if
             end do
          else if (ltype(l) == istwet) then
             nlevs = nlevgrnd
             do j = 1, nlevs
                if (j > nlevsoi) then
                   h2osoi_vol(c,j) = 0.0_r8
                else
                   h2osoi_vol(c,j) = 1.0_r8
                endif
             end do
          else if (ltype(l) == istice .or. ltype(l) == istice_mec) then
             nlevs = nlevgrnd 
             do j = 1, nlevs
                h2osoi_vol(c,j) = 1.0_r8
             end do
          endif
          do j = 1, nlevs
             h2osoi_vol(c,j) = min(h2osoi_vol(c,j),watsat(c,j))
        
             ! soil layers
             if (t_soisno(c,j) <= SHR_CONST_TKFRZ) then
                h2osoi_ice(c,j)  = dz(c,j)*denice*h2osoi_vol(c,j)
                h2osoi_liq(c,j) = 0._r8
             else
                h2osoi_ice(c,j) = 0._r8
                h2osoi_liq(c,j) = dz(c,j)*denh2o*h2osoi_vol(c,j)
             endif
          end do

#if (defined CN) || (defined CASA)
          ! soil water potential (added 10/21/03, PET)
          ! required for CN and CASA code
          if (ltype(l) == istsoil) then
             nlevs = nlevgrnd
             do j = 1, nlevs
                if (h2osoi_liq(c,j) > 0._r8) then
                   vwc = h2osoi_liq(c,j)/(dz(c,j)*denh2o)
                   psi = psisat(c,j) * (vwc/vwcsat(c,j))**bsw2(c,j)
                   soilpsi(c,j) = max(psi, -15.0_r8)
                   soilpsi(c,j) = min(soilpsi(c,j),0.0_r8)
                end if
             end do
          end if
#endif
       end if

    end do

    ! Set snow

    do j = -nlevsno+1, 0
       do c = begc,endc
          l = clandunit(c)
          if (.not. lakpoi(l)) then  !not lake
             if (j > snl(c)) then
                h2osoi_ice(c,j) = dz(c,j)*250._r8
                h2osoi_liq(c,j) = 0._r8
             end if
          end if
       end do
    end do


    ! initialize SNICAR fields:
    do c = begc,endc
       mss_bctot(c,:) = 0._r8
       mss_bcpho(c,:) = 0._r8
       mss_bcphi(c,:) = 0._r8
       mss_cnc_bcphi(c,:)=0._r8
       mss_cnc_bcpho(c,:)=0._r8

       mss_octot(c,:) = 0._r8
       mss_ocpho(c,:) = 0._r8
       mss_ocphi(c,:) = 0._r8
       mss_cnc_ocphi(c,:)=0._r8
       mss_cnc_ocpho(c,:)=0._r8
       
       mss_dst1(c,:) = 0._r8
       mss_dst2(c,:) = 0._r8
       mss_dst3(c,:) = 0._r8
       mss_dst4(c,:) = 0._r8
       mss_dsttot(c,:) = 0._r8
       mss_cnc_dst1(c,:)=0._r8
       mss_cnc_dst2(c,:)=0._r8
       mss_cnc_dst3(c,:)=0._r8
       mss_cnc_dst4(c,:)=0._r8
       
       if (snl(c) < 0) then
          snw_rds(c,snl(c)+1:0)        = snw_rds_min
          snw_rds(c,-nlevsno+1:snl(c)) = 0._r8
          snw_rds_top(c)               = snw_rds_min
          sno_liq_top(c) = h2osoi_liq(c,snl(c)+1) / (h2osoi_liq(c,snl(c)+1)+h2osoi_ice(c,snl(c)+1))
       elseif (h2osno(c) > 0._r8) then
          snw_rds(c,0)             = snw_rds_min
          snw_rds(c,-nlevsno+1:-1) = 0._r8
          snw_rds_top(c)           = spval
          sno_liq_top(c)           = spval
       else
          snw_rds(c,:)   = 0._r8
          snw_rds_top(c) = spval
          sno_liq_top(c) = spval
       endif
    enddo


  end subroutine mkarbinit

end module mkarbinitMod