program fill  !  Version 1.0
!
! This is a program takes two known CLM ASCII files and a intermediate RF ag file

! it adjusts the grass landand agriculture areas.  Years to process are
! 1870, 1890, 1910 and 1930
!
! by Johannes Feddema
!
  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
!
!  dimension PFT arrays to use
!
  real(r8) :: pft_1(0:numpft)      ! input PFT data 
  real(r8) :: pft_2(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 
!
  real(r8) :: landmask                      ! NCAR land mask value

  integer :: time_1                        ! time of earlier bounding condition

  integer :: time_2                        ! time interpolated to

  integer :: time_3                        ! time of later bounding condition

  real(r8) :: time_span                     ! time span between time 1 and time 2

  real(r8) :: time_gap                      ! time span between time 1 and out time

  real(r8) :: time_ratio                    ! percent of time gap to out year

!
  real(r8) :: veg_1                         ! total time 1 vegetation (excluding ag and bare)

  real(r8) :: veg_2                         ! total time 2 vegetation (excluding ag and bare)

  real(r8) :: veg_out                       ! total output vegetation (excluding ag and bare)

!
  real(r8) :: pft_1_pct                     ! pft percent of veg at time 1 

  real(r8) :: pft_2_pct                     ! pft percent of veg at time 2 

!
  integer :: i
  integer :: j
  integer :: k
  integer :: m
  integer :: itime
!
  real :: rdummy

  real(r8) :: lat

  real(r8) :: lon

  real(r8) :: sum_1

  real(r8) :: sum_2

  real(r8) :: sum_out

  real(r8) :: sum_pft

  real(r8) :: sumtest_pft

!
  real(r8) :: tree_out                 ! percent tree output
  real(r8) :: grass_out                ! percent grass output 
  real(r8) :: grazing                  ! grazing 
  real(r8) :: graz                     ! temp grazing value 
  real(r8) :: out_grazing              ! temp grazing value 
!  
  real(r8) :: pct_tree                 ! percent tree (tree = sum pft 1-8) 
  real(r8) :: pct_grass                ! percent grass (grass = sum pft 9-14) 
  real(r8) :: pct_grazing_1            ! percent of grass area used for grazing (grass = sum pft 9-14) 
  real(r8) :: pct_grazing_2            ! percent of grass area used for grazing (grass = sum pft 9-14) 
  real(r8) :: pct_grazing              ! percent of grass area used for grazing (grass = sum pft 9-14) 
  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 
!
  character*7 :: dummy
!
  Character*80 :: file10
  Character*80 :: file11
  Character*80 :: file12
  Character*80 :: file13
!
  character*40 :: rf1
  character*40 :: rf2
  character*40 :: out1
  character*40 :: out2
!
! Open the files for this time step
!
! NOTE the order of the files going back in time for HYDE
!
  rf1 = 'input/glcrop'
  rf2 = '_0_5.asc'
  out1 = 'output/pft-'
  out2 = '-hist.5x.5'
!
 do itime = 1, 4
  if (itime .eq. 1) then
   time_1 = 1850
   time_2 = 1870
   time_3 = 1900
  elseif(itime .eq. 2) then
   time_1 = 1870
   time_2 = 1890
   time_3 = 1900

  elseif(itime .eq. 3) then
   time_1 = 1900
   time_2 = 1910
   time_3 = 1950

  elseif(itime .eq. 4) then
   time_1 = 1910
   time_2 = 1930
   time_3 = 1950
  endif

!
  write(file10, '(a11,i4,a10)') out1,time_1,out2
  write(file11, '(a11,i4,a10)') out1,time_3,out2
  write(file12, '(a12,i4,a8)') rf1,time_2,rf2
  write(file13, '(a11,i4,a10)') out1,time_2,out2
!
  open(10, file = file10)          ! CLM PFT input file time 1 

  open(11, file = file11)          ! CLM input file time 3

  open(12, file = file12)          ! cropland input file for given year
  open(13, file = file13)          ! CLM PFT output file

!
  time_span = real(time_3 - time_1)

  time_gap = real(time_2 - time_1)

  time_ratio = time_gap / time_span

!
! read over the GIS header info 
!
  do i=1,6
    read(12,*)dummy
  enddo
!
! read the ag data
!
  do i = 1,nlat
    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 (10,*) lat,lon,landmask,(pft_1(m),m=0,numpft),dummy,dummy,  &
                  pct_grazing_1,dummy,dummy,dummy,dummy,dummy,dummy
     read (11,*) dummy,dummy,landmask,(pft_2(m),m=0,numpft),dummy,dummy,  &
                  pct_grazing_2,dummy,dummy,dummy,dummy,dummy,dummy
!

!   initialize 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 mask)

!

      if(rf_ag(i,j) .lt. 0.0) rf_ag(i,j) = 0.0

!
!  Set bare ground fraction
!
      pft_out(0) = pft_1(0)

!
!  Set AGRICULTURE based on R&F data 
!
      pft_out(15) = rf_ag(i,j) * 100.0
!	    

!  determine the vegetated areas for each period

!

      veg_1 = 100.0 - pft_1(0) - pft_1(15)
      veg_2 = 100.0 - pft_2(0) - pft_2(15) 

      veg_out = 100.0 - pft_out(0) - pft_out(15)

      if(veg_out .lt. 0.0)then
        pft_out(15) = 100.0 - pft_out(0)
        veg_out = pft_out(15)
      endif

!

!  correct the vegetation fraction for each pft

!

      do k = 1,14

!

!  calculate pft contribution to vegetated area at start

!

        if(veg_1 .lt. 0.01) then

          pft_1_pct = 0.0

        else

          pft_1_pct = pft_1(k) / veg_1

        endif

!

!  calculate pft contribution to grass area at end

!

        if(veg_2 .lt. 0.01) then

          pft_2_pct = 0.0

        else

          pft_2_pct = pft_2(k) / veg_2

        endif

!

!  calculate pft at out time

!

        if(veg_1 .lt. 0.001)then
!
!  can't use percent change because of 0 condition on one end of interval
!
          pft_out(k) = pft_2_pct * veg_out
        elseif(veg_2 .lt. 0.001)then
!
!  can't use percent change because of 0 condition on one end of interval
!
          pft_out(k) = pft_1_pct * veg_out
        else
         pft_out(k) = (pft_1_pct + (pft_2_pct - pft_1_pct) * time_ratio) * veg_out

        endif
!

      enddo
!
!   calculate the percent of grass area grazed assume linear change
!
      pct_grazing = pct_grazing_1 + (pct_grazing_2 - pct_grazing_1) * time_ratio  

!

!  DONE write out test info

!

     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
!
!  get grass and tree info
!
     tree_out = pft_out(1) + pft_out(2) + pft_out(3) + pft_out(4) + &
                pft_out(5) + pft_out(6) + pft_out(7) + pft_out(8)
     grass_out = pft_out(9) + pft_out(10) + pft_out(11) + &
                 pft_out(12) + pft_out(13) + pft_out(14)
!
!  initialize grazing values to 0._r8
!
     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
!
     out_grazing = pct_grazing / 100._r8 * grass_out
!
!  only work on land grid cells
!
     if (landmask == 100.) then
!
      grazing = out_grazing
!
!  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 (13,*) 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
     pct_grazing = 0._r8
!
    enddo
  enddo
!
!  end processing - close files and go to next case
!
  close (10)
  close (11)
  close (12)
  close (13)
!
 enddo
!
stop
end