! source file: /Users/jfidler/work/UVic_ESCM/2.9/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) = crops(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