SUBROUTINE SETUP

#include "com2d.h"

C     DECLINATION, ZENITH ANGLE, 
C     CHAPMAN FUNCTION AT LATITUDE LATP AND ELAPSED TIME T.
C         RESET DECLINATION WITH TIME IF REQUESTED. FORMULA AND
C         CONSTANTS ARE FROM THE ASTRONOMICAL ALMANAC, WITH
C         MEANL=MEAN SOLAR LONGITUDE (CORRECTED FOR ABERRATION),
C         MEANA=MEAN ANOMALY, ECLIPL=ECLIPTIC LONGITUDE.
          MEANL=DEC1+DEC2*DAY
          MEANA=DEC3+DEC4*DAY
          ECLIPL=MEANL+DEC5*SIN(MEANA)+DEC6*SIN(2.*MEANA)
          SIND=DEC7*SIN(ECLIPL)
          R=ASIN(SIND)
          COSD=COS(R)
          DECD=R*RTD
c          print *,' dec1,dec2,dec3,dec4,dec5,dec6,meanl,meana,
c     *  eclipl,dec7,dec8,rtd,r',dec1,dec2,dec3,dec4,dec5,
c     *  dec6,meanl,meana,eclipl,dec7,dec8,rtd,r
c          print *,' day,decd,sind,cosd=',day,decd,sind,cosd

C           CORRECT SOLAR FLUXES AT THE TOP OF THE ATMOSPHERE FOR
C           THE CURRENT ORBITAL POSITION OF THE EARTH.
            CORB=1.0/(1.0-DEC8*COS(MEANA))**2

            DO 200 JL=1,IL$
             SOLCYC=1.0
            FLUX(JL)=CORB*SOLCYC*FLUX0(JL)
200	continue
C
            DO 2100 IJ=1,L$
            DO 2100 IK=1,Z$
            GCR(IJ,IK)=GCRMAX(IJ,IK)
2100        CONTINUE

	DO 3000 IJ=1,L$
	CALL FILL(IJ,Z$)

      SINL=SIN(LATP*DTR)
      COSL=COS(LATP*DTR)+1.e-20

c      print *,' ij latp cosl sinl',ij,latp,cosl,sinl

C     SOLAR ZENITH ANGLE CHID.
      COSBC=SIND*SINL+COSD*COSL*COS(OMEGAD*DTR*(CLOCK+12.0))
      R=ACOS(COSBC)
      SINC=SIN(R)
      CHID(IJ)=R*RTD

C     GRAZING ALTITUDE ZGRZ. AT ALTITUDES FOR WHICH IT IS ZERO OR
C     LESS, THE SUN IS NOT UP; THE J COEFFICIENTS SHOULD BE SET TO
C     ZERO, AND THE CHAPMAN FUNCTION IS NOT REQUIRED.
      IF(CHID(IJ).LE.90.0)THEN
        ZGRZ(IJ)=1000.0
      ELSE
        ZGRZ(IJ)=SINC*(ZP+R0)-R0
      END IF
c	type *,'ij,zgrz(ij),zp,r0,chid(ij)',ij,zgrz(ij),zp,r0,chid(ij)
C     CHAPMAN FUNCTION. APPROXIMATE IT WITH THE SECANT(CHI) FOR
C     CHI LESS THAN 70 DEGREES.
      IF(CHID(IJ).LT.70.0)THEN       
        SFAC(IJ)=1.0/COSBC
c	print *,' ij=',ij,' cosbc=',cosbc,' sfac=',sfac(ij)
      ELSE IF(ZGRZ(IJ).GT.0.0)THEN
C       CHAPMAN FUNCTION SFAC FOR CHID GREATER THAN OR EQUAL TO 70
C       DEGREES WHEN ZGRZ IS POSITIVE.
        R=SQRT(0.50*R0/HBAR)
        S=R*ABS(COSBC)

        IF(S.LE.8.0)THEN
          S=(D1+D2*S)/(D3+D4*S+S**2)
        ELSE
          S=D5/(D6+S)
        END IF

        R=R*SQRT(PI)
        SFAC(IJ)=R*S
        IF(CHID(IJ).GT.90.0)SFAC(IJ)=2.0*R*EXP((R0+ZBAR)*(
     .  1.0-SINC)/HBAR) - SFAC(IJ)
      ELSE
        SFAC(IJ)=-9999.0
      END IF

C     ADJUST THE LENGTH OF DAY AND LENGTH OF NIGHT FOR THE
C     CURRENT LATITUDE. TOTAL DAYLIGHT=24*ACOS(-TAN(L)*TAN(D))/PI
C     FROM RUNDEL (1977), WITH L=LATITUDE, D=DECLINATION.
        R=SINL*SIND/(COSL*COSD)
        IF(R.GE.1.e0)THEN
          TDF=1.
        ELSE IF(R.LE.-1.)THEN
          TDF=0.
        ELSE
          TDF=ACOS(-R)*RPI
        END IF
        TNF=1.-TDF
       TFD(IJ)=TDF
c	print *,' ij,sfac, zgrz tfd=',ij,sfac(ij),zgrz(ij),
c     *  tfd(ij)
3000	CONTINUE

C	print *,' sfac=',sfac
C	print *,' zgrz=',zgrz

      SAVE
      RETURN
      END