! rawtogrim.f90 ! f90 main program (gridtogrim.f90) written on 23.03.01 by Tim Mitchell ! last modification on 17.12.01 ! now rewritten to allow for .glo loading (23.01.03) ! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz ! -o ./../grim/rawtogrim time.f90 filenames.f90 grimfiles.f90 ! glofiles.f90 ./../grim/rawtogrim.f90 2> /tyn1/tim/scratch/stderr.txt ! loads data file(s) into standard grim files program RawToGrim use Time use FileNames use GrimFiles use GloFiles implicit none real, pointer, dimension (:,:,:) :: Data real, pointer, dimension (:,:) :: MonthData real, pointer, dimension (:) :: RowReal integer, pointer, dimension (:) :: RowInt, YearAD, ReGrid integer, pointer, dimension (:,:) :: Grid, RefGrid character (len=80), pointer, dimension (:,:) :: RawFile character (len=80), pointer, dimension (:) :: OutFile character (len= 3), pointer, dimension (:) :: MonthsNames real, dimension (4) :: Bounds character (len=3), dimension (12),target :: MonthsFirstCap = & (/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'/) character (len=3), dimension (12),target :: MonthsNoCap = & (/'jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec'/) character (len=2), dimension (12) :: MonthsNumeral = & (/'01','02','03','04','05','06','07','08','09','10','11','12'/) real, parameter :: MissVal = -999.0 character (len=80), parameter :: Blank = "" real :: FileMissVal,Multiplier,Fraction integer :: LongN,LatN,DataN, RegN, FileN,RowN,ColN,CellN, BlockN, YearN,MonthN, RawFileN, RefBoxN integer :: FileHeadN,MonthHeadN,YearHeadN integer :: AllocStat,ReadStatus integer :: DateCol,Tot0,Tot1 integer :: XRow,XCol,XCell, XReg, XLong,XLat, XFile, XBlock, XFileHead,XMonthHead,XYearHead, XYear,XMonth integer :: XRawFile,XBound integer :: StringLen, PosChar, Col0Cell integer :: Cell0,Cell1, YearAD0,YearAD1 integer :: QIntReal,QDateGreen,QZip,QDataGrim,QNorthSouth,QSource integer :: SubLen, SubBeg, SuffixBeg,SuffixEnd,SuffixLen, GrimSubBeg, YearSubLen,YearSubBeg,YearSubEnd integer :: MonthSubLen,MonthSubBeg,MonthSubEnd integer :: LoadFileLen,CurrentPos,CurrentCod,TotMiss character (len=1) :: CharIntReal character (len=4) :: SaveSuffix, RefSuffix character (len=80) :: OrigSub, ThisSub, YearSub, SpecSub character (len=40) :: LineFormat character (len=80) :: GridFilePath,GivenFile,FilePath,AutoLoadPath,AutoSavePath,RefFile,GrimFile,LoadFile character (len=80) :: Info, InfoItem character (len=80) :: Trash character (len=20) :: TextFound logical, parameter :: True = .TRUE. !******************************************************************************* ! main program call Initialise if (QDataGrim.EQ.1) then call DomainFromData else call DomainFromGrim end if call GetInfoLine call GetDimensions if (QSource.EQ.1) call DataSpec call FirstExec if (FileN.GT.1) call AutoSpec if (QSource.EQ.1) call LoadDataFreeForm if (QSource.EQ.2) call LoadDataGlo call Finalise contains !******************************************************************************* ! initialise ! in this program we take raw data files organised as sequential spatial grids ! (one grid per month, with 12 months (Jan,Feb,...,Dec) * YearN per file, or one month per file) ! we store the data in standard grim files subroutine Initialise open (99,file="/tyn1/tim/scratch/log-gridtogrim.dat",status="replace",action="write") print* print*, " > ***** RawToGrim : store gridded monthly data *****" print* print*, " > Load free-form (=1) or from .glo (=2) ?" do read (*,*,iostat=ReadStatus), QSource if (ReadStatus.LE.0.AND.QSource.GE.1.AND.QSource.LE.2) exit end do if (QSource.EQ.1) then print*, " > Construct domain from data (=1) or from a grim file (=2) ? " do read (*,*,iostat=ReadStatus), QDataGrim if (ReadStatus.LE.0.AND.QDataGrim.GE.1.AND.QDataGrim.LE.2) exit end do else QDataGrim = 2 end if end subroutine Initialise !******************************************************************************* ! specify reference domain = RefGrid (LongN,LatN), with RefBoxN valid cells subroutine DomainFromData print*, " > Select the number of longitude, latitude cells: " do read (*,*,iostat=ReadStatus), LongN, LatN if (ReadStatus.LE.0.AND.LongN.GT.0.AND.LatN.GT.0) exit end do print*, " > Select the bounds of the grid (not cell centres). " print*, " > Enter min,max longitude, min,max latitude: " do read (*,*,iostat=ReadStatus), Bounds(1),Bounds(2),Bounds(3),Bounds(4) if (Bounds(1).LT.-182.OR.Bounds(2).GT.182.OR.Bounds(2).LE.Bounds(1)) then print*, " > Unacceptable longitudes. Try again." ReadStatus = 1 else if (Bounds(3).LT. -92.OR.Bounds(4).GT. 92.OR.Bounds(4).LE.Bounds(3)) then print*, " > Unacceptable latitudes. Try again." ReadStatus = 1 end if if (ReadStatus.LE.0) exit end do DataN = LongN * LatN RefBoxN = DataN Info = "" allocate (RefGrid (LongN,LatN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: DomainFromData: Allocation failure #####" RefGrid = MissVal XCell = 0 do XLong = 1, LongN do XLat = 1, LatN XCell = XCell + 1 RefGrid(XLong,XLat) = XCell end do end do end subroutine DomainFromData !******************************************************************************* ! specify reference domain = RefGrid (LongN,LatN), with RefBoxN valid cells subroutine DomainFromGrim print*, " > Enter the template grim filepath: " do read (*,*,iostat=ReadStatus), RefFile if (ReadStatus.LE.0.AND.RefFile.NE."") exit end do call LoadGrim (Data,RefGrid,YearAD,Bounds,Info,RefFile," ",RefSuffix) LongN = size (RefGrid,1) ; LatN = size (RefGrid,2) ; DataN = LongN * LatN RefBoxN = size (Data,3) Info = "" print "(a,3i8)", " > Grid dimensions and domain size: ", LongN,LatN,RefBoxN print "(a,4f8.2)", " > Grid bounds: ", (Bounds(XBound),XBound=1,4) deallocate (Data,YearAD,stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: DomainFromGrim: Deallocation failure #####" end subroutine DomainFromGrim !******************************************************************************* subroutine GetInfoLine print*, " > Name the grid: " do read (*,*,iostat=ReadStatus), InfoItem if (ReadStatus.LE.0.AND.InfoItem.NE."") exit end do Info = trim(InfoItem) print*, " > Name the domain: " do read (*,*,iostat=ReadStatus), InfoItem if (ReadStatus.LE.0.AND.InfoItem.NE."") exit end do Info = trim(Info) // " " // trim(InfoItem) print*, " > Name the scenario: " do read (*,*,iostat=ReadStatus), InfoItem if (ReadStatus.LE.0.AND.InfoItem.NE."") exit end do Info = trim(Info) // " " // trim(InfoItem) end subroutine GetInfoLine !******************************************************************************* subroutine GetDimensions print*, " > Select the first, last years to upload: " do read (*,*,iostat=ReadStatus), YearAD0, YearAD1 if (ReadStatus.LE.0.AND.YearAD1.GE.YearAD0) exit end do YearN = YearAD1 - YearAD0 + 1 MonthN = 12 print*, " > Select the number of upload executions: " do read (*,*,iostat=ReadStatus), FileN if (ReadStatus.LE.0) exit end do ! allocate Grid (LongN,LatN), & ! Grid = MissVal allocate (ReGrid(DataN), & YearAD(YearN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: GetDimensions: Allocation failure #####" ReGrid = MissVal do XYear = 1, YearN YearAD(XYear) = YearAD0 + XYear - 1 end do end subroutine GetDimensions !******************************************************************************* ! specify first execution subroutine FirstExec if (FileN.GT.1) print*, " > Specify the first execution. " if (QSource.EQ.1) then print "(a,i4,a4,i4,a2)", " > Enter the number of raw files (1 or ", YearN, " or ", (YearN*12), "):" do read (*,*,iostat=ReadStatus), RawFileN if (ReadStatus.LE.0.AND.RawFileN.NE.1.AND.RawFileN.NE.YearN.AND.RawFileN.NE.(YearN*12)) then print*, " > That number is not allowed. Try again." ReadStatus = 1 end if if (ReadStatus.LE.0) exit end do else RawFileN = (YearN*12) end if allocate (RawFile(FileN,RawFileN),& OutFile(FileN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: AutoSpec: Allocation failure #####" RawFile = "missing" ; OutFile = "missing" if (RawFileN.EQ.1) print*, " > Enter the filepath of the raw data file: " if (RawFileN.GT.1) print*, " > Enter the filepath of the first raw data file: " do read (*,*,iostat=ReadStatus), GivenFile if (ReadStatus.LE.0.AND.GivenFile.NE."") exit end do RawFile(1,1) = LoadPath (GivenFile," ") if (RawFileN.EQ.YearN) then TextFound = GetTextFromInt (YearAD0) YearSub = trim(adjustl(TextFound)) YearSubLen = len_trim(YearSub) YearSubBeg = index(RawFile(1,1),YearSub(1:YearSubLen)) YearSubEnd = YearSubBeg + YearSubLen - 1 if (YearSubBeg.GT.0) then do XRawFile = 2, RawFileN GivenFile = RawFile(1,1) TextFound = GetTextFromInt (YearAD0+XRawFile-1) GivenFile (YearSubBeg:YearSubEnd) = trim(adjustl(TextFound)) RawFile(1,XRawFile) = GivenFile end do end if else if (RawFileN.EQ.(YearN*12)) then TextFound = GetTextFromInt (YearAD0) YearSub = trim(adjustl(TextFound)) YearSubLen = len_trim(YearSub) if (YearN.EQ.1) then YearSubBeg = MissVal ; YearSubEnd = MissVal else YearSubBeg = index(RawFile(1,1),YearSub(1:YearSubLen)) YearSubEnd = YearSubBeg + YearSubLen - 1 end if MonthSubBeg = index(RawFile(1,1),'jan') if (MonthSubBeg.EQ.0) then MonthSubBeg = index(RawFile(1,1),'Jan') if (MonthSubBeg.EQ.0) then MonthSubBeg = index(RawFile(1,1),'01') if (MonthSubBeg.EQ.0) then MonthSubBeg = 0 else MonthSubLen = 2 end if else MonthSubLen = 3 MonthsNames => MonthsFirstCap end if else MonthSubLen = 3 MonthsNames => MonthsNoCap end if MonthSubEnd = MonthSubBeg + MonthSubLen - 1 if (YearSubBeg.NE.0.AND.MonthSubBeg.GT.0) then do XYear = 1, YearN do XMonth = 1, MonthN GivenFile = RawFile(1,1) if (YearSubBeg.NE.MissVal) then TextFound = GetTextFromInt (YearAD0+XYear-1) GivenFile (YearSubBeg:YearSubEnd) = trim(adjustl(TextFound)) end if if (MonthSubLen.EQ.2) then GivenFile (MonthSubBeg:MonthSubEnd) = MonthsNumeral (XMonth) else if (MonthSubLen.EQ.3) then GivenFile (MonthSubBeg:MonthSubEnd) = MonthsNames (XMonth) end if XRawFile = ((XYear-1)*12) + XMonth RawFile(1,XRawFile) = GivenFile end do end do end if end if print*, " > Enter the filepath of the grim file to save: " do read (*,*,iostat=ReadStatus), GivenFile if (ReadStatus.LE.0.AND.GivenFile.NE."") exit end do call ReviewCall (GivenFile," ",AutoSavePath,SaveSuffix,2) ! checks for file/suffix consistency OutFile(1) = AutoSavePath end subroutine FirstExec !******************************************************************************* ! specify automatics subroutine AutoSpec print*, " > Enter the substring to vary between executions: " do read (*,*,iostat=ReadStatus), OrigSub if (ReadStatus.GT.0) then print*, " > Bad format. Try again." else if (OrigSub.EQ."") then print*, " > Blank not permitted. Try again." end if if (ReadStatus.LE.0.AND.OrigSub.NE."") exit end do SubLen = len(trim(OrigSub)) print*, " > Enter the substring in each execution, starting with no.2:" do XFile = 2, FileN do read (*,*,iostat=ReadStatus), SpecSub if (ReadStatus.GT.0) then print*, " > Bad format. Try again." else if (SpecSub.EQ."") then print*, " > Blank not permitted. Try again." ReadStatus = 1 ! else if (len_trim(SpecSub).NE.SubLen) then ! print*, " > The substring has a wrong length. Try again." ! ReadStatus = 1 end if if (ReadStatus.LE.0) exit end do do XRawFile = 1, RawFileN GivenFile = RawFile(1,XRawFile) SubBeg = index(GivenFile,OrigSub(1:SubLen)) ! start of orig sub ! GivenFile (SubBeg:(SubBeg+SubLen-1)) = adjustl(trim(SpecSub)) ! replace orig sub RawFile(XFile,XRawFile) = GivenFile(1:(SubBeg-1)) // trim(SpecSub) // & GivenFile((SubBeg+SubLen):80) end do GivenFile = OutFile(1) SubBeg = index(GivenFile,OrigSub(1:SubLen)) ! start of orig sub ! GivenFile (SubBeg:(SubBeg+SubLen-1)) = adjustl(trim(SpecSub)) ! replace orig sub OutFile(XFile) = GivenFile(1:(SubBeg-1)) // trim(SpecSub) // & GivenFile((SubBeg+SubLen):80) end do end subroutine AutoSpec !******************************************************************************* ! specify data structure in original files subroutine DataSpec print*, " > Enter the number of headers per file, year, month: " do read (*,*,iostat=ReadStatus), FileHeadN,YearHeadN,MonthHeadN if (ReadStatus.LE.0.AND.FileHeadN.GE.0.AND.YearHeadN.GE.0.AND.MonthHeadN.GE.0) exit end do print*, " > Enter the number of columns, rows in a monthly grid: " do read (*,*,iostat=ReadStatus), ColN, RowN if (ReadStatus.LE.0.AND.RowN.GE.1.AND.ColN.GE.1) exit end do allocate (RowInt (ColN), & RowReal(ColN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: DataSpec: Allocation failure #####" RowInt = MissVal RowReal = MissVal CellN = RowN * ColN if (CellN.EQ.DataN) then if (RowN.EQ.LatN.AND.ColN.EQ.LongN) then print*, " > Define the raw data structure." print*, " > Top-bottom: N-S (=1) or S-N (=2) ? " do read (*,*,iostat=ReadStatus), QNorthSouth if (ReadStatus.LE.0.AND.QNorthSouth.GE.1.AND.QNorthSouth.LE.2) exit end do print*, " > Left-right: W-E" print "(a,i4,a)", " > Enter the long cell (1-", ColN, ") for the first column in a monthly grid: " do read (*,*,iostat=ReadStatus), Col0Cell if (ReadStatus.LE.0.AND.Col0Cell.GE.1.AND.Col0Cell.LE.ColN) exit end do if (QNorthSouth.EQ.1) XLat = LatN + 1 if (QNorthSouth.EQ.2) XLat = 0 XCell = 0 do XRow = 1, RowN if (QNorthSouth.EQ.1) XLat = XLat - 1 if (QNorthSouth.EQ.2) XLat = XLat + 1 do XLong = Col0Cell, ColN XCell = XCell + 1 ! Grid (XLong,XLat) = XCell ReGrid (XCell) = RefGrid (XLong,XLat) end do if (Col0Cell.GT.1) then do XLong = 1, (Col0Cell-1) XCell = XCell + 1 ! Grid (XLong,XLat) = XCell ReGrid (XCell) = RefGrid (XLong,XLat) end do end if end do BlockN = 1 else if (ColN.EQ.LatN.AND.RowN.EQ.LongN) then print*, " > ##### ERROR: rotated section X not written #####" else ! The grid is actually split by long in file if (Fraction.EQ.floor(Fraction)) then ! but we imagine it as a single grid print*, " > ##### ERROR: multiple blocks per month not written #####" else print*, " > ##### ERROR: difficult blocks not written #####" end if end if else print*, " > ##### ERROR: The cells in file do not match with grid. #####" end if print*, " > Enter the data line format, as multiples of i,f,e: " do read (*,*,iostat=ReadStatus), LineFormat StringLen = len (trim(LineFormat)) PosChar = max (scan(LineFormat,'I'),scan(LineFormat,'i'),scan(LineFormat,'F'),scan(LineFormat,'f')) PosChar = max (scan(LineFormat,'E'),scan(LineFormat,'e'),PosChar) CharIntReal = "" if (PosChar.GT.0) CharIntReal = LineFormat(PosChar:PosChar) if (LineFormat(1:1).NE."(" .OR. LineFormat(StringLen:StringLen).NE.")") then ReadStatus = 99 print*, " > Format is unacceptable. Retry." else if (CharIntReal.EQ.'I'.OR.CharIntReal.EQ.'i') then QIntReal = 1 else if (CharIntReal.EQ.'F'.OR.CharIntReal.EQ.'f') then QIntReal = 2 else if (CharIntReal.EQ.'F'.OR.CharIntReal.EQ.'f') then QIntReal = 2 else ReadStatus = 99 print*, " > Format is unacceptable. Retry." end if end if if (ReadStatus.LE.0) exit end do print*, " > Enter the multiplier by which to convert data: " do read (*,*,iostat=ReadStatus), Multiplier if (ReadStatus.LE.0) exit end do print*, " > Enter the file missing value: " do read (*,*,iostat=ReadStatus), FileMissVal if (ReadStatus.LE.0) exit end do end subroutine DataSpec !******************************************************************************* ! load data (from .glo files) subroutine LoadDataGlo allocate (Data(YearN,MonthN,RefBoxN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadDataFreeForm: Allocation failure #####" do XFile = 1, FileN print "(2a)", " > Constructing: ", trim(OutFile(XFile)) XRawFile = 0 ; Data = MissVal Tot0 = 0 ; Tot1 = 0 do XYear = 1, YearN do XMonth = 1, MonthN XRawFile = XRawFile + 1 ! load monthly data grid (longn,latn) call LoadGloGrid (RawFile(XFile,XRawFile),MonthData,Silent=1) do XLong = 1, LongN ! transfer to main data grid do XLat = 1, LatN if (RefGrid(XLong,XLat).NE.MissVal) & Data(XYear,XMonth,RefGrid(XLong,XLat)) = MonthData(XLong,XLat) if (MonthData(XLong,XLat).NE.MissVal) then Tot0 = Tot0 + 1 else Tot1 = Tot1 + 1 end if end do end do deallocate (MonthData, stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadDataGlo: Dealloc failure #####" end do end do print "(a,2i10)", " > Valid, missing data: ",Tot0,Tot1 call SaveGrim (Data,RefGrid,YearAD,Bounds,Info,OutFile(XFile)," ",SaveSuffix) end do end subroutine LoadDataGlo !******************************************************************************* ! load data (free form files) subroutine LoadDataFreeForm allocate (Data(YearN,MonthN,RefBoxN), stat=AllocStat) if (AllocStat.NE.0) print*, " > ##### ERROR: LoadDataFreeForm: Allocation failure #####" Data = MissVal do XFile = 1, FileN XRawFile = 0 if (RawFileN.EQ.1) call StartUpFile do XYear = 1, YearN if (RawFileN.EQ.YearN.AND.RawFileN.NE.1) call StartUpFile if (RawFileN.NE.(YearN*12).AND.YearHeadN.GT.0) then ! load year headers do XYearHead = 1, YearHeadN read (2,"(a80)"), Trash end do end if do XMonth = 1, MonthN if (RawFileN.EQ.(YearN*12)) call StartUpFile if (MonthHeadN.GT.0) then ! load month headers do XMonthHead = 1, MonthHeadN read (2,"(a80)"), Trash end do end if XCell = 0 do XRow = 1, RowN ! load data by row,col into long vector if (QIntReal.EQ.1) then read (2,fmt=LineFormat), (RowInt (XCol), XCol=1,ColN) do XCol = 1, ColN XCell = XCell + 1 if (ReGrid(XCell).NE.MissVal) Data(XYear,XMonth,ReGrid(XCell)) = real (RowInt(XCol)) end do else if (QIntReal.EQ.2) then read (2,fmt=LineFormat), (RowReal(XCol), XCol=1,ColN) do XCol = 1, ColN XCell = XCell + 1 if (ReGrid(XCell).NE.MissVal) Data(XYear,XMonth,ReGrid(XCell)) = RowReal(XCol) end do end if end do if (RawFileN.EQ.(YearN*12)) call ShutDownFile end do if (RawFileN.EQ.YearN.AND.RawFileN.NE.1) call ShutDownFile end do if (RawFileN.EQ.1) call ShutDownFile TotMiss = 0 do XYear = 1, YearN do XMonth = 1, MonthN do XCell = 1, RefBoxN if (Data(XYear,XMonth,XCell).NE.FileMissVal) then Data(XYear,XMonth,XCell) = Data(XYear,XMonth,XCell) * Multiplier else Data(XYear,XMonth,XCell) = MissVal TotMiss = TotMiss + 1 end if end do end do end do print "(2(a,i12))", " > Missing: ", TotMiss, " and Total: ", (YearN*MonthN*RefBoxN) print "(2a)", " > Saving: ", trim(OutFile(XFile)) call SaveGrim (Data,RefGrid,YearAD,Bounds,Info,OutFile(XFile)," ",SaveSuffix) Data = MissVal end do end subroutine LoadDataFreeForm !******************************************************************************* ! get file unzipped, open, and the headers read subroutine StartUpFile XRawFile = XRawFile + 1 LoadFile = LoadPath (RawFile(XFile,XRawFile)," ") print "(2a)", " > Loading: ", trim(LoadFile) LoadFileLen = len_trim(LoadFile) if (LoadFileLen.GT.1.AND.LoadFile((LoadFileLen-1):LoadFileLen).EQ.".Z") then QZip = 2 ! file is zipped call system ('uncompress ' // LoadFile) LoadFile ((LoadFileLen-1):LoadFileLen) = " " ! change filename to the unzipped file else QZip = 1 ! file not zipped end if open (2,file=LoadFile,status="old",access="sequential",action="read",form="formatted") if (FileHeadN.GT.0) then ! load file headers do XFileHead = 1, FileHeadN read (2,"(a80)"), Trash end do end if end subroutine StartUpFile !******************************************************************************* ! get file shut, rezipped subroutine ShutDownFile close (2) if (QZip.EQ.2) call system ('compress ' // LoadFile // ' &') end subroutine ShutDownFile !******************************************************************************* ! finalise subroutine Finalise close (99) print* end subroutine Finalise !******************************************************************************* end program RawToGrim