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

Contents 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 - (show annotations) (download)
Tue Mar 13 15:38:59 2007 UTC (17 years, 2 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 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 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 _RL getcon, fac1,fac2
464
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 parameter ( dayscy = 365*4 + 1 )
473
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 C DESCRIPTION OF INPUT VARIABLES
807 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 _RL getcon, cosz(irun), alat(irun), alon(irun), ra
911
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 parameter ( dayscy = 365*4 + 1 )
926
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 FUN(Y,PI,ECC,YRLEN,PER) = (TWO*PI/((ONE-ECC**2)**1.5))*(ONE/YRLEN)
937 . * (ONE - ECC*COS(Y-PER)) ** 2
938
939 C***********************************************************************
940 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 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 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 TH(KM) = (T1 + TWO*(T2 + T3) + T4) / SIX
977
978 DO 200 K=2,DAYSCY
979 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 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 _RL getcon
1139
1140 integer nsday
1141 PARAMETER ( NSDAY = 86400 )
1142 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