annfiles_pik.f000644 125236 002424 00000022616 07616514432 013472 0ustar00f709762cru000000 000000 ! annfiles.f90 ! module in which standard save to and load from .ann file routines are held ! contains: MakeANNColTitles, LoadANN, SaveANN ! 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 module ANNFiles use FileNames implicit none contains !******************************************************************************* ! make line format for .ann save subroutine MakeANNFormat (ColN,LineFormat,DecPlaceN) integer, intent (in) :: ColN integer, intent (in), optional :: DecPlaceN character (len=20), intent(out) :: LineFormat integer :: ReadStatus, PlaceN character (len=20) :: DeciTxtLong character (len= 4) :: ColNText character (len= 1) :: DeciTxtShort LineFormat = "" if (present(DecPlaceN)) then PlaceN = DecPlaceN if (PlaceN.LT.1) PlaceN = 1 if (PlaceN.GT.3) PlaceN = 3 DeciTxtLong = GetTextFromInt (PlaceN) DeciTxtLong = adjustl(DeciTxtLong) DeciTxtShort = DeciTxtLong(1:1) else print*, " > Enter the no. of decimal places to save (1...3): " do read (*,*,iostat=ReadStatus) DeciTxtShort if (ReadStatus.GT.0) print*, " > Not a string. Try again." if (DeciTxtShort.EQ."") print*, " > A null string. Try again." if (ReadStatus.LE.0.AND.DeciTxtShort.NE."") exit end do end if open (1,file="col-scratch.txt",status="replace") write (1,"(i4)") ColN rewind (1) read (1,"(a4)") ColNText close (1) LineFormat = '(i4,' // trim(adjustl(ColNText)) // 'f9.' // DeciTxtShort // ')' end subroutine MakeANNFormat !******************************************************************************* ! make headers for .ann files subroutine MakeANNHeaders (LineFormat,YearAD,ColTitles,Headers) integer, dimension (:), pointer :: YearAD character (len=9), dimension (:), pointer :: ColTitles character (len=200), dimension (4), intent (out) :: Headers character (len=20), intent(in) :: LineFormat real, parameter :: MissVal = -999.0 integer :: YearN, ColN, XCol character (len=12) :: Date, Time character (len=4 ) :: Year, Year0, Year1, ColNText character (len=2 ) :: Month, Day, Hour, Minute !*************************************** Headers = "" YearN = size (YearAD) ColN = size (ColTitles) open (1,file="year-scratch.txt",status="replace") write (1,"(3i4)") YearAD(1), YearAD(YearN), ColN rewind (1) read (1,"(3a4)") Year0, Year1, ColNText close (1) call date_and_time (Date, Time) Year = Date (1:4) Month = Date (5:6) Day = Date (7:8) Hour = Time (1:2) Minute= Time (3:4) Headers (1) = "Tyndall Centre file (www.tyndall.ac.uk) created on " & // Day // "." // Month // "." // Year // " at " // Hour // ":" // Minute & // " by Dr. Tim Mitchell" Headers (2) = "Annual (.ann) file format: year AD, with" // ColNText // " data columns" Headers (3) = "Period = " // Year0 // "-" // Year1 // ": missing value = -999.0 : format = " // LineFormat Headers (4) = 'YEAR' do XCol = 1, ColN Headers (4) = trim(Headers(4)) // adjustr(ColTitles(XCol)) end do end subroutine MakeANNHeaders !******************************************************************************* ! load headers for .ann data files subroutine LoadANNHeaders (FileName,LineFormat,ColTitles,YearAD) integer, dimension (:), pointer :: YearAD character (len=9), dimension (:), pointer :: ColTitles character (len=80), intent (in) :: FileName character (len=20), intent (out) :: LineFormat real, parameter :: MissVal = -999.0 integer :: Year0, Year1, YearN, AllocStat, XYear, XCol, ColN character (len=4) :: ColNText character (len=20) :: HeaderFormat character (len=80) :: LoadName, Trash !*************************************** open (1, file=FileName, status="old", action="read") read (1,*) Trash ! ownership and date stamp read (1,"(a40,i4)") Trash, ColN ! number of columns read (1,"(a9,i4,a1,i4,a36,a20)") Trash, Year0, Trash, Year1, Trash, LineFormat open (2,file="year-scratch.txt",status="replace") write (2,"(i4)") ColN rewind (2) read (2,"(a4)") ColNText close (2) allocate (ColTitles (ColN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadANNHeaders: Allocation failure #####" HeaderFormat = "(a4," // trim(adjustl(ColNText)) // "a9)" read (1,HeaderFormat) Trash, (ColTitles(XCol), XCol=1,ColN) close (1) YearN = Year1 - Year0 + 1 allocate (YearAD (YearN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadANNHeaders: Allocation failure #####" do XYear = 1, YearN YearAD (XYear) = XYear + Year0 - 1 end do end subroutine LoadANNHeaders !******************************************************************************* ! make titles to columns subroutine MakeANNColTitles (ColTitles) character (len=9), pointer, dimension (:) :: ColTitles integer :: XCol, ColN, ChosenCol integer :: QModify, ReadStatus print*, " > Modify column titles (1=no,2=yes) ? " do read (*,*,iostat=ReadStatus) QModify if (ReadStatus.LE.0.AND.QModify.GE.1.AND.QModify.LE.2) exit end do if (QModify.EQ.2) then ColN = size (ColTitles) print "(a16)", " COL TITLE" do XCol = 1, ColN print "(i7,a9)", XCol, ColTitles(XCol) end do do if (ColN.GT.1) then print*, " > Select column (0=all,-1=exit): " do read (*,*,iostat=ReadStatus) ChosenCol if (ReadStatus.LE.0.AND.ChosenCol.GE.-1.AND.ChosenCol.LE.ColN) exit end do else ChosenCol = 0 end if if (ChosenCol.EQ. 0) then do ChosenCol = 1, ColN print "(a44,i4)", " > Enter title (9 characters) for column: ", ChosenCol do read (*,*,iostat=ReadStatus) ColTitles(ChosenCol) if (ReadStatus.GT.0) print*, " > Not a string. Try again." if (ColTitles(ChosenCol).EQ."") print*, " > A blank. Try again." if (ColTitles(ChosenCol).NE."") ColTitles(ChosenCol) = adjustr(ColTitles(ChosenCol)) if (ReadStatus.LE.0.AND.ColTitles(ChosenCol).NE."") exit end do end do ChosenCol = -1 end if if (ChosenCol.GE. 1) then print*, " > Enter title (9 characters): " do read (*,*,iostat=ReadStatus) ColTitles(ChosenCol) if (ReadStatus.GT.0) print*, " > Not a string. Try again." if (ColTitles(ChosenCol).EQ."") print*, " > A blank. Try again." if (ColTitles(ChosenCol).NE."") ColTitles(ChosenCol) = adjustr(ColTitles(ChosenCol)) if (ReadStatus.LE.0.AND.ColTitles(ChosenCol).NE."") exit end do end if if (ChosenCol.EQ.-1) exit end do end if end subroutine MakeANNColTitles !******************************************************************************* ! load .ann files subroutine LoadANN (CallFile, YearAD, ColTitles, Data) real, pointer, dimension (:,:) :: Data ! YearN, ColN integer, pointer, dimension (:) :: YearAD character (len=9), pointer, dimension(:) :: ColTitles character (len=80), intent(in) :: CallFile ! can be blank real, parameter :: MissVal = -999.0 integer :: ReadStatus, AllocStat integer :: YearN, ColN integer :: XYear, XCol, XHeader integer :: Year character (len=99) :: Trash character (len=80) :: LoadFile character (len=20) :: LineFormat character (len= 4) :: Suffix !*************************************** Suffix = ".ann" LoadFile = LoadPath (CallFile,Suffix) call LoadANNHeaders (LoadFile,LineFormat,ColTitles,YearAD) YearN = size (YearAD) ColN = size (ColTitles) allocate (Data(YearN,ColN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadANN: Allocation failure #####" open (2, file=LoadFile, status="old", access="sequential", form="formatted", action="read") do XHeader = 1, 4 read (2, "(a99)") Trash end do do XYear = 1, YearN read (2 ,LineFormat) Year, (Data(XYear,XCol), XCol=1,ColN) end do close (2) end subroutine LoadANN !******************************************************************************* ! save .ann files: multiple columns with common precision subroutine SaveANN (CallFile, YearAD, ColTitles, Data, DecPlaceN) real, pointer, dimension (:,:) :: Data ! YearN, ColN integer, pointer, dimension (:) :: YearAD character (len=9), pointer, dimension(:) :: ColTitles ! can be 'X', no blanks integer, intent (in), optional :: DecPlaceN character (len=80), intent(in) :: CallFile ! can be blank character (len=200), dimension (4) :: Headers real, parameter :: MissVal = -999.0 integer :: ReadStatus, AllocStat integer :: YearN, ColN integer :: XYear, XCol, XHeader character (len=80) :: SaveFile character (len=20) :: LineFormat character (len= 4) :: Suffix !*************************************** Suffix = ".ann" YearN = size (YearAD) ColN = size (ColTitles) if (YearN.NE.size(Data,1)) then print*, " > SaveANN. Mismatch between YearAD and Data arrays. No save." else if (ColN .NE.size(Data,2)) then print*, " > SaveANN. Mismatch between ColTitles and Data arrays. No save." else if (present(DecPlaceN)) then call MakeANNFormat (ColN,LineFormat,DecPlaceN) else call MakeANNFormat (ColN,LineFormat) end if call MakeANNHeaders (LineFormat,YearAD,ColTitles,Headers) SaveFile = SavePath (CallFile,Suffix) open (2, file=SaveFile, status="replace", access="sequential", form="formatted", action="write") do XHeader = 1, 4 write (2,"(a)") trim(Headers(XHeader)) end do do XYear = 1, YearN write (2,LineFormat) YearAD(XYear), (Data(XYear,XCol), XCol=1,ColN) end do close (2) end if end subroutine SaveANN !******************************************************************************* end module ANNFiles basicfun_pik.f000644 125236 002424 00000013063 07616514431 013460 0ustar00f709762cru000000 000000 ! basicfun.f90 ! module procedure written by Tim Mitchell on 14.12.99 ! if you want the slice statistics routines, go look in SliceStats.f90 module BasicFun implicit none contains !******************************************************************************* ! subtract Two from One ! print*, " > Calc A-B (=1), 100*(A-B)/B (=2), A/B (=3) ?" function OpTwoFromOne (InputA,InputB,MenuChoice) real OpTwoFromOne real, parameter :: MissVal = -999.0 real :: InputA, InputB integer :: MenuChoice OpTwoFromOne = MissVal if (InputA.NE.MissVal.AND.InputB.NE.MissVal) then if (MenuChoice.EQ.1) then OpTwoFromOne = InputA - InputB else if (MenuChoice.EQ.2) then OpTwoFromOne = 100 * (InputA - InputB) / abs (InputB) else if (MenuChoice.EQ.3) then OpTwoFromOne = InputA / InputB end if end if end function OpTwoFromOne !******************************************************************************* ! calculate weighted mean ! if InputN is unknown it may be set to 'missing' function OpWeightedMean (Input, Weights, InputN, MissAccept) real OpWeightedMean real, pointer, dimension (:) :: Input, Weights real, parameter :: MissVal = -999.0 real :: InputSum, WeightSum, TotWeight, MissAccept, Thresh integer :: InputN, XInput if (InputN.EQ.MissVal) InputN = size (Input,1) InputSum = 0.0 WeightSum = 0.0 TotWeight = 0.0 do XInput = 1, InputN TotWeight = TotWeight + Weights(XInput) if (Input(XInput).NE.MissVal.AND.Weights(XInput).NE.MissVal) then InputSum = InputSum + (Input (XInput) * Weights (XInput)) WeightSum = WeightSum + Weights (XInput) end if end do Thresh = (100.0 - MissAccept) * TotWeight / 100.0 if (WeightSum.GE.Thresh) then OpWeightedMean = InputSum / WeightSum else OpWeightedMean = MissVal end if end function OpWeightedMean !******************************************************************************* ! calculate mean ! if InputN is unknown it may be set to 'missing' function OpCalcMean (Input, InputN, MissAccept, CallMissVal) real OpCalcMean real, pointer, dimension (:) :: Input real, optional :: CallMissVal real :: MissVal real :: InputSum, ValidN, MissAccept, Thresh integer :: InputN, XInput if (present(CallMissVal)) then MissVal = CallMissVal else MissVal = -999.0 end if if (InputN.EQ.MissVal) InputN = size (Input,1) InputSum = 0.0 ValidN = 0.0 Thresh = (100.0 - MissAccept) * InputN / 100.0 do XInput = 1, InputN if (Input(XInput).NE.MissVal) then InputSum = InputSum + Input(XInput) ValidN = ValidN + 1.0 end if end do if (ValidN.GE.Thresh) then OpCalcMean = InputSum / ValidN else OpCalcMean = MissVal end if end function OpCalcMean !******************************************************************************* ! calculate stdev ! QSampPop: 1=sample-sd, 2=population-sd, 3=sample-variance, 4=pop-variance ! if InputN is unknown it may be set to 'missing' function OpCalcStDev (Input, InputN, MissAccept, QSampPop, CallMissVal) real OpCalcStDev real, pointer, dimension (:) :: Input real, optional :: CallMissVal real :: MissVal real :: InputSum, InputSumSq, ValidN, MissAccept, Thresh, Result integer :: InputN, XInput, QSampPop if (present(CallMissVal)) then MissVal = CallMissVal else MissVal = -999.0 end if if (InputN.EQ.MissVal) InputN = size (Input,1) InputSum = 0.0 InputSumSq = 0.0 ValidN = 0.0 Thresh = (100.0 - MissAccept) * InputN / 100.0 do XInput = 1, InputN if (Input(XInput).NE.MissVal) then InputSumSq = InputSumSq + (Input(XInput) ** 2) InputSum = InputSum + Input(XInput) ValidN = ValidN + 1.0 end if end do if (ValidN.GE.Thresh.AND.ValidN.GE.2) then Result = (InputSumSq / ValidN) - ( (InputSum / ValidN) ** 2 ) if (Result.GE.0) then if (QSampPop.EQ.2.OR.QSampPop.EQ.4) Result = Result * ValidN / (ValidN - 1) if (QSampPop.EQ.1.OR.QSampPop.EQ.2) Result = sqrt (Result) OpCalcStDev = Result else OpCalcStDev = MissVal end if else OpCalcStDev = MissVal end if end function OpCalcStDev !******************************************************************************* ! operate on one value using another ! " > Divide =1, times =2, minus =3, add =4, sqroot =5, expon =6, abs =7" ! added 30.10.01: 8=log(InputA), 9=e(InputA), 10=ln(InputA) function OpAwithB (InputA,InputB,MenuChoice) real OpAwithB real, parameter :: MissVal = -999.0 real :: InputA, InputB integer :: MenuChoice, IntegerB OpAwithB = MissVal if (MenuChoice.EQ.1) then if (InputA.NE.MissVal.AND.InputB.NE.MissVal.AND.InputB.NE.0) OpAwithB = InputA / InputB else if (MenuChoice.EQ.2) then if (InputA.NE.MissVal.AND.InputB.NE.MissVal) OpAwithB = InputA * InputB else if (MenuChoice.EQ.3) then if (InputA.NE.MissVal.AND.InputB.NE.MissVal) OpAwithB = InputA - InputB else if (MenuChoice.EQ.4) then if (InputA.NE.MissVal.AND.InputB.NE.MissVal) OpAwithB = InputA + InputB else if (MenuChoice.EQ.5) then if (InputA.NE.MissVal.AND.InputA.GT.0) OpAwithB = sqrt (InputA) else if (MenuChoice.EQ.6) then IntegerB = nint (InputB) if (InputA.NE.MissVal.AND.InputB.NE.MissVal) OpAwithB = InputA ** IntegerB else if (MenuChoice.EQ.7) then if (InputA.NE.MissVal) OpAwithB = abs (InputA) else if (MenuChoice.EQ.8) then if (InputA.NE.MissVal.AND.InputA.GT.0) OpAwithB = log10 (InputA) else if (MenuChoice.EQ.9) then if (InputA.NE.MissVal) OpAwithB = exp (InputA) else if (MenuChoice.EQ.10) then if (InputA.NE.MissVal.AND.InputA.GT.0) OpAwithB = log (InputA) end if end function OpAwithB !******************************************************************************* end module BasicFun filenames_pik.f000644 125236 002424 00000026676 07616514432 013650 0ustar00f709762cru000000 000000 ! filenames.f90 ! module to hold standard routines for obtaining file/name/paths for save and load ! contains: ! LoadPath, SavePath for getting paths with(out) a specified suffix ! GetTextFromReal for getting text from a real ! GetIntFromText, GetTextFromInt for getting text from an integer, or vice versa ! GetBatch for getting a filtered selection of filepaths from the filesystem ! MakeBatch for making a set of new filepaths with common stem ! last modified 09.08.01 ! 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 module FileNames implicit none contains !******************************************************************************* ! when fed with the filter (a text string such as '/cru/mydir/*.glo'), which can be "blank" ! returns a string array (size FileN) with all the files that fit that filter subroutine GetBatch (CallFilter,Batch) character (len=80), pointer, dimension (:) :: Batch character (len=80), intent (in) :: CallFilter integer :: ReadStatus, AllocStat integer :: FileN integer :: XFile character (len=80), parameter :: Blank = "" character (len= 80) :: Filter, NamesFile, CountFile character (len=200) :: CommandLine !character (len= 80) :: Dummy1, Dummy2 !*************************************** do if (CallFilter.EQ.Blank) then print*, " > Enter the filter by which to identify the files:" do read (*,*,iostat=ReadStatus) Filter if (ReadStatus.LE.0.AND.Filter.NE."") exit end do else Filter = CallFilter end if NamesFile = 'deleteme.names.txt' CommandLine = 'ls -1 ' // trim(adjustl(Filter)) // ' > ' // trim(adjustl(NamesFile)) CommandLine = trim (CommandLine) call system (CommandLine) CountFile = 'deleteme.count.txt' CommandLine = 'wc -l ' // trim(adjustl(NamesFile)) // ' > ' // trim(adjustl(CountFile)) CommandLine = trim (CommandLine) call system (CommandLine) open (1, file=CountFile, status="old", access="sequential", form="formatted", action="read") read (1, fmt="(i8)") FileN ! read (1, fmt="(i8,x,a)") FileN, Dummy1, Dummy2 close (1) print "(a,i4)", " > Number of files found: ", FileN if (FileN.EQ.0) print*, " > Re-enter the filter." if (FileN.GT.0) exit end do allocate (Batch (FileN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: GetBatch: Allocation failure #####" Batch = "" open (1, file=NamesFile, status="old", access="sequential", form="formatted", action="read") do XFile = 1, FileN read (1,"(a80)") Batch (XFile) end do close (1) CommandLine = 'rm deleteme*' CommandLine = trim (CommandLine) call system (CommandLine) end subroutine GetBatch !******************************************************************************* ! obtain path for loading file (no .Z or .X required, although it is allowable) ! note that this will find any relevant zipped or binary file and return it ! it is permissible to call the routine with a blank suffix, in which case no suffix is required function LoadPath (CallFile,Suffix) character (len=80) :: CallFile, GivenFile, FoundFile, LoadPath character (len=4) :: Suffix integer :: ReadStatus, StrLen, CallFileLen, NullSuffix CallFileLen = len_trim(CallFile) ! remove any zip or binary expression if (CallFileLen.GE.2) then if (CallFile((CallFileLen-1):CallFileLen).EQ.".Z".OR.CallFile((CallFileLen-1):CallFileLen).EQ.".X") then GivenFile = CallFile(1:(CallFileLen-2)) // " " else GivenFile = CallFile end if else GivenFile = CallFile end if NullSuffix = 0 if (Suffix.EQ." ".OR.Suffix.EQ."") NullSuffix = 1 do if (GivenFile.EQ."") then print*, " > Enter the file, with suffix: ", Suffix do read (*,*,iostat=ReadStatus) GivenFile if (ReadStatus.GT.0) print*, " > Not a string. Try again." if (GivenFile.EQ."") print*, " > A null string. Try again." if (ReadStatus.LE.0.AND.GivenFile.NE."") exit end do end if StrLen = len(trim(GivenFile)) if (GivenFile((StrLen-1):StrLen).EQ.'.Z'.OR.GivenFile((StrLen-1):StrLen).EQ.'.X') then GivenFile((StrLen-1):StrLen) = ' ' ! remove any zip or binary StrLen = len(trim(GivenFile)) end if if (NullSuffix.EQ.0.AND.GivenFile((StrLen-3):StrLen).NE.Suffix) then ! file has wrong suffix print*, " > Suffix is not: ", Suffix ReadStatus = 1 ; GivenFile = "" else inquire (file=GivenFile, name=FoundFile) ! look for file FoundFile=GivenFile open (1, file=FoundFile, status="old", action="read", iostat=ReadStatus) if (ReadStatus .NE. 0) then FoundFile = trim(FoundFile) // ".Z" ! look for zipped file open (1, file=FoundFile, status="old", action="read", iostat=ReadStatus) if (ReadStatus .NE. 0) then inquire (file=GivenFile, name=FoundFile) ! look for file FoundFile=GivenFile FoundFile = trim(FoundFile) // ".X" ! look for binary file open (1, file=FoundFile, status="old", action="read", form="unformatted", iostat=ReadStatus) if (ReadStatus .NE. 0) then print*, " > Failed to find file." GivenFile = "" end if end if end if end if if (ReadStatus .EQ. 0) close (1) if (ReadStatus .EQ. 0) exit end do LoadPath = FoundFile end function LoadPath !******************************************************************************* ! obtain SavePath for save file ! note that this will check for any zipped or binary file with the same name and not permit it function SavePath (CallFile,Suffix) character (len=80) :: CallFile, GivenFile, FoundFile, SavePath, ZipFile, BinFile character (len=4) :: Suffix integer :: ReadStatus, StrLen, CallFileLen, NullSuffix CallFileLen = len_trim(CallFile) ! remove any zip or binary expression if (CallFileLen.GE.2) then if (CallFile((CallFileLen-1):CallFileLen).EQ.".Z".OR.CallFile((CallFileLen-1):CallFileLen).EQ.".X") then GivenFile = CallFile(1:(CallFileLen-2)) // " " else GivenFile = CallFile end if else GivenFile = CallFile end if NullSuffix = 0 if (Suffix.EQ." ") NullSuffix = 1 do if (GivenFile.EQ."") then print*, " > Enter the file, with suffix: ", Suffix do read (*,*,iostat=ReadStatus) GivenFile if (ReadStatus.GT.0) print*, " > Not a string. Try again." if (GivenFile.EQ."") print*, " > A null string. Try again." if (ReadStatus.LE.0.AND.GivenFile.NE."") exit end do end if StrLen = len(trim(GivenFile)) if (GivenFile((StrLen-1):StrLen).EQ.'.Z'.OR.GivenFile((StrLen-1):StrLen).EQ.'.X') then GivenFile((StrLen-1):StrLen) = ' ' ! remove any zip or binary StrLen = len(trim(GivenFile)) end if if (NullSuffix.EQ.0.AND.GivenFile((StrLen-3):StrLen).NE.Suffix) then ! file has wrong suffix print*, " > Suffix is not: ", Suffix GivenFile = "" ; GivenFile = "" else ! look for file inquire (file=GivenFile, name=FoundFile) FoundFile=GivenFile ! print *, GivenFile ! print *, FoundFile ! print *, Suffix ! print *, CallFile open (1, file=FoundFile, status="new", action="write", iostat=ReadStatus) if (ReadStatus .NE. 0) then print*, " > Failed to create file. Try again." GivenFile = "" else close (1) call system ('rm ' // FoundFile) ! remove successfully created file StrLen = len_trim (FoundFile) ZipFile = FoundFile(1:StrLen) // ".Z" open (1, file=ZipFile, status="new", action="write", iostat=ReadStatus) if (ReadStatus .NE. 0) then print*, " > A zipped file has this name. Try again." GivenFile = "" else close (1) call system ('rm ' // ZipFile) ! remove successfully created zip file BinFile = FoundFile(1:StrLen) // ".X" open (1, file=BinFile, status="new", action="write", iostat=ReadStatus) if (ReadStatus .NE. 0) then print*, " > A binary file has this name. Try again." GivenFile = "" else close (1) call system ('rm ' // BinFile) ! remove successfully created bin file end if end if end if end if if (GivenFile.NE."") exit end do SavePath = FoundFile end function SavePath !******************************************************************************* ! feed this function with an Real and it will return it as text function GetTextFromReal (Real) character (len=20) :: GetTextFromReal integer :: ReadStatus real :: Real open (98,file="temp.trash",status="replace",action="readwrite",iostat=ReadStatus) write (98,"(a15,f5.1)") " ", Real rewind(98) read (98,"(a20)") GetTextFromReal close (98) GetTextFromReal = adjustl(GetTextFromReal) end function GetTextFromReal !******************************************************************************* ! feed this function with an integer and it will return it as text function GetTextFromInt (Int) character (len=20) :: GetTextFromInt integer :: Int, ReadStatus open (98,file="temp.trash",status="replace",action="readwrite",iostat=ReadStatus) write (98,"(i20)") Int rewind(98) read (98,"(a20)") GetTextFromInt close (98) GetTextFromInt = adjustl(GetTextFromInt) end function GetTextFromInt !******************************************************************************* ! feed this function with text and it will return it as an integer function GetIntFromText (Text) character (len=20) :: Text integer :: GetIntFromText, ReadStatus open (98,file="temp.trash",status="replace",action="readwrite",iostat=ReadStatus) write (98,"(a20)") Text rewind(98) read (98,"(i20)") GetIntFromText close (98) end function GetIntFromText !******************************************************************************* ! feed it with the constant stem and tip, and the unique text strings ! to insert between the stem and tip in order to make the set of filepaths (Batch) subroutine MakeBatch (Stem,Tip,UniqueName,Batch) character (len=80), pointer, dimension (:) :: Batch character (len=20), pointer, dimension (:) :: UniqueName character (len=80), intent(in) :: Stem,Tip integer :: AllocStat integer :: XFile,XLetter,XFileCheck integer :: FileN,LetterN integer :: NameLen character (len=80) :: Name FileN = size (UniqueName,1) allocate (Batch(FileN),stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: MakeBatch: Allocation failure #####" Batch=" " ; Batch = trim(adjustl(Stem)) do XFile = 1, FileN ! iterate by File Name = trim(adjustl(UniqueName(XFile))) LetterN = len_trim(Name) if (LetterN.GE.1) then do XLetter = 1, LetterN ! for each letter in the Fileion name if (Name(XLetter:XLetter).NE." ".AND.Name(XLetter:XLetter).NE.".") then Batch(XFile) = trim(Batch(XFile)) // Name(XLetter:XLetter) ! add it to filepath else Batch(XFile) = trim(Batch(XFile)) // "_" ! or add _ to filepath end if end do else Batch(XFile) = trim(Stem) // "v" end if end do do XFile = 2, FileN ! check for duplicates XFileCheck = 0 do XFileCheck = XFileCheck + 1 if (Batch(XFile).EQ.Batch(XFileCheck)) then Name = GetTextFromInt(XFile) Batch(XFile) = trim(Batch(XFile)) // trim(Name) end if if (XFileCheck.EQ.(XFile-1)) exit end do end do do XFile = 1, FileN ! add the final suffixes Batch(XFile) = trim(Batch(XFile)) // trim(adjustl(Tip)) end do end subroutine MakeBatch !******************************************************************************* end module FileNames grimfiles_pik.f000644 125236 002424 00000104640 07616514432 013652 0ustar00f709762cru000000 000000 ! grimfiles.f90 ! module holding standard save and load GRIM (GRIdded Monthly) routines ! SaveGrim,LoadGrim,SaveGrimMonthPost,LoadGrimMonthPrep,SaveGrimMonth,LoadGrimMonth ! SaveGrip,LoadGrip,CheckGridAB ! written by Tim Mitchell on 22.03.01 ! last modified on 20.11.01 ! 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 module GrimFiles use FileNames use Time implicit none contains !******************************************************************************* ! this ensures that the data is saved in the order implicit in the grid dimensions: ! i.e. long (X) varies slowest (-180...180), and lat (Y) quickest (-90...90) subroutine SaveGrim (Data,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,SaveSuffix,NoZip,Silent) real, pointer, dimension (:,:,:) :: Data integer, pointer, dimension (:,:) :: Grid integer, pointer, dimension (:) :: YearAD real, dimension (4) :: Bounds integer, dimension (12) :: YearData character (len=80), dimension (5) :: Headers integer, optional, intent (in) :: NoZip,Silent character (len=80), intent (in) :: Info, CallFile character (len= 4), intent (in) :: CallSuffix character (len= 4), intent (out) :: SaveSuffix real, parameter :: MissVal = -999.0 real :: Divisor integer :: YearN,MonthN,BoxN,ExeN,WyeN integer :: XYear,XMonth,XBox,XHeader,XExe,XWye integer :: AllocStat, TotExceed character (len=80) :: SaveFile, Variable !*************************************** YearN = size (Data,1) ; MonthN = size (Data,2) ; BoxN = size (Data,3) ExeN = size (Grid,1) ; WyeN = size (Grid,2) ; TotExceed = 0 if (size(YearAD).NE.YearN) print*, " > ##### ERROR: SaveGrim: incompatible arrays #####" if (MonthN .NE. 12) print*, " > ##### ERROR: SaveGrim: weird calendar #####" !if (Bounds(1).LT.-180.OR.Bounds(2).GT.180.OR.Bounds(1).GT.Bounds(2)) & ! print*, " > ##### ERROR: SaveGrim: longitudes not within spherical geometry #####" ! !if (Bounds(3).LT.-90.OR.Bounds(4).GT.90.OR.Bounds(3).GT.Bounds(4)) & ! print*, " > ##### ERROR: SaveGrim: latitudes not within spherical geometry #####" if (Bounds(2).EQ.Bounds(1)) print*, " > ##### ERROR: SaveGrim: zero longitude range #####" if (Bounds(4).EQ.Bounds(3)) print*, " > ##### ERROR: SaveGrim: zero latitude range #####" call ReviewCall (CallFile,CallSuffix,SaveFile,SaveSuffix,2) ! checks for file/suffix consistency !*************************************** call CheckVariSuffix (SaveSuffix,Variable,Divisor) ! get Divisor if (.not.present(Silent)) print "(2a)", " > ", trim(Info) open (2,file=SaveFile,status="replace",access="sequential",form="formatted",action="write") call WriteHeaders (SaveSuffix,Info,Bounds,ExeN,WyeN,BoxN,YearAD(1),YearAD(YearN),"grim") do XExe = 1, ExeN do XWye = 1, WyeN if (Grid(XExe,XWye).NE.MissVal) then write (2,"(a9,i4,a1,i4)") "Grid-ref=", XExe, ",", XWye do XYear = 1, YearN do XMonth = 1, MonthN if (Data(XYear,XMonth,Grid(XExe,XWye)).NE.MissVal) then YearData(XMonth) = nint(Data(XYear,XMonth,Grid(XExe,XWye)) / Divisor) if (YearData(XMonth).GT.99999.OR.YearData(XMonth).LT.-9999) then YearData(XMonth) = nint(MissVal) TotExceed = TotExceed + 1 end if else YearData(XMonth) = nint(MissVal) end if end do write (2,"(12i5)") (YearData(XMonth), XMonth=1,MonthN) end do end if end do end do close (2) if (TotExceed.GT.0) print "(a,i9)", " > Box-months outside range: ", TotExceed !*************************************** if (present(NoZip)) then if (.not.present(Silent)) print*, " > The file is not being compressed." else if (.not.present(Silent)) print*, " > Compressing in background..." call system ('compress ' // SaveFile // ' &') end if end subroutine SaveGrim !******************************************************************************* ! when using this routine, do not assume that this Grid and your Grid are identical ! however, if the same cells are non-missing in the original saves to grim, ! the SaveGrim procedure ensures that the Grids obtained will be the same from the two loads from grim ! Data (A) may be defined (=Y) or not (=N) in the call; MasterYearAD (B) may (=Y) or may not (=N) be specified ! if A=Y,B=Y: loading requires sizes A=B, and common period between B and file ! if A=Y,B=N: loading requires sizes A=file ! if A=N,B=Y: loading requires common period between B and file, Data=(MasterYearN,12,BoxN) ! if A=N,B=N: loaded automatically, Data=(FileYearN,12,BoxN) ! if Data is defined in the call, it must also have already been initialised if that is necessary ! FileYearAD IS the genuine FileYearAD, and may differ from Data and MasterYearAD ! use Override with caution! it fills in FileYear0 (and hence FileYear1) subroutine LoadGrim (Data,Grid,FileYearAD,Bounds,Info,CallFile,CallSuffix,LoadSuffix,& MasterYearAD,Quiet,OverRide) real, pointer, dimension (:,:,:) :: Data integer, pointer, dimension (:,:) :: Grid integer, pointer, dimension (:) :: FileYearAD integer, pointer, optional, dimension (:) :: MasterYearAD ! may or may not be included: see top integer, dimension (12) :: YearData character (len=80), dimension (5) :: Headers real, dimension (4) :: Bounds integer, intent(in), optional :: Quiet,OverRide character (len=80), intent (out) :: Info character (len=80), intent (in) :: CallFile character (len= 4), intent (in) :: CallSuffix character (len= 4), intent (out) :: LoadSuffix real, parameter :: MissVal = -999.0 real :: Multiplier integer :: FileYearN,MasterYearN,MonthN,BoxN,ExeN,WyeN integer :: XFileYear,XMasterYear,XMonth,XBox,XHeader,XExe,XWye integer :: AllocStat integer :: LoadFileLen,FileYearAD0,FileYearAD1,FileYear0,FileYear1,MasterYear0,MasterYear1 integer :: QZip integer :: DataEn,DataMiss character (len=80) :: LoadFile, Variable, Trash !*************************************** call ReviewCall (CallFile,CallSuffix,LoadFile,LoadSuffix,1) ! checks for file/suffix consistency LoadFileLen = len_trim(LoadFile) if (LoadFileLen.GT.1.AND.LoadFile((LoadFileLen-1):LoadFileLen).EQ.".Z") then QZip = 2 ! file is zipped if (.not.present(Quiet)) print*, " > Uncompressing on the fly..." call system ('uncompress ' // LoadFile) LoadFile ((LoadFileLen-1):LoadFileLen) = " " ! change filename to the unzipped file else QZip = 1 ! file already unzipped end if call CheckVariSuffix (LoadSuffix,Variable,Multiplier) ! get multiplier !*************************************** open (2,file=LoadFile,status="old",access="sequential",form="formatted",action="read") call ReadHeaders (LoadSuffix,Info,Bounds,ExeN,WyeN,BoxN,FileYearAD0,FileYearAD1,"grim") if (Bounds(2).EQ.Bounds(1)) print*, " > ##### ERROR: LoadGrim: zero longitude range #####" if (Bounds(4).EQ.Bounds(3)) print*, " > ##### ERROR: LoadGrim: zero latitude range #####" FileYearN = FileYearAD1 - FileYearAD0 + 1 MonthN = 12 if (present(MasterYearAD).AND.associated(Data)) then if (size(Data,1).NE.size(MasterYearAD,1)) then print*, " > ##### ERROR: LoadGrim: call mismatch #####" MasterYearN = MissVal else MasterYearN = size(Data,1) end if else if (present(MasterYearAD)) then MasterYearN = size(MasterYearAD,1) else if (associated(Data)) then MasterYearN = size(Data,1) else MasterYearN = FileYearN end if if (associated(Data)) then if (size(Data,2).NE.MonthN) print*, " > ##### ERROR: LoadGrim: Months mismatch #####" if (size(Data,3).NE.BoxN) print*, " > ##### ERROR: LoadGrim: Domain mismatch #####" else allocate (Data (MasterYearN,MonthN,BoxN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadGRIM: Allocation failure: Data #####" Data = MissVal end if allocate (Grid (ExeN,WyeN), & FileYearAD (FileYearN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadGRIM: Allocation failure: All #####" Grid = MissVal if (present(OverRide)) then FileYearAD0=OverRide ; FileYearAD1=OverRide+FileYearN-1 end if do XFileYear = 1, FileYearN FileYearAD(XFileYear) = FileYearAD0 + XFileYear - 1 end do if (present(MasterYearAD)) then call CommonVecPer (FileYearAD,MasterYearAD,FileYear0,FileYear1,MasterYear0,MasterYear1) if (FileYear0.EQ.MissVal) then print "(a)", " > ##### ERROR: LoadGrim: The file and spec are mismatched in period #####" FileYear0 = MissVal ; FileYear1 = MissVal ; MasterYear0 = MissVal ; MasterYear1 = MissVal else if (.not.present(Quiet)) print "(a,i4,a1,i4))", " > Period loading: ", FileYearAD(FileYear0), & "-", FileYearAD(FileYear1) end if else if (size(Data,1).EQ.FileYearN) then FileYear0 = 1 ; FileYear1 = FileYearN ; MasterYear0 = 1 ; MasterYear1 = FileYearN else print*, " > ##### ERROR: LoadGrim: Period mismatch #####" FileYear0 = MissVal ; FileYear1 = MissVal ; MasterYear0 = MissVal ; MasterYear1 = MissVal end if end if if (FileYear0.NE.MissVal) then if (.not.present(Quiet)) print "(2a)", " > ", trim(Info) DataEn = 0 ; DataMiss = 0 do XBox = 1, BoxN read (2,"(a9,i4,a1,i4)") Trash, XExe, Trash, XWye Grid(XExe,XWye) = XBox do XFileYear = 1, FileYearN read (2,"(12i5)") (YearData(XMonth), XMonth=1,MonthN) if (XFileYear.GE.FileYear0.AND.XFileYear.LE.FileYear1) then XMasterYear = MasterYear0 + XFileYear - FileYear0 do XMonth = 1, MonthN DataEn = DataEn + 1 if (YearData(XMonth).EQ.MissVal) then Data(XMasterYear,XMonth,XBox) = MissVal DataMiss = DataMiss + 1 else Data(XMasterYear,XMonth,XBox) = real(YearData(XMonth)) * Multiplier end if end do end if end do end do if (.not.present(Quiet)) print "(a,2i12)", " > Total and missing values: ", DataEn, DataMiss else print*, " > ##### Nothing was loaded from the file. #####" end if close (2) !*************************************** if (QZip.EQ.2) then call system ('compress ' // LoadFile // ' &') if (.not.present(Quiet)) print*, " > Recompressing in the background..." end if end subroutine LoadGrim !******************************************************************************* ! saves GRIdded Seasonic files (12*month,4*season,1*annual) ! this ensures that the data is saved in the order implicit in the grid dimensions: ! i.e. long (X) varies slowest (-180...180), and lat (Y) quickest (-90...90) subroutine SaveGrip (Data,Grid,Bounds,Info,CallFile,CallSuffix,SaveSuffix) real, pointer, dimension (:,:) :: Data integer, pointer, dimension (:,:) :: Grid real, dimension (4) :: Bounds character (len=80), dimension (5) :: Headers character (len=80), intent (in) :: Info, CallFile character (len= 4), intent (in) :: CallSuffix character (len= 4), intent (out) :: SaveSuffix real, parameter :: MissVal = -999.0 real :: Divisor integer :: SeasonN,BoxN,ExeN,WyeN integer :: XSeason,XBox,XHeader,XExe,XWye integer :: AllocStat character (len=80) :: SaveFile, Variable, SaveBinFile !*************************************** SeasonN = size (Data,1) ; BoxN = size (Data,2) ExeN = size (Grid,1) ; WyeN = size (Grid,2) if (SeasonN .NE. 17) print*, " > ##### ERROR: SaveGrip: weird calendar #####" !if (Bounds(1).LT.-180.OR.Bounds(2).GT.180.OR.Bounds(1).GT.Bounds(2)) & ! print*, " > ##### ERROR: SaveGrip: longitudes not within spherical geometry #####" ! !if (Bounds(3).LT.-90.OR.Bounds(4).GT.90.OR.Bounds(3).GT.Bounds(4)) & ! print*, " > ##### ERROR: SaveGrip: latitudes not within spherical geometry #####" if (Bounds(2).EQ.Bounds(1)) print*, " > ##### ERROR: SaveGrip: zero longitude range #####" if (Bounds(4).EQ.Bounds(3)) print*, " > ##### ERROR: SaveGrip: zero latitude range #####" call ReviewCall (CallFile,CallSuffix,SaveFile,SaveSuffix,2) ! checks for file/suffix consistency SaveFile = trim(SaveFile) // ".X" !*************************************** print "(2a)", " > ", trim(Info) open (2,file=SaveFile,status="replace",access="sequential",form="unformatted",action="write") call WriteHeaders (SaveSuffix,Info,Bounds,ExeN,WyeN,BoxN,0,0,"grip") do XExe = 1, ExeN do XWye = 1, WyeN if (Grid(XExe,XWye).NE.MissVal) then write (2) XExe, XWye, (Data(XSeason,Grid(XExe,XWye)), XSeason=1,SeasonN) end if end do end do close (2) end subroutine SaveGrip !******************************************************************************* ! when using this routine, do not assume that this Grid and your Grid are identical ! however, if the same cells are non-missing in the original saves to grim, ! the SaveGrip procedure ensures that the Grids obtained will be the same from the two loads from grim ! Data (A) may be defined (=Y) or not (=N) in the call ! if Data is defined in the call, it must also have already been initialised if that is necessary subroutine LoadGrip (Data,Grid,Bounds,Info,CallFile,CallSuffix,LoadSuffix) real, pointer, dimension (:,:) :: Data real, dimension (17) :: Line integer, pointer, dimension (:,:) :: Grid character (len=80), dimension (5) :: Headers real, dimension (4) :: Bounds character (len=80), intent (out) :: Info character (len=80), intent (in) :: CallFile character (len= 4), intent (in) :: CallSuffix character (len= 4), intent (out) :: LoadSuffix real, parameter :: MissVal = -999.0 integer, parameter :: SeasonN = 17 integer :: BoxN,ExeN,WyeN integer :: XSeason,XBox,XHeader,XExe,XWye integer :: AllocStat integer :: LoadFileLen,QBin integer :: Year0,Year1 integer :: DataEn,DataMiss character (len=80) :: LoadFile, Variable, Trash !*************************************** call ReviewCall (CallFile,CallSuffix,LoadFile,LoadSuffix,1) ! checks for file/suffix consistency open (2,file=LoadFile,status="old",access="sequential",form="unformatted",action="read") call ReadHeaders (LoadSuffix,Info,Bounds,ExeN,WyeN,BoxN,Year0,Year1,"grip") if (Bounds(2).EQ.Bounds(1)) print*, " > ##### ERROR: LoadGrip: zero longitude range #####" if (Bounds(4).EQ.Bounds(3)) print*, " > ##### ERROR: LoadGrip: zero latitude range #####" if (associated(Data)) then if (size(Data,1).NE.SeasonN) print*, " > ##### ERROR: LoadGrip: SeasonN mismatch #####" if (size(Data,2).NE.BoxN) print*, " > ##### ERROR: LoadGrip: Domain mismatch #####" else allocate (Data (SeasonN,BoxN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadGrip: Allocation failure: Data #####" Data = MissVal end if allocate (Grid (ExeN,WyeN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadGrip: Allocation failure: All #####" Grid = MissVal if (Year1.GE.Year0) then print "(2a)", " > ", trim(Info) do XBox = 1, BoxN read (2) XExe, XWye, (Line(XSeason), XSeason=1,SeasonN) Grid (XExe,XWye) = XBox do XSeason = 1, SeasonN Data(XSeason,XBox) = Line(XSeason) end do end do else print*, " > Nothing was loaded from the file." end if close (2) end subroutine LoadGrip !******************************************************************************* ! this saves the data from individual month/year files into one big array ! so do work, save to year/month files, and call this ! hence Save Monthly Grim files Post-analysis subroutine SaveGrimMonthPost (DataFile,BoxN,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,SaveSuffix) real, pointer, dimension (:,:,:) :: Data real, pointer, dimension (:) :: DataMonth integer, pointer, dimension (:,:) :: Grid integer, pointer, dimension (:) :: YearAD character (len=80), pointer, dimension (:,:) :: DataFile character (len= 4), pointer, dimension (:) :: YearText character (len= 2), pointer, dimension (:) :: MonthText real, dimension (4) :: Bounds integer, intent (in) :: BoxN character (len=80), intent (in) :: Info, CallFile character (len= 4), intent (in) :: CallSuffix character (len= 4), intent (out) :: SaveSuffix real, parameter :: MissVal = -999.0 integer :: AllocStat integer :: YearN,MonthN,XYear,XMonth,XBox integer :: CurrentPos,CurrentCod character (len= 4) :: GivenText !*************************************** YearN = size (DataFile,1) MonthN = size (DataFile,2) allocate (Data(YearN,MonthN,BoxN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: SaveGrimMonthPost: Allocation failure #####" do XYear = 1, YearN ! get Data from individual time steps do XMonth = 1, MonthN call LoadGrimMonth (DataFile(XYear,XMonth),DataMonth) ! ... and dump to file do XBox = 1, BoxN Data(XYear,XMonth,XBox) = DataMonth(XBox) end do deallocate (DataMonth,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: SaveGrimMonthPost: Deallocation failure #####" end do end do call SaveGrim (Data,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,SaveSuffix) deallocate (Data,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: SaveGrimMonthPrep: Deallocation failure #####" end subroutine SaveGrimMonthPost !******************************************************************************* ! this loads the data into one big array ! saves it as individual month/year files subroutine LoadGrimMonthPrep (DataFile,BoxN,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,LoadSuffix) real, pointer, dimension (:,:,:) :: Data real, pointer, dimension (:) :: DataMonth integer, pointer, dimension (:,:) :: Grid integer, pointer, dimension (:) :: YearAD character (len=80), pointer, dimension (:,:) :: DataFile character (len= 4), pointer, dimension (:) :: YearText character (len= 2), pointer, dimension (:) :: MonthText real, dimension (4) :: Bounds integer, intent (out) :: BoxN character (len=80), intent (in) :: CallFile character (len=80), intent (out) :: Info character (len= 4), intent (in) :: CallSuffix character (len= 4), intent (out) :: LoadSuffix real, parameter :: MissVal = -999.0 integer :: AllocStat integer :: YearN,MonthN,XYear,XMonth,XBox integer :: CurrentPos,CurrentCod character (len= 4) :: GivenText !*************************************** call LoadGrim (Data,Grid,YearAD,Bounds,Info,CallFile,CallSuffix,LoadSuffix) YearN = size (Data,1) MonthN = size (Data,2) BoxN = size (Data,3) allocate (DataFile (YearN,MonthN), & YearText (YearN), & MonthText (MonthN), & DataMonth (BoxN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadGrimMonthPrep: Allocation failure #####" open (1,file="/cru/scratch2/f709762/year-scratch.txt",status="replace") write (1,"(i4)") YearAD(1) rewind (1) read (1,"(a4)") GivenText close (1) YearText(1) = GivenText ! get year AD in text form do XYear = 2, YearN GivenText = YearText (XYear-1) CurrentPos = 4 do CurrentCod = iachar(GivenText(CurrentPos:CurrentPos)) if (CurrentCod.LE.56) then GivenText(CurrentPos:CurrentPos) = achar(CurrentCod+1) CurrentPos = 0 else GivenText(CurrentPos:CurrentPos) = achar(48) CurrentPos = CurrentPos - 1 end if if (CurrentPos.EQ.0) exit end do YearText (XYear) = GivenText end do MonthText(1) = "01" ; MonthText(2) = "02" ; MonthText(3) = "03" ; MonthText(4) = "04" MonthText(5) = "05" ; MonthText(6) = "06" ; MonthText(7) = "07" ; MonthText(8) = "08" MonthText(9) = "09" ; MonthText(10) = "10" ; MonthText(11) = "11" ; MonthText(12) = "12" do XYear = 1, YearN ! split up Data into individual time steps do XMonth = 1, MonthN do XBox = 1, BoxN DataMonth(XBox) = Data(XYear,XMonth,XBox) end do DataFile(XYear,XMonth) = '/cru/scratch2/f709762/' // YearText(XYear) // '.' & // MonthText(XMonth) // LoadSuffix call SaveGrimMonth (DataFile(XYear,XMonth),DataMonth) ! ... and dump to file end do end do deallocate (Data,DataMonth,YearText,MonthText,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadGrimMonthPrep: Deallocation failure #####" end subroutine LoadGrimMonthPrep !******************************************************************************* ! checks the filename and suffix given in call ! returns the correct suffix, and correct filename with the correct suffix subroutine ReviewCall (CallFile,CallSuffix,CheckFile,CheckSuffix,LoadSave) character (len=80), intent (in) :: CallFile character (len= 4), intent (in) :: CallSuffix character (len=80), intent (out) :: CheckFile character (len= 4), intent (out) :: CheckSuffix integer, intent (in) :: LoadSave ! 1=load,2=save character (len=80), parameter :: Blank = "" character (len=80) :: Variable, WorkFile character (len= 4) :: FileSuffix real :: Factor integer :: SuffixStart,ReadStatus,ZipStart,CallFileLen !*************************************** CheckFile = "" ; CheckSuffix = "" ; Variable = "" CallFileLen = len_trim(CallFile) ! remove any zip or binary expression if (CallFileLen.GT.1) then if (CallFile((CallFileLen-1):CallFileLen).EQ.".Z".OR.CallFile((CallFileLen-1):CallFileLen).EQ.".X") then WorkFile = CallFile(1:(CallFileLen-2)) // " " else WorkFile = CallFile end if else WorkFile = CallFile end if if (CallSuffix.EQ."".OR.CallSuffix.EQ." ") then ! suffix is not present in call if (WorkFile.EQ."") then ! and there is no file in call do print*, " > Enter the grim filepath: " ! get filepath do read (*,*,iostat=ReadStatus) WorkFile if (ReadStatus.EQ.0) exit end do CheckSuffix = GetFileSuffix (WorkFile) ! get suffix from filepath if (CheckSuffix.NE."") exit end do ! check filepath if (LoadSave.EQ.1) CheckFile = LoadPath (WorkFile,CheckSuffix) if (LoadSave.EQ.2) CheckFile = SavePath (WorkFile,CheckSuffix) else ! and there is a file in call do CheckSuffix = GetFileSuffix (WorkFile) ! get suffix from filepath if (CheckSuffix.EQ."") then print*, " > Enter the grim filepath: " ! get filepath do read (*,*,iostat=ReadStatus) WorkFile if (ReadStatus.EQ.0) exit end do end if if (CheckSuffix.NE."") exit end do ! check filepath if (LoadSave.EQ.1) CheckFile = LoadPath (WorkFile,CheckSuffix) if (LoadSave.EQ.2) CheckFile = SavePath (WorkFile,CheckSuffix) end if else ! suffix present in call call CheckVariSuffix (CallSuffix,Variable,Factor) if (Variable.EQ."") then ! but unrecognised... print*, " > Supplied suffix is invalid. Supply another. " ! therefore obtain another do CheckSuffix = GetFreshSuffix () if (Variable.NE."") exit end do else ! and is recognised CheckSuffix = CallSuffix end if if (WorkFile.NE."") then ! gotta suffix, investigate file in call FileSuffix = GetFileSuffix (WorkFile) ! get suffix from file if (FileSuffix.EQ.CheckSuffix) then ! suffix good, so check filepath if (LoadSave.EQ.1) CheckFile = LoadPath (WorkFile,CheckSuffix) if (LoadSave.EQ.2) CheckFile = SavePath (WorkFile,CheckSuffix) else if (FileSuffix.EQ."") then ! suffix bad, so get filepath print*, " > File suffix is not recognised." else print*, " > File suffix is not the same as the supplied suffix." end if if (LoadSave.EQ.1) CheckFile = LoadPath (Blank,CheckSuffix) if (LoadSave.EQ.2) CheckFile = SavePath (Blank,CheckSuffix) end if else ! gotta suffix but no file in call if (LoadSave.EQ.1) CheckFile = LoadPath (Blank,CheckSuffix) if (LoadSave.EQ.2) CheckFile = SavePath (Blank,CheckSuffix) end if end if end subroutine ReviewCall !******************************************************************************* ! get any legitimate suffix function GetFreshSuffix (Verbose, NoPrompt) integer, intent (in), optional :: Verbose integer, optional :: NoPrompt real :: Factor integer :: ReadStatus character (len=80) :: Variable character (len= 4) :: GetFreshSuffix, Suffix if (.not.present(NoPrompt)) print*, " > Enter the grim file suffix: " Variable = "" do read (*,*,iostat=ReadStatus) Suffix if (Suffix.EQ."") then print*, " > Blank is unacceptable. Try again." else call CheckVariSuffix (Suffix,Variable,Factor) if (Variable.EQ."") print*, " > Suffix unrecognised. Try again." end if if (ReadStatus.EQ.0.AND.Variable.NE."") exit end do if (present(Verbose)) print "(2a)", " > Variable: ", trim(Variable) GetFreshSuffix = Suffix end function GetFreshSuffix !******************************************************************************* ! extract suffix from file function GetFileSuffix (CallFile,Verbose) real :: Factor integer, intent (in), optional :: Verbose integer :: ReadStatus, CallFileLen, SuffixStart character (len=80) :: Variable, CallFile, WorkFile character (len= 4) :: GetFileSuffix, Suffix CallFileLen = len_trim(CallFile) ! remove any zip or binary expression if (CallFileLen.GT.1) then if (CallFile((CallFileLen-1):CallFileLen).EQ.".Z".OR.CallFile((CallFileLen-1):CallFileLen).EQ.".X") then WorkFile = CallFile(1:(CallFileLen-2)) // " " else WorkFile = CallFile end if else WorkFile = CallFile end if Suffix = "" ; Variable = "" ; GetFileSuffix = " " SuffixStart = index (WorkFile,".",.TRUE.) ! identify suffix if (SuffixStart.GT.0) then ! check suffix Suffix = WorkFile(SuffixStart:(SuffixStart+3)) call CheckVariSuffix (Suffix,Variable,Factor) if (Variable.EQ."") print*, " > Suffix unrecognised." if (Variable.NE."") then if (present(Verbose)) print "(2a)", " > Variable: ", trim(Variable) GetFileSuffix = Suffix ! return suffix end if else print*, " > No identifiable suffix in filepath." end if end function GetFileSuffix !******************************************************************************* ! relies on Suffix being correct !! subroutine WriteHeaders (Suffix,Info,Bounds,ExeN,WyeN,BoxN,YearAD0,YearAD1,FileType) real, dimension (4) :: Bounds character (len=80), intent (in) :: Info character (len= 4), intent (in) :: Suffix, FileType integer, intent(in) :: BoxN,ExeN,WyeN,YearAD0,YearAD1 real :: Multi integer :: XBound character (len=80) :: Variable character (len=12) :: Date,Time character (len= 4) :: Year character (len= 2) :: Month,Day,Hour,Minute !*************************************** call CheckVariSuffix (Suffix,Variable,Multi) call date_and_time (Date, Time) Year = Date (1:4) Month = Date (5:6) Day = Date (7:8) Hour = Time (1:2) Minute= Time (3:4) if (FileType.EQ."grim") then write (2,"(a15,a4,a17,a2,a1,a2,a1,a4,a4,a2,a1,a2,a20)") "Tyndall Centre ", FileType, " file created on ", & Day, ".", Month, ".", Year, " at ", Hour, ":", Minute, " by Dr. Tim Mitchell" write (2,"(a4,a3,a)") Suffix, " = ", trim(Variable) write (2,"(a)") Info write (2,"(a6,f7.2,a1,f7.2,a8,f7.2,a1,f7.2,a12,i4,a1,i4,a1)") & "[Long=", Bounds(1), ",", Bounds(2), "] [Lati=", Bounds(3), ",", Bounds(4), & "] [Grid X,Y=", ExeN, ",", WyeN, "]" write (2,"(a7,i8,a9,i4,a1,i4,a9,f10.4,a16)") & "[Boxes=", BoxN, "] [Years=", YearAD0, "-", YearAD1, "] [Multi=", Multi, "] [Missing=-999]" else if (FileType.EQ."grip") then write (2) FileType, Day, Month, Year, Hour, Minute write (2) Suffix, Variable write (2) Info write (2) (Bounds(XBound),XBound=1,4), ExeN, WyeN write (2) BoxN, 0, 0, 1.0, -999.0 end if end subroutine WriteHeaders !******************************************************************************* ! checks filename suffix against header suffix subroutine ReadHeaders (FileNameSuffix,Info,Bounds,ExeN,WyeN,BoxN,YearAD0,YearAD1,FileType) real, dimension (4) :: Bounds character (len= 4), intent (in) :: FileNameSuffix,FileType character (len=80), intent (out) :: Info integer, intent(out) :: BoxN,ExeN,WyeN,YearAD0,YearAD1 real :: LoadMulti, LoadMissVal integer :: XBound character (len=80) :: Trash, LoadVariable character (len= 4) :: HeaderSuffix, LoadFileType character (len= 4) :: Year character (len= 2) :: Month,Day,Hour,Minute !*************************************** if (FileType.EQ."grim") then read (2,"(a15,a4)") Trash, LoadFileType read (2,"(a4)") HeaderSuffix read (2,"(a80)") Info read (2,"(a6,f7.2,a1,f7.2,a8,f7.2,a1,f7.2,a12,i4,a1,i4)") & Trash, Bounds(1), Trash, Bounds(2), Trash, Bounds(3), Trash, Bounds(4), & Trash, ExeN, Trash, WyeN read (2,"(a7,i8,a9,i4,a1,i4)") Trash, BoxN, Trash, YearAD0, Trash, YearAD1 else if (FileType.EQ."grip") then read (2) LoadFileType, Day, Month, Year, Hour, Minute ! read (2), HeaderSuffix, LoadVariable read (2) HeaderSuffix read (2) Info read (2) (Bounds(XBound),XBound=1,4), ExeN, WyeN read (2) BoxN, YearAD0,YearAD1, LoadMulti, LoadMissVal end if if (LoadFileType.EQ."file") LoadFileType = "grim" if (LoadFileType.NE.FileType) then print*, " > ##### ERROR: ReadHeaders: Mismatch between called and actual file type." YearAD1 = YearAD0 - 1 end if if (HeaderSuffix.NE.FileNameSuffix) then print*, " > ##### ERROR: ReadHeaders: Mismatch between filename and header suffixes." YearAD1 = YearAD0 - 1 end if end subroutine ReadHeaders !******************************************************************************* ! saves a single month's data in a temporary file subroutine SaveGrimMonth (CallFile,DataMonth) real, pointer, dimension (:) :: DataMonth character (len=80), intent (in) :: CallFile integer :: BoxN, XBox BoxN = size (DataMonth) open (3,file=CallFile,status="new",access="sequential",form="unformatted",action="write") write (3) BoxN do XBox = 1, BoxN write (3) DataMonth(XBox) end do close (3) end subroutine SaveGrimMonth !******************************************************************************* ! loads a single month's data from a temporary file subroutine LoadGrimMonth (CallFile,DataMonth) real, pointer, dimension (:) :: DataMonth character (len=80), intent (in) :: CallFile integer :: BoxN, XBox, AllocStat open (3,file=CallFile,status="old",access="sequential",form="unformatted",action="read") read (3) BoxN allocate (DataMonth(BoxN),stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadGrimMonth: Allocation failure #####" do XBox = 1, BoxN read (3) DataMonth(XBox) end do close (3) end subroutine LoadGrimMonth !******************************************************************************* ! BoxN : >0:grids are the same with this number of boxes, MissVal=grids are not the same subroutine CheckGridAB (GridA,GridB,BoxN) integer, pointer, dimension (:,:) :: GridA, GridB integer, intent(out) :: BoxN real, parameter :: MissVal = -999.0 integer :: ExeN,WyeN,XExe,XWye BoxN = 0 if (size(GridA,1).NE.size(GridB,1)) then BoxN = MissVal else ExeN = size(GridA,1) end if if (size(GridA,2).NE.size(GridB,2)) then BoxN = MissVal else WyeN = size(GridA,2) end if if (BoxN.NE.MissVal) then XExe = 0 do XExe = XExe + 1 ; XWye = 0 do XWye = XWye + 1 if (GridA(XExe,XWye).EQ.GridB(XExe,XWye)) then if (GridA(XExe,XWye).NE.MissVal) BoxN = BoxN + 1 else BoxN = MissVal end if if (XWye.EQ.WyeN) exit if (BoxN.EQ.MissVal) exit end do if (XExe.EQ.ExeN) exit if (BoxN.EQ.MissVal) exit end do end if end subroutine CheckGridAB !******************************************************************************* ! check given file suffix ! Variable is either blank (suffix unrecognised) or contains description of variable ! Factor is zero unless Variable is set subroutine CheckVariSuffix (Suffix,Variable,Factor) real, intent (out) :: Factor character (len=80), intent (out) :: Variable character (len= 4), intent (in) :: Suffix Variable = "undefined" if (Suffix.EQ.".alp") then Variable="precip (mm) gamma distribution alpha" Factor = 0.01 else if (Suffix.EQ.".bet") then Variable="precip (mm) gamma distribution beta" Factor = 0.0001 else if (Suffix.EQ.".cld") then Variable="cloud cover (percentage)" Factor = 0.1 else if (Suffix.EQ.".dew") then Variable="dew point temperature (degrees Celsius)" Factor = 0.1 else if (Suffix.EQ.".dir") then Variable="wind direction (degrees)" Factor = 0.1 else if (Suffix.EQ.".dpd") then Variable="dew point depression (degrees Celsius)" Factor = 0.1 else if (Suffix.EQ.".dsf") then Variable="downward shortwave flux (W/m2)" Factor = 0.1 else if (Suffix.EQ.".dtr") then Variable="diurnal temperature range (degrees Celsius)" Factor = 0.1 else if (Suffix.EQ.".elv") then Variable="elevation (km)" Factor = 0.001 else if (Suffix.EQ.".frs") then Variable="ground frost frequency (days)" Factor = 0.01 else if (Suffix.EQ.".oct") then Variable="cloud cover (octas)" Factor = 0.01 else if (Suffix.EQ.".pre") then Variable="precipitation (mm)" Factor = 0.1 else if (Suffix.EQ.".prc") then Variable="precip change (%)" Factor = 0.001 else if (Suffix.EQ.".prd") then Variable="precip change (%)" Factor = 0.1 else if (Suffix.EQ.".prp") then Variable="precipitation rate (mm/day)" Factor = 0.1 else if (Suffix.EQ.".prr") then Variable="precipitation rate (mm/day)" Factor = 0.001 else if (Suffix.EQ.".rad") then Variable="radiation (W/(m*m))" Factor = 0.1 else if (Suffix.EQ.".rhm") then Variable="relative humidity (%)" Factor = 0.1 else if (Suffix.EQ.".slp") then Variable="sea level pressure" Factor = 0.1 else if (Suffix.EQ.".spc") then Variable="sunshine percentage" Factor = 0.1 else if (Suffix.EQ.".spn") then Variable="sunshine (percentage of normal)" Factor = 0.1 else if (Suffix.EQ.".ssh") then Variable="sunshine hours" Factor = 0.01 else if (Suffix.EQ.".stp") then Variable="station pressure" Factor = 0.1 else if (Suffix.EQ.".tmp") then Variable="near-surface temperature (degrees Celsius)" Factor = 0.1 else if (Suffix.EQ.".tmk") then Variable="near-surface temperature (degrees Kelvin)" Factor = 0.1 else if (Suffix.EQ.".tmn") then Variable="near-surface temperature minimum (degrees Celsius)" Factor = 0.1 else if (Suffix.EQ.".tik") then Variable="near-surface temperature minimum (degrees Kelvin)" Factor = 0.1 else if (Suffix.EQ.".tmx") then Variable="near-surface temperature maximum (degrees Celsius)" Factor = 0.1 else if (Suffix.EQ.".twb") then Variable="wet-bulb temperature (degrees Celsius)" Factor = 0.1 else if (Suffix.EQ.".txk") then Variable="near-surface temperature maximum (degrees Kelvin)" Factor = 0.1 else if (Suffix.EQ.".tcl") then Variable="total cloud (fraction) from long-wave radiation" Factor = 0.0001 else if (Suffix.EQ.".tcs") then Variable="total cloud (fraction) from short-wave radiation" Factor = 0.0001 else if (Suffix.EQ.".vap") then Variable="vapour pressure (hPa)" Factor = 0.1 else if (Suffix.EQ.".wet") then Variable="wet day frequency (days)" Factor = 0.01 else if (Suffix.EQ.".wnd") then Variable="wind (m/s)" Factor = 0.1 else Variable = "" Factor = 0.0 end if end subroutine CheckVariSuffix !******************************************************************************* end module GRIMFiles time_pik.f000644 125236 002424 00000015452 07616516605 012635 0ustar00f709762cru000000 000000 ! time.f90 ! module to hold standard routines for manipulating temporal things ! contains: GaussSmooth,GetMonthLengths,CommonVecPer module Time implicit none contains !******************************************************************************* ! compare year vectors to find common period subroutine CommonVecPer (AYears,BYears,AStart,AEnd,BStart,BEnd) integer, pointer, dimension (:) :: AYears, BYears integer, intent (out) :: AStart,BStart,AEnd,BEnd integer :: AYearN,BYearN,ARest,BRest integer :: Offset,CommonLen real, parameter :: MissVal = -999.0 AYearN = size (AYears) BYearN = size (BYears) Offset = BYears(1) - AYears(1) + 1 if (Offset.GE.1) then AStart = Offset BStart = 1 else AStart = 1 BStart = 2 - Offset end if ARest = AYearN - AStart BRest = BYearN - BStart if (ARest.GE.0.AND.BRest.GE.0) then if (ARest.LT.BRest) then CommonLen = ARest + 1 else CommonLen = BRest + 1 end if AEnd = AStart + CommonLen - 1 BEnd = BStart + CommonLen - 1 else AStart=MissVal BStart=MissVal AEnd=MissVal BEnd=MissVal end if end subroutine CommonVecPer !******************************************************************************* ! smooths using Gaussian filtering subroutine GaussSmooth (TimeN,THalf,QRetainMiss,TSInVec,TSLowVec,TSHighVec,MissThresh) real, pointer, dimension (:) :: TSInVec,TSLowVec,TSHighVec ! def,alloc,init in call integer, intent (in) :: TimeN, THalf, QRetainMiss ! QRetainMiss:1=no,2=yes real, intent (in), optional :: MissThresh ! %missing=OK ; if>, output = MissVal real, allocatable, dimension (:) :: Weight, TSPad, TSPadExt real, parameter :: MissVal = -999.0 real :: WeightFactor, WeightRoot, WeightSum real :: LocalTot, LocalN real :: MeanStart, MeanEnd real :: Real0, Real1 real :: MissPercent integer :: AllocStat integer :: Integer0, Integer1 integer :: WeightN integer :: XWeight, XTime, XValue integer :: EndN integer :: MissTot, QSetMiss QSetMiss = 0 if (present(MissThresh)) then MissTot = 0 do XTime = 1, TimeN if (TSInVec(XTime).EQ.MissVal) MissTot = MissTot + 1 end do MissPercent = (real(MissTot) / real(TimeN)) * 100.0 if (MissPercent.GE.MissThresh) then TSInVec = 999 ; QSetMiss = 1 end if end if ! write (99,*), "##### computing no. weights" ! ################### WeightN = nint ((real(THalf)/2.5)+0.5) ! compute no. weights required if (modulo(WeightN,2).EQ.0) then WeightN = WeightN + 2 else WeightN = WeightN + 1 end if if (WeightN.LE.7) WeightN = 7 allocate (Weight(WeightN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Allocation failure #####" Weight = 0.0 ! compute weights ! write (99,*), "##### computing weights" ! ################### WeightFactor = -18.0 / (real(THalf*THalf)) WeightRoot = 1.0 / sqrt (2.0 * 3.1415927) Weight (1) = WeightRoot WeightSum = Weight (1) do XWeight = 2, WeightN Weight (XWeight) = WeightRoot * exp (WeightFactor * real(XWeight) * real(XWeight)) WeightSum = WeightSum + 2.0 * Weight (XWeight) end do do XWeight = 1, WeightN Weight (XWeight) = Weight (XWeight) / WeightSum end do ! pad the ts with local mean where missing ! write (99,*), "##### padding where missing" ! ################### allocate (TSPad(TimeN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Allocation failure #####" TSPad = 0.0 do XTime = 1, TimeN if (TSInVec(XTime).EQ.MissVal) then Real0 = max ((XTime-WeightN+1),1) Real1 = min ((XTime+WeightN-1),TimeN) Integer0 = nint (Real0) Integer1 = nint (Real1) LocalTot = 0.0 LocalN = 0.0 do XValue = Integer0, Integer1 if (TSInVec(XValue).NE.MissVal) then LocalTot = LocalTot + TSInVec(XValue) LocalN = LocalN + 1.0 end if end do if (LocalN.GT.0.AND.LocalTot.NE.0) TSPad (XTime) = LocalTot / LocalN else TSPad (XTime) = TSInVec (XTime) end if end do ! extend ends of ts by mean from each end ! write (99,*), "##### extending ends" ! ################### EndN = WeightN - 1 LocalTot = 0.0 LocalN = 0.0 do XTime = 1, EndN LocalTot = LocalTot + TSPad(XTime) LocalN = LocalN + 1.0 end do MeanStart = LocalTot / LocalN LocalTot = 0.0 LocalN = 0.0 do XTime = (TimeN-EndN+1), TimeN LocalTot = LocalTot + TSPad(XTime) LocalN = LocalN + 1.0 end do MeanEnd = LocalTot / LocalN allocate (TSPadExt(TimeN+2*EndN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Allocation failure #####" TSPadExt = 0.0 do XValue = 1, EndN TSPadExt (XValue) = MeanStart end do do XValue = (EndN+1), (EndN+TimeN) TSPadExt (XValue) = TSPad (XValue-EndN) end do do XValue = (EndN+TimeN+1), (TimeN+2*EndN) TSPadExt (XValue) = MeanEnd end do ! apply the filter ! write (99,*), "##### applying filter" ! ################### ! if it fails in the following section, it is probably because you haven't alloc all the arrays in call do XTime = 1, TimeN WeightSum = Weight (1) * TSPadExt (EndN+XTime) do XWeight = 2, WeightN WeightSum = WeightSum + Weight(XWeight) * & (TSPadExt(XTime+EndN-XWeight+1)+TSPadExt(XTime+EndN+XWeight-1)) end do TSLowVec (XTime) = WeightSum end do ! compute the high-residual ! write (99,*), "##### computing high residual" ! ################### do XTime = 1, TimeN TSHighVec (XTime) = TSInVec (XTime) - TSLowVec (XTime) end do ! reinsert MissVal unless called otherwise if (QRetainMiss.EQ.2) then do XTime = 1, TimeN if (TSInVec(XTime).EQ.MissVal) then TSLowVec (XTime) = MissVal TSHighVec(XTime) = MissVal end if end do end if deallocate (Weight, TSPad, TSPadExt) if (QSetMiss.EQ.1) then TSLowVec = MissVal ; TSHighVec = MissVal end if end subroutine GaussSmooth !******************************************************************************* ! get month lengths appropriate to years AD subroutine GetMonthLengths (YearAD,MonthLengths) integer, pointer, dimension (:,:) :: MonthLengths integer, pointer, dimension (:) :: YearAD real :: Rem004,Rem100,Rem400 integer :: XYear, XMonth, YearN, AllocStat YearN = size (YearAD) allocate ( MonthLengths(YearN,12), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: GetMonthLengths: Allocation failure #####" MonthLengths = 31 do XYear = 1, YearN MonthLengths(XYear,4 ) = 30 MonthLengths(XYear,6 ) = 30 MonthLengths(XYear,9 ) = 30 MonthLengths(XYear,11) = 30 Rem004 = mod (YearAD(XYear), 4) Rem100 = mod (YearAD(XYear),100) Rem400 = mod (YearAD(XYear),400) if (Rem004.EQ.0) then if (Rem100.EQ.0) then if (Rem400.EQ.0) then MonthLengths(XYear,2) = 29 else MonthLengths(XYear,2) = 28 end if else MonthLengths(XYear,2) = 29 end if else MonthLengths(XYear,2) = 28 end if end do end subroutine GetMonthLengths !******************************************************************************* end module Time tynsc20seqgrim_pik.f000600 125236 002424 00000047323 07617515336 014564 0ustar00f709762cru000000 000000 ! tynsc20seqgrim_pik.f90 ! f90 -o tynsc20seqgrim_pik time_pik.f90 filenames_pik.f90 grimfiles_pik.f90 annfiles_pik.f90 basicfun_pik.f90 tynsc20seqgrim_pik.f90 ! written by Tim Mitchell on 30.03.01 and subsequently modified ! it is designed to be run using Comms defined using a specification file from *unpack*.f90 ! 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 program TYNsc20SeqGrim use Time use FileNames use GrimFiles use AnnFiles use BasicFun implicit none real, pointer, dimension (:,:,:) :: OperA, OperB, OperC, OperD real, pointer, dimension (:,:,:) :: Grim1, Grim2, Grim3, Grim4, Grim5, Grim6 real, pointer, dimension (:,:,:) :: GrimXtra real, pointer, dimension (:,:) :: GridWeights, WeightsXtra, LinData,MonthData,SeasonData real, pointer, dimension (:) :: CommKay, VecA, VecB, VecC, VecD, GloSlice,AnnualData real, pointer, dimension (:) :: RegWeights,GloWeights real, dimension (4) :: Bounds, FileBounds integer, pointer, dimension (:,:) :: MapIDLReg ! shape of IDL mapping arrays integer, pointer, dimension (:) :: RegSizes ! region sizes + reg-->raw integer, pointer, dimension (:,:) :: CommOp, Grid, FileGrid, MonthLengths integer, pointer, dimension (:) :: YearAD, FileYearAD character (len=80), pointer, dimension (:,:) :: CommFile, CommInfo character (len=80), pointer, dimension (:) :: FileLineNames,RegFilePaths,OpsFile character (len=20), pointer, dimension (:) :: RegNames !names of individual regions character (len=09), pointer, dimension (:) :: ColTitles character (len= 2), dimension (12) :: MonthName = (/'01','02','03','04','05','06',& '07','08','09','10','11','12'/) character (len= 3), dimension (4) :: SeasonName = (/'MAM','JJA','SON','DJF'/) real, parameter :: MissVal = -999.0 real, parameter :: MissAccept = 10.0 integer, parameter :: SeasonN = 4 character (len=80), parameter :: Blank = "" real :: OpTot, OpEn, OpMiss, OpMean, OpThresh, OpValid, OpBelow, OpAbove real :: FileFactor, LSRAye,LSRBee,LSRPea integer :: AllocStat, CheckGrid integer :: YearN, MonthN, BoxN, ExeN, WyeN, ExecN, ArrayN, CommN, FileYearN, RegN integer :: XYear, XMonth, XBox, XExe, XWye, XExec, XArray, XComm, XElement, XPerYear, XBound, XSeason integer :: XFileYear, XMasterYear, XReg, XRegCheck, XLetter,XExeNear,XWyeNear integer :: AStart,AEnd,BStart,BEnd,XAyear,XBYear integer :: SuffixStart, FileSubBeg, InfoSubBeg, RegNameLen integer :: PerYear0, PerYear1, YearAD0, YearAD1, FileYear0, FileYear1, MasterYear0, MasterYear1 integer :: ThisMonth, ThisYear, NoZip, TotA,TotB,TotC,TotD, OpsFileN,Operate integer :: QMeanSum character (len=80) :: FileInfo,FilePath,FileVariable,GridFilePath,GloFile,GloInfo,GenericFile,SaveVariName character (len=80) :: Stem,Tip,ExtraStr,SpecFile,Filter,GridFile character (len=20) :: RegName character (len=40) :: RegTitle character (len= 4) :: FileSuffix, SaveSuffix !******************************************************************************* open (99,file="log.dat",status="replace",action="write") call FindLatestOpsFile call Initialise call Main call Finalise close (99) contains !******************************************************************************* subroutine FindLatestOpsFile Filter="*.ops" call GetBatch (Filter,OpsFile) OpsFileN = size(OpsFile) SpecFile = OpsFile(OpsFileN) end subroutine FindLatestOpsFile !******************************************************************************* subroutine Initialise nullify (OperA,OperB,OperC,OperD) nullify (Grim1,Grim2,Grim3,Grim4,Grim5,Grim6) nullify (GrimXtra,LinData,CommKay,VecA,VecB,VecC) nullify (CommOp,Grid,FileGrid,YearAD,FileYearAD) nullify (CommFile,CommInfo,FileLineNames) print "(2a)", " > Operating on: ", trim(SpecFile) print* call LoadSpec allocate (Grim1 (YearN,MonthN,BoxN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Initialise: Allocation failure: Grim1 #####" allocate (Grim2 (YearN,MonthN,BoxN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Initialise: Allocation failure: Grim2 #####" end subroutine Initialise !******************************************************************************* subroutine LoadSpec open (1,file=SpecFile,status="old",form="formatted",action="read") read (1,"(5i9)") YearN, MonthN, BoxN, ExeN, WyeN read (1,"(5i9)") ExecN, ArrayN, CommN, YearAD0, YearAD1 read (1,"(4f9.2)") (Bounds(XBound), XBound=1,4) print*, " > Basic specifications: " print "(a13,i4,a11,i2,a11,i6,a10,i4,a9,i4)", " > Years = ", YearN, ", Months = ", MonthN, & ", Domain = ", BoxN, ", Longs = ", ExeN, ", Lats = ", WyeN print "(a18,i2,a11,i1,a13,i2,a11,i4,a1,i4)", " > Executions = ", ExecN, ", Arrays = ", ArrayN, & ", Commands = ", CommN, ", Period = ", YearAD0, "-", YearAD1 allocate (CommOp (CommN,5), & CommKay (CommN), & CommFile (CommN,ExecN), & CommInfo (CommN,ExecN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadSpec: Allocation failure: Comm* #####" write (99,*) "SeqGrim : data found within comms file." do XComm = 1, CommN read (1,"(5i9,f9.2)") (CommOp(XComm,XElement), XElement=1,5), CommKay(XComm) write (99,"(6i6,f12.4)") XComm, (CommOp(XComm,XElement), XElement=1,5), CommKay(XComm) do XExec = 1, ExecN read (1,"(a)") CommFile(XComm,XExec) read (1,"(a)") CommInfo(XComm,XExec) write (99,"(a)") trim(adjustl(CommFile(XComm,XExec))) write (99,"(a)") trim(adjustl(CommInfo(XComm,XExec))) end do end do write (99,*) "SeqGrim : data found within comms file. END." close (1) allocate (YearAD (YearN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadSpec: Allocation failure: G+YAD #####" do XYear = 1, YearN YearAD(XYear) = YearAD0 + XYear - 1 end do GridFile = 'obs.clim6190.ateam.cld' call LoadGrim (GrimXtra,Grid,FileYearAD,FileBounds,FileInfo,GridFile,'.cld',FileSuffix) deallocate (GrimXtra,FileYearAD,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadSpec: Deallocation failure: Grid #####" end subroutine LoadSpec !******************************************************************************* subroutine Main do XExec = 1, ExecN print* call ResetAll call FindVariable do XComm = 1, CommN if (CommOp(XComm,1).EQ.-1) then print "(a,i1)", " > Resetting: ", CommOp(XComm,2) call ResetGrim else if (CommOp(XComm,1).EQ. 0) then print "(a,i1)", " > Loading into: ", CommOp(XComm,2) call GrabGrim else if (CommOp(XComm,1).EQ. 1) then print "(a,i1)", " > Loading single year into: ", CommOp(XComm,2) call GrabGrimYear else if (CommOp(XComm,1).EQ. 3) then call GrabAnn print "(a,i2,a,i1,a,i3,2(a,i4))", " > Loaded .ann (col ", CommOp(XComm,4), ") into: ", CommOp(XComm,2), & "; VALID:", TotA, "; COMMON: ",YearAD(AStart),"-",YearAD(AEnd) else if (CommOp(XComm,1).EQ.10) then print "(a,i1,a,i4,a1,i4,a,i1)", " > Climatology: ", CommOp(XComm,2), " = ", CommOp(XComm,4), "-", & CommOp(XComm,5), " mean from ", CommOp(XComm,3) call OpClim else if (CommOp(XComm,1).EQ.17) then call OpMaskRel print "(a,i1,a,i2,a,i1,a,i1,a,i1,a,i7)", " > Where ", CommOp(XComm,2), "[rel=", & nint(CommKay(XComm)), "]", CommOp(XComm,5), & ", ", CommOp(XComm,3), "=", CommOp(XComm,4), " ; changed: ", TotA else if (CommOp(XComm,1).EQ.30) then call OpGrimKay print "(3(a,i1),2(a,i10),a)", " > Operating: ", CommOp(XComm,2), " = f(", CommOp(XComm,3), & "): op ", CommOp(XComm,4), "; VALID: ",TotA,"=f(",TotB,")" else if (CommOp(XComm,1).EQ.31) then call OpGrimGrim print "(4(a,i1),3(a,i10),a)", " > Operating: ", CommOp(XComm,2), " = f(", CommOp(XComm,3), & ",", CommOp(XComm,5), "): op ", CommOp(XComm,4), & "; VALID: ",TotA,"=f(",TotB,",",TotC,")" else if (CommOp(XComm,1).EQ.43) then print "(a,i1)", " > Saving from: ", CommOp(XComm,2) call DumpGrim else if (CommOp(XComm,1).EQ.51) then call InterpolateMissing print "(2(a,i1),2(a,i9))", " > Interpol missing ", CommOp(XComm,2), "=", CommOp(XComm,3), & ": success=", TotA, ", failure=", TotB end if end do end do end subroutine Main !******************************************************************************* subroutine ResetAll if (associated(Grim1)) Grim1=MissVal if (associated(Grim2)) Grim2=MissVal end subroutine ResetAll !******************************************************************************* subroutine FindVariable XComm = 0 ; FileSuffix = "" ; FileVariable = "" ; FileFactor = MissVal do XComm = XComm + 1 if (CommFile(XComm,XExec).NE."") then SuffixStart = index(CommFile(XComm,XExec),'.',.TRUE.) if (SuffixStart.GT.0) then FilePath = CommFile(XComm,XExec) FileSuffix = FilePath(SuffixStart:(SuffixStart+3)) call CheckVariSuffix (FileSuffix,FileVariable,FileFactor) end if end if if (XComm.EQ.CommN.OR.FileVariable.NE."") exit end do if (FileSuffix.EQ.".pre".OR.FileSuffix.EQ.".frs".OR.FileSuffix.EQ.".wet") then QMeanSum = 2 else QMeanSum = 0 end if if (FileVariable.EQ."") FileVariable = "unknown" print "(2a)", " > Variable: ", trim(FileVariable) end subroutine FindVariable !******************************************************************************* ! command = -1 subroutine ResetGrim if (CommOp(XComm,2).EQ.1) then Grim1 = MissVal else if (CommOp(XComm,2).EQ.2) then Grim2 = MissVal else if (CommOp(XComm,2).EQ.3) then Grim3 = MissVal else if (CommOp(XComm,2).EQ.4) then Grim4 = MissVal else if (CommOp(XComm,2).EQ.5) then Grim5 = MissVal else if (CommOp(XComm,2).EQ.6) then Grim6 = MissVal end if end subroutine ResetGrim !******************************************************************************* ! command = 0 ! this has been altered to ignore the year structure and make it 2001-end subroutine GrabGrim call PointOperA print "(2a)", " > Seeking: ", trim(CommFile(XComm,XExec)) call LoadGrim (OperA,FileGrid,FileYearAD,FileBounds,FileInfo,CommFile(XComm,XExec)," ",FileSuffix,& MasterYearAD=YearAD,OverRide=2001) if (size(FileGrid,1).NE.size(Grid,1)) print*, " > ##### ERROR: loaded Grid mismatch: X" if (size(FileGrid,2).NE.size(Grid,2)) print*, " > ##### ERROR: loaded Grid mismatch: Y" CheckGrid = 0 do XExe = 1, ExeN do XWye = 1, WyeN if (FileGrid(XExe,XWye).NE.Grid(XExe,XWye)) CheckGrid = CheckGrid + 1 end do end do if (CheckGrid.NE.0) print*, " > ##### ERROR: Grid mismatches: ", CheckGrid nullify (OperA) deallocate (FileGrid,FileYearAD,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: GrabGrim: Deallocation failure: File* #####" end subroutine GrabGrim !******************************************************************************* ! command = 1 subroutine GrabGrimYear call PointOperA print "(2a)", " > Seeking: ", trim(CommFile(XComm,XExec)) call LoadGrim (GrimXtra,FileGrid,FileYearAD,FileBounds,FileInfo,CommFile(XComm,XExec)," ",FileSuffix) if (size(FileGrid,1).NE.size(Grid,1)) print*, " > ##### ERROR: loaded Grid mismatch: X" if (size(FileGrid,2).NE.size(Grid,2)) print*, " > ##### ERROR: loaded Grid mismatch: Y" CheckGrid = 0 do XExe = 1, ExeN do XWye = 1, WyeN if (FileGrid(XExe,XWye).NE.Grid(XExe,XWye)) CheckGrid = CheckGrid + 1 end do end do if (CheckGrid.NE.0) print*, " > ##### ERROR: Grid mismatches: ", CheckGrid if (size(FileYearAD).NE.1) print*, " > ##### ERROR: file period > 1 year #####" deallocate (FileGrid,FileYearAD,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: GrabGrimYear: Deallocation failure: File* #####" do XMonth = 1, MonthN do XBox = 1, BoxN do XYear = 1, YearN OperA(XYear,XMonth,XBox) = GrimXtra(1,XMonth,XBox) end do end do end do nullify (OperA) deallocate (GrimXtra,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: GrabGrimYear: Deallocation failure: GrimXtra #####" end subroutine GrabGrimYear !******************************************************************************* ! command = 3 subroutine GrabAnn call PointOperA print "(2a)", " > Seeking: ", trim(CommFile(XComm,XExec)) call LoadANN (CommFile(XComm,XExec),FileYearAD,ColTitles,LinData) call CommonVecPer (YearAD,FileYearAD,AStart,AEnd,BStart,BEnd) TotA = 0 if (AStart.EQ.-999.OR.BStart.EQ.-999) then print*, " > ##### ERROR: loaded YearAD wrong years #####" else if (size(LinData,2).LT.CommOp(XComm,4)) then print*, " > ##### ERROR: loaded array does not contain desired column #####" else do XAYear = AStart,AEnd XBYear = BStart + XAYear - 1 if (LinData(XBYear,CommOp(XComm,4)).NE.MissVal) TotA=TotA+1 do XMonth = 1, MonthN do XBox = 1, BoxN OperA(XAYear,XMonth,XBox) = LinData(XBYear,CommOp(XComm,4)) end do end do end do end if nullify (OperA) deallocate (FileYearAD,ColTitles,LinData,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: GrabAnn: Deallocation failure #####" end subroutine GrabAnn !******************************************************************************* ! command = 10 subroutine OpClim call PointOperA call PointOperB PerYear0 = 0 ; PerYear1 = 0 ; XYear = 0 do XYear = XYear + 1 if (YearAD(XYear).EQ.CommOp(XComm,4)) PerYear0 = XYear if (YearAD(XYear).EQ.CommOp(XComm,5)) PerYear1 = XYear if (PerYear1.GT.0.OR.XYear.EQ.YearN) exit end do OpThresh = real(PerYear1-PerYear0+1) * MissAccept if (PerYear1.GT.PerYear0) then do XBox = 1, BoxN do XMonth = 1, MonthN OpTot = 0.0 ; OpMiss = 0.0 ; OpEn = 0.0 do XYear = PerYear0, PerYear1 if (OperB(XYear,XMonth,XBox).NE.MissVal) then OpTot = OpTot + OperB(XYear,XMonth,XBox) OpEn = OpEn + 1.0 else OpMiss = OpMiss + 1.0 end if end do if (OpMiss.LE.OpThresh.AND.OpEn.GT.0) then OpMean = OpTot / OpEn else OpMean = MissVal end if do XYear = 1, YearN OperA(XYear,XMonth,XBox) = OpMean end do end do end do end if nullify (OperA,OperB) end subroutine OpClim !******************************************************************************* ! command = 17 subroutine OpMaskRel call PointOperA call PointOperB call PointOperC call PointOperD TotA=0 ; Operate=nint(CommKay(XComm)) do XMonth = 1, MonthN do XBox = 1, BoxN do XYear = 1, YearN if (OperA(XYear,XMonth,XBox).NE.MissVal.AND.OperC(XYear,XMonth,XBox).NE.MissVal) then if (Operate.EQ.-3.AND.OperA(XYear,XMonth,XBox).NE.OperC(XYear,XMonth,XBox)) then OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1 else if (Operate.EQ.-2.AND.OperA(XYear,XMonth,XBox).LT.OperC(XYear,XMonth,XBox)) then OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1 else if (Operate.EQ.-1.AND.OperA(XYear,XMonth,XBox).LE.OperC(XYear,XMonth,XBox)) then OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1 else if (Operate.EQ. 0.AND.OperA(XYear,XMonth,XBox).EQ.OperC(XYear,XMonth,XBox)) then OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1 else if (Operate.EQ. 1.AND.OperA(XYear,XMonth,XBox).GE.OperC(XYear,XMonth,XBox)) then OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1 else if (Operate.EQ. 2.AND.OperA(XYear,XMonth,XBox).GT.OperC(XYear,XMonth,XBox)) then OperB(XYear,XMonth,XBox)=OperD(XYear,XMonth,XBox) ; TotA=TotA+1 end if end if end do end do end do nullify (OperA,OperB,OperC,OperD) end subroutine OpMaskRel !******************************************************************************* ! command = 30 subroutine OpGrimKay call PointOperA call PointOperB TotA=0 ; TotB=0 ; TotC=0 do XYear = 1, YearN do XMonth = 1, MonthN do XBox = 1, BoxN if (OperB(XYear,XMonth,XBox).NE.MissVal) TotB=TotB+1 OperA(XYear,XMonth,XBox) = OpAwithB (OperB(XYear,XMonth,XBox),CommKay(XComm),CommOp(XComm,4)) if (OperA(XYear,XMonth,XBox).NE.MissVal) TotA=TotA+1 end do end do end do nullify (OperA,OperB) end subroutine OpGrimKay !******************************************************************************* ! command = 31 subroutine OpGrimGrim call PointOperA call PointOperB call PointOperC TotA=0 ; TotB=0 ; TotC=0 do XYear = 1, YearN do XMonth = 1, MonthN do XBox = 1, BoxN if (OperB(XYear,XMonth,XBox).NE.MissVal) TotB=TotB+1 if (OperC(XYear,XMonth,XBox).NE.MissVal) TotC=TotC+1 OperA(XYear,XMonth,XBox) = OpAwithB (OperB(XYear,XMonth,XBox),OperC(XYear,XMonth,XBox),CommOp(XComm,4)) if (OperA(XYear,XMonth,XBox).NE.MissVal) TotA=TotA+1 end do end do end do nullify (OperA,OperB,OperC) end subroutine OpGrimGrim !******************************************************************************* ! command = 43 subroutine DumpGrim call PointOperA call SaveGrim (OperA,Grid,YearAD,Bounds,CommInfo(XComm,XExec),CommFile(XComm,XExec)," ",FileSuffix,NoZip=1) nullify (OperA) end subroutine DumpGrim !******************************************************************************* ! command = 51 subroutine InterpolateMissing call PointOperA call PointOperB OperA = OperB TotA=0 ; TotB=0 do XYear = 1, YearN do XMonth = 1, MonthN do XExe = 1, ExeN do XWye = 1, WyeN if (Grid(XExe,XWye).NE.MissVal) then if (OperB(XYear,XMonth,Grid(XExe,XWye)).EQ.MissVal) then OpTot = 0.0 ; OpEn = 0.0 do XExeNear = XExe-1,XExe+1 do XWyeNear = XWye-1,XWye+1 if (XExeNear.GE.1.AND.XExeNear.LE.ExeN.AND.XWyeNear.GE.1.AND.XWyeNear.LE.WyeN) then if (Grid(XExeNear,XWyeNear).NE.MissVal) then if (OperB(XYear,XMonth,Grid(XExeNear,XWyeNear)).NE.MissVal) then OpTot = OpTot + OperB(XYear,XMonth,Grid(XExeNear,XWyeNear)) OpEn = OpEn + 1.0 end if end if end if end do end do if (OpEn.GE.1) then OperA(XYear,XMonth,Grid(XExe,XWye)) = OpTot / OpEn TotA = TotA + 1 else TotB = TotB + 1 end if end if end if end do end do end do end do nullify (OperA,OperB) end subroutine InterpolateMissing !******************************************************************************* subroutine PointOperA if (CommOp(XComm,2).EQ.1) OperA => Grim1 if (CommOp(XComm,2).EQ.2) OperA => Grim2 end subroutine PointOperA !******************************************************************************* subroutine PointOperB if (CommOp(XComm,3).EQ.1) OperB => Grim1 if (CommOp(XComm,3).EQ.2) OperB => Grim2 end subroutine PointOperB !******************************************************************************* subroutine PointOperC if (CommOp(XComm,5).EQ.1) OperC => Grim1 if (CommOp(XComm,5).EQ.2) OperC => Grim2 end subroutine PointOperC !******************************************************************************* subroutine PointOperD if (CommOp(XComm,4).EQ.1) OperD => Grim1 if (CommOp(XComm,4).EQ.2) OperD => Grim2 end subroutine PointOperD !******************************************************************************* subroutine Finalise deallocate (CommOp,CommKay,CommFile,CommInfo, stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Finalise: Deallocation failure: Comm* #####" deallocate (Grid,YearAD,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Finalise: Deallocation failure: Basic #####" if (associated(Grim1)) deallocate (Grim1, stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Finalise: Deallocation failure: Grim1 #####" if (associated(Grim2)) deallocate (Grim2, stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: Finalise: Deallocation failure: Grim2 #####" print* end subroutine Finalise !******************************************************************************* end program TYNsc20SeqGrim tynsc20unpack_pik.f000600 125236 002424 00000022472 07617514637 014377 0ustar00f709762cru000000 000000 ! 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 ! 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 ; Co