C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_clockstuff.F,v 1.31 2012/03/22 14:22:32 jmc Exp $
C $Name: $
#include "FIZHI_OPTIONS.h"
C-- File fizhi_clockstuff.F:
C-- Contents
C-- o SET_ALARM
C-- o GET_ALARM
C-- o ALARM (function)
C-- o ALARM2 (function)
C-- o ALARM2NEXT (function)
C-- o SET_TIME
C-- o GET_TIME
C-- o NSECF (function)
C-- o NHMSF (function)
C-- o NSECF2 (function)
C-- o FIXDATE
C-- o INTERP_TIME
C-- o TICK
C-- o TIC_TIME
C-- o NALARM (function)
C-- o NALARM2 (function)
C-- o INCYMD (function)
C-- o ASTRO
C-- o TIME_BOUND
C-- o TIME2FREQ2
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
subroutine SET_ALARM (tag,datein,timein,freq)
C***********************************************************************
C Purpose
C -------
C Utility to Set Internal Alarms
C
C Argument Description
C --------------------
C tag ....... Character String Tagging Alarm Process
C date ...... Begining Date for Alarm
C time ...... Begining Time for Alarm
C freq ...... Repeating Frequency Interval for Alarm
C
C***********************************************************************
implicit none
#include "EEPARAMS.h"
#include "chronos.h"
character*(*) tag
integer freq,datein,timein
C- functions:
INTEGER ILNBLNK
EXTERNAL
C- local variables:
integer myid
logical first,set
data first /.true./
integer n, iL
myid = 1
if(first) then
ntags = 1
tags(1) = tag
freqs(1) = freq
dates(1) = datein
times(1) = timein
iL = ILNBLNK(tag)
WRITE(standardMessageUnit,'(A,I8,A,I6.6,A,I10,2A)')
& ' Set Alarm for: ', datein, ' ', timein,
& ', with frequency: ', freq, ', and Tag: ',tag(1:iL)
else
set = .false.
do n=1,ntags
if(tag.eq.tags(n)) then
if( myid.eq.1 ) then
print *, 'Warning! Alarm has already been set for Tag: ',tag
print *, 'Changing Alarm Information:'
print *, 'Frequency: ',freqs(n),' (Old) ',freq,' (New)'
print *, ' Date0: ',dates(n),' (Old) ',datein,' (New)'
print *, ' Time0: ',times(n),' (Old) ',timein,' (New)'
endif
freqs(n) = freq
dates(n) = datein
times(n) = timein
set = .true.
endif
enddo
if(.not.set) then
ntags = ntags+1
if(ntags.gt.maxtag ) then
if( myid.eq.1 ) then
print *, 'Too many Alarms are Set!!'
print *, 'Maximum Number of Alarms = ',maxtag
endif
call MY_FINALIZE
call MY_EXIT (101)
endif
tags(ntags) = tag
freqs(ntags) = freq
dates(ntags) = datein
times(ntags) = timein
iL = ILNBLNK(tag)
WRITE(standardMessageUnit,'(A,I8,A,I6.6,A,I10,2A)')
& ' Set Alarm for: ', datein, ' ', timein,
& ', with frequency: ', freq, ', and Tag: ',tag(1:iL)
endif
endif
first = .false.
return
end
subroutine GET_ALARM (tag,datein,timein,freq,tleft)
C***********************************************************************
C Purpose
C -------
C Utility to Get Internal Alarm Information
C
C Input
C -----
C tag ....... Character String Tagging Alarm Process
C
C Output
C ------
C datein ...... Begining Date for Alarm
C timein ...... Begining Time for Alarm
C freq ........ Frequency Interval for Alarm
C tleft ....... Time Remaining (seconds) before Alarm is TRUE
C
C***********************************************************************
implicit none
character*(*) tag
integer freq,datein,timein,tleft
#include "chronos.h"
logical set,alarm
external
integer myid,n,nalarm,nsecf
myid = 1
set = .false.
do n=1,ntags
if (tag.eq.tags(n)) then
freq = freqs(n)
datein = dates(n)
timein = times(n)
if ( alarm(tag) ) then
tleft = 0
else
call GET_TIME (nymd,nhms)
tleft = nsecf(freq) - nalarm(freq,nymd,nhms,datein,timein )
endif
set = .true.
endif
enddo
if(.not.set) then
if( myid.eq.1 ) print *, 'Alarm has not been set for Tag: ',tag
freq = 0
datein = 0
timein = 0
tleft = 0
endif
return
end
LOGICAL FUNCTION ALARM (tag)
implicit none
character*(*) tag
#include "chronos.h"
integer datein,timein
integer n,nalarm
external
call GET_TIME (datein,timein)
alarm = .false.
do n=1,ntags
if( tags(n).eq.tag ) then
if( freqs(n).eq.0 ) then
alarm = (dates(n).eq.datein) .and. (times(n).eq.timein)
else
alarm = ( datein.gt.dates(n) .or.
& (datein.eq.dates(n) .and. timein.ge.times(n)) )
& .and. nalarm( freqs(n),datein,timein,dates(n),times(n) ).eq.0
endif
endif
enddo
return
end
LOGICAL FUNCTION ALARM2 (tag)
implicit none
character*(*) tag
#include "chronos.h"
integer datein,timein
integer n,nalarm2
external
call GET_TIME (datein,timein)
alarm2 = .false.
do n=1,ntags
if( tags(n).eq.tag ) then
if( freqs(n).eq.0 ) then
alarm2 = (dates(n).eq.datein) .and. (times(n).eq.timein)
else
alarm2 = ( datein.gt.dates(n) .or.
& (datein.eq.dates(n) .and. timein.ge.times(n)) )
& .and. nalarm2( freqs(n),datein,timein,dates(n),times(n) ).eq.0
endif
endif
enddo
return
end
LOGICAL FUNCTION ALARM2NEXT (tag,deltat)
implicit none
character*(*) tag
_RL deltat
#include "chronos.h"
integer datein,timein,ndt
integer dateminus,timeminus
integer n,nalarm2
external
ndt = int(deltat)
call GET_TIME (dateminus,timeminus)
datein = dateminus
timein = timeminus
call TICK(datein,timein,ndt)
alarm2next = .false.
do n=1,ntags
if( tags(n).eq.tag ) then
if( freqs(n).eq.0 ) then
alarm2next = (dates(n).eq.datein) .and. (times(n).eq.timein)
else
alarm2next = ( datein.gt.dates(n) .or.
& (datein.eq.dates(n) .and. timein.ge.times(n)) )
& .and. nalarm2( freqs(n),datein,timein,dates(n),times(n) ).eq.0
endif
endif
enddo
return
end
subroutine SET_TIME (datein,timein)
implicit none
integer datein,timein
#include "chronos.h"
integer myid
myid = 1
if( myid.eq.1 ) then
print *, 'Setting Clock'
print *, 'Date: ',datein
print *, 'Time: ',timein
endif
nymd = datein
nhms = timein
return
end
subroutine GET_TIME (datein,timein)
implicit none
integer datein,timein
#include "chronos.h"
datein = nymd
timein = nhms
return
end
function nsecf (nhms)
C***********************************************************************
C Purpose
C Converts NHMS format to Total Seconds
C
C***********************************************************************
implicit none
integer nhms, nsecf
nsecf = nhms/10000*3600 + mod(nhms,10000)/100*60 + mod(nhms,100)
return
end
function nhmsf (nsec)
C***********************************************************************
C Purpose
C Converts Total Seconds to NHMS format
C
C***********************************************************************
implicit none
integer nhmsf, nsec
nhmsf = nsec/3600*10000 + mod(nsec,3600)/60*100 + mod(nsec,60)
return
end
function nsecf2 (nhhmmss,nmmdd,nymd)
C***********************************************************************
C Purpose
C Computes the Total Number of seconds from NYMD using NHHMMSS & NMMDD
C
C Arguments Description
C NHHMMSS IntervaL Frequency (HHMMSS)
C NMMDD Interval Frequency (MMDD)
C NYMD Current Date (YYMMDD)
C
C NOTE:
C IF (NMMDD.ne.0), THEN HOUR FREQUENCY HH MUST BE < 24
C
C***********************************************************************
implicit none
integer nsecf2,nhhmmss,nmmdd,nymd
INTEGER NSDAY, NCYCLE
PARAMETER ( NSDAY = 86400 )
PARAMETER ( NCYCLE = 1461*24*3600 )
INTEGER YEAR, MONTH, DAY
c INTEGER MNDY(12,4)
INTEGER MNDY(12*4)
DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
& 397,34*0 /
integer nsecf,i,nsegm,nsegd,iday,iday2,nday
C***********************************************************************
C* COMPUTE # OF SECONDS FROM NHHMMSS *
C***********************************************************************
nsecf2 = nsecf( nhhmmss )
if( nmmdd.eq.0 ) return
C***********************************************************************
C* COMPUTE # OF DAYS IN A 4-YEAR CYCLE *
C***********************************************************************
DO I=15,48
c MNDY(I,1) = MNDY(I-12,1) + 365
MNDY(I) = MNDY(I-12) + 365
ENDDO
C***********************************************************************
C* COMPUTE # OF SECONDS FROM NMMDD *
C***********************************************************************
nsegm = nmmdd/100
nsegd = mod(nmmdd,100)
YEAR = NYMD / 10000
MONTH = MOD(NYMD,10000) / 100
DAY = MOD(NYMD,100)
c IDAY = MNDY( MONTH ,MOD(YEAR ,4)+1 )
IDAY = MNDY( MONTH +12*MOD(YEAR ,4) )
month = month + nsegm
If( month.gt.12 ) then
month = month - 12
year = year + 1
endif
c IDAY2 = MNDY( MONTH ,MOD(YEAR ,4)+1 )
IDAY2 = MNDY( MONTH +12*MOD(YEAR ,4) )
nday = iday2-iday
if(nday.lt.0) nday = nday + 1461
nday = nday + nsegd
nsecf2 = nsecf2 + nday*nsday
return
end
subroutine FIXDATE (nymd)
implicit none
integer nymd
C Modify 6-digit YYMMDD for dates between 1950-2050
C -------------------------------------------------
if (nymd .lt. 500101) then
nymd = 20000000 + nymd
else if (nymd .le. 991231) then
nymd = 19000000 + nymd
endif
return
end
subroutine INTERP_TIME ( nymd ,nhms ,
& nymd1,nhms1, nymd2,nhms2, fac1,fac2 )
C***********************************************************************
C
C PURPOSE:
C ========
C Compute interpolation factors, fac1 & fac2, to be used in the
C calculation of the instantanious boundary conditions, ie:
C
C q(i,j) = fac1*q1(i,j) + fac2*q2(i,j)
C where:
C q(i,j) => Boundary Data valid at (nymd , nhms )
C q1(i,j) => Boundary Data centered at (nymd1 , nhms1)
C q2(i,j) => Boundary Data centered at (nymd2 , nhms2)
C
C INPUT:
C ======
C nymd : Date (yymmdd) of Current Timestep
C nhms : Time (hhmmss) of Current Timestep
C nymd1 : Date (yymmdd) of Boundary Data 1
C nhms1 : Time (hhmmss) of Boundary Data 1
C nymd2 : Date (yymmdd) of Boundary Data 2
C nhms2 : Time (hhmmss) of Boundary Data 2
C
C OUTPUT:
C =======
C fac1 : Interpolation factor for Boundary Data 1
C fac2 : Interpolation factor for Boundary Data 2
C
C
C***********************************************************************
implicit none
integer nhms,nymd,nhms1,nymd1,nhms2,nymd2
_RL fac1,fac2
INTEGER YEAR , MONTH , DAY , SEC
INTEGER YEAR1, MONTH1, DAY1, SEC1
INTEGER YEAR2, MONTH2, DAY2, SEC2
_RL time00, time1, time2
INTEGER DAYSCY
parameter ( dayscy = 365*4 + 1 )
INTEGER MNDY(12*4)
LOGICAL FIRST
DATA FIRST/.TRUE./
DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
& 397,34*0 /
integer i,nsecf
C***********************************************************************
C* SET TIME BOUNDARIES *
C***********************************************************************
YEAR = NYMD / 10000
MONTH = MOD(NYMD,10000) / 100
DAY = MOD(NYMD,100)
SEC = NSECF(NHMS)
YEAR1 = NYMD1 / 10000
MONTH1 = MOD(NYMD1,10000) / 100
DAY1 = MOD(NYMD1,100)
SEC1 = NSECF(NHMS1)
YEAR2 = NYMD2 / 10000
MONTH2 = MOD(NYMD2,10000) / 100
DAY2 = MOD(NYMD2,100)
SEC2 = NSECF(NHMS2)
C***********************************************************************
C* COMPUTE DAYS IN 4-YEAR CYCLE *
C***********************************************************************
IF(FIRST) THEN
DO I=15,48
MNDY(I) = MNDY(I-12) + 365
ENDDO
FIRST=.FALSE.
ENDIF
C***********************************************************************
C* COMPUTE INTERPOLATION FACTORS *
C***********************************************************************
time00 = DAY + MNDY(MONTH +12*MOD(YEAR ,4)) + float(sec )/86400.
time1 = DAY1 + MNDY(MONTH1+12*MOD(YEAR1,4)) + float(sec1)/86400.
time2 = DAY2 + MNDY(MONTH2+12*MOD(YEAR2,4)) + float(sec2)/86400.
if( time00 .lt.time1 ) time00 = time00 + dayscy
if( time2.lt.time1 ) time2 = time2 + dayscy
fac1 = (time2-time00)/(time2-time1)
fac2 = (time00-time1)/(time2-time1)
RETURN
END
subroutine TICK (nymd,nhms,ndt)
C***********************************************************************
C Purpose
C Tick the Date (nymd) and Time (nhms) by NDT (seconds)
C
C***********************************************************************
implicit none
integer nymd,nhms,ndt
integer nsec,nsecf,incymd,nhmsf
IF(NDT.NE.0) THEN
NSEC = NSECF(NHMS) + NDT
IF (NSEC.GT.86400) THEN
DO WHILE (NSEC.GT.86400)
NSEC = NSEC - 86400
NYMD = INCYMD (NYMD,1)
ENDDO
ENDIF
IF (NSEC.EQ.86400) THEN
NSEC = 0
NYMD = INCYMD (NYMD,1)
ENDIF
IF (NSEC.LT.00000) THEN
DO WHILE (NSEC.LT.0)
NSEC = 86400 + NSEC
NYMD = INCYMD (NYMD,-1)
ENDDO
ENDIF
NHMS = NHMSF (NSEC)
ENDIF
#ifdef FIZHI_USE_FIXED_DAY
NYMD = 20040321
#endif
RETURN
END
subroutine TIC_TIME (mymd,mhms,ndt)
C***********************************************************************
C PURPOSE
C Tick the Clock by NDT (seconds)
C
C***********************************************************************
implicit none
#include "chronos.h"
integer mymd,mhms,ndt
integer nsec,nsecf,incymd,nhmsf
IF(NDT.NE.0) THEN
NSEC = NSECF(NHMS) + NDT
IF (NSEC.GT.86400) THEN
DO WHILE (NSEC.GT.86400)
NSEC = NSEC - 86400
NYMD = INCYMD (NYMD,1)
ENDDO
ENDIF
IF (NSEC.EQ.86400) THEN
NSEC = 0
NYMD = INCYMD (NYMD,1)
ENDIF
IF (NSEC.LT.00000) THEN
DO WHILE (NSEC.LT.0)
NSEC = 86400 + NSEC
NYMD = INCYMD (NYMD,-1)
ENDDO
ENDIF
NHMS = NHMSF (NSEC)
ENDIF
C Pass Back Current Updated Time
C ------------------------------
mymd = nymd
mhms = nhms
RETURN
END
FUNCTION NALARM (MHMS,NYMD,NHMS,NYMD0,NHMS0)
C***********************************************************************
C PURPOSE
C COMPUTES MODULO-FRACTION BETWEEN MHHS AND TOTAL TIME
C USAGE
C ARGUMENTS DESCRIPTION
C MHMS INTERVAL FREQUENCY (HHMMSS)
C NYMD CURRENT YYMMDD
C NHMS CURRENT HHMMSS
C NYMD0 BEGINNING YYMMDD
C NHMS0 BEGINNING HHMMSS
C
C***********************************************************************
implicit none
integer nalarm,MHMS,NYMD,NHMS,NYMD0,NHMS0
integer nsday, ncycle
PARAMETER ( NSDAY = 86400 )
PARAMETER ( NCYCLE = 1461*24*3600 )
INTEGER YEAR, MONTH, DAY, SEC, YEAR0, MONTH0, DAY0, SEC0
integer MNDY(12*4)
DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
& 397,34*0 /
integer i,nsecf,iday,iday0,nsec,nsec0,ntime
C***********************************************************************
C* COMPUTE # OF DAYS IN A 4-YEAR CYCLE *
C***********************************************************************
DO I=15,48
MNDY(I) = MNDY(I-12) + 365
ENDDO
C***********************************************************************
C* SET CURRENT AND BEGINNING TIMES *
C***********************************************************************
YEAR = NYMD / 10000
MONTH = MOD(NYMD,10000) / 100
DAY = MOD(NYMD,100)
SEC = NSECF(NHMS)
YEAR0 = NYMD0 / 10000
MONTH0 = MOD(NYMD0,10000) / 100
DAY0 = MOD(NYMD0,100)
SEC0 = NSECF(NHMS0)
C***********************************************************************
C* COMPUTE POSITIONS IN CYCLE FOR CURRENT AND BEGINNING TIMES *
C***********************************************************************
IDAY = (DAY -1) + MNDY( MONTH +12*MOD(YEAR ,4) )
IDAY0 = (DAY0-1) + MNDY( MONTH0+12*MOD(YEAR0,4) )
NSEC = IDAY *NSDAY + SEC
NSEC0 = IDAY0*NSDAY + SEC0
NTIME = NSEC-NSEC0
IF (NTIME.LT.0 ) NTIME = NTIME + NCYCLE
NALARM = NTIME
IF ( MHMS.NE.0 ) NALARM = MOD( NALARM,NSECF(MHMS) )
RETURN
END
FUNCTION NALARM2(MHMS,NYMD,NHMS,NYMD0,NHMS0)
C***********************************************************************
C PURPOSE
C COMPUTES MODULO-FRACTION BETWEEN MHHS AND TOTAL TIME
C USAGE
C ARGUMENTS DESCRIPTION
C MHMS INTERVAL FREQUENCY (MMDDHHMMSS)
C NYMD CURRENT YYMMDD
C NHMS CURRENT HHMMSS
C NYMD0 BEGINNING YYMMDD
C NHMS0 BEGINNING HHMMSS
C
C***********************************************************************
implicit none
integer nalarm2,MHMS,NYMD,NHMS,NYMD0,NHMS0
integer nsday, ncycle
PARAMETER ( NSDAY = 86400 )
PARAMETER ( NCYCLE = 1461*24*3600 )
INTEGER YEAR, MONTH, DAY, SEC, YEAR0, MONTH0, DAY0, SEC0
integer MNDY(12*4)
DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
& 397,34*0 /
INTEGER NDPM(12)
DATA NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
integer i,nsecf,iday,iday0,nsec,nsec0,ntime
integer NHMMSS,NMMDD
integer iloop
integer testnymd,testnhms,testndpm
C***********************************************************************
C* COMPUTE # OF DAYS IN A 4-YEAR CYCLE *
C***********************************************************************
DO I=15,48
MNDY(I) = MNDY(I-12) + 365
ENDDO
C***********************************************************************
C* SET CURRENT AND BEGINNING TIMES *
C***********************************************************************
YEAR = NYMD / 10000
MONTH = MOD(NYMD,10000) / 100
DAY = MOD(NYMD,100)
SEC = NSECF(NHMS)
YEAR0 = NYMD0 / 10000
MONTH0 = MOD(NYMD0,10000) / 100
DAY0 = MOD(NYMD0,100)
SEC0 = NSECF(NHMS0)
C***********************************************************************
C* COMPUTE POSITIONS IN CYCLE FOR CURRENT AND BEGINNING TIMES *
C***********************************************************************
IDAY = (DAY -1) + MNDY( MONTH +12*MOD(YEAR ,4) )
IDAY0 = (DAY0-1) + MNDY( MONTH0+12*MOD(YEAR0,4) )
NSEC = IDAY *NSDAY + SEC
NSEC0 = IDAY0*NSDAY + SEC0
NTIME = NSEC-NSEC0
IF(NTIME.LT.0) NTIME = NTIME + NCYCLE
NALARM2 = NTIME
IF(MHMS.NE.0)NALARM2 = MOD( NALARM2,NSECF(MHMS) )
IF(MHMS.GE.1000000) THEN
testnymd=nymd0
testnhms=nhms0
NMMDD = MHMS / 1000000
NHMMSS = MOD(MHMS,1000000)
do iloop=1,100000
testnymd=testnymd + nmmdd
testnhms=testnhms + nhmmss
year0=testnymd/10000
month0=mod(testnymd,10000)/100
day0 = mod(testnymd,100)
testndpm = ndpm(month0)
if( month0.eq.2 .and. mod(year0,4).eq.0) testndpm = 29
if(testnhms.ge.240000) then
testnhms = testnhms-240000
testnymd = testnymd + 1
day0 = day0 + 1
endif
if(day0.gt.testndpm) then
testnymd = testnymd - testndpm
testnymd = testnymd + 100
day0 = day0 - testndpm
month0 = month0 + 1
endif
if(month0.gt.12) then
month0 = month0 - 12
year0 = year0 + 1
testnymd = testnymd + 10000 - 1200
endif
sec0 = nsecf(testnhms)
iday0 = (day0-1) + MNDY(month0+12*mod(year0,4) )
nsec0 = iday0 *nsday + sec0
if( (testnymd.gt.nymd) .or.
& (testnymd.eq.testnymd) .and. (testnhms.gt.nhms) )
& go to 200
nalarm2 = nsec-nsec0
enddo
200 continue
ENDIF
RETURN
END
FUNCTION INCYMD (NYMD,M)
C***********************************************************************
C PURPOSE
C INCYMD: NYMD CHANGED BY ONE DAY
C MODYMD: NYMD CONVERTED TO JULIAN DATE
C DESCRIPTION OF INPUT VARIABLES
C NYMD CURRENT DATE IN YYMMDD FORMAT
C M +/- 1 (DAY ADJUSTMENT)
C
C***********************************************************************
implicit none
integer incymd,nymd,m
integer ny,nm,nd,ny00,modymd
INTEGER NDPM(12)
DATA NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
LOGICAL LEAP
DATA NY00 /1900 /
LEAP(NY) = MOD(NY,4).EQ.0 .AND. (NY.NE.0 .OR. MOD(NY00,400).EQ.0)
C***********************************************************************
C
NY = NYMD / 10000
NM = MOD(NYMD,10000) / 100
ND = MOD(NYMD,100) + M
IF (ND.EQ.0) THEN
NM = NM - 1
IF (NM.EQ.0) THEN
NM = 12
NY = NY - 1
ENDIF
ND = NDPM(NM)
IF (NM.EQ.2 .AND. LEAP(NY)) ND = 29
ENDIF
IF (ND.EQ.29 .AND. NM.EQ.2 .AND. LEAP(NY)) GO TO 20
IF (ND.GT.NDPM(NM)) THEN
ND = 1
NM = NM + 1
IF (NM.GT.12) THEN
NM = 1
NY = NY + 1
ENDIF
ENDIF
20 CONTINUE
INCYMD = NY*10000 + NM*100 + ND
RETURN
C***********************************************************************
C E N T R Y M O D Y M D
C***********************************************************************
ENTRY MODYMD (NYMD)
NY = NYMD / 10000
NM = MOD(NYMD,10000) / 100
ND = MOD(NYMD,100)
40 CONTINUE
IF (NM.LE.1) GO TO 60
NM = NM - 1
ND = ND + NDPM(NM)
IF (NM.EQ.2 .AND. LEAP(NY)) ND = ND + 1
GO TO 40
60 CONTINUE
MODYMD = ND
RETURN
END
SUBROUTINE ASTRO ( NYMD,NHMS,ALAT,ALON,IRUN,COSZ,RA )
C***********************************************************************
C
C INPUT:
C ======
C NYMD : CURRENT YYMMDD
C NHMS : CURRENT HHMMSS
C ALAT(IRUN):LATITUDES IN DEGREES.
C ALON(IRUN):LONGITUDES IN DEGREES. (0 = GREENWICH, + = EAST).
C IRUN : # OF POINTS TO CALCULATE
C
C OUTPUT:
C =======
C COSZ(IRUN) : COSINE OF ZENITH ANGLE.
C RA : EARTH-SUN DISTANCE IN UNITS OF
C THE ORBITS SEMI-MAJOR AXIS.
C
C NOTE:
C =====
C THE INSOLATION AT THE TOP OF THE ATMOSPHERE IS:
C
C S(I) = (SOLAR CONSTANT)*(1/RA**2)*COSZ(I),
C
C WHERE:
C RA AND COSZ(I) ARE THE TWO OUTPUTS OF THIS SUBROUTINE.
C
C***********************************************************************
implicit none
C Input Variables
C ---------------
integer nymd, nhms, irun
_RL getcon, cosz(irun), alat(irun), alon(irun), ra
C Local Variables
C ---------------
integer year, day, sec, month, iday, idayp1
integer dayscy
integer i,nsecf,k,km,kp
_RL hc
_RL pi, zero, one, two, six, dg2rd, yrlen, eqnx, ob, ecc, per
_RL daylen, fac, thm, thp, thnow, zs, zc, sj, cj
parameter ( zero = 0.0 )
parameter ( one = 1.0 )
parameter ( two = 2.0 )
parameter ( six = 6.0 )
parameter ( dayscy = 365*4 + 1 )
_RL TH(DAYSCY),T0,T1,T2,T3,T4,FUN,Y
INTEGER MNDY(12*4)
LOGICAL FIRST
DATA FIRST/.TRUE./
SAVE
DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
& 397,34*0 /
FUN(Y,PI,ECC,YRLEN,PER) = (TWO*PI/((ONE-ECC**2)**1.5))*(ONE/YRLEN)
& * (ONE - ECC*COS(Y-PER)) ** 2
C***********************************************************************
C* SET SOME CONSTANTS *
C***********************************************************************
pi = getcon('PI')
dg2rd = getcon('DEG2RAD')
yrlen = getcon('YRLEN')
ob = getcon('OBLDEG') * dg2rd
daylen = getcon('SDAY')
eqnx = getcon('VERNAL EQUINOX')
ecc = getcon('ECCENTRICITY')
per = getcon('PERIHELION') * dg2rd
C***********************************************************************
C* SET CURRENT TIME *
C***********************************************************************
YEAR = NYMD / 10000
MONTH = MOD(NYMD,10000) / 100
DAY = MOD(NYMD,100)
SEC = NSECF(NHMS)
C***********************************************************************
C* COMPUTE DAY-ANGLES FOR 4-YEAR CYCLE *
C***********************************************************************
IF(FIRST) THEN
DO 100 I=15,48
MNDY(I) = MNDY(I-12) + 365
100 CONTINUE
KM = INT(EQNX) + 1
FAC = KM-EQNX
T0 = ZERO
T1 = FUN(T0,PI,ECC,YRLEN,PER )*FAC
T2 = FUN(ZERO+T1/TWO,PI,ECC,YRLEN,PER)*FAC
T3 = FUN(ZERO+T2/TWO,PI,ECC,YRLEN,PER)*FAC
T4 = FUN(ZERO+T3,PI,ECC,YRLEN,PER )*FAC
TH(KM) = (T1 + TWO*(T2 + T3) + T4) / SIX
DO 200 K=2,DAYSCY
T1 = FUN(TH(KM),PI,ECC,YRLEN,PER )
T2 = FUN(TH(KM)+T1/TWO,PI,ECC,YRLEN,PER)
T3 = FUN(TH(KM)+T2/TWO,PI,ECC,YRLEN,PER)
T4 = FUN(TH(KM)+T3,PI,ECC,YRLEN,PER )
KP = MOD(KM,DAYSCY) + 1
TH(KP) = TH(KM) + (T1 + TWO*(T2 + T3) + T4) / SIX
KM = KP
200 CONTINUE
FIRST=.FALSE.
ENDIF
C***********************************************************************
C* COMPUTE EARTH-SUN DISTANCE TO CURRENT SECOND *
C***********************************************************************
IDAY = DAY + MNDY(MONTH+12*MOD(YEAR,4) )
IDAYP1 = MOD( IDAY,DAYSCY) + 1
THM = MOD( TH(IDAY) ,TWO*PI)
THP = MOD( TH(IDAYP1),TWO*PI)
IF(THP.LT.THM) THP = THP + TWO*PI
FAC = FLOAT(SEC)/DAYLEN
THNOW = THM*(ONE-FAC) + THP*FAC
ZS = SIN(THNOW) * SIN(OB)
ZC = SQRT(ONE-ZS*ZS)
RA = (1.-ECC*ECC) / ( ONE-ECC*COS(THNOW-PER) )
C***********************************************************************
C* COMPUTE COSINE OF THE ZENITH ANGLE *
C***********************************************************************
FAC = FAC*TWO*PI + PI
DO I = 1,IRUN
HC = COS( FAC+ALON(I)*DG2RD )
SJ = SIN(ALAT(I)*DG2RD)
CJ = SQRT(ONE-SJ*SJ)
COSZ(I) = SJ*ZS + CJ*ZC*HC
IF( COSZ(I).LT.ZERO ) COSZ(I) = ZERO
ENDDO
RETURN
END
subroutine TIME_BOUND(nymd,nhms,nymd1,nhms1,nymd2,nhms2,imnm,imnp)
C***********************************************************************
C PURPOSE
C Compute Date and Time boundaries.
C
C ARGUMENTS DESCRIPTION
C nymd .... Current Date
C nhms .... Current Time
C nymd1 ... Previous Date Boundary
C nhms1 ... Previous Time Boundary
C nymd2 ... Subsequent Date Boundary
C nhms2 ... Subsequent Time Boundary
C
C imnm .... Previous Time Index for Interpolation
C imnp .... Subsequent Time Index for Interpolation
C
C***********************************************************************
implicit none
integer nymd,nhms, nymd1,nhms1, nymd2,nhms2
C Local Variables
C ---------------
integer month,day,nyear,midmon1,midmon,midmon2
integer imnm,imnp
INTEGER DAYS(14), daysm, days0, daysp
DATA DAYS /31,31,28,31,30,31,30,31,31,30,31,30,31,31/
integer nmonf,ndayf,n
NMONF(N) = MOD(N,10000)/100
NDAYF(N) = MOD(N,100)
C*********************************************************************
C**** Find Proper Month and Time Boundaries for Climatological Data **
C*********************************************************************
MONTH = NMONF(NYMD)
DAY = NDAYF(NYMD)
daysm = days(month )
days0 = days(month+1)
daysp = days(month+2)
C Check for Leap Year
C -------------------
nyear = nymd/10000
if( 4*(nyear/4).eq.nyear ) then
if( month.eq.3 ) daysm = daysm+1
if( month.eq.2 ) days0 = days0+1
if( month.eq.1 ) daysp = daysp+1
endif
MIDMON1 = daysm/2 + 1
MIDMON = days0/2 + 1
MIDMON2 = daysp/2 + 1
IF(DAY.LT.MIDMON) THEN
imnm = month
imnp = month + 1
nymd2 = (nymd/10000)*10000 + month*100 + midmon
nhms2 = 000000
nymd1 = nymd2
nhms1 = nhms2
call TICK ( nymd1,nhms1, -midmon *86400 )
call TICK ( nymd1,nhms1,-(daysm-midmon1)*86400 )
ELSE
IMNM = MONTH + 1
IMNP = MONTH + 2
nymd1 = (nymd/10000)*10000 + month*100 + midmon
nhms1 = 000000
nymd2 = nymd1
nhms2 = nhms1
call TICK ( nymd2,nhms2,(days0-midmon)*86400 )
call TICK ( nymd2,nhms2, midmon2*86400 )
ENDIF
C -------------------------------------------------------------
C Note: At this point, imnm & imnp range between 01-14, where
C 01 -> Previous years December
C 02-13 -> Current years January-December
C 14 -> Next years January
C -------------------------------------------------------------
imnm = imnm-1
imnp = imnp-1
if( imnm.eq.0 ) imnm = 12
if( imnp.eq.0 ) imnp = 12
if( imnm.eq.13 ) imnm = 1
if( imnp.eq.13 ) imnp = 1
return
end
subroutine TIME2FREQ2(MMDD,NYMD,NHMS,timeleft)
C***********************************************************************
C PURPOSE
C COMPUTES TIME IN SECONDS UNTIL WE REACH THE NEXT MMDD
C (ASSUME that the target time is 0Z)
C
C ARGUMENTS DESCRIPTION
C MMDD FREQUENCY (MMDDHHMMSS)
C NYMD CURRENT YYMMDD
C NHMS CURRENT HHMMSS
C TIMELEFT TIME LEFT (SECONDS)
C
C NOTES - Only used when the frequency is in units of months
C Assumes that we always want to be at a month boundary
C***********************************************************************
implicit none
integer mmdd,nymd,nhms,timeleft,daysleft
integer nsday
PARAMETER ( NSDAY = 86400 )
integer year, month, day, sec
integer yearnext, monthnext, daynext
integer i,nsecf,iday,idaynext,nsec
integer testnymd
integer MNDY(12*4)
DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
& 397,34*0 /
C***********************************************************************
C* COMPUTE # OF DAYS IN A 4-YEAR CYCLE *
C***********************************************************************
DO I=15,48
MNDY(I) = MNDY(I-12) + 365
ENDDO
C***********************************************************************
C* SET CURRENT TIME ELEMENTS *
C***********************************************************************
YEAR = NYMD / 10000
MONTH = MOD(NYMD,10000) / 100
DAY = MOD(NYMD,100)
SEC = NSECF(NHMS)
C***********************************************************************
C* COMPUTE POSITIONS IN CYCLE FOR CURRENT AND BEGINNING TIMES *
C***********************************************************************
IDAY = (DAY -1) + MNDY( MONTH +12*MOD(YEAR ,4) )
NSEC = IDAY *NSDAY + SEC
testnymd=nymd + mmdd
yearnext=testnymd/10000
monthnext=mod(testnymd,10000)/100
daynext = 1
if(monthnext.gt.12) then
monthnext = monthnext - 12
yearnext = yearnext + 1
endif
testnymd = yearnext*10000 + monthnext*100 + daynext
idaynext = MNDY(monthnext+12*mod(yearnext,4) )
daysleft = idaynext - iday
if(daysleft.lt.0) daysleft = daysleft + 1461
timeleft = daysleft * nsday - sec
RETURN
END