subroutine setmtlm (is, ie, js, je)

#if defined O_mtlm
!-----------------------------------------------------------------------
!     Initialize the land surface and vegetation model
!-----------------------------------------------------------------------

      implicit none

      include "size.h"
      include "calendar.h"
      include "cembm.h"
      include "csbc.h"
      include "grdvar.h"
      include "coord.h"
      include "switch.h"
      include "levind.h"
      include "atm.h"
# if defined O_crop_data || defined O_pasture_data || defined O_agric_data
      include "veg.h"
# endif
# if defined O_ice && defined O_landice_data
#  if defined O_ice_cpts
      include "cpts.h"
#  endif
      include "ice.h"
# endif
      include "tmngr.h"
      include "mtlm.h"
      include "mtlm_data.h"

      character(120) :: fname, new_file_name

      integer i, ie, iem1, is, isp1, iou, j, je, jem1, js, jsp1, k, L, n

      logical exists, inqvardef

      real LAI_BAL, epsln
      parameter (epsln=1.0e-20)

# if defined O_mtlm_nsoil 
!------------------------------------------------------------------------------
! Soil and land surface variables
!------------------------------------------------------------------------------
! albedos   = Soil albedo
! cdensity  = Areal density of carbon in soil point (kg/m2)
! CDENS     = As above, indexed in terms of land points
! CL30      = Percentage of clay sized mineral particles from 0 - 30 cm in soil (%)
! CLAY      = As above, indexed in terms of land points
! CL150     = Percentage of clay sized mineral particles from 0 - 150 cm in soil(%) (not used)
! heatflow  = Geothermal heat flow W/m2
! iceexcess = Excess ice content (%)
! SA30      = Percentage of sand sized mineral particles from 0 - 30 cm in soil (%)
! SAND      = As above, indexed in terms of land points
! SA150     = Percentage of sand sized mineral particles from 0 - 150 cm in soil(%) (not used)
! SI30      = Percentage of silt sized mineral particles from 0 - 30 cm in soil (%)
! SILT      = As above, indexed in terms of land points
! SI150     = Percentage of silt sized mineral particles from 0 - 150 cm in soil(%) (not used)
! soildepth = Depth to bedrock (m)
! SOIL_DZ   = Depth to bedrock, indexed on land point grid (m)
! spintemp  = Temperature to initialize subsurface points to (K) (quickspinup option)
! s_map     = Fraction of a grid cell having a slope less than a set value
! s_val     = Slope value for the slope map
 
      character(120) :: text

      integer ib(10), ic(10), iem2, jem2, Z
      real tmpij(ie-2, je-2)
      real tmpi(NSLOPE), tmpijk(ie-2,je-2,NSLOPE)
      real s_val(NSLOPE), s_map(ie,je,NSLOPE)
      real SA030(ie,je), SA0150(ie,je), SA30150(ie,je)
      real SI030(ie,je), SI0150(ie,je), SI30150(ie,je)
      real CL030(ie,je), CL0150(ie,je), CL30150(ie,je) 
      real albedos(ie,je)
      real iceexcess(ie,je), heatflow(ie,je), spintemp(ie,je)
      real cdensity(ie,je), CDENS(POINTS) 
      real CLAY030(POINTS), SAND030(POINTS), SILT030(POINTS)
      real CLAY0150(POINTS), SAND0150(POINTS), SILT0150(POINTS)
      real CLAY30150(POINTS), SAND30150(POINTS), SILT30150(POINTS)
      real soildepth(ie,je), SOIL_DZ(POINTS)
      
      real c0, c1
      parameter (c0 = 0.0, c1 = 1.0)
# endif


      data MAF / 1.0, 1.0, 0.95, 0.95, 0.97, 0.95 /

      isp1 = is+1
      iem1 = ie-1
      jsp1 = js+1
      jem1 = je-1
      atlnd = 1.

!-----------------------------------------------------------------------
!     Default parameters
!-----------------------------------------------------------------------
# if defined O_mtlm_nsoil 
      iem2 = ie-2
      jem2 = je-2
#  if defined O_mtlm_gheatvari
      fname = new_file_name("L_soildata.nc")
      ib(:) = 1
      ic(1) = iem2
      ic(2) = jem2
      call openfile(fname, iou)
      call getvara ('heatflow', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      heatflow(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
      call closefile(iou)

#  endif
#  if defined O_mtlm_quickspin
      fname = new_file_name("L_spinuptemp.nc")
      ib(:) = 1
      ic(1) = iem2
      ic(2) = jem2
      call openfile(fname, iou)
      call getvara ('spintemp', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      spintemp(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
     
      call closefile(iou)

#  endif
#  if defined O_mtlm_exstruct || O_mtlm_exslab
      fname = new_file_name("L_soildata.nc")
      ib(:) = 1
      ic(1) = iem2
      ic(2) = jem2
      call openfile(fname, iou)
      call getvara ('groundice', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      iceexcess(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
     
      call closefile(iou)

#  endif
      fname = new_file_name("L_soildata.nc")
      ib(:) = 1
      ic(1) = NSLOPE
      ic(2) = iem2
      ic(3) = jem2
      call openfile(fname, iou)
      call getvara ('fracslope', iou, NSLOPE*iem2*jem2, ib, ic, tmpijk
     &,              c1, c0)

      s_map(2:iem1, 2:jem1, 1:NSLOPE) = tmpijk(1:iem2, 1:jem2, 1:NSLOPE)
      call getvara ('slope', iou, NSLOPE, ib, ic, tmpi, c1, c0)
      SLOPEVAL(1:NSLOPE) = tmpi(1:NSLOPE)
      
      ib(:) = 1
      ic(1) = iem2
      ic(2) = jem2
     
      call getvara ('SA030', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      SA030(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
      call getvara ('SA0150', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      SA0150(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
      call getvara ('SA30150', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      SA30150(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)

      call getvara ('SI030', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      SI030(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
      call getvara ('SI0150', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      SI0150(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
      call getvara ('SI30150', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      SI30150(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)

      call getvara ('CL030', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      CL030(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
      call getvara ('CL0150', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      CL0150(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)

      call getvara ('CL30150', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      CL30150(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)

      call getvara ('soil_albedo', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      albedos(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
     
      ib(:) = 1
      ic(1) = iem2
      ic(2) = jem2
     
      call getvara ('cdensity', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      cdensity(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
      
      call closefile(iou)

# if defined O_mtlm_soildepth
      fname = new_file_name("L_soildata.nc")
      ib(:) = 1
      ic(1) = iem2
      ic(2) = jem2
      call openfile(fname, iou)
      call getvara ('soildepth', iou, iem2*jem2, ib, ic, tmpij
     &,              c1, c0)
      soildepth(2:iem1, 2:jem1) = tmpij(1:iem2, 1:jem2)
     
      call closefile(iou)

#  endif
      SCHON = 1
      SCH_COUNT = 0
      YEARCOUNT = 0.
      H2OHEAT = 0.

# endif
# if defined O_embm
      LHC = vlocn*1.e-4
      LHF = flice*1.e-4
# else
      LHC = 2.501E6
      LHF = 0.334E6
# endif
      SIGMA = 5.67E-8
      DAY_YEAR = yrlen
      SEC_DAY = daylen
      SEC_YEAR = DAY_YEAR*SEC_DAY
      STEP_DAY = INT(SEC_DAY/TIMESTEP)
!     set longitude to half the STEP_DAY interval to line up with noon.
!     longitude is set to a constant to treat all latitudes equally.
      LONG(:) = 360./STEP_DAY/2.
      LAND_COUNTER = 0
      dtlnd = TIMESTEP

!----------------------------------------------------------------------
!     Initialization of arrays
!----------------------------------------------------------------------
      PSTAR(:) = 1.e5
      DTEMP_DAY(:) = 0.
      LYING_SNOW(:) = 0.
      TSTAR(:,:) = 280.
      TSOIL(:) = 280.0
      TS1(:) = 280.0
# if !defined O_mtlm_nsoil 
      CS(:) = 10.0
      M(:) = 242.
# endif
      MNEG(:) = 0.
      FSMC(:) = 1.0
      RESP_S_DR(:) = 0.0
      ALBSOIL(:) = 0.3
      ALBSNOW(:) = 0.6
      Z0S(:) = 0.0003
      FRACA(:) = 0.0
      FRAC(:,:) = 0.1
      LAI(:,1:2) = 6.
      LAI(:,3:5) = 2.
      HT(:,1:2) = 21.46
      HT(:,3:4) = 0.794
      HT(:,5) = 1.587
      NPP_DR(:,:) = 0.0
      G_LEAF_DR(:,:) = 0.0
      RESP_W_DR(:,:) = 0.0
# if defined O_mtlm_nsoil 
      ISEXSLAB(:) = .false.
      FROZENMAP(:,:) = 0.
      TSTAR_P(:,:) = 280.
      TSTAR_S(:) = 280.0
      TGND(:,:) = 263.15
      CS(:) = 15.0
      CLAY030(:) = 0.
      SAND030(:) = 0.
      SILT030(:) = 0.
      CLAY0150(:) = 0.
      SAND0150(:) = 0.
      SILT0150(:) = 0.
      CLAY30150(:) = 0.
      SAND30150(:) = 0.
      SILT30150(:) = 0. 
      SOIL_DZ(:) = 10.
      DFROZEN_C(:,:) = 730.
      DFROZEN_A(:,:) = 0.
      DUNFROZEN_C(:,:) = 730.
      M(:,:) = 0.
      SU(:,:) = 0.
      SF(:,:) = 0.
      SLOPEMAP(:,:) = 0.
      FSMC_PFT(:,:) = 1.0
      RESP_S_DD(:,:) = 0.0
      GEOHEAT(:) = 0.0
      EXCESSI(:) = 0.0
      EXCOL(:) = 0.0
      ANN_DWET(:) = 0.0
      SCHDMIN(:) = 0.0
      DMIN_HOLD(:) = 9999
      TRACK_PFC(:) = 0.
      SQU_PFC(:) = 0.
      FZN_SL(:) = 0.
#  if defined O_mtlm_kovdiff
      YEARKOV_SW = 0
      CS_P(:,:)=0.0
      AFp(:,:)= 0.
      AF_DTa(:,:) = 0.0
      RESP_S_DP(:,:) = 0.0
      PREG_C(:) = 0.0
#  endif       
      
!----------------------------------------------------------------------
! Define soil carbon with depth
!----------------------------------------------------------------------
! CS set to 10 kg/m2 need this to match put in layer
      CS_D(:,:)=0.0
      CS_D(:,1)=150.0
# endif

!-----------------------------------------------------------------------
!     Define externally dependent arrays
!-----------------------------------------------------------------------
      L = 0
      do j=jsp1,jem1
        do i=isp1,iem1
          if (kmt(i,j) .le. klmax) then
            L = L + 1
# if defined O_mtlm_nsoil 
            SLOPEMAP(L,:) = s_map(i,j,:)
#  if defined O_mtlm_gheatvari
            GEOHEAT(L) = heatflow(i,j)
#  endif           

#  if defined O_mtlm_gheatconst
            GEOHEAT(L) = 50./1000.
#  endif

#  if defined O_mtlm_quickspin
            SPINUPTEMP(L) = spintemp(i,j)
#  endif

#  if defined O_mtlm_exstruct
#   if defined O_mtlm_initPF
                if(iceexcess(i,j) .gt. 0.) then
                   EXCESSI(L) = iceexcess(i,j)
                endif
#   endif
#  endif

#  if defined O_mtlm_exslab
#    if defined O_mtlm_initPF              
                if(iceexcess(i,j) .gt. 0.) then
                   EXCESSI(L) = iceexcess(i,j)
                   ISEXSLAB(L) = .true.    
                endif
#    endif         
#  endif
# endif
            GAREA(L) = dxt(i)*dyt(j)*cst(j)*1e-4
# if defined O_mtlm_pressure
            PSTAR(L) = PSTAR(L)*exp((-1.)*elev(i,j)/7.3e5)
# endif
# if defined O_crop_data || defined O_pasture_data || defined O_agric_data
            FRACA(L) = agric(i,j,2)
# endif
            LAT(L) = tlat(i,j)
            sbc(i,j,isca) = 1. - ALBSOIL(L)
            sbc(i,j,ievap) = 0.
            sbc(i,j,isens) = 0.
            sbc(i,j,ilwr) = 0.
          endif
        enddo
      enddo

      if (L .gt. POINTS) then
        print*, "==> Error: Number of land points is inconsistent"
        print*, "==>        set POINTS in size.h to: ", L
        stop
      endif

!----------------------------------------------------------------------
!     Initialize the non-vegetation fractions
!----------------------------------------------------------------------
      FRAC(:,SOIL) = 1.0
      do N=1,NPFT
        FRAC(:,SOIL) = FRAC(:,SOIL) - FRAC(:,N)
      enddo

!----------------------------------------------------------------------
!     Initialize the vegetation carbon contents
!----------------------------------------------------------------------
      CV(:) = 0.
      G_LEAF_PHEN(:,:) = 0.0
      do N=1,NPFT
        LAI_BAL = (A_WS(N)*ETA_SL(N)*HT(1,N)/A_WL(N))
     &                 **(1.0/(B_WL(N)-1))
        C_VEG(:,N) = 2*SIGL(N)*LAI_BAL
     &             + A_WS(N)*ETA_SL(N)*HT(:,N)*LAI_BAL
        CV(:) = CV(:) + C_VEG(:,N)*FRAC(:,N)
      enddo

!----------------------------------------------------------------------
!     Derive vegetation parameters from the areal fractions and the
!     structural properties.
!----------------------------------------------------------------------
      L = 0
      LAND_INDEX(:) = 0
      do j=jsp1,jem1
        do i=isp1,iem1
          if (kmt(i,j) .le. klmax) then
            L = L + 1
            LAND_INDEX(L) = L
          endif
        enddo
      enddo
# if !defined O_mtlm_nsoil 
      do N=1,NPFT
        call PFT_SPARM (L, LAND_INDEX, N, ALBSOIL, HT(1,N), LAI(1,N)
     &,                 ALBSNC(1,N), ALBSNF(1,N), CATCH(1,N), Z0(1,N))
      enddo
# endif

!----------------------------------------------------------------------
!     Define other vegetation parameters
!----------------------------------------------------------------------
      VEG_FRAC(:) = 0.0
      do N=1,NPFT
        VEG_FRAC(:) = VEG_FRAC(:) + FRAC(:,N)
      enddo
      FRAC_VS(:) = VEG_FRAC(:) + FRAC(:,SOIL)

!----------------------------------------------------------------------
!       Reading a restart file
!----------------------------------------------------------------------
      if (.not. init) then
        fname = new_file_name ('restart_mtlm.nc')
        inquire (file=trim(fname), exist=exists)
        if (exists) call mtlm_rest_in (fname, is, ie, js, je)
# if defined O_restart_2
        fname = new_file_name ("restart_2_mtlm.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) call mtlm_rest_in (fname, is, ie, js, je)
# endif
      endif

!----------------------------------------------------------------------
!     Create the VEG_INDEX array of land points with each type
!----------------------------------------------------------------------
      L = 0
      land_map(:,:) = 0
      LAND_PTS = 0
      LAND_INDEX(:) = 0
      do j=jsp1,jem1
        do i=isp1,iem1
          if (kmt(i,j) .le. klmax) then
            L = L + 1
# if defined O_ice && defined O_landice_data
            if (aicel(i,j,2) .lt. 0.5 .and. tmsk(i,j) .lt. 0.5) then
              land_map(i,j) = L
              LAND_PTS = LAND_PTS + 1
              LAND_INDEX(LAND_PTS) = L
            endif
# else
            if (tmsk(i,j) .lt. 0.5) then
              land_map(i,j) = L
              LAND_PTS = LAND_PTS + 1
              LAND_INDEX(LAND_PTS) = L
            endif
# endif
# if defined O_mtlm_nsoil            
            CLAY030(L) = CL030(i,j)
            SAND030(L) = SA030(i,j)
            SILT030(L) = SI030(i,j)
            CLAY0150(L) = CL0150(i,j)
            SAND0150(L) = SA0150(i,j)
            SILT0150(L) = SI0150(i,j)
            CLAY30150(L) = CL30150(i,j)
            SAND30150(L) = SA30150(i,j)
            SILT30150(L) = SI30150(i,j)
            CDENS(L) = cdensity(i,j)
            SOIL_DZ(L) = soildepth(i,j)           
#  if defined O_mtlm_soildepth
            SOIL_DZ(L) = soildepth(i,j)
#  endif
# endif
          endif
        enddo
      enddo
# if defined O_mtlm_nsoil            

      call  SETSOIL(init, SAND030, SILT030, CLAY030
     &,             SAND0150, SILT0150, CLAY0150
     &,             SAND30150, SILT30150, CLAY30150, CDENS, SOIL_DZ)
# endif     
      do N=1,NPFT
        VEG_PTS(N) = 0
        do J=1,LAND_PTS
          L = LAND_INDEX(J)
          if (FRAC(L,N) .gt. FRAC_MIN + epsln) then
            VEG_PTS(N) = VEG_PTS(N) + 1
            VEG_INDEX(VEG_PTS(N),N) = L
          endif
        enddo
      enddo

# if defined O_mtlm_segday
      if (DAY_TRIF .lt. segtim) then
        print*, ""
        print*, "==> Warning: DAY_TRIF is set to be less than segtim."
        print*, "             with option mtlm_segday, triffid will"
        print*, "             only be done once every coupling time."
      endif
      if (DAY_PHEN .lt. segtim) then
        print*, ""
        print*, "==> Warning: DAY_PHEN is set to be less than segtim."
        print*, "             with option mtlm_segday, phenology will"
        print*, "             only be done once every coupling time."
      endif
      if (segtim .lt. 1.) then
        print*, ""
        print*, "==> Error: segtim must be greater than one when using"
        print*, "           the option mtlm_segday. Turn off this"
        print*, "           option if segtim is less than one."
        stop
      endif
# endif
      if (STEP_DAY .gt. STEPSM) then
        print*, ""
        print*, "==> Error: STEPSM is too small. Increase TIMESTEP or "
        print*, "           set STEPSM in size.h to: ", STEP_DAY
        stop
      endif
# if defined O_mtlm_segday
      if (mod(SEC_DAY*segtim,TIMESTEP) .gt. 1.e-6) then
        print*, ""
        print*, "==> Error: there must be an integral number of mtlm "
        print*, "           timesteps in a coupling time."
        stop
      endif
      if (DAY_TRIF .gt. 1) then
        if (mod(FLOAT(DAY_TRIF),segtim) .gt. 1.e-6) then
          print*, '==> Error: there must be an integral number coupling'
     &,     ' timesteps within DAY_TRIF when using O_mtlm_segday.'
          stop
        else
          DAY_TRIF = int(float(DAY_TRIF)/segtim)
        endif
      endif
      if (DAY_PHEN .gt. 1) then
        if (mod(FLOAT(DAY_PHEN),segtim) .gt. 1.e-6) then
          print*, '==> Error: there must be an integral number coupling'
     &,     ' timesteps within DAY_PHEN when using O_mtlm_segday.'
          stop
        else
          DAY_PHEN = int(float(DAY_PHEN)/segtim)
        endif
      endif
# else
      if (mod(SEC_DAY,TIMESTEP) .gt. 1.e-6) then
        print*, ""
        print*, "==> Error: there must be an integral number of mtlm "
        print*, "           timesteps in a day."
        stop
      endif
# endif
# if defined O_time_averages

!-----------------------------------------------------------------------
!     zero time average accumulators
!-----------------------------------------------------------------------

      call ta_mtlm_tavg (is, ie, js, je, 0)
# endif
# if defined O_time_step_monitor

!-----------------------------------------------------------------------
!     zero integrated time average accumulators
!-----------------------------------------------------------------------

      call ta_mtlm_tsi (is, ie, js, je, 0)
# endif
#endif

      return
      end