program future ! Version 1.0 ! ! This is a program to create land cover conditions for future ! SRES scenarios based on existing present day information and ! IMAGE 2.2 SRES projections ! ! 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_ag = 1 !ag land cover type to be evaluated -- IMAGE integer, parameter :: lc_grass = 2 !grass land cover type to be evaluated -- IMAGE integer, parameter :: water = 0 !water value - cells ignored in calculations -- IMAGE ! ! 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 ! integer :: svegin1(259200) ! input IMAGE 1D landcover array time 1 integer :: svegin2(259200) ! input IMAGE 1D landcover array time 2 real :: sresx(259200) ! input IMAGE x coordinates to read data file real :: sresy(259200) ! input IMAGE y coordinates to read data file integer :: sres_lc_in(nlon,nlat) ! input IMAGE landcover time 1 integer :: sres_lc_out(nlon,nlat) ! input IMAGE landcover time 2 real(r8) :: gradient(nlon,nlat) ! input slope information (gradient units) ! 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) :: in_image_ag ! initial IMAGE input agriculture real(r8) :: out_image_ag ! initial IMAGE output agriculture real(r8) :: in_ag ! input agriculture real(r8) :: out_ag ! output agriculture real(r8) :: slope_scale ! ag scaling factor based on slope real(r8) :: max_veg ! total potential vegetation for ag calculation real(r8) :: max_ag ! max ag allowed in cell based on slope real(r8) :: scale ! scale veriable to adjust IMAGE/input difference real(r8) :: ag_im_diff ! change in Image ag area real(r8) :: ag_change ! scaled change in Image ag area real(r8) :: ag_test ! initial IMAGE agriculture estimate real(r8) :: ag_output ! Final agriculture output ! 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_test ! preliminary input IMAGE grazing real(r8) :: out_grazing_test ! preliminary output IMAGE grazing real(r8) :: in_grazing ! input IMAGE grazing real(r8) :: out_grazing ! output IMAGE grazing real(r8) :: clm_graz ! CLM grass with IMAGE 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) :: pct_grazing_in ! percent of grass area input 1990 (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 :: npts integer :: iyr integer :: iskip integer :: nx integer :: ny integer :: i integer :: j integer :: k integer :: m integer :: ix integer :: iy integer :: scen ! real :: rdummy real(r8) :: lat real(r8) :: 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 :: fname1_2 ! A2, B1 etc - 2 letter character*20 :: fname2_2 ! A2, B1 etc - 2 letter character*80 :: fname_in_2 ! A2, B1 etc - 2 letter character*80 :: fname_out_2 ! A2, B1 etc - 2 letter character*20 :: fname1_3 ! A1B - 3 letter character*20 :: fname2_3 ! A1B - 3 letter character*80 :: fname_in_3 ! A1B - 3 letter character*80 :: fname_out_3 ! A1B - 3 letter integer itime, year_in, year_out ! real(r8) :: lc_frac ! function ! ! Start loop to do all files to 2100 ! ! write(6,*) ' enter the scenario to simulate' write(6,*) ' 1 = A2' write(6,*) ' 2 = B1' write(6,*) ' 3 = B2' write(6,*) ' 4 = A1B' read(5,*)scen ! if(scen.eq.1)then fname1_2 = 'output/pft-' fname2_2 = '-A2.5x.5' elseif(scen.eq.2)then fname1_2 = 'output/pft-' fname2_2 = '-B1.5x.5' elseif(scen.eq.3)then fname1_2 = 'output/pft-' fname2_2 = '-B2.5x.5' elseif(scen.eq.4)then fname1_3 = 'output/pft-' fname2_3 = '-A1B.5x.5' endif ! ! go through the time loop to process each decade year from 2000 to 2100 ! do itime = 1,11 ! ! Open the files for this time step ! ! write file names ! year_in = 2000 + (itime - 2) * 10 year_out = 2000 + (itime - 1) * 10 if(scen.lt.4)then write(fname_in_2,'(a11,i4,a8)')fname1_2,year_in,fname2_2 ! A2, B1 etc (2 letters) write(fname_out_2,'(a11,i4,a8)')fname1_2,year_out,fname2_2 ! A2, B1 etc (2 letters) else write(fname_in_3,'(a11,i4,a9)')fname1_3,year_in,fname2_3 ! A1B (3 letters) write(fname_out_3,'(a11,i4,a9)')fname1_3,year_out,fname2_3 ! A1B (3 letters) endif write(6,*)'processing ',year_out ! ! NOTE the order of the files going back in time for HYDE ! open(10, file ='input/worldnc.grd') ! IMAGE land mask info ! if(scen.eq.1)then open(11, file ='input/landcovera2.dat') ! IMAGE landuse info elseif(scen.eq.2)then open(11, file ='input/landcoverb1.dat') ! IMAGE landuse info elseif(scen.eq.3)then open(11, file ='input/landcoverb2.dat') ! IMAGE landuse info elseif(scen.eq.4)then open(11, file ='input/landcoverA1B.dat') ! IMAGE landuse info endif ! open(12, file ='input/gradient05.asc') ! slope file in gradients if(itime .eq. 1)then open(13, file = 'output/pft-1990-hist.5x.5') ! CLM PFT input file time 1 else if(scen.lt.4)then open(13, file = fname_in_2) ! CLM PFT input file time 1 else open(13, file = fname_in_3) ! CLM PFT input file time 1 endif endif open(14, file = 'output/pft-pot.5x.5') ! CLM PFT potential file if(scen.lt.4)then open(15, file = fname_out_2) ! CLM output file time 2 else open(15, file = fname_out_3) ! CLM output file time 2 endif ! ! read in SRES data SET iyr to read the appropriate year ! iyr = itime + 2 iskip=2*iyr+1 ! ! initialize SRES veg array to ocean value (0) ! this array is from 90,-180 across and then ! down to -90,180 ! do i=1,nlat do j=1,nlon sres_lc_in(j,i) = -9999 sres_lc_out(j,i) = -9999 enddo enddo ! ! read the IMAGE 2.0 SRES data ! skip headers and read lat long locations ! do i=1,15 read(10,*)dummy enddo ! do i=1,259200 read(10,*,end=2)sresx(i),sresy(i) npts=i enddo 2 write(6,*)'done read sres grid read: npts = ',npts ! ! skip lines to the year for which you want to read the lu data ! do i=1,iskip read(11, *) dummy enddo read(11,*) (svegin1(k),k=1,npts) read(11, *) dummy read(11,*) (svegin2(k),k=1,npts) ! write(6,*)'read sres vegetation and lu data' ! ! overwrite the locations where there are real values ! do i=1,npts nx=int((sresx(i)+180)/0.5)+1 ny=int((89.5-sresy(i))/0.5)+1 ! change to 90 if offset by 0.5 deg ! ! create the SRES agriculture grid ! sres_lc_in(nx,ny)=svegin1(i) sres_lc_out(nx,ny)=svegin2(i) ! enddo ! ! Read GIS slope data ! do i=1,6 read(12,*)dummy enddo ! do i = 1,360 read(12,*) (gradient(j,i),j=1,720) enddo ! ! 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, & pct_grazing_in,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 ! ! Get IMAGE agriculture information ! ! get slope information and calculate slope factor ! if (gradient(i,j) .lt. 0.0)gradient(i,j) = 0.0 if (gradient(i,j) .ge. 15._r8) then slope_scale = 0.5_r8 else slope_scale = 1._r8 - (gradient (i,j) / 15._r8 * 0.5_r8) endif ! ! calculate the fraction of grazing land for each time period (hyde data) ! assume grazing can occupy no more than 60 percent i.e. scale to max 75% ! in_image_ag = lc_frac(sres_lc_in,i,j,lc_ag,water) * 100._r8 out_image_ag = lc_frac(sres_lc_out,i,j,lc_ag,water) * 100._r8 ! ! Apply a weighting scheme for 2000, 2010 and 2020 to make a smoother transition ! if(itime .eq.1)then in_ag = pft_in(15) out_ag = 0.75 * pft_in(15) + 0.25 * out_image_ag elseif(itime .eq.2)then in_ag = pft_in(15) out_ag = 0.5 * pft_in(15) + 0.5 * out_image_ag elseif(itime .eq.3)then in_ag = pft_in(15) out_ag = 0.25 * pft_in(15) + 0.75 * out_image_ag else in_ag = in_image_ag out_ag = out_image_ag endif ! ! scale ag in and ag out to reflect the effects of slope and accounting for bare ground ! max_veg = 100._r8 - pft_in(0) max_ag = max_veg * slope_scale ! ! Scale the values so that they don't completely scale to ! 100 percent. Scaling is dependent on the size change in ! the Image data. The larger the one time change the greater ! greater the difference allowed. This is to accomodate e.g. the ! Amazon where there is 0 ag today but where high intensity ag ! can be sustained if required (80-90%). If decreasing the change is ! scaled by 50% so that an area that presently has lots of Ag does ! not reduce to 0% ag at the end e.g. China ! ! calculate scaling factor and the area difference based on IMAGE ! if(in_ag.lt. 0.01)then scale = 1.0 else if(in_image_ag.lt. 0.01 .and. out_image_ag .lt. 0.01)then scale = 0.0 else scale = pft_in(15)/in_ag endif endif if(scale .gt. 1.0) scale = 1.0 ag_im_diff = (out_ag - in_ag) * scale ! ! check for gain or loss ! if (ag_im_diff .ge. 0.0) then ! ! gain in agricultural area from IMAGE data ! funcion results in a max value of 100% if a grid cell instantaneously ! goes from 0 to 100%. Otherwise max values by increments get to about 85% ! ag_change = ((100._r8 - pft_in(15))/100._r8)**0.25 * ag_im_diff * slope_scale ag_test = pft_in(15) + ag_change ! else ! ! loss of ag in IMAGE -- reduce ag so that IMAGE = 0 does not remove all AG ! -- remove only 50% of the IMAGE change at a time which results in a ! remainder of about 35% Ag when Image shows no agriculture ! ag_test = pft_in(15) + (ag_im_diff * 0.5) ! endif ! ! test to see if more than max possible or less than 0.0 ! if (ag_test .gt. max_ag) ag_test = max_ag if (ag_test .lt. 0.0) ag_test = 0._r8 ! ag_output = ag_test ! ! Get IMAGE grazing information ! ! calculate the fraction of grazing land for each time period (hyde data) ! assume grazing can occupy no more than 60 percent i.e. scale to max 75% ! in_grazing_test = lc_frac(sres_lc_in,i,j,lc_grass,water) * 75.0 out_grazing_test = lc_frac(sres_lc_out,i,j,lc_grass,water) * 75.0 ! ! Initialize tree and grass 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) ! ! Modify grazing for the first three decade values to smoothen the transition ! from historical to IMAGE 2.2 data ! if(itime.eq.1)then in_grazing = 0.25_r8 * in_grazing_test + & 0.75_r8 * pct_grazing_in / 100._r8 * grass_in out_grazing = 0.25_r8 * out_grazing_test + & 0.75_r8 * pct_grazing_in / 100._r8 * grass_in elseif(itime.eq.2)then in_grazing = 0.5_r8 * in_grazing_test + & 0.5_r8 * pct_grazing_in / 100._r8 * grass_in out_grazing = 0.5_r8 * out_grazing_test + & 0.5_r8 * pct_grazing_in / 100._r8 * grass_in elseif(itime.eq.3)then in_grazing = 0.75_r8 * in_grazing_test + & 0.25_r8 * pct_grazing_in / 100._r8 * grass_in out_grazing = 0.75_r8 * out_grazing_test + & 0.25_r8 * pct_grazing_in / 100._r8 * grass_in else in_grazing = in_grazing_test out_grazing = out_grazing_test endif ! ! Only process grid cells where there is some ag or grazing information ! if (pft_in(15).gt.0.0 .or. ag_output.gt.0.0 .or. & in_grazing.gt.0.0.or.out_grazing.gt.0.0)then ! ! Set AGRICULTURE and calculate ag stats ! pft_out(15) = ag_output ! ! get the change in agriculture ! delta_ag = pft_out(15) - 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 ! ! Write testing stats ! 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._r8 - pft_out(0) - pft_out(15) sum_pft = 0._r8 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 ! ! end of individual processes, close time loop ! close (10) close (11) close (12) close (13) close (14) close (15) 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