! source file: /Users/csomes/Research/Models/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