! tynsc20unpack_pik.f90
! f90 -o tynsc20unpack_pik tynsc20unpack_pik.f90
! written by Tim Mitchell on 28.06.02
! this program constructs a specification file for use in TYNsc20SeqGrim
! Originally written to construct a specification file for ATEAM data-set
! This file has been modified by Soenke Zaehle (PIK) on 30.07.02 for use with XLF compiler
! (PIK). changes are: 
! - removed "," after read and write statements
! - replaced SCRATCH statements with REPLACE statements
! plus some others (noticed by TDM)
! subsequently modified again by TDM to suit the TYN SC 2.0 data-set

program TynSC20UnPack

!*******************************************************************************
! use statements

implicit none

!*******************************************************************************
! specifications

real, allocatable, dimension (:,:)			:: Constants
real, allocatable, dimension (:)			:: CommKay
real, dimension (4)					:: RefBounds

integer, allocatable, dimension (:,:)			:: CommOp
integer, allocatable, dimension (:) 			:: Variable

character (len=80), allocatable, dimension (:,:)	:: CommFile, CommInfo
character (len=40), allocatable, dimension (:,:) 	:: FilePatt,FileGloT,File2085
character (len=10), allocatable, dimension (:) 		:: Emissions,Model

character (len=04), dimension (5) :: Suffix = (/".cld",".dtr",".pre",".tmp",".vap"/)

real, parameter :: MissVal = -999.0

integer :: ReadStatus, AllocStat
integer :: NYear,NVariable,NModel,NEmissions,NBound,NElement,NComm
integer :: XYear,XVariable,XModel,XEmissions,XBound,XElement,XComm
integer :: YearAD0,YearAD1
integer :: QModel,QEmissions
integer :: NMonth,NBox,NExe,NWye,NArray

character (len=20) :: OpsFile
character (len=12) :: Date,Time
character (len=04) :: YearAD0Txt,YearAD1Txt

!*******************************************************************************
! main commands

call Initialise
call GetScenarios
call Selections
print*, "  > Preparing specification..."
call CalcOutputs
print*, "  > Writing .ops file..."
call WriteOpsFile
print*, "  > Now unpacking..."
call system ('tynsc20seqgrim_pik')
call Finalise

!*******************************************************************************
! subroutines

contains

!*******************************************************************************
! set initial things

subroutine Initialise

print*
print*, "  > ******** Program to unpack a TYN SC 2.0 scenario ********"
print*, "  >       Tim Mitchell, www.tyndall.ac.uk, January 2003"
print*, "  >          modified by Soenke Zaehle (PIK) for XLF"
print*

end subroutine Initialise

!*******************************************************************************
! load scenario details

subroutine GetScenarios

open  (1,file="scenarios.txt",status="old",action="read")

read  (1,"(2i6)") NModel, NEmissions

allocate (Emissions(NEmissions), &
	  Model    (NModel),         &
	  FilePatt (NModel,NEmissions), &
	  FileGloT (NModel,NEmissions), &
	  Constants(NModel,NEmissions), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GetScenarios: Allocation failure #####"

do XEmissions = 1, NEmissions
  read (1,"(a10)") Emissions(XEmissions)
end do

do XModel = 1, NModel
  read (1,"(a10)") Model(XModel)
  
  do XEmissions = 1, NEmissions
    read (1,"(a40)") FilePatt(XModel,XEmissions)
    read (1,"(a40)") FileGloT(XModel,XEmissions)
    read (1,"(f6.3)") Constants(XModel,XEmissions)
  end do
end do

close (1)

end subroutine GetScenarios

!*******************************************************************************
! make selections

subroutine Selections

print*, "  > Specify the details of the 21st century file to construct."

print*, "  > Enter the first, last years AD (in range 2001-2100): "
do
	read (*,*,iostat=ReadStatus) YearAD0, YearAD1
	if (YearAD1.LT.YearAD0.OR.YearAD0.LT.2001.OR.YearAD1.GT.2100) then
	  ReadStatus = 1
	  print*, "  > Unacceptable range. Try again."
	end if
	if (ReadStatus.LE.0) exit
end do
NYear  = YearAD1 - YearAD0 + 1

print*, "  > Enter the number of climate variables required (1-5): "
do
	read (*,*,iostat=ReadStatus) NVariable
	if (NVariable.LT.1.OR.NVariable.GT.5) then
	  ReadStatus = 1
	  print*, "  > Unacceptable number. Try again."
	end if
	if (ReadStatus.LE.0) exit
end do

allocate (Variable(NVariable), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: UnPack: Allocation failure #####"

print "(a,i1,a)", "   > Specify the variables by listing ", NVariable, " integer(s): "
print*, "  >   (cld=1,dtr=2,pre=3,tmp=4,vap=5)"
do
	read (*,*,iostat=ReadStatus) (Variable(XVariable),XVariable=1,NVariable)
	do XVariable = 1, NVariable
	  if (Variable(XVariable).LT.1.OR.Variable(XVariable).GT.5) ReadStatus = 1
	end do
	if (ReadStatus.GT.0) print*, "  > Unacceptable sequence. Try again."
	if (ReadStatus.LE.0) exit
end do

print*, "  > Select the GCM to use (0=list): "
do
	read (*,*,iostat=ReadStatus) QModel
	if (QModel.EQ.0) then
	  do XModel = 1, NModel
	    print "(a,i2,2a)", "   >   ", XModel, ": ", trim(Model(XModel))
	  end do
	else if (QModel.LT.0.OR.QModel.GT.NModel) then
	  ReadStatus = 1
	end if
	if (ReadStatus.GT.0) print*, "  > Out of range. Try again."
	if (ReadStatus.LE.0.AND.QModel.NE.0) exit
end do

print*, "  > Select the SRES emissions scenario to use (0=list): "
do
	read (*,*,iostat=ReadStatus) QEmissions
	if (QEmissions.EQ.0) then
	  do XEmissions = 1, NEmissions
	    if (FilePatt(QModel,XEmissions).NE."none") &
	    		print "(a,i2,2a)", "   >   ", XEmissions, ": ", trim(Emissions(XEmissions))
	  end do
	else if (QEmissions.LT.0.OR.QEmissions.GT.NEmissions) then
	  ReadStatus = 1
	else if (FilePatt(QModel,QEmissions).EQ."none") then
	  ReadStatus = 1
	end if
	if (ReadStatus.GT.0) print*, "  > Out of range. Try again."
	if (ReadStatus.LE.0.AND.QEmissions.NE.0) exit
end do

end subroutine Selections

!*******************************************************************************
! calculate the outputs to place in the .ops file

subroutine CalcOutputs

NMonth = 12 ; NBox = 67420 ; NExe = 720 ; NWye = 360
NArray = 2 ; NComm = 15 ; NElement = 5
RefBounds = (/-180.0,180.0,-90.0,90.0/)

open  (98,file='trash',status="replace",action="readwrite",iostat=ReadStatus)
write (98,"(2i4)") YearAD0,YearAD1
rewind(98)
read  (98,"(2a4)") YearAD0Txt,YearAD1Txt
close (98)

allocate (CommOp  (NComm,NElement), &
	  CommFile(NComm,NVariable), &
	  CommInfo(NComm,NVariable), &
	  CommKay (NComm), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcOutputs: Allocation failure #####"
CommOp   = MissVal ; CommKay  = MissVal ; CommFile = "" ; CommInfo = ""

do XVariable = 1, NVariable
  CommFile  (2,XVariable) = trim(FilePatt(QModel,QEmissions)) // Suffix(Variable(XVariable))
  CommFile  (5,XVariable) = FileGloT (QModel,QEmissions)
  CommFile  (7,XVariable) = "iavar.1901-2000" // Suffix(Variable(XVariable))
  CommFile  (9,XVariable) = "obs.clim6190.globe.lan" // Suffix(Variable(XVariable))
  CommFile (11,XVariable) = "min" // Suffix(Variable(XVariable)) // ".ann"
  CommFile (13,XVariable) = "max" // Suffix(Variable(XVariable)) // ".ann"
  CommFile (15,XVariable) = trim(Emissions(QEmissions)) // "." // trim(Model(QModel)) // "." // &
  				YearAD0Txt // "-" // YearAD1Txt // Suffix(Variable(XVariable))
  CommInfo (15,XVariable) = "SRES=" // trim(Emissions(QEmissions)) // " GCM=" // trim(Model(QModel)) // &
  				" Period=" // YearAD0Txt // "-" // YearAD1Txt // &
  				" Variable=" // Suffix(Variable(XVariable))
end do

CommOp(1,1)=-1 ; CommOp(1,2)=1
CommOp(2,1)=1 ; CommOp(2,2)=1				! load 2080s anomalies
CommOp(3,1)=30 ; CommOp(3,2)=1 ; CommOp(3,3)=1 ; CommOp(3,4)=1 ; CommKay(3)=Constants(QModel,QEmissions)
CommOp(4,1)=51 ; CommOp(4,2)=1 ; CommOp(4,3)=1		! interpolate gaps in response patterns
CommOp(5,1)=3 ; CommOp(5,2)=2 ; CommOp(5,4)=1		! load global mean T time-series
CommOp(6,1)=31 ; CommOp(6,2)=1 ; CommOp(6,3)=1 ; CommOp(6,4)=2 ; CommOp(6,5)=2
CommOp(7,1)=0 ; CommOp(7,2)=2
CommOp(8,1)=31 ; CommOp(8,2)=1 ; CommOp(8,3)=1 ; CommOp(8,4)=4 ; CommOp(8,5)=2
CommOp(9,1)=1 ; CommOp(9,2)=2
CommOp(10,1)=31 ; CommOp(10,2)=1 ; CommOp(10,3)=1 ; CommOp(10,4)=4 ; CommOp(10,5)=2
CommOp(11,1)=3 ; CommOp(11,2)=2 ; CommOp(11,4)=1
CommOp(12,1)=17 ; CommOp(12,2)=1 ; CommOp(12,3)=1 ; CommOp(12,4)=2 ; CommOp(12,5)=2 ; CommKay(12)=-2
CommOp(13,1)=3 ; CommOp(13,2)=2 ; CommOp(13,4)=1
CommOp(14,1)=17 ; CommOp(14,2)=1 ; CommOp(14,3)=1 ; CommOp(14,4)=2 ; CommOp(14,5)=2 ; CommKay(14)=2
CommOp(15,1)=43 ; CommOp(15,2)=1


end subroutine CalcOutputs

!*******************************************************************************
! write .ops file

subroutine WriteOpsFile

call date_and_time (Date, Time)
OpsFile = Date(3:8) // Time(1:4) // ".ops"

open  (1,file=OpsFile,status="new",form="formatted",action="write")

write (1,"(5i9)") NYear, NMonth, NBox, NExe, NWye
write (1,"(5i9)") NVariable, NArray, NComm, YearAD0, YearAD1
write (1,"(4f9.2)") (RefBounds(XBound),XBound=1,4)

do XComm = 1, NComm
  write (1,"(5i9,f9.2)") (CommOp(XComm,XElement), XElement=1,5), CommKay(XComm)

  do XVariable = 1, NVariable
    write (1,"(a)") trim(CommFile(XComm,XVariable))
    write (1,"(a)") trim(CommInfo(XComm,XVariable))
  end do
end do

close (1)

print*, "  > Specifications saved as: ", trim(OpsFile)

end subroutine WriteOpsFile

!*******************************************************************************
! set final things

subroutine Finalise

print*

end subroutine Finalise

!*******************************************************************************

end program ynSC20UnPack