! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/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) !------------------------------------------------------------------------------ ! 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) 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 !----------------------------------------------------------------------- iem2 = ie-2 jem2 = je-2 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) 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) SCHON = 1 SCH_COUNT = 0 YEARCOUNT = 0. H2OHEAT = 0. 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 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 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. YEARKOV_SW = 0 CS_P(:,:)=0.0 AFp(:,:)= 0. AF_DTa(:,:) = 0.0 RESP_S_DP(:,:) = 0.0 PREG_C(:) = 0.0 !---------------------------------------------------------------------- ! 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 !----------------------------------------------------------------------- ! Define externally dependent arrays !----------------------------------------------------------------------- L = 0 do j=jsp1,jem1 do i=isp1,iem1 if (kmt(i,j) .le. klmax) then L = L + 1 SLOPEMAP(L,:) = s_map(i,j,:) GEOHEAT(L) = heatflow(i,j) 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 !---------------------------------------------------------------------- ! 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 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) endif enddo enddo call SETSOIL(init, SAND030, SILT030, CLAY030 &, SAND0150, SILT0150, CLAY0150 &, SAND30150, SILT30150, CLAY30150, CDENS, SOIL_DZ) 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