Program readandtransfer
!==============================================================       
!    Purpose: 
!       1. Read annual mean SSI for 850-1850
!       2. Transfer SSI to user-provided spectral grid 
!==============================================================       
       implicit none

integer, parameter :: nyear = 1001 ! number of years in PMIP data file
integer year(nyear)
integer i, j, iy
real*8 tyear

! Variables for input spectral grid
! DO NOT CHANGE
integer, parameter :: n1   = 3780 ! number of spectral bins in the input SSI dataset	   
real*8 flux_in(nyear,n1) ! input SSI (W/m**2/nm) 
real*8 tsi_in(nyear) ! input TSI (W/m**2) 
real*8, dimension (n1) :: sps, spe, spl, spld, wgt

! Variables for output spectral grid (example from ECHAM-5)
! CHANGE AS NECESSARY
integer, parameter :: n2 = 6          ! number of spectral bins in the output 	   

! DO NOT CHANGE
real*8 flux_out(nyear,n2)             ! output SSI (W/m**2) 
real*8 flux_out_spectral(nyear,n2)    ! output SSI (W/m**2/nm) 
real*8, dimension (n2) :: spsm, spem, splm, spdm

!1. Define output grid (nm) example from ECHAM-5

! CHANGE AS NECESSARY
data spsm /185.2,250.0,440.0,690.0,1190.0,2380.0/ ! lower bound (nm) of the output spectral bins
data spem /250.0,440.0,690.0,1190.0,2380.0,4000.0/! upper bound (nm) of the output spectral bins 

do i = 1, n2
spdm(i) = (spem(i)-spsm(i))  ! length of the output spectral bins
enddo


!2. read input data

open(21,file='pmip4_corr_c14.dat',status='old') 

do i = 1,11
read(21,*)
enddo

read(21,'(10(e12.4))') sps 
read(21,'(10(e12.4))') spld 

do i = 1, n1
spe(i) = sps(i)+spld(i)
enddo

do iy = 1,nyear
read(21,*) year(iy),tsi_in(iy) 
read(21,'(10(e12.4))') flux_in(iy,:) 
enddo


do j = 1, n2 

wgt = 0.0

do i = 1, n1 

if (sps(i) > spem(j)) goto 700
if (spe(i) < spsm(j)) goto 700
if ( (sps(i) >= spsm(j)) .and. (spe(i) <= spem(j))) wgt(i) = 1.0
if (sps(i) < spsm(j) .and. spe(i) > spem(j)) &
wgt(i) = 1.0/spld(i)*(spe(i)-spsm(j))
if (sps(i) > spsm(j) .and. spe(i) > spem(j) .and. sps(i) < spem(j)) &
wgt(i) = 1.0/spld(i-1)*(spe(i-1)-spem(j-1))
if (sps(i) < spsm(j) .and. spe(i) > spem(j)) &
wgt(i) = 1.0/spld(i)*(spem(j)-spsm(j))

700 continue 

enddo ! i = 1,n1

do iy = 1, nyear
flux_out(iy,j) = 0.
do i = 1,n1 
flux_out(iy,j) = flux_out(iy,j) + flux_in(iy,i)*wgt(i)*spld(i)
enddo
enddo

enddo  ! j=1,n2

! Correct total incoming solar irradiance if far infrared is missing

do iy = 1, nyear
flux_out(iy,n2) = flux_out(iy,n2) + (tsi_in(iy) - sum(flux_out(iy,:)))
enddo

do iy = 1, nyear
do j = 1,n2
  flux_out_spectral(iy,j) = flux_out(iy,j)/spdm(j)
enddo
enddo

open(21,file='ssi_output.dat',status='unknown')
open(22,file='tsi_output.dat',status='unknown')

do iy = 1, nyear
write(21,'(i4,1x,6(1pe11.4,1x))') year(iy),flux_out(iy,:)
tyear=year(iy)+0.5
write(22,*) tyear,tsi_in(iy)
enddo

END PROGRAM