! opnormals.f90
! f90 program written by Tim Mitchell on 04.10.02
! last modified on 04.10.02
! program to manipulate CRU ts normals files (.nrm)
! f90 -o ./../obs/opnormals filenames.f90 time.f90 grimfiles.f90 crutsfiles.f90 convert.f90 
!	basicfun.f90 ./../obs/opnormals.f90

program OpNormals

use FileNames
use Time
use GrimFiles
use CRUtsFiles
use Convert
use BasicFun

implicit none

real, pointer, dimension (:)			:: ALat,ALon,AElv,BLat,BLon,BElv,CLat,CLon,CElv

integer, pointer, dimension (:,:)		:: NormalsA,NormalsB,NormalsC
integer, pointer, dimension (:,:)		:: StnInfoA,StnInfoB,StnInfoC
integer, pointer, dimension (:)			:: AStn,BStn,CStn
integer, pointer, dimension (:)			:: CfromA,CfromB
integer, pointer, dimension (:)			:: ANmlYr0,ANmlYr1,ANmlSrc,ANmlInc
integer, pointer, dimension (:)			:: BNmlYr0,BNmlYr1,BNmlSrc,BNmlInc
integer, pointer, dimension (:)			:: CNmlYr0,CNmlYr1,CNmlSrc,CNmlInc

character (len=20), pointer, dimension (:)	:: StnNameA,StnNameB,StnNameC
character (len=13), pointer, dimension (:)	:: StnCtyA,StnCtyB,StnCtyC
character (len=09), pointer, dimension (:)	:: StnLocalA,StnLocalB,StnLocalC

real, parameter 	:: MissVal = -999.0
integer, parameter 	:: DataMissVal = -9999

real :: FactorA,FactorB,Value,Constant

integer :: ReadStatus,AllocStat
integer :: QMenu,QConversion,QRestrict,QOperation
integer :: XMonth,XAStn,XBStn,XCStn
integer :: NMonth,NAStn,NBStn,NCStn
integer :: VarSuffBeg,Log1,Log2,Log3,Log4,NormalBeg,NormalEnd,AScore,BScore
integer :: Min,MinImpose,Max,MaxImpose

character (len=80) :: GivenFile,LoadFileA,LoadFileB,SaveFileA,SaveFileB,SaveFileC
character (len=04) :: VarSuffA,VarSuffB

call Intro
call Menu
if (QMenu.EQ. 1) call ConvertHdr
if (QMenu.EQ.11) call ConvertVar
if (QMenu.EQ.12) call MergeTwo
if (QMenu.EQ.18) call AopConstant
if (QMenu.EQ.80) call FilterYrs
call Finish

contains

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

subroutine Intro

open (99,file="/cru/mikeh1/f709762/scratch/log/log-normals.dat",status="replace",action="write")

print*
print*, "  > ***** OpNormals: manipulates .nrm files *****"
print*

NMonth = 12

end subroutine Intro

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

subroutine Menu

print*, "  > Select an operation to perform (0=list): "
do
	read (*,*,iostat=ReadStatus), QMenu
	
	if (QMenu.EQ.0) then
	  print*, "  >  1. Convert into a .hdr file."
	  print*, "  > 11. Convert variable."
	  print*, "  > 12. Merge two .nrm files."
	  print*, "  > 18. A.k --> B"
	  print*, "  > 80. Filter unlikely years. [ANTIQUE]"
	end if
	
	if (ReadStatus.LE.0.AND.QMenu.NE.0) exit
end do    

end subroutine Menu

!*******************************************************************************
!  1. convert .nrm to .hdr

subroutine ConvertHdr

print*, "  > Select the .nrm file to load: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".nrm")

print*, "  > Select the .hdr file to save: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,".hdr")

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,Lat=ALat,Lon=ALon,Elv=AElv, &
		    CallFile=LoadFileA,Normals=NormalsA,NmlYr0=ANmlYr0,NmlYr1=ANmlYr1, &
		    NmlSrc=ANmlSrc,NmlInc=ANmlInc)
call MakeStnInfo (StnInfoB,AStn,ALat,ALon,AElv,1,YearAD0=ANmlYr0,YearAD1=ANmlYr1)
call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileA,HeadOnly=1)

end subroutine ConvertHdr

!*******************************************************************************
! 11. convert variable

subroutine ConvertVar

print*, "  > Select the conversion (0=list): "
do
	read (*,*,iostat=ReadStatus), QConversion
	
	if (QConversion.EQ.0) then
	  print*, "  >  1. .snh to .spc"
	  print*, "  >  2. .spc to .cld"
	end if
	
	if (ReadStatus.LE.0.AND.QConversion.NE.0) exit
end do    

if      (QConversion.EQ.1) then
  VarSuffA=".snh" ; VarSuffB=".spc" ; FactorA=0.1 ; FactorB=0.1
else if (QConversion.EQ.2) then
  VarSuffA=".spc" ; VarSuffB=".cld" ; FactorA=0.1 ; FactorB=0.1
end if

print*, "  > Select the .nrm file to load: "
do
	VarSuffBeg = 0
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0) VarSuffBeg = index(GivenFile,VarSuffA)
	if (VarSuffBeg.GT.0) exit
end do
LoadFileA = LoadPath (GivenFile,".nrm")

print*, "  > Select the .nrm file to save: "
do
	VarSuffBeg = 0
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0) VarSuffBeg = index(GivenFile,VarSuffB)
	if (VarSuffBeg.GT.0) exit
end do
SaveFileB = SavePath (GivenFile,".nrm")

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,Lat=ALat,Lon=ALon,Elv=AElv, &
		    CallFile=LoadFileA,Normals=NormalsA,NmlYr0=ANmlYr0,NmlYr1=ANmlYr1, &
		    NmlSrc=ANmlSrc,NmlInc=ANmlInc)

NAStn=size(AStn,1) ; Log1=0 ; Log2=0
do XMonth = 1, NMonth
  do XAStn = 1, NAStn
    if (NormalsA(XMonth,XAStn).NE.DataMissVal) then
      Value = real(NormalsA(XMonth,XAStn))*FactorA      
      
      if (QConversion.EQ.1) Value=SunPCfromHours(Value,ALat(XAStn),XMonth)
      if (QConversion.EQ.2) Value=100.0-Value
      
      if (Value.NE.DataMissVal) then
      	NormalsA(XMonth,XAStn)=nint(Value/FactorB)
      else
        Log2=Log2+1
      end if
    else
      Log1=Log1+1
    end if
  end do
end do
Log3=(NMonth*NAStn) ; Log4=((NMonth*NAStn)-Log1-Log2)
print "(a)", "   > Totals: checked, original missing, didn't convert, and valid: "
print "(4(a,i6))", "              ",Log3,"            ",Log1,"          ",Log2, &
			"     ",Log4

call MakeStnInfo (StnInfoB,AStn,ALat,ALon,AElv,1,YearAD0=ANmlYr0,YearAD1=ANmlYr1)
call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileB, &
              Normals=NormalsA,NmlSrc=ANmlSrc,NmlInc=ANmlInc)

end subroutine ConvertVar

!*******************************************************************************
! 12. merge 2 .nrm files

subroutine MergeTwo

print*, "  > Select the MASTER .nrm file to load: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".nrm")

print*, "  > Select the SUPPLEMENTARY .nrm file to load: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileB = LoadPath (GivenFile,".nrm")

print*, "  > Select the .nrm file to save: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileC = SavePath (GivenFile,".nrm")

print*, "  > Enter the normals period to target (beg,end): "
do
	read (*,*,iostat=ReadStatus), NormalBeg,NormalEnd
	if (ReadStatus.LE.0.AND.NormalEnd.GE.NormalBeg) exit
end do

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,Lat=ALat,Lon=ALon,Elv=AElv, &
		    CallFile=LoadFileA,Normals=NormalsA,NmlYr0=ANmlYr0,NmlYr1=ANmlYr1, &
		    NmlSrc=ANmlSrc,NmlInc=ANmlInc)
call LoadCTS (StnInfoB,StnLocalB,StnNameB,StnCtyB,Code=BStn,Lat=BLat,Lon=BLon,Elv=BElv, &
		    CallFile=LoadFileB,Normals=NormalsB,NmlYr0=BNmlYr0,NmlYr1=BNmlYr1, &
		    NmlSrc=BNmlSrc,NmlInc=BNmlInc)
call CombineStn (AStn,BStn,CStn,CfromA,CfromB)		! see CRUtsFiles
NAStn=size(AStn) ; NBStn=size(BStn) ; NCStn=size(CStn)

allocate (CStn(NCStn), CLat(NCStn), CLon(NCStn), CElv(NCStn), &
	  CNmlYr0(NCStn), CNmlYr1(NCStn), StnLocalC(NCStn), StnNameC(NCStn), StnCtyC(NCStn), &
	  NormalsC(NMonth,NCStn), CNmlSrc(NCStn), CNmlInc(NCStn), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: MergeTwo: Allocation failure: C* #####"

write (99,*), "STATIONS OVERWRITTEN:"
Log1=0 ; Log2=0
do XCStn = 1, NCStn
  if      (CfromA(XCStn).NE.-999.AND.CfromB(XCStn).NE.-999) then
    if (BNmlYr0(CfromB(XCStn)).NE.9999.AND.BNmlYr1(CfromB(XCStn)).NE.9999) then
      if (BNmlYr0(CfromB(XCStn)).GE.NormalBeg.AND.BNmlYr1(CfromB(XCStn)).LE.NormalEnd) then
        if (ANmlYr0(CfromA(XCStn)).GE.NormalBeg.AND.ANmlYr1(CfromA(XCStn)).LE.NormalEnd) then
          AScore = (ANmlYr1(CfromA(XCStn))-ANmlYr0(CfromA(XCStn))+1)*ANmlInc(CfromA(XCStn))
          BScore = (BNmlYr1(CfromB(XCStn))-BNmlYr0(CfromB(XCStn))+1)*BNmlInc(CfromB(XCStn))
          if (AScore.GE.BScore) then
            CfromB(XCStn) = -999
          else
            write (99,"(i7,x,a3,2i5,i4,12i5)"), CStn(XCStn), "Asc",ANmlYr0(CfromA(XCStn)),&
            	ANmlYr1(CfromA(XCStn)),ANmlInc(CfromA(XCStn)),(NormalsA(XMonth,CfromA(XCStn)),XMonth=1,NMonth)
            write (99,"(7x,x,a3,2i5,i4,12i5)"), "Bsc",BNmlYr0(CfromB(XCStn)),&
            	BNmlYr1(CfromB(XCStn)),BNmlInc(CfromB(XCStn)),(NormalsB(XMonth,CfromB(XCStn)),XMonth=1,NMonth)
            CfromA(XCStn) = -999
          end if
        else
          write (99,"(i7,x,a3,2i5,i4,12i5)"), CStn(XCStn), "Ayr",ANmlYr0(CfromA(XCStn)),&
            	ANmlYr1(CfromA(XCStn)),ANmlInc(CfromA(XCStn)),(NormalsA(XMonth,CfromA(XCStn)),XMonth=1,NMonth)
          write (99,"(7x,x,a3,2i5,i4,12i5)"), "Byr",BNmlYr0(CfromB(XCStn)),&
            	BNmlYr1(CfromB(XCStn)),BNmlInc(CfromB(XCStn)),(NormalsB(XMonth,CfromB(XCStn)),XMonth=1,NMonth)
          CfromA(XCStn) = -999
        end if
      else
        CfromB(XCStn) = -999
      end if
    else
      CfromB(XCStn) = -999
    end if
  end if
  
  if (CfromA(XCStn).NE.-999.AND.CfromB(XCStn).EQ.-999) then
    Log1=Log1+1
    do XMonth = 1, NMonth
      NormalsC(XMonth,XCStn) = NormalsA(XMonth,CfromA(XCStn))
    end do    
  else if (CfromA(XCStn).EQ.-999.AND.CfromB(XCStn).NE.-999) then
    Log2=Log2+1
    do XMonth = 1, NMonth
      NormalsC(XMonth,XCStn) = NormalsB(XMonth,CfromB(XCStn))
    end do
  end if
  
  if      (CfromA(XCStn).NE.-999.AND.CfromB(XCStn).EQ.-999) then
    CStn(XCStn)=AStn(CfromA(XCStn)) ; CLat(XCStn)=ALat(CfromA(XCStn))
    CLon(XCStn)=ALon(CfromA(XCStn)) ; CElv(XCStn)=AElv(CfromA(XCStn))
    CNmlYr0(XCStn)=ANmlYr0(CfromA(XCStn)) ; CNmlYr1(XCStn)=ANmlYr1(CfromA(XCStn))
    StnLocalC(XCStn)=StnLocalA(CfromA(XCStn)) ; StnNameC(XCStn)=StnNameA(CfromA(XCStn))
    StnCtyC(XCStn)=StnCtyA(CfromA(XCStn)) ; CNmlSrc(XCStn)=ANmlSrc(CfromA(XCStn))
    CNmlInc(XCStn)=ANmlInc(CfromA(XCStn))    
  else if (CfromA(XCStn).EQ.-999.AND.CfromB(XCStn).NE.-999) then
    CStn(XCStn)=BStn(CfromB(XCStn)) ; CLat(XCStn)=BLat(CfromB(XCStn))
    CLon(XCStn)=BLon(CfromB(XCStn)) ; CElv(XCStn)=BElv(CfromB(XCStn))
    CNmlYr0(XCStn)=BNmlYr0(CfromB(XCStn)) ; CNmlYr1(XCStn)=BNmlYr1(CfromB(XCStn))
    StnLocalC(XCStn)=StnLocalB(CfromB(XCStn)) ; StnNameC(XCStn)=StnNameB(CfromB(XCStn))
    StnCtyC(XCStn)=StnCtyB(CfromB(XCStn)) ; CNmlSrc(XCStn)=BNmlSrc(CfromB(XCStn))
    CNmlInc(XCStn)=BNmlInc(CfromB(XCStn))    
  end if
end do
print*, "  > Overwritten stations dumped to log."
print "(a,2i7)", "   > Stns: kept from master, added from supplement: ", Log1,Log2

call MakeStnInfo (StnInfoC,CStn,CLat,CLon,CElv,1,YearAD0=CNmlYr0,YearAD1=CNmlYr1)
call SaveCTS (StnInfoC,StnLocalC,StnNameC,StnCtyC,CallFile=SaveFileC, &
              Normals=NormalsC,NmlSrc=CNmlSrc,NmlInc=CNmlInc)

end subroutine MergeTwo

!*******************************************************************************
! 18. A.k --> B

subroutine AopConstant

print*, "  > Select the .nrm file to load: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".nrm")

print*, "  > Select the .nrm file to save: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,".nrm")

print*, "  > Select the operation (A/k=1, A*k=2, A-k=3, A+k=4, sqrt=5, **k=6, abs=7): "
do
	read (*,*,iostat=ReadStatus), QOperation
	if (ReadStatus.LE.0 .AND. QOperation.GE.1 .AND. QOperation.LE.7) exit
end do

Constant = -999.0
if (QOperation.NE.5.AND.QOperation.NE.7) then
  print*, "  > Select the constant: "
  do
	read (*,*,iostat=ReadStatus), Constant
	if (ReadStatus.LE.0) exit
  end do
end if

call EnquireRestrict

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,Lat=ALat,Lon=ALon,Elv=AElv, &
		    CallFile=LoadFileA,Normals=NormalsA,NmlYr0=ANmlYr0,NmlYr1=ANmlYr1, &
		    NmlSrc=ANmlSrc,NmlInc=ANmlInc)
NAStn = size(NormalsA,2)

do XAStn = 1, NAStn
    do XMonth = 1, 12
      if (NormalsA(XMonth,XAStn).NE.-9999) NormalsA(XMonth,XAStn) = &
      		OpAwithB (real(NormalsA(XMonth,XAStn)),Constant,QOperation)
      if (QRestrict.Eq.1.OR.QRestrict.Eq.3) then
        if (NormalsA(XMonth,XAStn).LT.Min) NormalsA(XMonth,XAStn) = MinImpose
      end if
      if (QRestrict.Eq.2.OR.QRestrict.Eq.3) then
        if (NormalsA(XMonth,XAStn).GT.Max) NormalsA(XMonth,XAStn) = MaxImpose
      end if
    end do
end do

call SaveCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileA, &
              Normals=NormalsA,NmlSrc=ANmlSrc,NmlInc=ANmlInc)

end subroutine AopConstant

!*******************************************************************************
! 80. Filter unlikely years. [ANTIQUE]

subroutine FilterYrs

print*, "  > Select the .nrm file to load: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
LoadFileA = LoadPath (GivenFile,".nrm")

print*, "  > Select the .nrm file to save: "
do
	read (*,*,iostat=ReadStatus), GivenFile
	if (ReadStatus.LE.0.AND.GivenFile.NE."") exit
end do
SaveFileA = SavePath (GivenFile,".nrm")

print*, "  > Operating..."
call LoadCTS (StnInfoA,StnLocalA,StnNameA,StnCtyA,Code=AStn,Lat=ALat,Lon=ALon,Elv=AElv, &
		    CallFile=LoadFileA,Normals=NormalsA,NmlYr0=ANmlYr0,NmlYr1=ANmlYr1, &
		    NmlSrc=ANmlSrc,NmlInc=ANmlInc)

NAStn = size(NormalsA,2)
do XAStn = 1, NAStn
  if      (ANmlYr0(XAStn).LT.1600.OR.ANmlYr0(XAStn).GT.2100) then
    ANmlYr0(XAStn) = 9999 ; ANmlYr1(XAStn) = 9999
  else if (ANmlYr1(XAStn).LT.1600.OR.ANmlYr1(XAStn).GT.2100) then
    ANmlYr0(XAStn) = 9999 ; ANmlYr1(XAStn) = 9999
  end if
end do

call MakeStnInfo (StnInfoB,AStn,ALat,ALon,AElv,1,YearAD0=ANmlYr0,YearAD1=ANmlYr1)
call SaveCTS (StnInfoB,StnLocalA,StnNameA,StnCtyA,CallFile=SaveFileA, &
              Normals=NormalsA,NmlSrc=ANmlSrc,NmlInc=ANmlInc)

end subroutine FilterYrs

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

subroutine EnquireRestrict

print*, "  > Impose min (=1), max (=2), both (=3), neither (=0):"
do
	read (*,*,iostat=ReadStatus), QRestrict
	if (ReadStatus.LE.0.AND.QRestrict.GE.0.AND.QRestrict.LE.3) exit
end do
if (QRestrict.EQ.1.OR.QRestrict.EQ.3) then
  print*, "  > Enter the minimum value to accept, and value to impose:"
  do
	read (*,*,iostat=ReadStatus), Min,MinImpose
	if (ReadStatus.LE.0) exit
  end do
end if
if (QRestrict.EQ.2.OR.QRestrict.EQ.3) then
  print*, "  > Enter the maximum value to accept, and value to impose:"
  do
	read (*,*,iostat=ReadStatus), Max,MaxImpose
	if (ReadStatus.LE.0) exit
  end do
end if

end subroutine EnquireRestrict

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

subroutine Finish

print*

close (99)

end subroutine Finish

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

end program OpNormals