! addnorm.f90
! f90 program written by Tim Mitchell on 19.09.03
! program to add normals to existing dtb using ref ts
! converted to use on beo1 with Portland Group compiler 
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
!	-o ./../cruts/addnorm time.f90 filenames.f90 grimfiles.f90 
!	crutsfiles.f90 gridops.f90 sortmod.f90 plainperm.f90 regress.f90 
!	ghcnrefiter.f90 ./../cruts/addnorm.f90 2> /tyn1/tim/scratch/stderr.txt

program AddNorm

use Time
use FileNames
use GrimFiles
use CRUtsFiles
use GridOps
use SortMod
use PlainPerm
use Regress
use GHCNrefIter

implicit none

logical, pointer, dimension (:,:)		:: Trust
logical, pointer, dimension (:)			:: GotRef,GetRef

real, pointer, dimension (:)			:: Lat,Lon,Elv

integer, pointer, dimension (:,:,:)		:: Data,DataSrc,RefTS
integer, pointer, dimension (:,:)		:: Norm,NormSrc
integer, pointer, dimension (:,:)		:: Info
integer, pointer, dimension (:)			:: CodeNew,YearAD
integer, dimension (12)				:: Min,Max

character (len=20), pointer, dimension (:) 	:: Name
character (len=13), pointer, dimension (:) 	:: Cty
character (len=09), pointer, dimension (:) 	:: Local

real, parameter		:: MissVal = -999.0

integer, parameter	:: DataMissVal = -9999
integer, parameter 	:: MultiSource 	= -888

logical :: Differ

real :: OpTot,OpEn,Factor,DecayDistance,MissThresh

integer :: ReadStatus,AllocStat
integer :: NYear,NMonth,NStn
integer :: XYear,XMonth,XStn
integer :: StampInt,SuffixBeg,SubBeg, ClimBeg,ClimEnd, YearAD0,YearAD1, NormMin

character (len=80) :: LoadFileDtb,LoadFileDts,SavefileDtb,SaveFileDts
character (len=80) :: GivenFile,Variable
character (len=20) :: Stamp20
character (len=10) :: CurrentTime,FileTime,Stamp10
character (len=08) :: CurrentDate,FileDate
character (len=04) :: VariableSuffix

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

call Intro
call Specify
call SelectVariable
call LoadUp
call GetNorm
call SaveDown

contains

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

subroutine GetNorm

allocate (RefTS(NYear,NMonth,NStn),Trust(NYear,NStn), &
	  GotRef(NStn),GetRef(NStn),stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### GetNorm: alloc #####"
RefTS=DataMissVal ; Trust=.TRUE. ; GotRef=.FALSE. ; GetRef=.TRUE.
ClimBeg=YearAD0-YearAD(1)+1 ; ClimEnd=YearAD1-YearAD(1)+1

print*, "  > Calculating normals..."
call MakeRefTS (RefTS,Data,Lat,Lon,GetRef,Trust,DecayDistance,Differ, &
		      .TRUE.,(ClimEnd-ClimBeg+1-NormMin),GotRef, &
		      NeedBeg=ClimBeg,NeedEnd=ClimEnd,Verbose=1)
print "(a,21x,a,i7,a)", "   > Normals are here","for ", count(GotRef), " stns"

do XStn=1,NStn
  if (GotRef(XStn)) then
    do XMonth=1,NMonth
      if (Norm(XMonth,XStn).EQ.DataMissVal) then
        OpTot=0 ; OpEn=0
        do XYear=ClimBeg,ClimEnd
          if (RefTS(XYear,XMonth,XStn).NE.DataMissVal) then
	    OpTot=OpTot+RefTS(XYear,XMonth,XStn)
	    OpEn=OpEn+1
	  end if
        end do
	
!	write (99,"(2i7,e16.5,i7)"), CodeNew(XStn),XMonth,OpTot,nint(OpEn)
	Norm   (XMonth,XStn)=nint(OpTot/OpEn)
        NormSrc(XMonth,XStn)=MultiSource
      end if
    end do
  end if
end do

deallocate (RefTS,Trust,GetRef,GotRef,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### GetNormals: dealloc #####"

end subroutine GetNorm

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

subroutine Intro

open (99,file="./../../../scratch/log-addnorm.dat",status="replace",action="write")

print*
print*, "  > ***** AddNorm: add normals to .dtb *****"
print*

NMonth=12 ; YearAD0=1961 ; YearAD1=1990

call date_and_time (CurrentDate, CurrentTime)
Stamp10  = CurrentDate(3:8) // CurrentTime(1:4) 
Stamp20  = adjustr(Stamp10)
StampInt = GetIntFromText(Stamp20)

end subroutine Intro

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

subroutine LoadUp

SubBeg = index(LoadFileDtb,"/",.TRUE.)
FileDate=LoadFileDtb((SubBeg+5):(SubBeg+10)) 
FileTime=LoadFileDtb((SubBeg+11):(SubBeg+14))

call LoadCTS (Info,Local,Name,Cty,Data=Data,YearAD=YearAD, &
		DtbNormals=Norm,Silent=1,CallFile=LoadFileDtb, &
		YearADMin=1961,YearADMax=1990)
NStn = size(Data,3) ; NYear = size(YearAD,1)
deallocate (Info,Local,Name,Cty,YearAD,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: LoadUp: Dealloc #####"

print "(a,2(a2,a1),a2,a4,a2,a1,a2,a,i7,a,i10,a)", "   > Database   file on ", &
		FileDate(5:6),".",FileDate(3:4),".",FileDate(1:2)," at ", &
		FileTime(1:2),":",FileTime(3:4)," has ",NStn," stns;", &
		count(Data.NE.DataMissVal), " data"

call LoadCTS (Info,Local,Name,Cty,Code=CodeNew, &
		Lat=Lat,Lon=Lon,Elv=Elv,DtbNormals=NormSrc,Data=DataSrc, &
		YearAD=YearAD,CallFile=LoadFileDts,Silent=1, &
		YearADMin=1961,YearADMax=1990)

end subroutine LoadUp

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

subroutine SaveDown

call SaveCTS (Info,Local,Name,Cty,CallFile=SaveFileDtb, &
		Data=Data,DtbNormals=Norm,YearAD=YearAD,Silent=1)
call SaveCTS (Info,Local,Name,Cty,CallFile=SaveFileDts, &
		Data=DataSrc,DtbNormals=NormSrc,YearAD=YearAD,Silent=1)
print*

end subroutine SaveDown

!*******************************************************************************
! select the climate variable, by suffix

subroutine SelectVariable

SubBeg = index(LoadFileDtb,"/",.TRUE.)
VariableSuffix = "." // LoadFileDtb((SubBeg+1):(SubBeg+3))

do
   if (VariableSuffix.EQ."    ") then
    do
      print*, "  > Identify the variable suffix (form:'.???'): "
      read (*,*,iostat=ReadStatus), VariableSuffix
      if (ReadStatus.LE.0) exit
    end do
   end if
  
   call CheckVariSuffix (VariableSuffix,Variable,Factor)
   if (Variable.EQ."") VariableSuffix = "    "
   if (Variable.NE."") print "(2a)", "   > The variable is: ", trim(Variable)
   if (Variable.NE."") exit
end do

if      (VariableSuffix.EQ.".pre") then
  DecayDistance = 450  ; Min=0 ; Max=99900 ; Differ=.FALSE.
else if (VariableSuffix.EQ.".tmp") then
  DecayDistance = 1200 ; Min=-2730 ; Max=9000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".dtr") then 
  DecayDistance = 750  ; Min=0 ; Max=10000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".cld") then 
  DecayDistance = 750  ; Min=0 ; Max=1000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".spc") then	  	! assumed, after cld
  DecayDistance = 750  ; Min=0 ; Max=1000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".vap") then 
  DecayDistance = 1200 ; Min=0 ; Max=10000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".rhm") then 	! assumed, after vap, tmp
  DecayDistance = 1200 ; Min=0 ; Max=1000 ; Differ=.TRUE.
else if (VariableSuffix.EQ.".wet") then 
  DecayDistance = 450  ; Min=0  ; Differ=.FALSE.
!  DecayDistance = 450  ; Min=0  ; Differ=.TRUE.
  Max=(/3100,2900,3100,3000,3100,3000,3100,3100,3000,3100,3000,3100/)
else if (VariableSuffix.EQ.".rd0") then 
  DecayDistance = 450  ; Min=0  ; Differ=.FALSE.
  Max=(/310,290,310,300,310,300,310,310,300,310,300,310/)
else
  print*, "  > @@@@@ ERROR: SelectVariable: no var @@@@@"
end if

end subroutine SelectVariable

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

subroutine Specify

print*, "  > Select the .dtb file to load:"
do
	read (*,*,iostat=ReadStatus), GivenFile
        SuffixBeg = index(GivenFile,".dtb")
	if (ReadStatus.LE.0.AND.SuffixBeg.GT.0) exit
end do
LoadFileDtb  = LoadPath (GivenFile,".dtb")
GivenFile=LoadFileDtb ; SuffixBeg = index(GivenFile,".dtb")
GivenFile = trim(GivenFile(1:(SuffixBeg-1))//".dts"//GivenFile((SuffixBeg+4):80))
LoadFileDts  = LoadPath (GivenFile,".dts")

print*, "  > Select the .dtb file to save:"
do
	read (*,*,iostat=ReadStatus), GivenFile
        SuffixBeg = index(GivenFile,".dtb")
	if (ReadStatus.LE.0.AND.SuffixBeg.GT.0) exit
end do
SaveFileDtb  = SavePath (GivenFile,".dtb")
GivenFile=SaveFileDtb ; SuffixBeg = index(GivenFile,".dtb")
GivenFile = trim(GivenFile(1:(SuffixBeg-1))//".dts"//GivenFile((SuffixBeg+4):80))
SaveFileDts  = SavePath (GivenFile,".dts")

print*, "  > Specify the missing percentage permitted: "
do
	read (*,*,iostat=ReadStatus), MissThresh
	if (ReadStatus.LE.0.AND.MissThresh.GE.0.AND.MissThresh.LE.100) exit
end do
NormMin = ceiling((100.0-MissThresh)*(YearAD1-YearAD0+1)/100.0)
print*, "  > Data required for a normal: ", NormMin	! @@@@@@@@@@@@@@@@@@@@@@@

end subroutine Specify

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

end program AddNorm