%GCL
// -----------------------------------------------------------------------------
// Merge La2004 files for the Rupelian Stage (ca. -33.9 ... -28.4 Myr)
// la2004_merge.gcl - created 2005-Dec-26, last updated 2006-Jan-15
// (c) Erik De Sonville, http://users.skynet.be/erik.desonville/
// Script language: disgcl (H. Michels), http://www.dislin.de/
// -----------------------------------------------------------------------------
// Background: http://users.skynet.be/erik.desonville/en/astro.paleo.htm
// Input data is obtained from the La2004 model, see
// http://www.imcce.fr/Equipes/ASD/insola/earth/earth.html
// This script merges and reformats input from various La2004 data files into
// one output file suitable for further processing (i.e. plotting).
// -----------------------------------------------------------------------------
// Input files used:
// - INSOLN.LA2004.BTL.ASC
// file with orbital data, from -51000 to 0 kyr, per kyr
// (generated by IMCCE/ASD, by numerical integration)
// - Output file from INSOLA.EXE, choice 4
// file with mean annual insolation (global),
// from -34000 to -28000 kyr, per kyr
// - Output files from INSOLA.EXE, choice 2
// twelve files with mean daily insolation at 51°N,
// for the mean longitudes 0°, 30°, ... 330°,
// corresponding to conventional dates Mar 21, Apr 21, ... Feb 21,
// from -34000 to -28000 kyr, per kyr
// Note: "D" exponent indicators must be replaced with "E" before running script
// -----------------------------------------------------------------------------
// Output file with merged data runs from -34000 to -28000 kyr, per kyr
// Output data format:
// Year
// column 01-06 (%6i) kiloyear from J2000.0
// Orbital elements
// column 09-16 (%8.6f) eccentricity
// column 18-23 (%6.3f) obliquity (deg)
// column 25-31 (%7.3f) longitude of perihelion (deg)
// Global insolation (W/m^2)
// column 34-40 (%7.3f) mean annual (global)
// Insolation (W/m^2) @ 51°N
// column 42-48 (%7.3f) mean annual (@ 51°N) = average calculated from the
// twelve mean daily insolation files for mean longitudes 0°, 30°, ... 330°
// column 50-56 (%7.3f) mean daily @ mean longitude 270° (winter solstice)
// column 58-64 (%7.3f) mean daily @ mean longitude 0° (vernal equinox)
// column 66-72 (%7.3f) mean daily @ mean longitude 90° (summer solstice)
// column 74-80 (%7.3f) mean daily @ mean longitude 180° (autumnal equinox)
// -----------------------------------------------------------------------------
// output format strings
fmt = '%6i %8.6f %6.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'
fmtavg = 'Mean %8.6f %6.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'
fmtmin = 'Min. %8.6f %6.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'
fmtmax = 'Max. %8.6f %6.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n'
// debug switches
debug = 0
debu9 = 0
// Variables
pi = 3.14159265359
todeg = 180./pi
stat = 0 // return status of functions
yr0 = -34000 // first, last, number of kyears
yr9 = -28000
nyr = yr9-yr0+1
eccy = falloc(nyr) // "n.btl" orbital element arrays
obly = falloc(nyr)
lphy = falloc(nyr)
rglo = falloc(nyr) // "yearly mean" insolation array (global)
ravg = falloc(nyr) // "yearly mean" insolation array (@ 51°N)
r000 = falloc(nyr) // "daily mean" insolation arrays (@ 51°N),
r030 = falloc(nyr) // for mean longitudes 000, 030 etc.
r060 = falloc(nyr) // i.e. for "Mar 21", "Apr 21" etc.
r090 = falloc(nyr)
r120 = falloc(nyr)
r150 = falloc(nyr)
r180 = falloc(nyr)
r210 = falloc(nyr)
r240 = falloc(nyr)
r270 = falloc(nyr)
r300 = falloc(nyr)
r330 = falloc(nyr)
hbtl = -1 // file handles
hglo = -1
hout = -1
h000 = -1
h030 = -1
h060 = -1
h090 = -1
h120 = -1
h150 = -1
h180 = -1
h210 = -1
h240 = -1
h270 = -1
h300 = -1
h330 = -1
kyear = 0. // work variables
eccen = 0.
obliq = 0.
lonph = 0.
insol = 0.
delta = 0.
delt2 = 0.
// Read arrays for eccentricity, obliquity, longitude of perihelion from file
hbtl = fopen('.\data\btl\la2004n.btl.asc','r')
if (hbtl==-1) goto ebtl
// skip the non-relevant part
do i=0,-yr9-1
stat = fscanf(hbtl,'%f%f%f%f',kyear,eccen,obliq,lonph)
if (stat==-1) goto ebtl
end do
if (debu9!=0) print 'skip done'
// read orbital element arrays
do i=0,nyr-1
j=nyr-i-1
stat = fscanf(hbtl,'%f%f%f%f',kyear,eccen,obliq,lonph)
if (stat==-1) goto ebtl
delta = (kyear+34000.)-j
delt2 = delta*delta
if (delt2>=1.e-8) goto ebtl
eccy[j] = eccen
obly[j] = obliq
lphy[j] = lonph
end do
stat = fclose(hbtl)
if (stat==-1) goto ebtl
if (debu9!=0) print 'btl read'
if (debug!=0) print 'first/last eccy',eccy[0],eccy[nyr-1]
if (debug!=0) print 'first/last obly',obly[0],obly[nyr-1]
if (debug!=0) print 'first/last lphy',lphy[0],lphy[nyr-1]
obly = obly * todeg
lphy = lphy * todeg
// Read insolation array for yearly mean (global) from file
hglo = fopen('.\data\nxx\nxx.year.asc','r')
if (hglo==-1) goto eglo
do i=0,nyr-1
stat = fscanf(hglo,'%f%f',kyear,insol)
if (stat==-1) goto eglo
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto eglo
rglo[i] = insol
end do
stat = fclose(hglo)
if (stat==-1) goto eglo
if (debug!=0) print 'first/last rglo',rglo[0],rglo[nyr-1]
// Read insolation array for mean longitude 000 ("Mar 21") from file
h000 = fopen('.\data\n51\n51.ml000.asc','r')
if (h000==-1) goto e000
do i=0,nyr-1
stat = fscanf(h000,'%f%f',kyear,insol)
if (stat==-1) goto e000
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e000
r000[i] = insol
end do
stat = fclose(h000)
if (stat==-1) goto e000
if (debug!=0) print 'first/last r000',r000[0],r000[nyr-1]
// Read insolation array for mean longitude 030 ("Apr 21") from file
h030 = fopen('.\data\n51\n51.ml030.asc','r')
if (h030==-1) goto e030
do i=0,nyr-1
stat = fscanf(h030,'%f%f',kyear,insol)
if (stat==-1) goto e030
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e030
r030[i] = insol
end do
stat = fclose(h030)
if (stat==-1) goto e030
if (debug!=0) print 'first/last r030',r030[0],r030[nyr-1]
// Read insolation array for mean longitude 060 ("May 21") from file
h060 = fopen('.\data\n51\n51.ml060.asc','r')
if (h060==-1) goto e060
do i=0,nyr-1
stat = fscanf(h060,'%f%f',kyear,insol)
if (stat==-1) goto e060
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e060
r060[i] = insol
end do
stat = fclose(h060)
if (stat==-1) goto e060
if (debug!=0) print 'first/last r060',r060[0],r060[nyr-1]
// Read insolation array for mean longitude 090 ("Jun 21") from file
h090 = fopen('.\data\n51\n51.ml090.asc','r')
if (h090==-1) goto e090
do i=0,nyr-1
stat = fscanf(h090,'%f%f',kyear,insol)
if (stat==-1) goto e090
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e090
r090[i] = insol
end do
stat = fclose(h090)
if (stat==-1) goto e090
if (debug!=0) print 'first/last r090',r090[0],r090[nyr-1]
// Read insolation array for mean longitude 120 ("Jul 21") from file
h120 = fopen('.\data\n51\n51.ml120.asc','r')
if (h120==-1) goto e120
do i=0,nyr-1
stat = fscanf(h120,'%f%f',kyear,insol)
if (stat==-1) goto e120
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e120
r120[i] = insol
end do
stat = fclose(h120)
if (stat==-1) goto e120
if (debug!=0) print 'first/last r120',r120[0],r120[nyr-1]
// Read insolation array for mean longitude 150 ("Aug 21") from file
h150 = fopen('.\data\n51\n51.ml150.asc','r')
if (h150==-1) goto e150
do i=0,nyr-1
stat = fscanf(h150,'%f%f',kyear,insol)
if (stat==-1) goto e150
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e150
r150[i] = insol
end do
stat = fclose(h150)
if (stat==-1) goto e150
if (debug!=0) print 'first/last r150',r150[0],r150[nyr-1]
// Read insolation array for mean longitude 180 ("Sep 21") from file
h180 = fopen('.\data\n51\n51.ml180.asc','r')
if (h180==-1) goto e180
do i=0,nyr-1
stat = fscanf(h180,'%f%f',kyear,insol)
if (stat==-1) goto e180
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e180
r180[i] = insol
end do
stat = fclose(h180)
if (stat==-1) goto e180
if (debug!=0) print 'first/last r180',r180[0],r180[nyr-1]
// Read insolation array for mean longitude 210 ("Oct 21") from file
h210 = fopen('.\data\n51\n51.ml210.asc','r')
if (h210==-1) goto e210
do i=0,nyr-1
stat = fscanf(h210,'%f%f',kyear,insol)
if (stat==-1) goto e210
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e210
r210[i] = insol
end do
stat = fclose(h210)
if (stat==-1) goto e210
if (debug!=0) print 'first/last r210',r210[0],r210[nyr-1]
// Read insolation array for mean longitude 240 ("Nov 21") from file
h240 = fopen('.\data\n51\n51.ml240.asc','r')
if (h240==-1) goto e240
do i=0,nyr-1
stat = fscanf(h240,'%f%f',kyear,insol)
if (stat==-1) goto e240
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e240
r240[i] = insol
end do
stat = fclose(h240)
if (stat==-1) goto e240
if (debug!=0) print 'first/last r240',r240[0],r240[nyr-1]
// Read insolation array for mean longitude 270 ("Dec 21") from file
h270 = fopen('.\data\n51\n51.ml270.asc','r')
if (h270==-1) goto e270
do i=0,nyr-1
stat = fscanf(h270,'%f%f',kyear,insol)
if (stat==-1) goto e270
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e270
r270[i] = insol
end do
stat = fclose(h270)
if (stat==-1) goto e270
if (debug!=0) print 'first/last r270',r270[0],r270[nyr-1]
// Read insolation array for mean longitude 300 ("Jan 21") from file
h300 = fopen('.\data\n51\n51.ml300.asc','r')
if (h300==-1) goto e300
do i=0,nyr-1
stat = fscanf(h300,'%f%f',kyear,insol)
if (stat==-1) goto e300
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e300
r300[i] = insol
end do
stat = fclose(h300)
if (stat==-1) goto e300
if (debug!=0) print 'first/last r300',r300[0],r300[nyr-1]
// Read insolation array for mean longitude 330 ("Feb 21") from file
h330 = fopen('.\data\n51\n51.ml330.asc','r')
if (h330==-1) goto e330
do i=0,nyr-1
stat = fscanf(h330,'%f%f',kyear,insol)
if (stat==-1) goto e330
delta = (kyear+34000.)-i
delt2 = delta*delta
if (delt2>=1.e-8) goto e3300
r330[i] = insol
end do
stat = fclose(h330)
if (stat==-1) goto e330
if (debug!=0) print 'first/last r330',r330[0],r330[nyr-1]
// Calculate mean yearly insolation (@ 51°N)
ravg = (r000+r030+r060+r090+r120+r150+r180+r210+r240+r270+r300+r330)/12.
if (debug!=0) print 'first/last ravg',ravg[0],ravg[nyr-1]
// Calculate mean,min,max values for eccy,obly,rglo,ravg,r270,r000,r090,r180
avgecc = 0.
minecc = eccy[0]
maxecc = eccy[0]
avgobl = 0.
minobl = obly[0]
maxobl = obly[0]
avgglo = 0.
minglo = rglo[0]
maxglo = rglo[0]
avgavg = 0.
minavg = ravg[0]
maxavg = ravg[0]
avg270 = 0.
min270 = r270[0]
max270 = r270[0]
avg000 = 0.
min000 = r000[0]
max000 = r000[0]
avg090 = 0.
min090 = r090[0]
max090 = r090[0]
avg180 = 0.
min180 = r180[0]
max180 = r180[0]
do i=0,nyr-1
eccen = eccy[i]
avgecc = avgecc+eccen
if (eccenmaxecc) maxecc=eccen
obliq = obly[i]
avgobl = avgobl+obliq
if (obliqmaxobl) maxobl=obliq
insol = rglo[i]
avgglo = avgglo+insol
if (insolmaxglo) maxglo=insol
insol = ravg[i]
avgavg = avgavg+insol
if (insolmaxavg) maxavg=insol
insol = r270[i]
avg270 = avg270+insol
if (insolmax270) max270=insol
insol = r000[i]
avg000 = avg000+insol
if (insolmax000) max000=insol
insol = r090[i]
avg090 = avg090+insol
if (insolmax090) max090=insol
insol = r180[i]
avg180 = avg180+insol
if (insolmax180) max180=insol
end do
avgecc = avgecc/nyr
avgobl = avgobl/nyr
avgglo = avgglo/nyr
avgavg = avgavg/nyr
avg270 = avg270/nyr
avg000 = avg000/nyr
avg090 = avg090/nyr
avg180 = avg180/nyr
// Write output file
hout = fopen('.\data\n51out\insola51.asc','w')
if (hout==-1) goto eout
stat = fprintf(hout,fmtavg,avgecc,avgobl,avgglo,avgavg,avg270,avg000,avg090,avg180)
if (stat==-1) goto eout
stat = fprintf(hout,fmtmin,minecc,minobl,minglo,minavg,min270,min000,min090,min180)
if (stat==-1) goto eout
stat = fprintf(hout,fmtmax,maxecc,maxobl,maxglo,maxavg,max270,max000,max090,max180)
if (stat==-1) goto eout
stat = fprintf(hout,'// #\n') // # for end-of-header
if (stat==-1) goto eout
do i=0,nyr-1
j=yr0+i
stat = fprintf(hout,fmt,j,eccy[i],obly[i],lphy[i],rglo[i],ravg[i],r270[i],r000[i],r090[i],r180[i])
if (stat==-1) goto eout
end do
stat = fclose(hout)
if (stat==-1) goto eout
// Exits
exitok:
print 'Done okay'
goto endpgm
ebtl:
print 'Error btl'
goto endpgm
eglo:
print 'Error glo'
goto endpgm
e000:
print 'Error 000'
goto endpgm
e030:
print 'Error 030'
goto endpgm
e060:
print 'Error 060'
goto endpgm
e090:
print 'Error 090'
goto endpgm
e120:
print 'Error 120'
goto endpgm
e150:
print 'Error 150'
goto endpgm
e180:
print 'Error 180'
goto endpgm
e210:
print 'Error 210'
goto endpgm
e240:
print 'Error 240'
goto endpgm
e270:
print 'Error 270'
goto endpgm
e300:
print 'Error 300'
goto endpgm
e330:
print 'Error 330'
goto endpgm
eout:
print 'Error out'
goto endpgm
endpgm: