program hist_vcf75_graz ! Version 1.0 ! ! This is a program to create land cover conditions for known years of R&F agricultural and HYDE ! grazing data from present back to 1870 ! ! by Johannes Feddema ! Department of Geography ! University of Kansas ! Lawrence, KS 66045 ! implicit none ! ! dimensions - basic parameters first ! integer, parameter :: r8 = selected_real_kind(12) integer, parameter :: nlon = 720 !input grid : longitude points integer, parameter :: nlat = 360 !input grid : latitude points integer, parameter :: numpft = 16 !input grid : latitude points integer, parameter :: lc_class = 2 !land cover type to be evaluated -- HYDE only integer, parameter :: water = -9999 !land cover type to be ignored -- HYDE only ! ! dimension PFT arrays to use ! real(r8) :: pft_in(0:numpft) ! input PFT data real(r8) :: pft_pot(0:numpft) ! input potential PFT data real(r8) :: pft_out(0:numpft) ! output PFT data ! ! 0 = bare ! 1 = needle leaf evergreen - temperate ! 2 = needle leaf evergreen - boreal ! 3 = needle leaf deciduous - boreal ! 4 = broadleaf evergreen tree - tropical ! 5 = broadleaf evergreen tree - temperate ! 6 = broadleaf deciduous tree - tropical ! 7 = broadleaf deciduous tree - temperate ! 8 = broadleaf deciduous tree - boreal ! 9 = broadleaf evergreen shrub - temperate ! 10= broadleaf deciduous shrub - temperate ! 11= broadleaf deciduous shrub - boreal ! 12= C3 arctic grass ! 13= C3 grass ! 14= C4 grass ! 15= crop ! 16= crop 2 !! Not used ! ! dimension other arrays ! real(r8) :: rf_ag(nlon,nlat) ! input R&F ag for time 1 integer :: hyde_lc_in(nlon,nlat) ! input HYDE landcover time 1 integer :: hyde_lc_out(nlon,nlat) ! input HYDE landcover time 2 ! real(r8) :: landmask ! NCAR land mask value ! real(r8) :: delta_ag ! change in ag area real(r8) :: del_per_ag ! percent change in ag real(r8) :: in_veg ! total input vegetation (excluding ag and bare) real(r8) :: out_veg ! total output vegetation (excluding ag and bare) real(r8) :: pot_veg ! total potential vegetation (excluding ag and bare) ! real(r8) :: grass_in ! input total grass real(r8) :: grass_pot ! potential total grass real(r8) :: grass_out ! output total grass real(r8) :: in_grazing ! input HYDE grazing real(r8) :: out_grazing ! output HYDE grazing real(r8) :: clm_graz ! CLM grass with HYDE percent change grazing real(r8) :: grass_diff ! in/out grass difference real(r8) :: grass_ipot_diff ! in/pot grass difference real(r8) :: grass_opot_diff ! out/pot grass difference real(r8) :: pct_grass_change ! percent change in total grass real(r8) :: grass_in_pct ! input pft to total grass percentage real(r8) :: grass_pot_pct ! potential pft to total grass percentage ! real(r8) :: tree_in ! input total tree real(r8) :: tree_pot ! input total tree real(r8) :: tree1_diff ! in/pot grass difference real(r8) :: tree_out ! input total tree real(r8) :: tree_diff ! in/out grass difference real(r8) :: tree_ipot_diff ! in/pot grass difference real(r8) :: tree_opot_diff ! out/pot grass difference real(r8) :: pct_tree_change ! percent change in total tree real(r8) :: tree_in_pct ! input pft to total grass percentage real(r8) :: tree_pot_pct ! potential pft to total grass percentage ! real(r8) :: grazing ! total grazing area (percent of total area in cell real(r8) :: pct_grazing ! percent of grass area used for grazing (grass = sum pft 9-14) real(r8) :: graz ! grazing area tracking variable real(r8) :: graz_9 ! percent of pft 9 grazed real(r8) :: graz_10 ! percent of pft 10 grazed real(r8) :: graz_11 ! percent of pft 11 grazed real(r8) :: graz_12 ! percent of pft 12 grazed real(r8) :: graz_13 ! percent of pft 13 grazed real(r8) :: graz_14 ! percent of pft 14 grazed ! integer :: i integer :: j integer :: k integer :: m integer :: ix integer :: iy integer :: scen integer :: time_1 integer :: time_2 ! real :: rdummy real :: lat real :: lon real(r8) :: sum_in real(r8) :: sum_out real(r8) :: sum_pot real(r8) :: sum_pft real(r8) :: sumtest_pft ! character*7 :: dummy character*20 :: gisheader(6) Character*80 :: file10 Character*80 :: file11 Character*80 :: file12 Character*80 :: file13 Character*80 :: file15 ! character*40 :: hyde1 character*40 :: hyde2 character*40 :: rf1 character*40 :: rf2 character*40 :: out1 character*40 :: out2 ! real(r8) :: lc_frac ! function ! ! Open the files for this time step ! ! The files have to be processed in reverse chronological order ! ! get file name information ! hyde1 = 'input/h' hyde2 = '.asc' rf1 = 'input/glcrop' rf2 = '_0_5.asc' out1 = 'output/pft-' out2 = '-hist.5x.5' ! ! begin processing do all 7 time steps back in time ! do scen = 1,7 ! ! get the file name information ! if (scen .eq. 1) then time_1 = 1990 time_2 = 1970 elseif (scen .eq. 2) then time_1 = 1970 time_2 = 1950 elseif (scen .eq. 3) then time_1 = 1950 time_2 = 1900 elseif (scen .eq. 4) then time_1 = 1900 time_2 = 1850 elseif (scen .eq. 5) then time_1 = 1850 time_2 = 1800 elseif (scen .eq. 6) then time_1 = 1800 time_2 = 1750 elseif (scen .eq. 7) then time_1 = 1750 time_2 = 1700 endif ! write(file10, '(a7,i4,a4)') hyde1,time_1,hyde2 write(file11, '(a7,i4,a4)') hyde1,time_2,hyde2 write(file12, '(a12,i4,a8)') rf1,time_2,rf2 write(file13, '(a11,i4,a10)') out1,time_1,out2 write(file15, '(a11,i4,a10)') out1,time_2,out2 ! open(10, file = file10) ! file for time 1 open(11, file = file11) ! CLM PFT input file time 1 open(12, file = file12) ! file for time 2 open(13, file = file13) ! file for decadal change time 2 open(14, file = 'output/pft-pot.5x.5') ! CLM PFT potential file open(15, file = file15) ! file for decadal change time 2 ! ! read over the GIS header info and write GIS header for output file ! do i=1,6 read(10,'(A20)')gisheader(i) read(11,*)dummy read(12,*)dummy enddo ! ! read the data ! do i = 1,nlat read(10,*) (hyde_lc_in(j,i),j=1,nlon) read(11,*) (hyde_lc_out(j,i),j=1,nlon) read(12,*) (rf_ag(j,i),j=1,nlon) enddo write(6,*)'done GIS file reads' ! ! Read in NCAR data ! do j=nlat,1,-1 write(6,*)'processing lat ', j do i=1,nlon ! ! Read in NCAR data: landmask and pft info ! read (13,*) rdummy,rdummy,landmask,(pft_in(m),m=0,numpft),rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy,rdummy,rdummy read (14,*) rdummy,(pft_pot(m),m=0,numpft) ! ! set pft_out to 0 ! do k=0,numpft pft_out(k) = 0.0 enddo ! ! only work on land grid cells ! if (landmask == 100.) then ! ! set gis missing values to 0.0 (different land masks) ! if(rf_ag(i,j) .lt. 0.0) rf_ag(i,j) = 0.0 ! ! Get HYDE grazing information ! ! calculate the fraction of grazing land for each time period (hyde data) ! assume grazing can occupy no more than 75 percent i.e. scale to max 75% ! in_grazing = lc_frac(hyde_lc_in,i,j,lc_class,water) * 75.0 out_grazing = lc_frac(hyde_lc_out,i,j,lc_class,water) * 75.0 ! ! Initialize tree and grass values and set lon and lat values ! lat = -1.0 * real(j) / 2.0 + 90.25 lon = real(i) / 2.0 - 180.25 ! ! Determine the grass proportion of the grid cell ! grass_in = pft_in(9) + pft_in(10) + pft_in(11) + & pft_in(12) + pft_in(13) + pft_in(14) tree_in = 100.0 - grass_in - pft_in(0) - pft_in(15) grass_pot = pft_pot(9) + pft_pot(10) + pft_pot(11) + & pft_pot(12) + pft_pot(13) + pft_pot(14) tree_pot = 100.0 - grass_pot - pft_in(0) ! ! set out values equal to in values for non changed cells and as initial conditions ! grass_out = grass_in tree_out = tree_in ! ! Set bare ground fraction ! pft_out(0) = pft_in(0) ! ! Only process grid cells where there is some ag or grazing information ! if (pft_in(15).gt.0.0 .or. rf_ag(i,j).gt.0.0 .or. & in_grazing.gt.0.0.or.out_grazing.gt.0.0)then ! ! Set AGRICULTURE based on R&F data and calculate ag stats ! pft_out(15) = rf_ag(i,j) * 100.0 ! ! get the change in agriculture ! delta_ag = rf_ag(i,j) * 100.0 - pft_in(15) ! out_veg = 100.0 - pft_out(0) - pft_out(15) in_veg = 100.0 - pft_in(0) - pft_in(15) pot_veg = 100.0 - pft_in(0) ! if (out_veg .lt. 0.0) then out_veg = 0.0 pft_out(15) = 100.0 - pft_out(0) endif ! ! determine the change in tree cover based on Ag change -- provisional tree_out ! ! calculate the provisional tree cover change straight ! if (delta_ag .lt. 0.0) then ! ! Agriculture area is decreasing back in time -- ! allocate portion of ag land to be tree ! tree_out = tree_in - tree_in / in_veg * delta_ag ! ! Make adjustment for tree level to go to potential level as ! agriculture area approaches 0 ! tree_out = tree_in - (tree_pot - tree_in) * delta_ag / pft_in(15) ! elseif(delta_ag .gt. 0.0)then ! ! Agriculture area is increasing back in time -- ! assume regrowth to present is based on potential vegetation ratios ! tree_out = tree_in - tree_pot / pot_veg * delta_ag ! else ! ! No change in delta ag -- keep tree area the same for now ! tree_out = tree_in ! ! Make adjustment for tree level to go to potential level as ! agriculture area approaches 0 ! if(pft_out(15) .eq. 0.0) tree_out = tree_pot ! endif ! ! make sure tree amount is in appropriate range ! if(tree_out .gt. out_veg) tree_out = out_veg if(tree_out .lt. 0.0) tree_out = 0.0 ! ! END Initial Tree calculation ! ! START grass ! ! get initial grass percentage ! grass_out = out_veg - tree_out ! ! If grass is insufficient for grazing adjust grass and tree proportions ! if (out_grazing.gt.grass_out) then ! ! see if the change in tree percent is enough to make the grazing portion ! if (out_grazing .lt. (grass_out + tree_out - tree_in)) then ! ! keeping the trees at incoming level is sufficient to allow grazing ! grass_out = out_grazing tree_out = out_veg - grass_out ! else ! ! need to reduce the tree percentage from in values to allow grazing ! scale this according the trees present and prior grazing values ! ! first do a check that compares the CLM based grass coverage change ! make the change based on the percent difference between HYDE and CLM data ! if (in_grazing .gt. 0.01)then clm_graz = grass_in * out_grazing/in_grazing else clm_graz = out_grazing endif if (clm_graz .ge. out_grazing) clm_graz = out_grazing ! if (clm_graz .lt. out_veg) then ! ! If clm percentage adjusted grazing is less than remnant veg then ! expand grazing if not things remain as is ! if(clm_graz .gt. (out_veg - tree_in)) then ! ! expand grazing based on this value if it is less than the total remnant veg ! grass_out = clm_graz tree_out = out_veg - grass_out ! endif ! else ! ! CLM adjusted value is still greater than the grass + tree areas -- remove trees ! grass_out = out_veg tree_out = 0.0 ! endif endif endif ! ! END grass now check to make sure all proportions do not exceed limits ! ! check to see if totals add up and get tree out ! if ((pft_out(0) + pft_out(15) + grass_out) .gt. 100.0) then grass_out = 100.0 - pft_out(0) - pft_out(15) tree_out = 0.0 else tree_out = 100.0 - pft_out(0) - pft_out(15) - grass_out endif ! ! Adjust the specific TREE PFT proportions to move towards ! the potential distribution as appropriate ! ! use the difference between in tree and potential tree percentages to do the calibration ! tree_ipot_diff = tree_pot - tree_in tree_opot_diff = tree_pot - tree_out tree_diff = tree_out - tree_in ! ! check to see if the potential tree proportion is reached ! if (tree_ipot_diff .eq. 0.0)then ! ! correct area of tree at outset, but change relative ratios of PFTs ! to go to potential distribution (just in case) ! pct_tree_change = 0.0 ! elseif((tree_diff * tree_ipot_diff) .le. 0.0)then ! ! change is moving further away from potential target ! assume the distribution of tree species remains as is ! pct_tree_change = 0.0 ! elseif ((tree_ipot_diff * tree_opot_diff) .lt. 0.0) then ! ! change in tree area overshoots the tree level at potential ! change individual tree ratios to potential level ! pct_tree_change = 1.0 ! else ! ! change is towards potential distribution - change the ! relatave abundance of species by this ratio to the ! ratio at potential level ! if(abs(tree_diff) .lt. 0.01) then pct_tree_change = 0.0 else pct_tree_change = tree_diff / tree_ipot_diff endif endif ! ! Use this percentage change to change the proportions of individual tree PFTs ! do k = 1,8 ! ! calculate pft contribution to tree area at start ! if(tree_in .lt. 0.001) then tree_in_pct = 0.0 else tree_in_pct = pft_in(k) / tree_in endif ! ! calculate pft contribution to tree area at potential ! if(tree_pot .lt. 0.001) then tree_pot_pct = 0.0 else tree_pot_pct = pft_pot(k) / tree_pot endif ! ! calculate the percent value for the out time period based on linear change ! to the potential type (overshoots not adjusted) ! if(tree_pot .lt. 0.01) then pft_out(k) = tree_out * tree_in_pct elseif(tree_in .lt. 0.01)then pft_out(k) = tree_out * tree_pot_pct else pft_out(k) = tree_out * & (tree_in_pct - (tree_in_pct - tree_pot_pct) * pct_tree_change) endif ! enddo ! ! DONE tree adjustments ! ! ! Adjust the specific GRASS PFT proportions to move towards ! the potential distribution as appropriate ! ! use the difference between in grass and potential grass percentages to do the calibration ! grass_ipot_diff = grass_pot - grass_in grass_opot_diff = grass_pot - grass_out grass_diff = grass_out - grass_in ! ! check to see if the potential tree proportion is reached ! if (grass_ipot_diff .eq. 0.0)then ! ! correct area of grass at outset, but change relative ratios of PFTs ! to go to potential distribution (just in case) ! pct_grass_change = 0.0 ! elseif((grass_diff * grass_ipot_diff) .lt. 0.0)then ! ! change is moving further away from potential target ! assume the distribution of grass species remains as is ! pct_grass_change = 0.0 ! elseif ((grass_ipot_diff * grass_opot_diff) .lt. 0.0) then ! ! change in grass area overshoots the grass level at potential ! change individual grass ratios to potential level ! pct_grass_change = 1.0 else ! ! change is towards potential distribution - change the ! relatave abundance of species by this ratio to the ! ratio at potential level ! if(grass_diff .lt. 0.001) then pct_grass_change = 0.0 else pct_grass_change = grass_diff / grass_ipot_diff endif endif ! ! Use this perecentage change to change the proportions of individual grass PFTs ! do k = 9,14 ! ! calculate pft contribution to grass area at start ! if(grass_in .lt. 0.001) then grass_in_pct = 0.0 else grass_in_pct = pft_in(k) / grass_in endif ! ! calculate pft contribution to grass area at potential ! if(grass_pot .lt. 0.001) then grass_pot_pct = 0.0 else grass_pot_pct = pft_pot(k) / grass_pot endif ! ! calculate the percent value for the out time period based on linear change ! to the potential type (overshoots not adjusted) ! if(grass_pot .lt. 0.01) then pft_out(k) = grass_out * grass_in_pct elseif(grass_in .lt. 0.01)then pft_out(k) = grass_out * grass_pot_pct else pft_out(k) = grass_out * & (grass_in_pct - (grass_in_pct - grass_pot_pct) * pct_grass_change) endif ! enddo ! ! DONE grass adjustments ! else do k = 0,16 pft_out(k) = pft_in(k) enddo endif ! end ag/grazing check endif ! end if landmask if ! ! make sure the data is precise to the 6th decimal place ! sumtest_pft = 100.0 - pft_out(0) - pft_out(15) sum_pft = 0.0 do k = 1, 14 sum_pft = sum_pft + pft_out(k) enddo ! do k = 1, 14 if(sum_pft .gt. 0.00000000000001)then pft_out(k) = pft_out(k) * sumtest_pft / sum_pft endif enddo ! ! ! initialize grazing values to 0._r8 ! grazing = 0._r8 pct_grazing = 0._r8 graz = 0._r8 graz_9 = 0._r8 graz_10 = 0._r8 graz_11 = 0._r8 graz_12 = 0._r8 graz_13 = 0._r8 graz_14 = 0._r8 ! ! only work on land grid cells ! if (landmask == 100.) then ! if(out_grazing .lt. grass_out)then grazing = out_grazing else grazing = grass_out endif ! ! determine the percentage of all grasses grazed ! if(grass_out.lt. 0.01)then pct_grazing = 0.0 else pct_grazing = grazing / grass_out * 100._r8 endif ! ! determine how much of each grass is to be grazed ! ! set grazing area ! graz = grazing if (graz .lt. 0.01) graz = 0._r8 ! ! assume that C3 and C4 grasses will be consumed 1st ! if (graz .gt. 0.01) then ! only process grazing if there is any ! if((pft_out(14) + pft_out(13)) .gt. 0.01)then if(graz .gt. (pft_out(14) + pft_out(13)))then if(pft_out(13).gt.0.01)graz_13 = 100._r8 if(pft_out(14).gt.0.01)graz_14 = 100._r8 graz = graz - pft_out(14) - pft_out(13) else if(pft_out(13) .gt. 0.01) then graz_13 = graz / (pft_out(14) + pft_out(13)) * 100._r8 endif ! if(pft_out(14) .gt. 0.01) then graz_14 = graz / (pft_out(14) + pft_out(13)) * 100._r8 endif graz = 0._r8 endif endif if (graz .lt. 0.01) graz = 0._r8 ! ! now test for C3 arctic grass ! if(pft_out(12) .gt. 0.01)then if(graz .gt. pft_out(12))then graz_12 = 100._r8 else graz_12 = graz / pft_out(12) * 100.0_r8 endif endif graz = graz - graz_12/100._r8 * pft_out(12) if (graz .lt. 0.01) graz = 0._r8 ! ! now test for shrubs temp evergreen first ! if(pft_out(9) .gt. 0.01)then if(graz .gt. pft_out(9))then graz_9 = 100._r8 else graz_9 = graz / pft_out(9) * 100.0_r8 endif endif graz = graz - graz_9/100._r8 * pft_out(9) if (graz .lt. 0.01) graz = 0._r8 ! ! now test for decid temperate/tropical shrub ! if(pft_out(10) .gt. 0.01)then if(graz .gt. pft_out(10))then graz_10 = 100._r8 else graz_10 = graz / pft_out(10) * 100.0_r8 endif endif graz = graz - graz_10/100._r8 * pft_out(10) if (graz .lt. 0.01) graz = 0._r8 ! ! now test for boreal shrub ! if(pft_out(11) .gt. 0.01)then if(graz .gt. pft_out(11))then graz_11 = 100._r8 else graz_11 = graz / pft_out(11) * 100.0_r8 endif endif graz = graz - graz_11/100._r8 * pft_out(11) if (graz .gt. 0.01) write(6,*)'not all grazing allocated',graz ! endif ! end no grazing block ! if(pct_grazing .lt. 0.01)pct_grazing = 0._r8 if(graz_9 .lt. 0.01)graz_9 = 0._r8 if(graz_10 .lt. 0.01)graz_10 = 0._r8 if(graz_11 .lt. 0.01)graz_11 = 0._r8 if(graz_12 .lt. 0.01)graz_12 = 0._r8 if(graz_13 .lt. 0.01)graz_13 = 0._r8 if(graz_14 .lt. 0.01)graz_14 = 0._r8 ! endif ! landmask if block ! ! write out the results ! write (15,*) lat,lon,landmask,(pft_out(m),m=0,numpft),tree_out,grass_out, & pct_grazing,graz_9,graz_10,graz_11,graz_12,graz_13,graz_14 ! ! reinitialize tree and grass percentages ! tree_out = 0._r8 grass_out = 0._r8 ! enddo enddo ! ! Close files and return the scen loop ! close(10) close(13) close(11) close(12) close(15) close(14) ! enddo ! stop end ! !===================================================================== ! function lc_frac(lc,ix,iy,lc_class,water) ! ! Function lc_extent to calculate the proportion of a specific landcover ! type in a square area surrounding (and including) a grid cell, not ! including water in the averaging ! integer, parameter :: r8 = selected_real_kind(12) integer :: lc(720,360) ! landcover data array integer :: ix ! x or longitude index integer :: iy ! y or latitude index integer :: lc_class ! integer value of land cover class to be assessed integer :: water ! integer value of water (not counted) integer :: start_x ! lower longitude index for search grid integer :: end_x ! upper longitude index for search grid integer :: start_y ! lower latitude index for search grid integer :: end_y ! upper latitude index for search grid integer :: index ! index value for counting integer :: sum_core ! central cell with lc values integer :: sum_rad1 ! sum of cells with lc values in 1st ring integer :: sum_rad2 ! sum of cells with lc values in 2nd ring integer :: nrad1 ! sum of all cells in 1st ring integer :: nrad2 ! sum of all cells in 2nd ring real(r8) :: delta_core ! relative change in center cell real(r8) :: delta_rad1 ! relative change in 1st ring real(r8) :: delta_rad2 ! relative change in 2nd ring integer :: indy ! index for latitude real(r8) :: lc_frac ! result ! ! set the search grid dimensions ! start_x = ix-2 end_x = ix+2 start_y = iy-2 end_y = iy+2 ! ! initialize the sums for the values and for the total land grid cell count ! sum_core=0 sum_rad1=0 sum_rad2=0 n_core=0 n_rad1=0 n_rad2=0 ! ! Go through each cell in the 25 cell grid and sum areas ! do i = start_y,end_y do j = start_x,end_x ! ! make sure to wrap around at the date line ! and identify the adjacent grid cell at the poles (over 360 or 180 degrees) ! indx = j indy = i if(indy.lt.1)then indy = 1 - indy indx = indx + 360 endif if(indy.gt.360)then indy = 360 - (indy-361) indx = indx + 360 endif if(indx.lt.1)indx = indx + 720 if(indx.gt.720)indx = indx - 720 ! ! check to see if there is water grid cell ! if(lc(indx,indy).ne.water)then ! ! sum all non water grid cells in area ! ! first evaluate the outer 2 rows and add to totals for outer ring ! if (j.eq.start_x.or.j.eq.end_x)then n_rad2 = n_rad2 + 1 if(lc(indx,indy) .eq. lc_class) sum_rad2 = sum_rad2 + 1 endif ! ! evaluate the middle 2 rows and add to totals for outer and middle rings ! ! outer cross rows ! if(j.eq.(start_x + 1).or.j.eq.(end_x - 1))then if(i.eq.start_y.or.i.eq.end_y)then n_rad2 = n_rad2 + 1 if(lc(indx,indy) .eq. lc_class) sum_rad2 = sum_rad2 + 1 endif ! ! middle cross rows ! if(i.eq.(start_y + 1).or.i.eq.(end_y - 1))then n_rad1 = n_rad1 + 1 if(lc(indx,indy) .eq. lc_class) sum_rad1 = sum_rad1 + 1 endif ! ! center cross row ! if(i.eq.iy)then n_rad1 = n_rad1 + 1 if(lc(indx,indy) .eq. lc_class) sum_rad1 = sum_rad1 + 1 endif endif ! ! evaluate the middle row and add to totals for outer, middle and center cell ! if(j.eq.(ix))then ! ! outer rows ! if(i.eq.start_y.or.i.eq.end_y)then n_rad2 = n_rad2 + 1 if(lc(indx,indy) .eq. lc_class) sum_rad2 = sum_rad2 + 1 endif ! ! middle rows ! if(i.eq.(start_y + 1).or.i.eq.(end_y - 1))then n_rad1 = n_rad1 + 1 if(lc(indx,indy) .eq. lc_class) sum_rad1 = sum_rad1 + 1 endif ! ! center row ! if(i.eq.iy)then if(lc(indx,indy) .eq. lc_class) sum_core = sum_core + 1 endif endif ! endif ! end no value/water if check ! enddo enddo ! ! calculate the fraction of land cover in search area ! delta_core = real(sum_core) ! if(n_rad1.gt.0)then delta_rad1 = real(sum_rad1) / real(n_rad1) else delta_rad1 = 0.0 endif ! if(n_rad2.gt.0)then delta_rad2 = real(sum_rad2) / real(n_rad2) else delta_rad2 = 0.0 endif ! ! calculate the relative change using weights for each section ! lc_frac = 0.5 * delta_core + 0.333 * delta_rad1 + 0.167 * delta_rad2 ! return end