pkg/radtrans/radtrans_declination_spencer.F

b55e95f1ff Oliv*0001 #include "RADTRANS_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: RADTRANS_RDECLINATION_SPENCER
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE RADTRANS_DECLINATION_SPENCER(
                0009      O                            delta,
                0010      I                            myTime, myIter, myThid )
                0011
                0012 C     !DESCRIPTION:
                0013 C     Compute declination for current date using Spencer's Fourier series
                0014
                0015 C     !USES:
                0016       IMPLICIT NONE
                0017 #include "SIZE.h"
                0018 #include "EEPARAMS.h"
                0019 #include "PARAMS.h"
                0020 #include "cal.h"
                0021
                0022 C     !INPUT PARAMETERS:
                0023       _RL myTime
                0024       INTEGER myIter, myThid
                0025
                0026 C     !OUTPUT PARAMETERS:
                0027 C     delta :: declination in radians
                0028       _RL delta
                0029 CEOP
                0030
                0031 #ifdef ALLOW_RADTRANS
                0032
                0033 C     !LOCAL VARIABLES:
                0034       INTEGER year0,myDate(4),yearStartDate(4),difftime(4)
                0035       _RL secondsInYear, myDateSeconds, yearangle
                0036
                0037 C find how many seconds have passed in the current year
                0038
                0039       CALL CAL_GETDATE( myIter, myTime, myDate, myThid )
                0040
                0041       year0            = int(myDate(1)/10000.)
                0042       yearStartDate(1) = year0 * 10000 + 101
                0043       yearStartDate(2) = 0
                0044       yearStartDate(3) = myDate(3)
                0045       yearStartDate(4) = myDate(4)
                0046       CALL cal_TimePassed(yearStartDate,myDate,difftime,myThid)
                0047       CALL cal_ToSeconds (difftime,myDateSeconds,myThid)
                0048
                0049       IF ( myDate(3) .EQ. 2 ) THEN
                0050         secondsInYear = ndaysleap * secondsperday
                0051       ELSE
                0052         secondsInYear = ndaysnoleap * secondsperday
                0053       ENDIF
                0054
                0055 C convert to radians
                0056       yearangle = 2.0 _d 0*PI*myDateSeconds/secondsInYear
                0057
                0058 C approximate declination according to Spencer 1971
                0059
                0060       delta = 0.006918 _d 0
                0061      &      - 0.399912 _d 0*cos(yearangle)
                0062      &      + 0.070257 _d 0*sin(yearangle)
                0063      &      - 0.006758 _d 0*cos(2.0 _d 0*yearangle)
                0064      &      + 0.000907 _d 0*sin(2.0 _d 0*yearangle)
                0065      &      - 0.002697 _d 0*cos(3.0 _d 0*yearangle)
                0066      &      + 0.001480 _d 0*sin(3.0 _d 0*yearangle)
                0067
                0068 #endif /* ALLOW_RADTRANS */
                0069
                0070       RETURN
                0071       END
                0072