! stangeneral.f90
! module in which various routines with general applicability are held
! contains: MakeSeasonKey, GetMonthLengths, CommonVecPer

module StanGeneral

implicit none

contains

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

subroutine MakeSeasonKey (YearAD,SeasonKey)

integer, pointer, dimension (:,:,:) 	:: SeasonKey 
integer, pointer, dimension (:,:) 	:: MonthLengths
integer, pointer, dimension (:) 	:: YearAD

real, parameter :: MissVal = -999.0

integer :: XYear, XMonth, XDay
integer :: Day0, Day1, Month0, Month1, DayBeg, DayEnd
integer :: YearN, MonthN, DayN
integer :: AllocStat, ReadStatus
integer :: Overlap

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

call GetMonthLengths (YearAD,MonthLengths)

DayN = 31 ; MonthN = 12

Day0 = 0 ; Month0 = 0
print*, "  > Enter the first day, month of the season: "
do
	read (*,*,iostat=ReadStatus), Day0, Month0
	if (ReadStatus.GT.0) print*, "  > Not integers. Try again."
	if (Month0.LT.1.OR.Month0.GT.MonthN) then
		print*, "  > Month out of range. Try again."
		Month0 = 0
	end if
	if (Day0.LT.1.OR.Day0.GT.DayN) then
		print*, "  > Day out of range. Try again."
		Day0 = 0
	end if
	if (ReadStatus.LE.0.AND.Day0.NE.0.AND.Month0.NE.0) exit
end do
      
Day1 = 0 ; Month1 = 0
print*, "  > Enter the last day, month of the season (overlaps OK): "
do
	read (*,*,iostat=ReadStatus), Day1, Month1
	if (ReadStatus.GT.0) print*, "  > Not integers. Try again."
	if (Month1.LT.1.OR.Month1.GT.MonthN) then
		print*, "  > Month out of range. Try again."
		Month1 = 0
	end if
	if (Day1.LT.1.OR.Day1.GT.DayN) then
		print*, "  > Day out of range. Try again."
		Day1 = 0
	end if
	if (ReadStatus.LE.0.AND.Day1.NE.0.AND.Month1.NE.0) exit
end do

Overlap = 0      
if      (Month1.LT.Month0) then
  Overlap = 1
else if (Month1.EQ.Month0) then
  if (Day1.LT.Day0) Overlap = 1
end if

YearN = size (YearAD)

allocate (SeasonKey(YearN,MonthN,DayN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SeasonKey: Allocation failure #####"
SeasonKey = MissVal

if (Overlap.EQ.0) then		! no overlap
  do XYear = 1, YearN
    do XMonth = Month0, Month1
      DayBeg = 1 ; DayEnd = DayN
      if (XMonth.EQ.Month0) DayBeg = Day0
      if (XMonth.EQ.Month1) DayEnd = Day1
      
      do XDay = DayBeg, DayEnd
        if (MonthLengths(XYear,XMonth).GE.XDay) SeasonKey(XYear,XMonth,XDay) = XYear
      end do
    end do
  end do
end if

if (Overlap.EQ.1) then		! overlap
  do XYear = 1, (YearN-1)
    do XMonth = Month0, MonthN
      DayBeg = 1 ; DayEnd = DayN
      if (XMonth.EQ.Month0) DayBeg = Day0
      
      do XDay = DayBeg, DayEnd
        if (MonthLengths(XYear,XMonth).GE.XDay) SeasonKey(XYear,XMonth,XDay) = XYear
      end do
    end do

    do XMonth = 1, Month1
      DayBeg = 1 ; DayEnd = DayN
      if (XMonth.EQ.Month1) DayEnd = Day1
      
      do XDay = DayBeg, DayEnd
        if (MonthLengths((XYear+1),XMonth).GE.XDay) SeasonKey((XYear+1),XMonth,XDay) = XYear
      end do
    end do
  end do
end if

deallocate (MonthLengths, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SeasonKey: Deallocation failure #####"

end subroutine MakeSeasonKey

!*******************************************************************************
! 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

!*******************************************************************************
! 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

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

end module StanGeneral