! source file: /sfs/fs6/home-geomar/smomw258/UVic_ESCM/2.9/updates/02/source/mtlm/setmtlm.F
      subroutine setmtlm (is, ie, js, je)

!-----------------------------------------------------------------------
!     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"

      include "veg.h"

      include "ice.h"

      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)

      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
!-----------------------------------------------------------------------

      LHC = vlocn*1.e-4
      LHF = flice*1.e-4

      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
      CS(:) = 10.0
      M(:) = 242.
      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

!-----------------------------------------------------------------------
!     Define externally dependent arrays
!-----------------------------------------------------------------------
      L = 0
      do j=jsp1,jem1
        do i=isp1,iem1
          if (kmt(i,j) .le. klmax) then
            L = L + 1
            GAREA(L) = dxt(i)*dyt(j)*cst(j)*1e-4

            FRACA(L) = agric(i,j,2)

            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
      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

!----------------------------------------------------------------------
!     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)

      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 (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

          endif
        enddo
      enddo
      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 (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

      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 (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

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

      call ta_mtlm_tavg (is, ie, js, je, 0)

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

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

      return
      end