astro.for

天文辐射计算

** $Id: astro.for 1.2 1997/10/02 15:20:17 LEM release $
      SUBROUTINE ASTRO (IDAY,LAT,AVRAD,
     &                  DAYL,DAYLP,SINLD,COSLD,DIFPP,ATMTR,DSINBE)

*  Purpose: This subroutine calculates astronomic daylength,
*           diurnal radiation characteristics such as the atmospheric
*           transmission, diffuse radiation etc.. This routine has
*           been modified so that it uses arrays to hold some input
*           output variables for faster processing 

*  FORMAL PARAMETERS:  (I=input,O=output,C=control,IN=init,T=time)
*  name   type meaning                                    units  class
*  ----   ---- -------                                    -----  -----
*  IDAY    I4  Day number (Jan 1st = 1)                      -      I
*  LAT     R4  Latitude of the site                       degrees   I
*  AVRAD   R4  Daily shortwave radiation                  J m-2 d-1 I
*  DAYL    R4  Astronomical daylength (base = 0 degrees)     h      O
*  DAYLP   R4  Astronomical daylength (base =-4 degrees)     h      O
*  SINLD   R4  Seasonal offset of sine of solar height       -      O
*  COSLD   R4  Amplitude of sine of solar height             -      O
*  DIFPP   R4  Diffuse irradiation perpendicular to direction of     
*              light                                      J m-2 s-1 O
*  ATMTR   R4  Daily atmospheric transmission                -      O
*  DSINBE  R4  Daily total of effective solar height         s      O

*  FATAL ERROR CHECKS: none
*  SUBROUTINES and FUNCTIONS called : none
*  FILE usage : none

*  Authors: Daniel van Kraalingen
*  Date   : April 1991


*     formal parameters
      IMPLICIT REAL(A-Z)
      INTEGER IDAY

**    local parameters
      REAL DECA(366),SCA(366)
      REAL LATA(366),AVRADA(366),DAYLA(366),DAYLPA(366)
      REAL SINLDA(366),COSLDA(366),DIFPPA(366),DSNBE(366),ATMTRA(366)

      PARAMETER (PI=3.1415926, ANGLE=-4., RAD=0.0174533)
      SAVE

*     initialize arrays
      DATA SCA /366*0./, LATA /366*-99./, AVRADA /366*0./

*     dependents of daynumber only
      IF (SCA(IDAY).EQ.0.) THEN
*        if solar constant for that day has not been calculated
         DEC = -ASIN(SIN(23.45*RAD)*COS(2.*PI*(REAL(IDAY)+10.)/365.))
         SC  = 1370.*(1.+0.033*COS(2.*PI*REAL(IDAY)/365.))
         DECA(IDAY) = DEC
         SCA(IDAY)  = SC
      END IF

*     dependents of latitude and global shortwave irradiation
      IF (LATA(IDAY).NE.LAT.OR.AVRADA(IDAY).NE.AVRAD) THEN
*        if latitude and irradiation at IDAY are different from previous
*        calls
         LATA(IDAY)   = LAT
         AVRADA(IDAY) = AVRAD
         SINLD = SIN(RAD*LAT)*SIN(DECA(IDAY))
         COSLD = COS(RAD*LAT)*COS(DECA(IDAY))
         AOB   = SINLD/COSLD

*        daylengths
         DAYL  = 12.0*(1.+2.*ASIN(AOB)/PI)
         DAYLP = 12.0*(1.+2.*ASIN((-SIN(ANGLE*RAD)+SINLD)/COSLD)/PI)

*        integrals of sine of solar height
         DSINB  = 3600.*(DAYL*SINLD+24.*COSLD*SQRT(1.-AOB**2)/PI)
         DSINBE = 3600.*(DAYL*(SINLD+0.4*(SINLD**2+COSLD**2*0.5))+
     &             12.*COSLD*(2.+3.*0.4*SINLD)*SQRT(1.-AOB**2)/PI)

*        extraterrestrial radiation and atmospheric transmission
         ANGOT  = SCA(IDAY)*DSINB
         ATMTR  = AVRAD/ANGOT

*        estimate fraction diffuse irradiation
         IF (ATMTR.GT.0.75) FRDIF = 0.23
         IF (ATMTR.LE.0.75.AND.ATMTR.GT.0.35) FRDIF = 1.33-1.46*ATMTR
         IF (ATMTR.LE.0.35.AND.ATMTR.GT.0.07)
     &       FRDIF = 1.-2.3*(ATMTR-0.07)**2
         IF (ATMTR.LE.0.07) FRDIF = 1.

         DIFPP = FRDIF*ATMTR*0.5*SCA(IDAY)

*        store in arrays for future use
         SINLDA(IDAY) = SINLD
         COSLDA(IDAY) = COSLD
         DAYLA(IDAY)  = DAYL
         DAYLPA(IDAY) = DAYLP
         DSNBE(IDAY)  = DSINBE
         ATMTRA(IDAY) = ATMTR
         DIFPPA(IDAY) = DIFPP
      ELSE
*        values have been calculated during a previous call to ASTRO
         SINLD  = SINLDA(IDAY)
         COSLD  = COSLDA(IDAY)
         DAYL   = DAYLA(IDAY)
         DAYLP  = DAYLPA(IDAY)
         DSINBE = DSNBE(IDAY)
         ATMTR  = ATMTRA(IDAY)
         DIFPP  = DIFPPA(IDAY)
      END IF

      RETURN
      END