- #1
tch2
- 4
- 0
In the following code, highlighted in bold, I am using a nested DO loop to create a running sum of the yearly emissions ("emissions") read out of histem.txt, to create a cumulative emissions curve.
The code generates an "Unexpected STATEMENT FUNCTION statement for sum in cem(i) = sum, however?
The code generates an "Unexpected STATEMENT FUNCTION statement for sum in cem(i) = sum, however?
Code:
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 modeled 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)
[B]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 [/B]
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