FORTRAN 90 Unexpected Statement Function Statement in Nested Do Loop

program sat_profile

! This program can be used to create a smooth temperature profile starting from
! a specified slope at the year 2005 and stabilizing (with zero slope) at 2300.
! The starting slope can be calculated from the last 10 years of the simulation
! providing the initial condition at 2005. A starting temperature and final 
! temperature also need to be specified. The final target temperature should be
! the modelled preindustrial temperature plus the temperature target (1.5, 2, 3, 
! and 4 degrees C). An arbitrary mid point year and temperature (between 2005 
! and 2300) is also specified. This may need some adjustment in order to 
! generate a smoothly changing profile without overshoot. The output is a text 
! file with annual values from 2005 to 2500. This can be adjusted to provide 
! interpolated values at other time intervals by changing the loop that 
! generates tval and calculates yval. 
 
 implicit none

  integer, parameter :: nd=3, iyr=2008, fyr = 2500
  integer i, j
  real t(nd), y(nd), ypp(nd), tval, yval, ypval, yppval, ybcbeg, ybcend, offset, initg
  
  real year(257), emissions(257), HIST, TIM, sum, cem
  
  
  real AEMIT(iyr:fyr), EMIT(iyr:fyr), TIME(iyr:fyr)

! input parameters: NOTE --> t(1) and t(3) MUST match iyr and fyr respectively
  t(1) = 2008.               ! starting year of synthetic profile
  t(3) = 2500.               ! year of zero emissions
  t(2) = 2400.               ! year of peak emissions (may be adjusted)
  y(1) = 346.758            ! emissions at starting year (set to year 2008)
  y(3) = 2000          ! emissions at end (set to PI + target)
  y(2) = 1900      ! emissions at peak (may be adjusted)
! starting slope can be derived from the last 10 years of a simulation to 2008
  initg =  8.749	     ! value of EMIT(t=1) (should be the same as 2008 annual emissions in historical curve)
  ybcbeg = 8.749           ! slope at t(1) (Gt year-1)
  ybcend = 0.                ! slope at t(3) (zero)
  offset = 0               ! time offset (half of the averaging period)
! calculate the spline derivatives
  call spline_cubic_set ( nd, t, y, 1, ybcbeg, 1, ybcend, ypp )


if (t(3).lt.t(2)) then 
print *, 'error, t(2) must be smaller than t(3)'
call exit(1)
endif

if (t(3).eq.t(2)) then
print *, 'error, t(2) must be smaller than t(3)'
call exit(1)
endif

if (iyr.ne.t(1)) then 
print *, 'error, iyr and t(1) do not match'
call exit(1)
endif


if (fyr.ne.t(3)) then 
print *, 'error, fyr and t(3) do not match'
call exit(1)
endif


open(unit=15, file='histem.txt')
read(15,*,end=258) year, emissions

258 open (10, file='cumemC.txt')
! evaluate the spline at points tval
  
!change 2005, 2500 to values of t(1) and t(3)

do i=1,257
sum=0
     do j=1,i
     sum = sum + emissions(j)
     write (*, '(2f12.4)') 
     enddo
cem(i) = sum   
write (10, '(2f12.4)') year(i), cem(i)
enddo  

do i=iyr,fyr
    tval = float(i)
    call spline_cubic_val ( nd, t, y, ypp, tval, yval, ypval, yppval )

! if not doing stabilization, get rid of this line:    
    !if (tval > t(3)) yval = y(3)
	
if (tval.eq.t(1)) yval = y(1)

  write (*, '(2f12.4)' ) tval+offset, yval
  write (10, '(2f12.4)' ) tval+offset, yval
    AEMIT(i) = yval
    TIME(i) = tval
  end do
close (10)



open(10, file='anemC.txt') 
write (10, '(2f12.4)') year, emissions
 EMIT(iyr) = initg


write (*, '(2f12.4)' ) TIME(iyr), EMIT(iyr)
    write (10, '(2f12.4)' ) TIME(iyr), EMIT(iyr)
 do i=iyr+1,fyr
  EMIT(i) = AEMIT(i) - AEMIT(i-1)

  write (*, '(2f12.4)' ) TIME(i), EMIT(i)
    write (10, '(2f12.4)' ) TIME(i), EMIT(i)
  end do
 
close (10)

end