pkg/sun/sun_sfcsolz.F

27f9df093b Oliv*0001 #include "SUN_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: SUN_SFCSOLZ
                0005 C     !INTERFACE: ======================================================
                0006       SUBROUTINE SUN_SFCSOLZ(
                0007      O                        solz,
                0008      I                        isec,
                0009      I                        bi, bj, iMin, iMax, jMin, jMax,
                0010      I                        myTime, myIter, myThid )
                0011
                0012 C     !DESCRIPTION:
                0013 C     Computes solar zenith angle above sea surface
                0014
                0015 C     !USES: ===========================================================
                0016       IMPLICIT NONE
                0017 #include "SIZE.h"
                0018 #include "EEPARAMS.h"
                0019 #include "PARAMS.h"
6d4cc12339 Oliv*0020 #include "GRID.h"
27f9df093b Oliv*0021 #include "SUN_FIELDS.h"
                0022
                0023 C     !INPUT PARAMETERS: ===============================================
                0024 C     isec :: time of day in seconds when to compute the zenith angle
                0025 C             .LT. 0 means current time
                0026 C             .GE. 0 is the time in seconds since midnight
                0027 C             so 43200 is noon
                0028 C     myTime :: time at end of current (sub)timestep
                0029       _RL myTime
                0030       INTEGER isec, bi, bj, iMin, iMax, jMin, jMax, myIter, myThid
                0031
                0032 C     !OUTPUT PARAMETERS: ==============================================
                0033 C     solz :: solar zenith angle above sea surface in degrees
                0034       _RL solz(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
                0035 CEOP
                0036
                0037 #ifdef ALLOW_SUN
                0038
                0039 C     !FUNCTIONS: ======================================================
                0040       INTEGER SUN_JD
                0041       EXTERNAL SUN_JD
                0042
                0043 C     !LOCAL VARIABLES: ================================================
                0044       INTEGER i,j,l
                0045       INTEGER iyr,imon,iday,isecnow,lp,wd,myDate(4)
6d4cc12339 Oliv*0046       _RL sec, suni(3), sung(3), sunv, sunn, sune, rs, gha, rlon
27f9df093b Oliv*0047       _RL rjd, t
                0048       _RL xls, gs, xlm, asc, dpsi, eps
                0049
                0050 C  Get current date and time of day: iyr/imon/iday+isecnow
                0051       CALL CAL_GETDATE( myIter,myTime,myDate,myThid )
                0052       CALL CAL_CONVDATE( mydate,iyr,imon,iday,isecnow,lp,wd,myThid )
                0053       IF ( isec .GE. 0 )THEN
                0054 C       overwrite time of day as requested
                0055         sec = isec
                0056       ELSE
                0057         sec = isecnow
                0058       ENDIF
                0059
                0060 C  Compute floating point days since Jan 1.5, 2000
                0061 C  Note that the Julian day starts at noon on the specified date
                0062       rjd = SUN_JD(iyr,imon,iday)
                0063       t = rjd - 2451545 _d 0 + (sec-43200 _d 0)/86400 _d 0
                0064
                0065 C  Compute solar ephemeris parameters
                0066       CALL SUN_EPHPARMS (t,xls,gs,xlm,asc)
                0067
                0068 C  Compute nutation corrections
                0069       CALL SUN_NUTATE (t,xls,gs,xlm,asc,dpsi,eps)
                0070
                0071 C  Compute unit sun vector in geocentric inertial coordinates
                0072       CALL SUN_SUN2000 (t,xls,gs,xlm,asc,dpsi,eps,suni,rs)
                0073
                0074 C  Get Greenwich mean sidereal angle
                0075       CALL SUN_GHA2000 (t,dpsi,eps,gha)
                0076       gha = gha*deg2rad
                0077
6d4cc12339 Oliv*0078       IF ( isec .GE. 0 )THEN
                0079
                0080        DO j=jMin,jMax
                0081         DO i=iMin,iMax
                0082
                0083 C        isec is in local time, need to transform gha
                0084          rlon = XC(i,j,bi,bj)*deg2rad
                0085
27f9df093b Oliv*0086 C  Transform Sun vector into geocentric rotating frame
6d4cc12339 Oliv*0087          sung(1) = suni(1)*COS(gha-rlon) + suni(2)*SIN(gha-rlon)
                0088          sung(2) = suni(2)*COS(gha-rlon) - suni(1)*SIN(gha-rlon)
                0089          sung(3) = suni(3)
27f9df093b Oliv*0090
                0091 C  Compute components of spacecraft and sun vector in the
                0092 C  vertical (up), North (no), and East (ea) vectors frame
6d4cc12339 Oliv*0093          sunv = 0 _d 0
                0094          sunn = 0 _d 0
                0095          sune = 0 _d 0
                0096          DO l=1,3
                0097           sunv = sunv + sung(l)*SUN_up(i,j,bi,bj,l)
                0098           sunn = sunn + sung(l)*SUN_no(i,j,bi,bj,l)
                0099           sune = sune + sung(l)*SUN_ea(i,j,bi,bj,l)
                0100          ENDDO
                0101
                0102 C  Compute the solar zenith
                0103          solz(i,j) = ATAN2(SQRT(sunn*sunn+sune*sune), sunv)/deg2rad
27f9df093b Oliv*0104         ENDDO
6d4cc12339 Oliv*0105        ENDDO
                0106
                0107       ELSE
                0108
                0109 C  Transform Sun vector into geocentric rotating frame
                0110        sung(1) = suni(1)*COS(gha) + suni(2)*SIN(gha)
                0111        sung(2) = suni(2)*COS(gha) - suni(1)*SIN(gha)
                0112        sung(3) = suni(3)
                0113
                0114        DO j=jMin,jMax
                0115         DO i=iMin,iMax
                0116 C  Compute components of spacecraft and sun vector in the
                0117 C  vertical (up), North (no), and East (ea) vectors frame
                0118          sunv = 0 _d 0
                0119          sunn = 0 _d 0
                0120          sune = 0 _d 0
                0121          DO l=1,3
                0122           sunv = sunv + sung(l)*SUN_up(i,j,bi,bj,l)
                0123           sunn = sunn + sung(l)*SUN_no(i,j,bi,bj,l)
                0124           sune = sune + sung(l)*SUN_ea(i,j,bi,bj,l)
                0125          ENDDO
27f9df093b Oliv*0126
                0127 C  Compute the solar zenith
6d4cc12339 Oliv*0128          solz(i,j) = ATAN2(SQRT(sunn*sunn+sune*sune), sunv)/deg2rad
                0129         ENDDO
27f9df093b Oliv*0130        ENDDO
6d4cc12339 Oliv*0131
                0132       ENDIF
27f9df093b Oliv*0133
                0134 #endif /* ALLOW_SUN */
                0135
                0136       RETURN
                0137       END