pkg/chronos/chronos.F

22cc8d854e Andr*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003       subroutine set_alarm (tag,date,time,freq)
                0004 C***********************************************************************
                0005 C  Purpose
                0006 C  -------
                0007 C     Utility to Set Internal Alarms
                0008 C
                0009 C  Argument Description
                0010 C  --------------------
                0011 C     tag ....... Character String Tagging Alarm Process
                0012 C     date ...... Begining Date for Alarm
                0013 C     time ...... Begining Time for Alarm
                0014 C     freq ...... Repeating Frequency Interval for Alarm
                0015 C
                0016 C***********************************************************************
                0017
                0018       implicit none
                0019       character*(*) tag
                0020       integer       freq,date,time
                0021
                0022 #ifdef ALLOW_USE_MPI
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "EESUPPORT.h"
                0026 #endif
                0027
                0028 #include "chronos.h"
                0029
                0030 #ifdef ALLOW_USE_MPI
                0031 c MPI Utilities
                0032 c -------------
                0033 #include "mpif.h"
                0034       integer  mpi_comm_model,ierror
                0035 #endif
                0036
                0037       integer myid
                0038       logical first,set
                0039       data          first /.true./
                0040
                0041       integer n
                0042 #ifdef ALLOW_USE_MPI
                0043       call mpi_comm_rank ( mpi_comm_model,myid,ierror )
                0044 #else
                0045       myid = 1
                0046 #endif
                0047
                0048       if(first) then
                0049          ntags    = 1
                0050           tags(1) = tag
                0051          freqs(1) = freq
                0052          dates(1) = date
                0053          times(1) = time
                0054          if( myid.eq.1 ) write(6,100) date,time,freq,tags(1)
                0055       else
                0056
                0057       set = .false.
                0058       do n=1,ntags
                0059        if(tag.eq.tags(n)) then
                0060         if( myid.eq.1 ) then
                0061          print *, 'Warning!  Alarm has already been set for Tag: ',tag
                0062          print *, 'Changing  Alarm Information:'
                0063          print *, 'Frequency: ',freqs(n),' (Old) ',freq,' (New)'
                0064          print *, '    Date0: ',dates(n),' (Old) ',date,' (New)'
                0065          print *, '    Time0: ',times(n),' (Old) ',time,' (New)'
                0066         endif
                0067         freqs(n) = freq
                0068         dates(n) = date
                0069         times(n) = time
                0070         set = .true.
                0071        endif
                0072       enddo
                0073       if(.not.set) then
                0074             ntags = ntags+1
                0075          if(ntags.gt.maxtag ) then
                0076             if( myid.eq.1 ) then
                0077             print *, 'Too many Alarms are Set!!'
                0078             print *, 'Maximum Number of Alarms = ',maxtag
                0079             endif
                0080          call my_finalize
                0081          call my_exit (101)
                0082          endif
                0083           tags(ntags) = tag
                0084          freqs(ntags) = freq
                0085          dates(ntags) = date
                0086          times(ntags) = time
                0087          if( myid.eq.1 ) write(6,100) date,time,freq,tags(ntags)
                0088       endif
                0089       endif
                0090
                0091       first = .false.
                0092   100 format(1x,'Setting Alarm for: ',i8,2x,i6.6,',  with frequency: ',
                0093      .       i8,', and Tag: ',a80)
                0094       return
                0095       end
                0096
                0097       subroutine get_alarm (tag,date,time,freq,tleft)
                0098 C***********************************************************************
                0099 C  Purpose
                0100 C  -------
                0101 C     Utility to Get Internal Alarm Information
                0102 C
                0103 C  Input
                0104 C  -----
                0105 C     tag ....... Character String Tagging Alarm Process
                0106 C
                0107 C  Output
                0108 C  ------
                0109 C     date ...... Begining  Date for Alarm
                0110 C     time ...... Begining  Time for Alarm
                0111 C     freq ...... Frequency Interval for Alarm
                0112 C     tleft ..... Time Remaining (seconds) before Alarm is TRUE
                0113 C
                0114 C***********************************************************************
                0115
                0116       implicit none
                0117       character*(*) tag
                0118       integer freq,date,time,tleft
                0119
                0120 #ifdef ALLOW_USE_MPI
                0121 #include "SIZE.h"
                0122 #include "EEPARAMS.h"
                0123 #include "EESUPPORT.h"
                0124 #endif
                0125
                0126 #include "chronos.h"
                0127
                0128 #ifdef ALLOW_USE_MPI
                0129 c MPI Utilities
                0130 c -------------
                0131 #include "mpif.h"
                0132       integer  mpi_comm_model,ierror
                0133 #endif
                0134
                0135       logical set,alarm
                0136       external alarm
                0137       integer myid,n,nalarm,nsecf
                0138
                0139 #ifdef ALLOW_USE_MPI
                0140       call mpi_comm_rank ( mpi_comm_model,myid,ierror )
                0141 #else
                0142       myid = 1
                0143 #endif
                0144
                0145       set = .false.
                0146       do n=1,ntags
                0147       if(tag.eq.tags(n)) then
                0148        freq  = freqs(n)
                0149        date  = dates(n)
                0150        time  = times(n)
                0151
                0152        if( alarm(tag) ) then
                0153        tleft = 0
                0154        else
                0155        call get_time (nymd,nhms)
                0156        tleft = nsecf(freq) - nalarm(freq,nymd,nhms,date,time )
                0157        endif
                0158
                0159        set = .true.
                0160       endif
                0161       enddo
                0162
                0163       if(.not.set) then
                0164       if( myid.eq.1 ) print *, 'Alarm has not been set for Tag: ',tag
                0165       freq  = 0
                0166       date  = 0
                0167       time  = 0
                0168       tleft = 0
                0169       endif
                0170
                0171       return
                0172       end
                0173
                0174       function alarm (tag)
                0175       implicit none
                0176       character*(*) tag
                0177       integer date,time
                0178       logical alarm
                0179 #include "chronos.h"
                0180
                0181       integer n,modalarm,nalarm,freq,date0,time0
                0182       modalarm(freq,date0,time0) = nalarm (freq,date,time,date0,time0 )
                0183
                0184       call get_time (date,time)
                0185
                0186       alarm = .false.
                0187       do n=1,ntags
                0188       if( tags(n).eq.tag  ) then
                0189           if( freqs(n).eq.0 ) then
                0190           alarm = (dates(n).eq.date) .and. (times(n).eq.time)
                0191           else
                0192           alarm = ( date.gt.dates(n) .or.
                0193      .             (date.eq.dates(n) .and. time.ge.times(n)) ) .and.
                0194      .              modalarm( freqs(n),dates(n),times(n) ).eq.0
                0195           endif
                0196       endif
                0197       enddo
                0198
                0199       return
                0200       end
                0201
                0202       subroutine set_time (date,time)
                0203       implicit none
                0204       integer  date,time
                0205
                0206 #ifdef ALLOW_USE_MPI
                0207 #include "SIZE.h"
                0208 #include "EEPARAMS.h"
                0209 #include "EESUPPORT.h"
                0210 #endif
                0211
                0212 #include "chronos.h"
                0213
                0214 #ifdef ALLOW_USE_MPI
                0215 c MPI Utilities
                0216 c -------------
                0217 #include "mpif.h"
                0218       integer  mpi_comm_model,ierror
                0219 #endif
                0220       integer myid
                0221
                0222 #ifdef ALLOW_USE_MPI
                0223       call mpi_comm_rank ( mpi_comm_model,myid,ierror )
                0224 #else
                0225       myid = 1
                0226 #endif
                0227       if(  myid.eq.1 ) then
                0228       print *, 'Setting Clock'
                0229       print *, 'Date: ',date
                0230       print *, 'Time: ',time
                0231       endif
                0232
                0233       nymd = date
                0234       nhms = time
                0235       return
                0236       end
                0237
                0238       subroutine get_time (date,time)
                0239       implicit none
                0240       integer date,time
                0241
                0242 #include "chronos.h"
                0243
                0244       date = nymd
                0245       time = nhms
                0246       return
                0247       end
                0248
                0249       function nsecf (nhms)
                0250 C***********************************************************************
                0251 C  Purpose
                0252 C     Converts NHMS format to Total Seconds
                0253 C
                0254 C***********************************************************************
                0255       implicit none
                0256       integer  nhms, nsecf
                0257       nsecf =  nhms/10000*3600 + mod(nhms,10000)/100*60 + mod(nhms,100)
                0258       return
                0259       end
                0260
                0261       function nhmsf (nsec)
                0262 C***********************************************************************
                0263 C  Purpose
                0264 C     Converts Total Seconds to NHMS format
                0265 C
                0266 C***********************************************************************
                0267       implicit none
                0268       integer  nhmsf, nsec
                0269       nhmsf =  nsec/3600*10000 + mod(nsec,3600)/60*100 + mod(nsec,60)
                0270       return
                0271       end
                0272
                0273       function nsecf2 (nhhmmss,nmmdd,nymd)
                0274 C***********************************************************************
                0275 C  Purpose
                0276 C     Computes the Total Number of seconds from NYMD using NHHMMSS & NMMDD
                0277 C
                0278 C  Arguments   Description
                0279 C     NHHMMSS  IntervaL Frequency (HHMMSS)
                0280 C     NMMDD    Interval Frequency (MMDD)
                0281 C     NYMD     Current  Date      (YYMMDD)
                0282 C
                0283 C  NOTE:
                0284 C     IF (NMMDD.ne.0), THEN HOUR FREQUENCY HH MUST BE < 24
                0285 C
                0286 C***********************************************************************
                0287       implicit none
                0288
                0289       integer nsecf2,nhhmmss,nmmdd,nymd
                0290
                0291       INTEGER NSDAY, NCYCLE
                0292       PARAMETER ( NSDAY  = 86400 )
                0293       PARAMETER ( NCYCLE = 1461*24*3600 )
                0294
                0295       INTEGER YEAR, MONTH, DAY
                0296
                0297       INTEGER  MNDY(12,4)
                0298       DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
                0299      .           397,34*0 /
                0300
                0301       integer nsecf,i,nsegm,nsegd,iday,iday2,nday
                0302
                0303 C***********************************************************************
                0304 C*                 COMPUTE # OF SECONDS FROM NHHMMSS                   *
                0305 C***********************************************************************
                0306
                0307       nsecf2 = nsecf( nhhmmss )
                0308
                0309       if( nmmdd.eq.0 ) return
                0310
                0311 C***********************************************************************
                0312 C*                 COMPUTE # OF DAYS IN A 4-YEAR CYCLE                 *
                0313 C***********************************************************************
                0314
                0315       DO I=15,48
                0316       MNDY(I,1) = MNDY(I-12,1) + 365
                0317       ENDDO
                0318
                0319 C***********************************************************************
                0320 C*                 COMPUTE # OF SECONDS FROM NMMDD                     *
                0321 C***********************************************************************
                0322
                0323       nsegm =     nmmdd/100
                0324       nsegd = mod(nmmdd,100)
                0325
                0326       YEAR   = NYMD / 10000
                0327       MONTH  = MOD(NYMD,10000) / 100
                0328       DAY    = MOD(NYMD,100)
                0329
                0330       IDAY   = MNDY( MONTH ,MOD(YEAR ,4)+1 )
                0331       month = month + nsegm
                0332       If( month.gt.12 ) then
                0333       month = month - 12
                0334       year = year + 1
                0335       endif
                0336       IDAY2  = MNDY( MONTH ,MOD(YEAR ,4)+1 )
                0337
                0338                     nday = iday2-iday
                0339       if(nday.lt.0) nday = nday + 1461
                0340                     nday = nday + nsegd
                0341
                0342       nsecf2 = nsecf2 + nday*nsday
                0343
                0344       return
                0345       end
                0346
                0347       subroutine fixdate (nymd)
                0348       implicit none
                0349       integer nymd
                0350
                0351 c Modify 6-digit YYMMDD for dates between 1950-2050
                0352 c -------------------------------------------------
                0353       if (nymd .lt. 500101) then
                0354         nymd = 20000000 + nymd
                0355       else if (nymd .le. 991231) then
                0356         nymd = 19000000 + nymd
                0357       endif
                0358
                0359       return
                0360       end
                0361
                0362       subroutine interp_time ( nymd ,nhms ,
                0363      .                         nymd1,nhms1, nymd2,nhms2, fac1,fac2 )
                0364 C***********************************************************************
                0365 C
                0366 C  PURPOSE:
                0367 C  ========
                0368 C    Compute interpolation factors, fac1 & fac2, to be used in the
                0369 C    calculation of the instantanious boundary conditions, ie:
                0370 C
                0371 C               q(i,j) = fac1*q1(i,j) + fac2*q2(i,j)
                0372 C    where:
                0373 C               q(i,j) => Boundary Data valid    at (nymd  , nhms )
                0374 C              q1(i,j) => Boundary Data centered at (nymd1 , nhms1)
                0375 C              q2(i,j) => Boundary Data centered at (nymd2 , nhms2)
                0376 C
                0377 C  INPUT:
                0378 C  ======
                0379 C    nymd     : Date (yymmdd) of Current Timestep
                0380 C    nhms     : Time (hhmmss) of Current Timestep
                0381 C    nymd1    : Date (yymmdd) of Boundary Data 1
                0382 C    nhms1    : Time (hhmmss) of Boundary Data 1
                0383 C    nymd2    : Date (yymmdd) of Boundary Data 2
                0384 C    nhms2    : Time (hhmmss) of Boundary Data 2
                0385 C
                0386 C  OUTPUT:
                0387 C  =======
                0388 C    fac1     : Interpolation factor for Boundary Data 1
                0389 C    fac2     : Interpolation factor for Boundary Data 2
                0390 C
                0391 C
                0392 C***********************************************************************
                0393       implicit none
                0394
                0395       integer nhms,nymd,nhms1,nymd1,nhms2,nymd2
                0396       _RL fac1,fac2
                0397
                0398       INTEGER  YEAR , MONTH , DAY , SEC
                0399       INTEGER  YEAR1, MONTH1, DAY1, SEC1
                0400       INTEGER  YEAR2, MONTH2, DAY2, SEC2
                0401
                0402       _RL time, time1, time2
                0403
                0404       INTEGER    DAYSCY
                0405       PARAMETER (DAYSCY = 365*4+1)
                0406
                0407       INTEGER MNDY(12,4)
                0408
                0409       LOGICAL FIRST
                0410       DATA    FIRST/.TRUE./
                0411
                0412       DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
                0413      .           397,34*0 /
                0414
                0415       integer i,nsecf
                0416
                0417 C***********************************************************************
                0418 C*                         SET TIME BOUNDARIES                         *
                0419 C***********************************************************************
                0420
                0421       YEAR   = NYMD / 10000
                0422       MONTH  = MOD(NYMD,10000) / 100
                0423       DAY    = MOD(NYMD,100)
                0424       SEC    = NSECF(NHMS)
                0425
                0426       YEAR1  = NYMD1 / 10000
                0427       MONTH1 = MOD(NYMD1,10000) / 100
                0428       DAY1   = MOD(NYMD1,100)
                0429       SEC1   = NSECF(NHMS1)
                0430
                0431       YEAR2  = NYMD2 / 10000
                0432       MONTH2 = MOD(NYMD2,10000) / 100
                0433       DAY2   = MOD(NYMD2,100)
                0434       SEC2   = NSECF(NHMS2)
                0435
                0436 C***********************************************************************
                0437 C*                    COMPUTE DAYS IN 4-YEAR CYCLE                     *
                0438 C***********************************************************************
                0439
                0440       IF(FIRST) THEN
                0441       DO I=15,48
                0442       MNDY(I,1) = MNDY(I-12,1) + 365
                0443       ENDDO
                0444       FIRST=.FALSE.
                0445       ENDIF
                0446
                0447 C***********************************************************************
                0448 C*                     COMPUTE INTERPOLATION FACTORS                   *
                0449 C***********************************************************************
                0450
                0451       time  = DAY  + MNDY(MONTH ,MOD(YEAR ,4)+1) + float(sec )/86400.
                0452       time1 = DAY1 + MNDY(MONTH1,MOD(YEAR1,4)+1) + float(sec1)/86400.
                0453       time2 = DAY2 + MNDY(MONTH2,MOD(YEAR2,4)+1) + float(sec2)/86400.
                0454
                0455       if( time .lt.time1 ) time  = time  + dayscy
                0456       if( time2.lt.time1 ) time2 = time2 + dayscy
                0457
                0458       fac1  = (time2-time)/(time2-time1)
                0459       fac2  = (time-time1)/(time2-time1)
                0460
                0461       RETURN
                0462       END
                0463
                0464       subroutine tick (nymd,nhms,ndt)
                0465 C***********************************************************************
                0466 C  Purpose
                0467 C     Tick the Date (nymd) and Time (nhms) by NDT (seconds)
                0468 C
                0469 C***********************************************************************
                0470       implicit none
                0471
                0472       integer nymd,nhms,ndt
                0473
                0474       integer nsec,nsecf,incymd,nhmsf
                0475
                0476       IF(NDT.NE.0) THEN
                0477       NSEC = NSECF(NHMS) + NDT
                0478
                0479       IF (NSEC.GT.86400)  THEN
                0480       DO WHILE (NSEC.GT.86400)
                0481       NSEC = NSEC - 86400
                0482       NYMD = INCYMD (NYMD,1)
                0483       ENDDO
                0484       ENDIF
                0485
                0486       IF (NSEC.EQ.86400)  THEN
                0487       NSEC = 0
                0488       NYMD = INCYMD (NYMD,1)
                0489       ENDIF
                0490
                0491       IF (NSEC.LT.00000)  THEN
                0492       DO WHILE (NSEC.LT.0)
                0493       NSEC = 86400 + NSEC
                0494       NYMD = INCYMD (NYMD,-1)
                0495       ENDDO
                0496       ENDIF
                0497
                0498       NHMS = NHMSF (NSEC)
                0499       ENDIF
                0500
                0501       RETURN
                0502       END
                0503
                0504       subroutine tic_time (mymd,mhms,ndt)
                0505 C***********************************************************************
                0506 C  PURPOSE
                0507 C     Tick the Clock by NDT (seconds)
                0508 C
                0509 C***********************************************************************
                0510       implicit none
                0511 #include "chronos.h"
                0512
                0513       integer mymd,mhms,ndt
                0514
                0515       integer nsec,nsecf,incymd,nhmsf
                0516
                0517       IF(NDT.NE.0) THEN
                0518       NSEC = NSECF(NHMS) + NDT
                0519
                0520       IF (NSEC.GT.86400)  THEN
                0521       DO WHILE (NSEC.GT.86400)
                0522       NSEC = NSEC - 86400
                0523       NYMD = INCYMD (NYMD,1)
                0524       ENDDO
                0525       ENDIF
                0526
                0527       IF (NSEC.EQ.86400)  THEN
                0528       NSEC = 0
                0529       NYMD = INCYMD (NYMD,1)
                0530       ENDIF
                0531
                0532       IF (NSEC.LT.00000)  THEN
                0533       DO WHILE (NSEC.LT.0)
                0534       NSEC = 86400 + NSEC
                0535       NYMD = INCYMD (NYMD,-1)
                0536       ENDDO
                0537       ENDIF
                0538
                0539       NHMS = NHMSF (NSEC)
                0540       ENDIF
                0541
                0542 c Pass Back Current Updated Time
                0543 c ------------------------------
                0544       mymd = nymd
                0545       mhms = nhms
                0546
                0547       RETURN
                0548       END
                0549
                0550       FUNCTION NALARM (MHMS,NYMD,NHMS,NYMD0,NHMS0)
                0551 C***********************************************************************
                0552 C  PURPOSE
                0553 C     COMPUTES MODULO-FRACTION BETWEEN MHHS AND TOTAL TIME
                0554 C  USAGE
                0555 C  ARGUMENTS   DESCRIPTION
                0556 C     MHMS     INTERVAL FREQUENCY (HHMMSS)
                0557 C     NYMD     CURRENT   YYMMDD
                0558 C     NHMS     CURRENT   HHMMSS
                0559 C     NYMD0    BEGINNING YYMMDD
                0560 C     NHMS0    BEGINNING HHMMSS
                0561 C
                0562 C***********************************************************************
                0563       implicit none
                0564
                0565       integer nalarm,MHMS,NYMD,NHMS,NYMD0,NHMS0
                0566
                0567       integer nsday, ncycle
                0568       PARAMETER ( NSDAY  = 86400 )
                0569       PARAMETER ( NCYCLE = 1461*24*3600 )
                0570
                0571       INTEGER YEAR, MONTH, DAY, SEC, YEAR0, MONTH0, DAY0, SEC0
                0572
                0573       integer MNDY(12,4)
                0574       DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
                0575      .           397,34*0 /
                0576
                0577       integer i,nsecf,iday,iday0,nsec,nsec0,ntime
                0578
                0579 C***********************************************************************
                0580 C*                 COMPUTE # OF DAYS IN A 4-YEAR CYCLE                 *
                0581 C***********************************************************************
                0582
                0583       DO I=15,48
                0584       MNDY(I,1) = MNDY(I-12,1) + 365
                0585       ENDDO
                0586
                0587 C***********************************************************************
                0588 C*                   SET CURRENT AND BEGINNING TIMES                   *
                0589 C***********************************************************************
                0590
                0591       YEAR   = NYMD / 10000
                0592       MONTH  = MOD(NYMD,10000) / 100
                0593       DAY    = MOD(NYMD,100)
                0594       SEC    = NSECF(NHMS)
                0595
                0596       YEAR0  = NYMD0 / 10000
                0597       MONTH0 = MOD(NYMD0,10000) / 100
                0598       DAY0   = MOD(NYMD0,100)
                0599       SEC0   = NSECF(NHMS0)
                0600
                0601 C***********************************************************************
                0602 C*      COMPUTE POSITIONS IN CYCLE FOR CURRENT AND BEGINNING TIMES     *
                0603 C***********************************************************************
                0604
                0605       IDAY   = (DAY -1) + MNDY( MONTH ,MOD(YEAR ,4)+1 )
                0606       IDAY0  = (DAY0-1) + MNDY( MONTH0,MOD(YEAR0,4)+1 )
                0607
                0608       NSEC   = IDAY *NSDAY + SEC
                0609       NSEC0  = IDAY0*NSDAY + SEC0
                0610
                0611                        NTIME  = NSEC-NSEC0
                0612       IF (NTIME.LT.0 ) NTIME  = NTIME + NCYCLE
                0613                        NALARM = NTIME
                0614       IF ( MHMS.NE.0 ) NALARM = MOD( NALARM,NSECF(MHMS) )
                0615
                0616       RETURN
                0617       END
                0618
                0619       FUNCTION INCYMD (NYMD,M)
                0620 C***********************************************************************
                0621 C  PURPOSE
                0622 C     INCYMD:  NYMD CHANGED BY ONE DAY
                0623 C     MODYMD:  NYMD CONVERTED TO JULIAN DATE
                0624 C  DESCRIPTION OF PARAMETERS
                0625 C     NYMD     CURRENT DATE IN YYMMDD FORMAT
                0626 C     M        +/- 1 (DAY ADJUSTMENT)
                0627 C
                0628 C***********************************************************************
                0629       implicit none
                0630       integer incymd,nymd,m
                0631
                0632       integer ny,nm,nd,ny00,modymd
                0633
                0634       INTEGER NDPM(12)
                0635       DATA NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
                0636       LOGICAL LEAP
                0637       DATA NY00 /1900 /
                0638       LEAP(NY) = MOD(NY,4).EQ.0 .AND. (NY.NE.0 .OR. MOD(NY00,400).EQ.0)
                0639
                0640 C***********************************************************************
                0641 C
                0642       NY = NYMD / 10000
                0643       NM = MOD(NYMD,10000) / 100
                0644       ND = MOD(NYMD,100) + M
                0645
                0646       IF (ND.EQ.0) THEN
                0647       NM = NM - 1
                0648       IF (NM.EQ.0) THEN
                0649           NM = 12
                0650           NY = NY - 1
                0651       ENDIF
                0652       ND = NDPM(NM)
                0653       IF (NM.EQ.2 .AND. LEAP(NY))  ND = 29
                0654       ENDIF
                0655
                0656       IF (ND.EQ.29 .AND. NM.EQ.2 .AND. LEAP(NY))  GO TO 20
                0657
                0658       IF (ND.GT.NDPM(NM)) THEN
                0659       ND = 1
                0660       NM = NM + 1
                0661       IF (NM.GT.12) THEN
                0662           NM = 1
                0663           NY = NY + 1
                0664       ENDIF
                0665       ENDIF
                0666
                0667    20 CONTINUE
                0668       INCYMD = NY*10000 + NM*100 + ND
                0669
                0670       RETURN
                0671
                0672 C***********************************************************************
                0673 C                      E N T R Y    M O D Y M D
                0674 C***********************************************************************
                0675
                0676       ENTRY MODYMD (NYMD)
                0677
                0678       NY = NYMD / 10000
                0679       NM = MOD(NYMD,10000) / 100
                0680       ND = MOD(NYMD,100)
                0681
                0682    40 CONTINUE
                0683       IF (NM.LE.1)  GO TO 60
                0684       NM = NM - 1
                0685       ND = ND + NDPM(NM)
                0686       IF (NM.EQ.2 .AND. LEAP(NY))  ND = ND + 1
                0687       GO TO 40
                0688
                0689    60 CONTINUE
                0690       MODYMD = ND
                0691
                0692       RETURN
                0693       END
                0694
                0695       SUBROUTINE ASTRO ( NYMD,NHMS,ALAT,ALON,IRUN,COSZ,RA )
                0696 C***********************************************************************
                0697 C
                0698 C  INPUT:
                0699 C  ======
                0700 C    NYMD      : CURRENT YYMMDD
                0701 C    NHMS      : CURRENT HHMMSS
                0702 C    ALAT(IRUN):LATITUDES  IN DEGREES.
                0703 C    ALON(IRUN):LONGITUDES IN DEGREES. (0 = GREENWICH, + = EAST).
                0704 C    IRUN      : # OF POINTS TO CALCULATE
                0705 C
                0706 C  OUTPUT:
                0707 C  =======
                0708 C    COSZ(IRUN)  : COSINE OF ZENITH ANGLE.
                0709 C    RA          : EARTH-SUN DISTANCE IN UNITS OF
                0710 C                  THE ORBITS SEMI-MAJOR AXIS.
                0711 C
                0712 C  NOTE:
                0713 C  =====
                0714 C  THE INSOLATION AT THE TOP OF THE ATMOSPHERE IS:
                0715 C
                0716 C  S(I) = (SOLAR CONSTANT)*(1/RA**2)*COSZ(I),
                0717 C
                0718 C  WHERE:
                0719 C  RA AND COSZ(I) ARE THE TWO OUTPUTS OF THIS SUBROUTINE.
                0720 C
                0721 C***********************************************************************
                0722
                0723       implicit none
                0724
                0725 c Input Variables
                0726 c ---------------
                0727       integer nymd, nhms, irun
                0728       _RL    cosz(irun), alat(irun), alon(irun), ra
                0729
                0730 c Local Variables
                0731 c ---------------
                0732       integer year, day, sec, month, iday, idayp1
                0733       integer dayscy
                0734       integer i,nsecf,k,km,kp
                0735
                0736       _RL hc
                0737       _RL pi, zero, one, two, six, dg2rd, yrlen, eqnx, ob, ecc, per
                0738       _RL daylen, fac, thm, thp, thnow, zs, zc, sj, cj
                0739
                0740       parameter ( pi    = 3.1415926535898)
                0741       parameter ( zero  = 0.0 )
                0742       parameter ( one   = 1.0 )
                0743       parameter ( two   = 2.0 )
                0744       parameter ( six   = 6.0 )
                0745       parameter ( dg2rd = pi/180. )
                0746
                0747       parameter ( yrlen  = 365.25  )
                0748       parameter ( dayscy = 365*4+1 )
                0749       parameter ( eqnx   =  80.9028)
                0750       parameter ( ob     =  23.45*dg2rd )
                0751       parameter ( ecc    =   0.0167 )
                0752       parameter ( per    = 102.0*dg2rd)
                0753       parameter ( daylen = 86400.)
                0754
                0755       _RL      TH(DAYSCY),T0,T1,T2,T3,T4,FUN,Y,MNDY(12,4)
                0756
                0757       LOGICAL FIRST
                0758       DATA    FIRST/.TRUE./
                0759       SAVE
                0760
                0761       DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,
                0762      .           397,34*0 /
                0763
                0764       FUN(Y) = (TWO*PI/((ONE-ECC**2)**1.5))*(ONE/YRLEN)
                0765      .       * (ONE - ECC*COS(Y-PER)) ** 2
                0766
                0767 C***********************************************************************
                0768 C*                          SET CURRENT TIME                           *
                0769 C***********************************************************************
                0770
                0771       YEAR  = NYMD / 10000
                0772       MONTH = MOD(NYMD,10000) / 100
                0773       DAY   = MOD(NYMD,100)
                0774       SEC   = NSECF(NHMS)
                0775
                0776 C***********************************************************************
                0777 C*                 COMPUTE DAY-ANGLES FOR 4-YEAR CYCLE                 *
                0778 C***********************************************************************
                0779
                0780       IF(FIRST) THEN
                0781            DO 100 I=15,48
                0782            MNDY(I,1) = MNDY(I-12,1) + 365
                0783 100        CONTINUE
                0784
                0785            KM  = INT(EQNX) + 1
                0786            FAC = KM-EQNX
                0787            T0 = ZERO
                0788            T1 = FUN(T0         )*FAC
                0789            T2 = FUN(ZERO+T1/TWO)*FAC
                0790            T3 = FUN(ZERO+T2/TWO)*FAC
                0791            T4 = FUN(ZERO+T3    )*FAC
                0792            TH(KM) = (T1 + TWO*(T2 + T3) + T4) / SIX
                0793
                0794            DO 200 K=2,DAYSCY
                0795            T1 = FUN(TH(KM)       )
                0796            T2 = FUN(TH(KM)+T1/TWO)
                0797            T3 = FUN(TH(KM)+T2/TWO)
                0798            T4 = FUN(TH(KM)+T3    )
                0799            KP = MOD(KM,DAYSCY) + 1
                0800            TH(KP) = TH(KM) + (T1 + TWO*(T2 + T3) + T4) / SIX
                0801            KM = KP
                0802  200       CONTINUE
                0803
                0804            FIRST=.FALSE.
                0805       ENDIF
                0806
                0807 C***********************************************************************
                0808 C*            COMPUTE EARTH-SUN DISTANCE TO CURRENT SECOND             *
                0809 C***********************************************************************
                0810
                0811       IDAY   = DAY + MNDY(MONTH,MOD(YEAR,4)+1)
                0812       IDAYP1 = MOD( IDAY,DAYSCY) + 1
                0813       THM    = MOD( TH(IDAY)  ,TWO*PI)
                0814       THP    = MOD( TH(IDAYP1),TWO*PI)
                0815
                0816       IF(THP.LT.THM) THP = THP + TWO*PI
                0817       FAC   = FLOAT(SEC)/DAYLEN
                0818       THNOW = THM*(ONE-FAC) + THP*FAC
                0819
                0820       ZS = SIN(THNOW) * SIN(OB)
                0821       ZC = SQRT(ONE-ZS*ZS)
                0822       RA = (1.-ECC*ECC) / ( ONE-ECC*COS(THNOW-PER) )
                0823
                0824 C***********************************************************************
                0825 C*                 COMPUTE COSINE OF THE ZENITH ANGLE                  *
                0826 C***********************************************************************
                0827
                0828       FAC  = FAC*TWO*PI + PI
                0829       DO I = 1,IRUN
                0830
                0831       HC = COS( FAC+ALON(I)*DG2RD )
                0832       SJ = SIN(ALAT(I)*DG2RD)
                0833       CJ = SQRT(ONE-SJ*SJ)
                0834
                0835           COSZ(I) = SJ*ZS + CJ*ZC*HC
                0836       IF( COSZ(I).LT.ZERO ) COSZ(I) = ZERO
                0837       ENDDO
                0838
                0839       RETURN
                0840       END
                0841
                0842       subroutine time_bound(nymd,nhms,nymd1,nhms1,nymd2,nhms2,imnm,imnp)
                0843 C***********************************************************************
                0844 C  PURPOSE
                0845 C     Compute Date and Time boundaries.
                0846 C
                0847 C  ARGUMENTS   DESCRIPTION
                0848 C     nymd .... Current    Date
                0849 C     nhms .... Current    Time
                0850 C     nymd1 ... Previous   Date Boundary
                0851 C     nhms1 ... Previous   Time Boundary
                0852 C     nymd2 ... Subsequent Date Boundary
                0853 C     nhms2 ... Subsequent Time Boundary
                0854 C
                0855 C     imnm .... Previous   Time Index for Interpolation
                0856 C     imnp .... Subsequent Time Index for Interpolation
                0857 C
                0858 C***********************************************************************
                0859
                0860       implicit none
                0861       integer  nymd,nhms, nymd1,nhms1, nymd2,nhms2
                0862
                0863 c Local Variables
                0864 c ---------------
                0865       integer  month,day,nyear,midmon1,midmon,midmon2
                0866       integer  imnm,imnp
                0867       INTEGER  DAYS(14), daysm, days0, daysp
                0868       DATA     DAYS /31,31,28,31,30,31,30,31,31,30,31,30,31,31/
                0869
                0870       integer nmonf,ndayf,n
                0871       NMONF(N) = MOD(N,10000)/100
                0872       NDAYF(N) = MOD(N,100)
                0873
                0874 C*********************************************************************
                0875 C**** Find Proper Month and Time Boundaries for Climatological Data **
                0876 C*********************************************************************
                0877
                0878       MONTH  = NMONF(NYMD)
                0879       DAY    = NDAYF(NYMD)
                0880
                0881       daysm  = days(month  )
                0882       days0  = days(month+1)
                0883       daysp  = days(month+2)
                0884
                0885 c Check for Leap Year
                0886 c -------------------
                0887       nyear = nymd/10000
                0888       if( 4*(nyear/4).eq.nyear ) then
                0889       if( month.eq.3 ) daysm = daysm+1
                0890       if( month.eq.2 ) days0 = days0+1
                0891       if( month.eq.1 ) daysp = daysp+1
                0892       endif
                0893
                0894       MIDMON1 = daysm/2 + 1
                0895       MIDMON  = days0/2 + 1
                0896       MIDMON2 = daysp/2 + 1
                0897
                0898
                0899       IF(DAY.LT.MIDMON) THEN
                0900          imnm = month
                0901          imnp = month + 1
                0902          nymd2 = (nymd/10000)*10000 + month*100 + midmon
                0903          nhms2 = 000000
                0904          nymd1 = nymd2
                0905          nhms1 = nhms2
                0906          call tick ( nymd1,nhms1,       -midmon  *86400 )
                0907          call tick ( nymd1,nhms1,-(daysm-midmon1)*86400 )
                0908       ELSE
                0909          IMNM = MONTH + 1
                0910          IMNP = MONTH + 2
                0911          nymd1 = (nymd/10000)*10000 + month*100 + midmon
                0912          nhms1 = 000000
                0913          nymd2 = nymd1
                0914          nhms2 = nhms1
                0915          call tick ( nymd2,nhms2,(days0-midmon)*86400 )
                0916          call tick ( nymd2,nhms2,       midmon2*86400 )
                0917       ENDIF
                0918
                0919 c -------------------------------------------------------------
                0920 c Note:  At this point, imnm & imnp range between 01-14, where
                0921 c        01    -> Previous years December
                0922 c        02-13 -> Current  years January-December
                0923 c        14    -> Next     years January
                0924 c -------------------------------------------------------------
                0925
                0926       imnm = imnm-1
                0927       imnp = imnp-1
                0928
                0929       if( imnm.eq.0  ) imnm = 12
                0930       if( imnp.eq.0  ) imnp = 12
                0931       if( imnm.eq.13 ) imnm = 1
                0932       if( imnp.eq.13 ) imnp = 1
                0933
                0934       return
                0935       end