/[MITgcm]/MITgcm/verification/fizhi-cs-aqualev20/code/fizhi_clockstuff.F
ViewVC logotype

Annotation of /MITgcm/verification/fizhi-cs-aqualev20/code/fizhi_clockstuff.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (hide annotations) (download)
Tue Mar 13 15:38:59 2007 UTC (17 years, 3 months ago) by molod
Branch: MAIN
CVS Tags: checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint59, checkpoint58x_post, checkpoint58y_post, checkpoint58w_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint60, checkpoint61
Changes since 1.1: +29 -26 lines
Bring this experiment code up to date with background fizhi code

1 molod 1.2 C $Header: /u/gcmpack/MITgcm/verification/fizhi-cs-aqualev20/code/fizhi_clockstuff.F,v 1.1 2006/04/03 20:55:14 molod Exp $
2 molod 1.1 C $Name: $
3    
4     #include "FIZHI_OPTIONS.h"
5     subroutine set_alarm (tag,datein,timein,freq)
6     C***********************************************************************
7     C Purpose
8     C -------
9     C Utility to Set Internal Alarms
10     C
11     C Argument Description
12     C --------------------
13     C tag ....... Character String Tagging Alarm Process
14     C date ...... Begining Date for Alarm
15     C time ...... Begining Time for Alarm
16     C freq ...... Repeating Frequency Interval for Alarm
17     C
18     C***********************************************************************
19    
20     implicit none
21     character*(*) tag
22     integer freq,datein,timein
23    
24     #ifdef ALLOW_USE_MPI
25     #include "SIZE.h"
26     #include "EEPARAMS.h"
27     #include "EESUPPORT.h"
28     #endif
29    
30     #include "chronos.h"
31    
32     #ifdef ALLOW_USE_MPI
33     c MPI Utilities
34     c -------------
35     c#include "mpif.h"
36     c integer mpi_comm_model,ierror
37     integer ierror
38     #endif
39    
40     integer myid
41     logical first,set
42     data first /.true./
43    
44     integer n
45     #ifdef ALLOW_USE_MPI
46     call mpi_comm_rank ( mpi_comm_model,myid,ierror )
47     #else
48     myid = 1
49     #endif
50    
51     if(first) then
52     ntags = 1
53     tags(1) = tag
54     freqs(1) = freq
55     dates(1) = datein
56     times(1) = timein
57     if( myid.eq.1 ) write(6,100) datein,timein,freq,tags(1)
58     else
59    
60     set = .false.
61     do n=1,ntags
62     if(tag.eq.tags(n)) then
63     if( myid.eq.1 ) then
64     print *, 'Warning! Alarm has already been set for Tag: ',tag
65     print *, 'Changing Alarm Information:'
66     print *, 'Frequency: ',freqs(n),' (Old) ',freq,' (New)'
67     print *, ' Date0: ',dates(n),' (Old) ',datein,' (New)'
68     print *, ' Time0: ',times(n),' (Old) ',timein,' (New)'
69     endif
70     freqs(n) = freq
71     dates(n) = datein
72     times(n) = timein
73     set = .true.
74     endif
75     enddo
76     if(.not.set) then
77     ntags = ntags+1
78     if(ntags.gt.maxtag ) then
79     if( myid.eq.1 ) then
80     print *, 'Too many Alarms are Set!!'
81     print *, 'Maximum Number of Alarms = ',maxtag
82     endif
83     call my_finalize
84     call my_exit (101)
85     endif
86     tags(ntags) = tag
87     freqs(ntags) = freq
88     dates(ntags) = datein
89     times(ntags) = timein
90     if( myid.eq.1 ) write(6,100) datein,timein,freq,tags(ntags)
91     endif
92     endif
93    
94     first = .false.
95     100 format(1x,'Setting Alarm for: ',i8,2x,i6.6,', with frequency: ',
96     . i10,', and Tag: ',a80)
97     return
98     end
99    
100     subroutine get_alarm (tag,datein,timein,freq,tleft)
101     C***********************************************************************
102     C Purpose
103     C -------
104     C Utility to Get Internal Alarm Information
105     C
106     C Input
107     C -----
108     C tag ....... Character String Tagging Alarm Process
109     C
110     C Output
111     C ------
112     C datein ...... Begining Date for Alarm
113     C timein ...... Begining Time for Alarm
114     C freq ........ Frequency Interval for Alarm
115     C tleft ....... Time Remaining (seconds) before Alarm is TRUE
116     C
117     C***********************************************************************
118    
119     implicit none
120     character*(*) tag
121     integer freq,datein,timein,tleft
122    
123     #ifdef ALLOW_USE_MPI
124     #include "SIZE.h"
125     #include "EEPARAMS.h"
126     #include "EESUPPORT.h"
127     #endif
128    
129     #include "chronos.h"
130    
131     #ifdef ALLOW_USE_MPI
132     c MPI Utilities
133     c -------------
134     c#include "mpif.h"
135     c integer mpi_comm_model,ierror
136     integer ierror
137     #endif
138    
139     logical set,alarm
140     external alarm
141     integer myid,n,nalarm,nsecf
142    
143     #ifdef ALLOW_USE_MPI
144     call mpi_comm_rank ( mpi_comm_model,myid,ierror )
145     #else
146     myid = 1
147     #endif
148    
149     set = .false.
150     do n=1,ntags
151     if(tag.eq.tags(n)) then
152     freq = freqs(n)
153     datein = dates(n)
154     timein = times(n)
155    
156     if( alarm(tag) ) then
157     tleft = 0
158     else
159     call get_time (nymd,nhms)
160     tleft = nsecf(freq) - nalarm(freq,nymd,nhms,datein,timein )
161     endif
162    
163     set = .true.
164     endif
165     enddo
166    
167     if(.not.set) then
168     if( myid.eq.1 ) print *, 'Alarm has not been set for Tag: ',tag
169     freq = 0
170     datein = 0
171     timein = 0
172     tleft = 0
173     endif
174    
175     return
176     end
177    
178     function alarm (tag)
179     implicit none
180     character*(*) tag
181     integer datein,timein
182     logical alarm
183     #include "chronos.h"
184    
185     integer n,modalarm,nalarm,freq,date0,time0
186     modalarm(freq,date0,time0)=nalarm(freq,datein,timein,date0,time0)
187    
188     call get_time (datein,timein)
189    
190     alarm = .false.
191     do n=1,ntags
192     if( tags(n).eq.tag ) then
193     if( freqs(n).eq.0 ) then
194     alarm = (dates(n).eq.datein) .and. (times(n).eq.timein)
195     else
196     alarm = ( datein.gt.dates(n) .or.
197     . (datein.eq.dates(n) .and. timein.ge.times(n)) ) .and.
198     . modalarm( freqs(n),dates(n),times(n) ).eq.0
199     endif
200     endif
201     enddo
202    
203     return
204     end
205    
206     function alarm2 (tag)
207     implicit none
208     character*(*) tag
209     integer datein,timein
210     logical alarm2
211     #include "chronos.h"
212    
213     integer n,modalarm,nalarm,nalarm2,freq,date0,time0
214     modalarm(freq,date0,time0)=nalarm2(freq,datein,timein,date0,time0)
215    
216     call get_time (datein,timein)
217    
218     alarm2 = .false.
219     do n=1,ntags
220     if( tags(n).eq.tag ) then
221     if( freqs(n).eq.0 ) then
222     alarm2 = (dates(n).eq.datein) .and. (times(n).eq.timein)
223     else
224     alarm2 = ( datein.gt.dates(n) .or.
225     . (datein.eq.dates(n) .and. timein.ge.times(n)) ) .and.
226     . modalarm( freqs(n),dates(n),times(n) ).eq.0
227     endif
228     endif
229     enddo
230    
231     return
232     end
233    
234     function alarm2next (tag,deltat)
235     implicit none
236     character*(*) tag
237     _RL deltat
238     integer datein,timein,ndt
239     integer dateminus,timeminus
240     logical alarm2next
241     #include "chronos.h"
242    
243     integer n,modalarm,nalarm,nalarm2,freq,date0,time0
244     modalarm(freq,date0,time0)=nalarm2(freq,datein,timein,date0,time0)
245    
246     ndt = int(deltat)
247     call get_time (dateminus,timeminus)
248     datein = dateminus
249     timein = timeminus
250     call tick(datein,timein,ndt)
251    
252     alarm2next = .false.
253     do n=1,ntags
254     if( tags(n).eq.tag ) then
255     if( freqs(n).eq.0 ) then
256     alarm2next = (dates(n).eq.datein) .and. (times(n).eq.timein)
257     else
258     alarm2next = ( datein.gt.dates(n) .or.
259     . (datein.eq.dates(n) .and. timein.ge.times(n)) ) .and.
260     . modalarm( freqs(n),dates(n),times(n) ).eq.0
261     endif
262     endif
263     enddo
264    
265     return
266     end
267    
268     subroutine set_time (datein,timein)
269     implicit none
270     integer datein,timein
271    
272     #ifdef ALLOW_USE_MPI
273     #include "SIZE.h"
274     #include "EEPARAMS.h"
275     #include "EESUPPORT.h"
276     #endif
277    
278     #include "chronos.h"
279    
280     #ifdef ALLOW_USE_MPI
281     c MPI Utilities
282     c -------------
283     c#include "mpif.h"
284     c integer mpi_comm_model,ierror
285     integer ierror
286     #endif
287     integer myid
288    
289     #ifdef ALLOW_USE_MPI
290     call mpi_comm_rank ( mpi_comm_model,myid,ierror )
291     #else
292     myid = 1
293     #endif
294     if( myid.eq.1 ) then
295     print *, 'Setting Clock'
296     print *, 'Date: ',datein
297     print *, 'Time: ',timein
298     endif
299    
300     nymd = datein
301     nhms = timein
302     return
303     end
304    
305     subroutine get_time (datein,timein)
306     implicit none
307     integer datein,timein
308    
309     #include "chronos.h"
310    
311     datein = nymd
312     timein = nhms
313     return
314     end
315    
316     function nsecf (nhms)
317     C***********************************************************************
318     C Purpose
319     C Converts NHMS format to Total Seconds
320     C
321     C***********************************************************************
322     implicit none
323     integer nhms, nsecf
324     nsecf = nhms/10000*3600 + mod(nhms,10000)/100*60 + mod(nhms,100)
325     return
326     end
327    
328     function nhmsf (nsec)
329     C***********************************************************************
330     C Purpose
331     C Converts Total Seconds to NHMS format
332     C
333     C***********************************************************************
334     implicit none
335     integer nhmsf, nsec
336     nhmsf = nsec/3600*10000 + mod(nsec,3600)/60*100 + mod(nsec,60)
337     return
338     end
339    
340     function nsecf2 (nhhmmss,nmmdd,nymd)
341     C***********************************************************************
342     C Purpose
343     C Computes the Total Number of seconds from NYMD using NHHMMSS & NMMDD
344     C
345     C Arguments Description
346     C NHHMMSS IntervaL Frequency (HHMMSS)
347     C NMMDD Interval Frequency (MMDD)
348     C NYMD Current Date (YYMMDD)
349     C
350     C NOTE:
351     C IF (NMMDD.ne.0), THEN HOUR FREQUENCY HH MUST BE < 24
352     C
353     C***********************************************************************
354     implicit none
355    
356     integer nsecf2,nhhmmss,nmmdd,nymd
357    
358     INTEGER NSDAY, NCYCLE
359     PARAMETER ( NSDAY = 86400 )
360     PARAMETER ( NCYCLE = 1461*24*3600 )
361    
362     INTEGER YEAR, MONTH, DAY
363    
364     INTEGER MNDY(12,4)
365     DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
366     . 397,34*0 /
367    
368     integer nsecf,i,nsegm,nsegd,iday,iday2,nday
369    
370     C***********************************************************************
371     C* COMPUTE # OF SECONDS FROM NHHMMSS *
372     C***********************************************************************
373    
374     nsecf2 = nsecf( nhhmmss )
375    
376     if( nmmdd.eq.0 ) return
377    
378     C***********************************************************************
379     C* COMPUTE # OF DAYS IN A 4-YEAR CYCLE *
380     C***********************************************************************
381    
382     DO I=15,48
383     MNDY(I,1) = MNDY(I-12,1) + 365
384     ENDDO
385    
386     C***********************************************************************
387     C* COMPUTE # OF SECONDS FROM NMMDD *
388     C***********************************************************************
389    
390     nsegm = nmmdd/100
391     nsegd = mod(nmmdd,100)
392    
393     YEAR = NYMD / 10000
394     MONTH = MOD(NYMD,10000) / 100
395     DAY = MOD(NYMD,100)
396    
397     IDAY = MNDY( MONTH ,MOD(YEAR ,4)+1 )
398     month = month + nsegm
399     If( month.gt.12 ) then
400     month = month - 12
401     year = year + 1
402     endif
403     IDAY2 = MNDY( MONTH ,MOD(YEAR ,4)+1 )
404    
405     nday = iday2-iday
406     if(nday.lt.0) nday = nday + 1461
407     nday = nday + nsegd
408    
409     nsecf2 = nsecf2 + nday*nsday
410    
411     return
412     end
413    
414     subroutine fixdate (nymd)
415     implicit none
416     integer nymd
417    
418     c Modify 6-digit YYMMDD for dates between 1950-2050
419     c -------------------------------------------------
420     if (nymd .lt. 500101) then
421     nymd = 20000000 + nymd
422     else if (nymd .le. 991231) then
423     nymd = 19000000 + nymd
424     endif
425    
426     return
427     end
428    
429     subroutine interp_time ( nymd ,nhms ,
430     . nymd1,nhms1, nymd2,nhms2, fac1,fac2 )
431     C***********************************************************************
432     C
433     C PURPOSE:
434     C ========
435     C Compute interpolation factors, fac1 & fac2, to be used in the
436     C calculation of the instantanious boundary conditions, ie:
437     C
438     C q(i,j) = fac1*q1(i,j) + fac2*q2(i,j)
439     C where:
440     C q(i,j) => Boundary Data valid at (nymd , nhms )
441     C q1(i,j) => Boundary Data centered at (nymd1 , nhms1)
442     C q2(i,j) => Boundary Data centered at (nymd2 , nhms2)
443     C
444     C INPUT:
445     C ======
446     C nymd : Date (yymmdd) of Current Timestep
447     C nhms : Time (hhmmss) of Current Timestep
448     C nymd1 : Date (yymmdd) of Boundary Data 1
449     C nhms1 : Time (hhmmss) of Boundary Data 1
450     C nymd2 : Date (yymmdd) of Boundary Data 2
451     C nhms2 : Time (hhmmss) of Boundary Data 2
452     C
453     C OUTPUT:
454     C =======
455     C fac1 : Interpolation factor for Boundary Data 1
456     C fac2 : Interpolation factor for Boundary Data 2
457     C
458     C
459     C***********************************************************************
460     implicit none
461    
462     integer nhms,nymd,nhms1,nymd1,nhms2,nymd2
463 molod 1.2 _RL getcon, fac1,fac2
464 molod 1.1
465     INTEGER YEAR , MONTH , DAY , SEC
466     INTEGER YEAR1, MONTH1, DAY1, SEC1
467     INTEGER YEAR2, MONTH2, DAY2, SEC2
468    
469     _RL time00, time1, time2
470    
471     INTEGER DAYSCY
472 molod 1.2 parameter ( dayscy = 365*4 + 1 )
473 molod 1.1
474     INTEGER MNDY(12,4)
475    
476     LOGICAL FIRST
477     DATA FIRST/.TRUE./
478    
479     DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
480     . 397,34*0 /
481    
482     integer i,nsecf
483    
484     C***********************************************************************
485     C* SET TIME BOUNDARIES *
486     C***********************************************************************
487    
488     YEAR = NYMD / 10000
489     MONTH = MOD(NYMD,10000) / 100
490     DAY = MOD(NYMD,100)
491     SEC = NSECF(NHMS)
492    
493     YEAR1 = NYMD1 / 10000
494     MONTH1 = MOD(NYMD1,10000) / 100
495     DAY1 = MOD(NYMD1,100)
496     SEC1 = NSECF(NHMS1)
497    
498     YEAR2 = NYMD2 / 10000
499     MONTH2 = MOD(NYMD2,10000) / 100
500     DAY2 = MOD(NYMD2,100)
501     SEC2 = NSECF(NHMS2)
502    
503     C***********************************************************************
504     C* COMPUTE DAYS IN 4-YEAR CYCLE *
505     C***********************************************************************
506    
507     IF(FIRST) THEN
508     DO I=15,48
509     MNDY(I,1) = MNDY(I-12,1) + 365
510     ENDDO
511     FIRST=.FALSE.
512     ENDIF
513    
514     C***********************************************************************
515     C* COMPUTE INTERPOLATION FACTORS *
516     C***********************************************************************
517    
518     time00 = DAY + MNDY(MONTH ,MOD(YEAR ,4)+1) + float(sec )/86400.
519     time1 = DAY1 + MNDY(MONTH1,MOD(YEAR1,4)+1) + float(sec1)/86400.
520     time2 = DAY2 + MNDY(MONTH2,MOD(YEAR2,4)+1) + float(sec2)/86400.
521    
522     if( time00 .lt.time1 ) time00 = time00 + dayscy
523     if( time2.lt.time1 ) time2 = time2 + dayscy
524    
525     fac1 = (time2-time00)/(time2-time1)
526     fac2 = (time00-time1)/(time2-time1)
527    
528     RETURN
529     END
530    
531     subroutine tick (nymd,nhms,ndt)
532     C***********************************************************************
533     C Purpose
534     C Tick the Date (nymd) and Time (nhms) by NDT (seconds)
535     C
536     C***********************************************************************
537     implicit none
538    
539     integer nymd,nhms,ndt
540    
541     integer nsec,nsecf,incymd,nhmsf
542    
543     IF(NDT.NE.0) THEN
544     NSEC = NSECF(NHMS) + NDT
545    
546     IF (NSEC.GT.86400) THEN
547     DO WHILE (NSEC.GT.86400)
548     NSEC = NSEC - 86400
549     NYMD = INCYMD (NYMD,1)
550     ENDDO
551     ENDIF
552    
553     IF (NSEC.EQ.86400) THEN
554     NSEC = 0
555     NYMD = INCYMD (NYMD,1)
556     ENDIF
557    
558     IF (NSEC.LT.00000) THEN
559     DO WHILE (NSEC.LT.0)
560     NSEC = 86400 + NSEC
561     NYMD = INCYMD (NYMD,-1)
562     ENDDO
563     ENDIF
564    
565     NHMS = NHMSF (NSEC)
566     ENDIF
567    
568     NYMD = 20040321
569    
570     RETURN
571     END
572    
573     subroutine tic_time (mymd,mhms,ndt)
574     C***********************************************************************
575     C PURPOSE
576     C Tick the Clock by NDT (seconds)
577     C
578     C***********************************************************************
579     implicit none
580     #include "chronos.h"
581    
582     integer mymd,mhms,ndt
583    
584     integer nsec,nsecf,incymd,nhmsf
585    
586     IF(NDT.NE.0) THEN
587     NSEC = NSECF(NHMS) + NDT
588    
589     IF (NSEC.GT.86400) THEN
590     DO WHILE (NSEC.GT.86400)
591     NSEC = NSEC - 86400
592     NYMD = INCYMD (NYMD,1)
593     ENDDO
594     ENDIF
595    
596     IF (NSEC.EQ.86400) THEN
597     NSEC = 0
598     NYMD = INCYMD (NYMD,1)
599     ENDIF
600    
601     IF (NSEC.LT.00000) THEN
602     DO WHILE (NSEC.LT.0)
603     NSEC = 86400 + NSEC
604     NYMD = INCYMD (NYMD,-1)
605     ENDDO
606     ENDIF
607    
608     NHMS = NHMSF (NSEC)
609     ENDIF
610    
611     c Pass Back Current Updated Time
612     c ------------------------------
613     mymd = nymd
614     mhms = nhms
615    
616     RETURN
617     END
618    
619     FUNCTION NALARM (MHMS,NYMD,NHMS,NYMD0,NHMS0)
620     C***********************************************************************
621     C PURPOSE
622     C COMPUTES MODULO-FRACTION BETWEEN MHHS AND TOTAL TIME
623     C USAGE
624     C ARGUMENTS DESCRIPTION
625     C MHMS INTERVAL FREQUENCY (HHMMSS)
626     C NYMD CURRENT YYMMDD
627     C NHMS CURRENT HHMMSS
628     C NYMD0 BEGINNING YYMMDD
629     C NHMS0 BEGINNING HHMMSS
630     C
631     C***********************************************************************
632     implicit none
633    
634     integer nalarm,MHMS,NYMD,NHMS,NYMD0,NHMS0
635    
636     integer nsday, ncycle
637     PARAMETER ( NSDAY = 86400 )
638     PARAMETER ( NCYCLE = 1461*24*3600 )
639    
640     INTEGER YEAR, MONTH, DAY, SEC, YEAR0, MONTH0, DAY0, SEC0
641    
642     integer MNDY(12,4)
643     DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
644     . 397,34*0 /
645    
646     integer i,nsecf,iday,iday0,nsec,nsec0,ntime
647    
648     C***********************************************************************
649     C* COMPUTE # OF DAYS IN A 4-YEAR CYCLE *
650     C***********************************************************************
651    
652     DO I=15,48
653     MNDY(I,1) = MNDY(I-12,1) + 365
654     ENDDO
655    
656     C***********************************************************************
657     C* SET CURRENT AND BEGINNING TIMES *
658     C***********************************************************************
659    
660     YEAR = NYMD / 10000
661     MONTH = MOD(NYMD,10000) / 100
662     DAY = MOD(NYMD,100)
663     SEC = NSECF(NHMS)
664    
665     YEAR0 = NYMD0 / 10000
666     MONTH0 = MOD(NYMD0,10000) / 100
667     DAY0 = MOD(NYMD0,100)
668     SEC0 = NSECF(NHMS0)
669    
670     C***********************************************************************
671     C* COMPUTE POSITIONS IN CYCLE FOR CURRENT AND BEGINNING TIMES *
672     C***********************************************************************
673    
674     IDAY = (DAY -1) + MNDY( MONTH ,MOD(YEAR ,4)+1 )
675     IDAY0 = (DAY0-1) + MNDY( MONTH0,MOD(YEAR0,4)+1 )
676    
677     NSEC = IDAY *NSDAY + SEC
678     NSEC0 = IDAY0*NSDAY + SEC0
679    
680     NTIME = NSEC-NSEC0
681     IF (NTIME.LT.0 ) NTIME = NTIME + NCYCLE
682     NALARM = NTIME
683     IF ( MHMS.NE.0 ) NALARM = MOD( NALARM,NSECF(MHMS) )
684    
685     RETURN
686     END
687    
688     FUNCTION NALARM2(MHMS,NYMD,NHMS,NYMD0,NHMS0)
689     C***********************************************************************
690     C PURPOSE
691     C COMPUTES MODULO-FRACTION BETWEEN MHHS AND TOTAL TIME
692     C USAGE
693     C ARGUMENTS DESCRIPTION
694     C MHMS INTERVAL FREQUENCY (MMDDHHMMSS)
695     C NYMD CURRENT YYMMDD
696     C NHMS CURRENT HHMMSS
697     C NYMD0 BEGINNING YYMMDD
698     C NHMS0 BEGINNING HHMMSS
699     C
700     C***********************************************************************
701     implicit none
702    
703     integer nalarm2,MHMS,NYMD,NHMS,NYMD0,NHMS0
704    
705     integer nsday, ncycle
706     PARAMETER ( NSDAY = 86400 )
707     PARAMETER ( NCYCLE = 1461*24*3600 )
708    
709     INTEGER YEAR, MONTH, DAY, SEC, YEAR0, MONTH0, DAY0, SEC0
710    
711     integer MNDY(12,4)
712     DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
713     . 397,34*0 /
714     INTEGER NDPM(12)
715     DATA NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
716    
717     integer i,nsecf2,nsecf,iday,iday0,nsec,nsec0,ntime
718     integer NHMMSS,NMMDD
719     integer iloop
720     integer testnymd,testnhms,testndpm
721    
722     C***********************************************************************
723     C* COMPUTE # OF DAYS IN A 4-YEAR CYCLE *
724     C***********************************************************************
725    
726     DO I=15,48
727     MNDY(I,1) = MNDY(I-12,1) + 365
728     ENDDO
729    
730     C***********************************************************************
731     C* SET CURRENT AND BEGINNING TIMES *
732     C***********************************************************************
733    
734     YEAR = NYMD / 10000
735     MONTH = MOD(NYMD,10000) / 100
736     DAY = MOD(NYMD,100)
737     SEC = NSECF(NHMS)
738    
739     YEAR0 = NYMD0 / 10000
740     MONTH0 = MOD(NYMD0,10000) / 100
741     DAY0 = MOD(NYMD0,100)
742     SEC0 = NSECF(NHMS0)
743    
744     C***********************************************************************
745     C* COMPUTE POSITIONS IN CYCLE FOR CURRENT AND BEGINNING TIMES *
746     C***********************************************************************
747    
748     IDAY = (DAY -1) + MNDY( MONTH ,MOD(YEAR ,4)+1 )
749     IDAY0 = (DAY0-1) + MNDY( MONTH0,MOD(YEAR0,4)+1 )
750    
751     NSEC = IDAY *NSDAY + SEC
752     NSEC0 = IDAY0*NSDAY + SEC0
753    
754     NTIME = NSEC-NSEC0
755     IF(NTIME.LT.0) NTIME = NTIME + NCYCLE
756     NALARM2 = NTIME
757     IF(MHMS.NE.0)NALARM2 = MOD( NALARM2,NSECF(MHMS) )
758     IF(MHMS.GE.1000000) THEN
759     testnymd=nymd0
760     testnhms=nhms0
761     NMMDD = MHMS / 1000000
762     NHMMSS = MOD(MHMS,1000000)
763     do iloop=1,100000
764     testnymd=testnymd + nmmdd
765     testnhms=testnhms + nhmmss
766     year0=testnymd/10000
767     month0=mod(testnymd,10000)/100
768     day0 = mod(testnymd,100)
769     testndpm = ndpm(month0)
770     if( month0.eq.2 .and. mod(year0,4).eq.0) testndpm = 29
771     if(testnhms.ge.240000) then
772     testnhms = testnhms-240000
773     testnymd = testnymd + 1
774     day0 = day0 + 1
775     endif
776     if(day0.gt.testndpm) then
777     testnymd = testnymd - testndpm
778     testnymd = testnymd + 100
779     day0 = day0 - testndpm
780     month0 = month0 + 1
781     endif
782     if(month0.gt.12) then
783     month0 = month0 - 12
784     year0 = year0 + 1
785     testnymd = testnymd + 10000 - 1200
786     endif
787     sec0 = nsecf(testnhms)
788     iday0 = (day0-1) + mndy(month0,mod(year0,4)+1)
789     nsec0 = iday0 *nsday + sec0
790     if( (testnymd.gt.nymd) .or.
791     . (testnymd.eq.testnymd) .and. (testnhms.gt.nhms) )
792     . go to 200
793     nalarm2 = nsec-nsec0
794     enddo
795     200 continue
796     ENDIF
797    
798     RETURN
799     END
800    
801     FUNCTION INCYMD (NYMD,M)
802     C***********************************************************************
803     C PURPOSE
804     C INCYMD: NYMD CHANGED BY ONE DAY
805     C MODYMD: NYMD CONVERTED TO JULIAN DATE
806 molod 1.2 C DESCRIPTION OF INPUT VARIABLES
807 molod 1.1 C NYMD CURRENT DATE IN YYMMDD FORMAT
808     C M +/- 1 (DAY ADJUSTMENT)
809     C
810     C***********************************************************************
811     implicit none
812     integer incymd,nymd,m
813    
814     integer ny,nm,nd,ny00,modymd
815    
816     INTEGER NDPM(12)
817     DATA NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
818     LOGICAL LEAP
819     DATA NY00 /1900 /
820     LEAP(NY) = MOD(NY,4).EQ.0 .AND. (NY.NE.0 .OR. MOD(NY00,400).EQ.0)
821    
822     C***********************************************************************
823     C
824     NY = NYMD / 10000
825     NM = MOD(NYMD,10000) / 100
826     ND = MOD(NYMD,100) + M
827    
828     IF (ND.EQ.0) THEN
829     NM = NM - 1
830     IF (NM.EQ.0) THEN
831     NM = 12
832     NY = NY - 1
833     ENDIF
834     ND = NDPM(NM)
835     IF (NM.EQ.2 .AND. LEAP(NY)) ND = 29
836     ENDIF
837    
838     IF (ND.EQ.29 .AND. NM.EQ.2 .AND. LEAP(NY)) GO TO 20
839    
840     IF (ND.GT.NDPM(NM)) THEN
841     ND = 1
842     NM = NM + 1
843     IF (NM.GT.12) THEN
844     NM = 1
845     NY = NY + 1
846     ENDIF
847     ENDIF
848    
849     20 CONTINUE
850     INCYMD = NY*10000 + NM*100 + ND
851    
852     RETURN
853    
854     C***********************************************************************
855     C E N T R Y M O D Y M D
856     C***********************************************************************
857    
858     ENTRY MODYMD (NYMD)
859    
860     NY = NYMD / 10000
861     NM = MOD(NYMD,10000) / 100
862     ND = MOD(NYMD,100)
863    
864     40 CONTINUE
865     IF (NM.LE.1) GO TO 60
866     NM = NM - 1
867     ND = ND + NDPM(NM)
868     IF (NM.EQ.2 .AND. LEAP(NY)) ND = ND + 1
869     GO TO 40
870    
871     60 CONTINUE
872     MODYMD = ND
873    
874     RETURN
875     END
876    
877     SUBROUTINE ASTRO ( NYMD,NHMS,ALAT,ALON,IRUN,COSZ,RA )
878     C***********************************************************************
879     C
880     C INPUT:
881     C ======
882     C NYMD : CURRENT YYMMDD
883     C NHMS : CURRENT HHMMSS
884     C ALAT(IRUN):LATITUDES IN DEGREES.
885     C ALON(IRUN):LONGITUDES IN DEGREES. (0 = GREENWICH, + = EAST).
886     C IRUN : # OF POINTS TO CALCULATE
887     C
888     C OUTPUT:
889     C =======
890     C COSZ(IRUN) : COSINE OF ZENITH ANGLE.
891     C RA : EARTH-SUN DISTANCE IN UNITS OF
892     C THE ORBITS SEMI-MAJOR AXIS.
893     C
894     C NOTE:
895     C =====
896     C THE INSOLATION AT THE TOP OF THE ATMOSPHERE IS:
897     C
898     C S(I) = (SOLAR CONSTANT)*(1/RA**2)*COSZ(I),
899     C
900     C WHERE:
901     C RA AND COSZ(I) ARE THE TWO OUTPUTS OF THIS SUBROUTINE.
902     C
903     C***********************************************************************
904    
905     implicit none
906    
907     c Input Variables
908     c ---------------
909     integer nymd, nhms, irun
910 molod 1.2 _RL getcon, cosz(irun), alat(irun), alon(irun), ra
911 molod 1.1
912     c Local Variables
913     c ---------------
914     integer year, day, sec, month, iday, idayp1
915     integer dayscy
916     integer i,nsecf,k,km,kp
917    
918     _RL hc
919     _RL pi, zero, one, two, six, dg2rd, yrlen, eqnx, ob, ecc, per
920     _RL daylen, fac, thm, thp, thnow, zs, zc, sj, cj
921    
922     parameter ( one = 1.0 )
923     parameter ( two = 2.0 )
924     parameter ( six = 6.0 )
925 molod 1.2 parameter ( dayscy = 365*4 + 1 )
926 molod 1.1
927     _RL TH(DAYSCY),T0,T1,T2,T3,T4,FUN,Y,MNDY(12,4)
928    
929     LOGICAL FIRST
930     DATA FIRST/.TRUE./
931     SAVE
932    
933     DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
934     . 397,34*0 /
935    
936 molod 1.2 FUN(Y,PI,ECC,YRLEN,PER) = (TWO*PI/((ONE-ECC**2)**1.5))*(ONE/YRLEN)
937 molod 1.1 . * (ONE - ECC*COS(Y-PER)) ** 2
938    
939     C***********************************************************************
940 molod 1.2 C* SET SOME CONSTANTS *
941     C***********************************************************************
942     pi = getcon('PI')
943     dg2rd = getcon('DEG2RAD')
944     yrlen = getcon('YRLEN')
945     ob = getcon('OBLDEG') * dg2rd
946     daylen = getcon('SDAY')
947     eqnx = getcon('VERNAL EQUINOX')
948     ecc = getcon('ECCENTRICITY')
949     per = getcon('PERIHELION') * dg2rd
950    
951     C***********************************************************************
952 molod 1.1 C* SET CURRENT TIME *
953     C***********************************************************************
954    
955     YEAR = NYMD / 10000
956     MONTH = MOD(NYMD,10000) / 100
957     DAY = MOD(NYMD,100)
958     SEC = NSECF(NHMS)
959    
960     C***********************************************************************
961     C* COMPUTE DAY-ANGLES FOR 4-YEAR CYCLE *
962     C***********************************************************************
963    
964     IF(FIRST) THEN
965     DO 100 I=15,48
966     MNDY(I,1) = MNDY(I-12,1) + 365
967     100 CONTINUE
968    
969     KM = INT(EQNX) + 1
970     FAC = KM-EQNX
971     T0 = ZERO
972 molod 1.2 T1 = FUN(T0,PI,ECC,YRLEN,PER )*FAC
973     T2 = FUN(ZERO+T1/TWO,PI,ECC,YRLEN,PER)*FAC
974     T3 = FUN(ZERO+T2/TWO,PI,ECC,YRLEN,PER)*FAC
975     T4 = FUN(ZERO+T3,PI,ECC,YRLEN,PER )*FAC
976 molod 1.1 TH(KM) = (T1 + TWO*(T2 + T3) + T4) / SIX
977    
978     DO 200 K=2,DAYSCY
979 molod 1.2 T1 = FUN(TH(KM),PI,ECC,YRLEN,PER )
980     T2 = FUN(TH(KM)+T1/TWO,PI,ECC,YRLEN,PER)
981     T3 = FUN(TH(KM)+T2/TWO,PI,ECC,YRLEN,PER)
982     T4 = FUN(TH(KM)+T3,PI,ECC,YRLEN,PER )
983 molod 1.1 KP = MOD(KM,DAYSCY) + 1
984     TH(KP) = TH(KM) + (T1 + TWO*(T2 + T3) + T4) / SIX
985     KM = KP
986     200 CONTINUE
987    
988     FIRST=.FALSE.
989     ENDIF
990    
991     C***********************************************************************
992     C* COMPUTE EARTH-SUN DISTANCE TO CURRENT SECOND *
993     C***********************************************************************
994    
995     IDAY = DAY + MNDY(MONTH,MOD(YEAR,4)+1)
996     IDAYP1 = MOD( IDAY,DAYSCY) + 1
997     THM = MOD( TH(IDAY) ,TWO*PI)
998     THP = MOD( TH(IDAYP1),TWO*PI)
999    
1000     IF(THP.LT.THM) THP = THP + TWO*PI
1001     FAC = FLOAT(SEC)/DAYLEN
1002     THNOW = THM*(ONE-FAC) + THP*FAC
1003    
1004     ZS = SIN(THNOW) * SIN(OB)
1005     ZC = SQRT(ONE-ZS*ZS)
1006     RA = (1.-ECC*ECC) / ( ONE-ECC*COS(THNOW-PER) )
1007    
1008     C***********************************************************************
1009     C* COMPUTE COSINE OF THE ZENITH ANGLE *
1010     C***********************************************************************
1011    
1012     FAC = FAC*TWO*PI + PI
1013     DO I = 1,IRUN
1014    
1015     HC = COS( FAC+ALON(I)*DG2RD )
1016     SJ = SIN(ALAT(I)*DG2RD)
1017     CJ = SQRT(ONE-SJ*SJ)
1018    
1019     COSZ(I) = SJ*ZS + CJ*ZC*HC
1020     IF( COSZ(I).LT.ZERO ) COSZ(I) = ZERO
1021     ENDDO
1022    
1023     RETURN
1024     END
1025    
1026     subroutine time_bound(nymd,nhms,nymd1,nhms1,nymd2,nhms2,imnm,imnp)
1027     C***********************************************************************
1028     C PURPOSE
1029     C Compute Date and Time boundaries.
1030     C
1031     C ARGUMENTS DESCRIPTION
1032     C nymd .... Current Date
1033     C nhms .... Current Time
1034     C nymd1 ... Previous Date Boundary
1035     C nhms1 ... Previous Time Boundary
1036     C nymd2 ... Subsequent Date Boundary
1037     C nhms2 ... Subsequent Time Boundary
1038     C
1039     C imnm .... Previous Time Index for Interpolation
1040     C imnp .... Subsequent Time Index for Interpolation
1041     C
1042     C***********************************************************************
1043    
1044     implicit none
1045     integer nymd,nhms, nymd1,nhms1, nymd2,nhms2
1046    
1047     c Local Variables
1048     c ---------------
1049     integer month,day,nyear,midmon1,midmon,midmon2
1050     integer imnm,imnp
1051     INTEGER DAYS(14), daysm, days0, daysp
1052     DATA DAYS /31,31,28,31,30,31,30,31,31,30,31,30,31,31/
1053    
1054     integer nmonf,ndayf,n
1055     NMONF(N) = MOD(N,10000)/100
1056     NDAYF(N) = MOD(N,100)
1057    
1058     C*********************************************************************
1059     C**** Find Proper Month and Time Boundaries for Climatological Data **
1060     C*********************************************************************
1061    
1062     MONTH = NMONF(NYMD)
1063     DAY = NDAYF(NYMD)
1064    
1065     daysm = days(month )
1066     days0 = days(month+1)
1067     daysp = days(month+2)
1068    
1069     c Check for Leap Year
1070     c -------------------
1071     nyear = nymd/10000
1072     if( 4*(nyear/4).eq.nyear ) then
1073     if( month.eq.3 ) daysm = daysm+1
1074     if( month.eq.2 ) days0 = days0+1
1075     if( month.eq.1 ) daysp = daysp+1
1076     endif
1077    
1078     MIDMON1 = daysm/2 + 1
1079     MIDMON = days0/2 + 1
1080     MIDMON2 = daysp/2 + 1
1081    
1082    
1083     IF(DAY.LT.MIDMON) THEN
1084     imnm = month
1085     imnp = month + 1
1086     nymd2 = (nymd/10000)*10000 + month*100 + midmon
1087     nhms2 = 000000
1088     nymd1 = nymd2
1089     nhms1 = nhms2
1090     call tick ( nymd1,nhms1, -midmon *86400 )
1091     call tick ( nymd1,nhms1,-(daysm-midmon1)*86400 )
1092     ELSE
1093     IMNM = MONTH + 1
1094     IMNP = MONTH + 2
1095     nymd1 = (nymd/10000)*10000 + month*100 + midmon
1096     nhms1 = 000000
1097     nymd2 = nymd1
1098     nhms2 = nhms1
1099     call tick ( nymd2,nhms2,(days0-midmon)*86400 )
1100     call tick ( nymd2,nhms2, midmon2*86400 )
1101     ENDIF
1102    
1103     c -------------------------------------------------------------
1104     c Note: At this point, imnm & imnp range between 01-14, where
1105     c 01 -> Previous years December
1106     c 02-13 -> Current years January-December
1107     c 14 -> Next years January
1108     c -------------------------------------------------------------
1109    
1110     imnm = imnm-1
1111     imnp = imnp-1
1112    
1113     if( imnm.eq.0 ) imnm = 12
1114     if( imnp.eq.0 ) imnp = 12
1115     if( imnm.eq.13 ) imnm = 1
1116     if( imnp.eq.13 ) imnp = 1
1117    
1118     return
1119     end
1120     subroutine time2freq2(MMDD,NYMD,NHMS,timeleft)
1121     C***********************************************************************
1122     C PURPOSE
1123     C COMPUTES TIME IN SECONDS UNTIL WE REACH THE NEXT MMDD
1124     C (ASSUME that the target time is 0Z)
1125     C
1126     C ARGUMENTS DESCRIPTION
1127     C MMDD FREQUENCY (MMDDHHMMSS)
1128     C NYMD CURRENT YYMMDD
1129     C NHMS CURRENT HHMMSS
1130     C TIMELEFT TIME LEFT (SECONDS)
1131     C
1132     C NOTES - Only used when the frequency is in units of months
1133     C Assumes that we always want to be at a month boundary
1134     C***********************************************************************
1135     implicit none
1136    
1137     integer mmdd,nymd,nhms,timeleft,daysleft
1138 molod 1.2 _RL getcon
1139 molod 1.1
1140     integer nsday
1141 molod 1.2 PARAMETER ( NSDAY = 86400 )
1142 molod 1.1 integer year, month, day, sec
1143     integer yearnext, monthnext, daynext, secnext
1144     integer i,nsecf2,nsecf,iday,idaynext,nsec,nsecnext,ntime
1145     integer testnymd
1146     integer MNDY(12,4)
1147     DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
1148     . 397,34*0 /
1149    
1150     C***********************************************************************
1151     C* COMPUTE # OF DAYS IN A 4-YEAR CYCLE *
1152     C***********************************************************************
1153     DO I=15,48
1154     MNDY(I,1) = MNDY(I-12,1) + 365
1155     ENDDO
1156     C***********************************************************************
1157     C* SET CURRENT TIME ELEMENTS *
1158     C***********************************************************************
1159     YEAR = NYMD / 10000
1160     MONTH = MOD(NYMD,10000) / 100
1161     DAY = MOD(NYMD,100)
1162     SEC = NSECF(NHMS)
1163     C***********************************************************************
1164     C* COMPUTE POSITIONS IN CYCLE FOR CURRENT AND BEGINNING TIMES *
1165     C***********************************************************************
1166     IDAY = (DAY -1) + MNDY( MONTH ,MOD(YEAR ,4)+1 )
1167     NSEC = IDAY *NSDAY + SEC
1168    
1169     testnymd=nymd + mmdd
1170     yearnext=testnymd/10000
1171     monthnext=mod(testnymd,10000)/100
1172     daynext = 1
1173     if(monthnext.gt.12) then
1174     monthnext = monthnext - 12
1175     yearnext = yearnext + 1
1176     endif
1177     testnymd = yearnext*10000 + monthnext*100 + daynext
1178     idaynext = mndy(monthnext,mod(yearnext,4)+1)
1179     daysleft = idaynext - iday
1180     if(daysleft.lt.0) daysleft = daysleft + 1461
1181    
1182     timeleft = daysleft * nsday - sec
1183    
1184     RETURN
1185     END

  ViewVC Help
Powered by ViewVC 1.1.22