      SUBROUTINE GFACAVG

      INCLUDE 'PARAMETERS.INC'

      INCLUDE 'PLANET.INC'

      INCLUDE 'PARAMS_GFAC.INC'

      DOUBLE PRECISION COSSZA(NTRAP,NLAT)

      COMMON

     &      /BLOCK6/  MCROSS(MB)

     &      /BLOCK12/ CROSS(NWAVE,MB)

     &      /BLOCK14/ MTAU(MB)

     &      /BLOCK16/ COSSZA(NTRAP,NLAT)  

     &      /BLOCK17/ NANGLE(NALT,NLAT)

     &      /BLOCK18/ CONC(MB,NALT,NLAT) , OLDCONC(MB,NALT,NLAT)

     &      /BLOCK19/ TRANS(NWAVE,NALT,NLAT)

     &      /BLOCK30/ NPI2(NLAT) 

     &      /BLOCK41/ S(NTRAP,NWAVE,NALT,NLAT)

     &      /END/ NSTEPS , NMOL , NCALC , IDELT

     &      /PHOTOE/  GFACTOR(NALT,NLAT,6)

      REAL SZA,F10P7,GFACTEMP,GFACDIU(NTRAP,NALT,NLAT)
      INTEGER FEATURE



C THIS SUBROUTINE COMPUTES THE DIURNAL AVERAGE OF EXP(-TAU). IT CAN BE RUN 
C IN TWO MODES. SETTING ISPHERE = 0 PRODUCES A PLANE-PARALLEL ATMOSPHERE. THIS
C MODE DIFFERS FROM THE SECANT APPROXIMATION IN THAT THE DIFFERENTIAL
C SLANT PATH IS APPROXIMATED BY A MORE ACCURATE FORMULA FOR HIGH SOLAR 
C ZENITH ANGLES ( THE SECANT BLOWS UP WHEN CHI = PI/2 ). SETTING ISPHERE = 1
C PRODUCES THE FULL RAY-TRACE CALCULATION FOR A SPHERICAL ATMOSPHERE.  
 
 
C COMPUTE G-FACTORS FOR SZA < 90 DEGREES

C CALCULATE DIFFERENTIAL OPTICAL DEPTH AT KTH LAYER
      F10P7 = 150.              !CONTROL LEVEL OF SOLAR ACTIVITY
      DO IFEAT=1,3              ! LOOP THROUGH 3 FEATURES
      DO L = 1 , NLAT
        DO K = 36 , NALT
          COLDEN = CONC(NMOL,K,L)*4E5   !ASSUME 4 KM SCALE HT. (A KLUGE)
          DO NANG = 1 , NPI2(L)
           SZA = ACOS(COSSZA(NANG,L))*180./3.1416     !SZA IN DEGREES
           CALL GFACDRIVER ('GFAC.DAT',SZA,F10P7,
     &                GFACTEMP,COLDEN,IFEAT,1)
          GFACDIU(NANG,K,L) = GFACTEMP               !FOR DIURNAL AVERAGING
          ENDDO
        ENDDO

      ENDDO     ! END LOOP OVER LATITUDE            

C############################################################################


C CALCULATE DIURNAL AVERAGE OF EACH GFACTOR
C USE TRAPEZOIDAL RULE TO DO TIME INTEGRATION

      DO L = 1 , NLAT
      NEND = NPI2(L) - 1

        DO K = 36 , NALT
            SUM = 0.0D0
            DO NANG = 1 , NEND
              SUM = SUM + 0.5D0*(GFACDIU(NANG,K,L) + 
     &                           GFACDIU(NANG+1,K,L))
            ENDDO

            IF( NEND .EQ. 0 ) THEN  ! SPECIAL CASE
              SUM = 0.5D0*GFACDIU(1,K,L)
            ENDIF

            GFACTOR(K,L,IFEAT) = SUM/(NTRAP-1.0)
    
        ENDDO 

      ENDDO     ! END LOOP OVER LATITUDE
      ENDDO     ! END LOOP OVER IFEAT  

C******** NOW DO O2 PHOTOIONIZATION (FROM DIFFERENT FILE) ******
C THIS MOSTLY LY-BETA, EXCEPT AT HIGH SUN (LOW SZA) WHEN OTHER THINGS ADD IN
C

      IFEAT = 4
      DO L = 1 , NLAT
        DO K = 36 , NALT
          COLDEN = CONC(NMOL,K,L)*6E5     !ASSUME 6 KM SCALE HT.
          DO NANG = 1 , NPI2(L)
           SZA = ACOS(COSSZA(NANG,L))*180./3.1416     !SZA IN DEGREES
           CALL GFACDRIVER ('PHOION.DAT',SZA,F10P7,
     &                GFACTEMP,COLDEN,IFEAT,1)
           GFACDIU(NANG,K,L) = GFACTEMP              !FOR DIURNAL AVERAGING
          ENDDO
        ENDDO

      ENDDO     ! END LOOP OVER LATITUDE   
C
C AS ABOVE, USE TRAPEZOIDAL RULE AND CALC DIURNAL AVERAGE
C
      DO L = 1 , NLAT
      NEND = NPI2(L) - 1

        DO K = 36 , NALT
            SUM = 0.0D0
            DO NANG = 1 , NEND
              SUM = SUM + 0.5D0*(GFACDIU(NANG,K,L) + 
     &                           GFACDIU(NANG+1,K,L))
            ENDDO

            IF( NEND .EQ. 0 ) THEN  ! SPECIAL CASE
              SUM = 0.5D0*GFACDIU(1,K,L)
            ENDIF

            GFACTOR(K,L,IFEAT) = SUM/(NTRAP-1.0)
    
        ENDDO 

      ENDDO     ! END LOOP OVER LATITUDE
  
      RETURN
      END

