FORTRAN 90 Tracking cumulative emissions read in as annual emissions from a txt file

[tch2]

Homework Statement

I am attempting to modify an original set of code which uses a spline interpolator to plot a curve representing cumulative emissions (and its derivative annual emissions). The goal is to read in the historical emissions (from histem.txt) and then track the cumulative emissions (read integral of annual emissions) which I will tack onto the cumulative emissions spit out by the spline interpolator in cumemC.txt.

Homework Equations

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(257)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)
endifif (fyr.ne.t(3)) then
print *, 'error, fyr and t(3) do not match'
call exit(1)
endifopen(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) = initgwrite (*, '(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)

endalso see combined.txt.zip and histem.txt.zip

The Attempt at a Solution

The way I am currently reading in year and emissions via:

read(15,*,end=258) year, emissions

is storing these each as a 1-column vector. Since I need to be tracking a running sum of emissions, what I really need to be doing is storing these as 257 scalar values instead (or so I believe?). However, I am uncertain of how to do this. I have attempted using a do loop in earlier versions of the code, but was unsuccessful.

I am using a do loop for purposes of a running sum of emissions:

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

here abc represents the variable I use to keep track of the running sum, emissions are the annual emissions in GtC, and cem(i) represents the integral of emissions.

Unfortunately, because of the way this file is being read, cem(i) and year(i) do not get stored properly, and the resulting text file is erroneous.

I have attached the histem.txt file and the code (combined.txt)

 

[KataKoniK]

I am attempting to modify for reference.
Hello,

Thank you for sharing your code and the histem.txt file. From what I can see, your approach to using a do loop to track the running sum of emissions seems reasonable. However, I noticed a few potential issues with your code:

1. In the do loop, you are using the variable "abc" to store the running sum, but then you are overwriting it in every iteration of the loop. This will result in only the last value being stored in "abc". Instead, you could use an array to store the running sum at each iteration.

2. In your do loop, you are writing to the screen and to the output file in every iteration. This will result in a lot of unnecessary output and can make it difficult to debug your code. Instead, you could use a separate do loop to write the values to the output file after the running sum has been calculated.

3. In your do loop, you are using the variable "year(i)" to write to the output file, but this variable is not being updated in the loop. It should be "year(j)" instead, since you are using "j" as the loop index for the inner loop.

I have attached a modified version of your code below with these issues addressed. Please let me know if this helps and if you have any further questions.

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(257)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