! program to generate area weights
! written by Tim Mitchell
! last modified on 24.10.01
! pgf90 -Mstandard -Minfo -fast -Mscalarsse -Mvect=sse -Mflushz 
! 	-o ./../goglo/makeweighted filenames.f90 glofiles.f90 initialmod.f90 
!	./../goglo/makeweighted.f90

program MakeWeighted

use FileNames
use GloFiles
use InitialMod

implicit none

real, pointer, dimension (:)   :: LatWeights, Weights

integer, pointer, dimension (:,:) :: WorkMapIDLRaw, WorkMapIDLReg		
						! shape of IDL mapping array
integer, pointer, dimension (:)   :: WorkRegSizes, WorkMapRawReg	
						! region sizes + reg-->raw

character (len=20), pointer, dimension (:) :: WorkRegNames
						! names of individual regions

real, parameter :: MissVal = -999.0
real, parameter :: Pi = 3.1415927

real :: NorthLat, SouthLat, Weight, Increment, Initial, MidLatDeg, MidLatRad

integer :: WorkGridChosen,WorkLongN,WorkLatN,WorkDataN, workRegN
integer :: XLat, XLong, XReg, ReadStatus			

character (len=10) :: WorkGridTitle			
character (len=40) :: WorkRegTitle			
character (len=80) :: WorkGridFilePath, GloTitle, Blank			

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

Blank = ""

print*
print*, "  > ***** MakeWeighted : constructs a weights .glo file *****"
print*

call GridSelect (WorkGridChosen,WorkGridTitle,WorkLongN,WorkLatN,WorkDataN,WorkGridFilePath)
call RegSelect  (WorkGridChosen, WorkLongN, WorkLatN, WorkDataN, WorkMapIDLReg, WorkRegSizes, WorkRegNames, &
			WorkRegTitle, WorkRegN)

allocate (Weights (WorkRegN))

print "(a,i4)", "   > Total latitude rows: ", WorkLatN
print*, "  > Enter the southerly and northerly extremes of the range: "
do
	read (*,*,iostat=ReadStatus), SouthLat, NorthLat
	if (ReadStatus.LE.0.AND.SouthLat.LE.NorthLat) exit
end do

Increment = (NorthLat-SouthLat)/real(WorkLatN)
Initial   = SouthLat + (Increment/2.0)
print "(a,f10.4)", "   > Degrees per latitude row: ", Increment

do XLat = 1, WorkLatN
  MidLatDeg = Initial + (real(XLat-1)*Increment)
  MidLatRad = (MidLatDeg/90.0) * (Pi/2.0)  
  Weight = cos (MidLatRad)
  
  do XLong = 1, WorkLongN
    Weights (WorkMapIDLReg (XLong,XLat)) = Weight
  end do
end do

call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,Blank,Blank,Weights,WorkMapIDLReg)

deallocate (Weights)

print*
close (99)

end program MakeWeighted