#include "dims.h"
      SUBROUTINE OPFLUX
      use init_module,only: secs,sin_sundec,cos_sundec
      use cons_module,only: len1,dlamda,dphi,pi
      implicit none
C     ****
C     ****     CALCULATES UPWARD O+ NUMBER FLUX IN T7
C     ****
#include "params.h"
#include "vscr.h"
#include "phys.h"
#include "trgm.h"
!
! Local:
      real :: phid,phin,ppolar,rtd,rlat,coslat,sinlat
      integer :: i
C     ****
C     ****     SET PHID AND PHIN
C     ****
      DATA PHID,PHIN/+2.0E+8,-2.0E+8/
C     DATA PHID,PHIN/0.,0./
C
C     ****     SET PPOLAR (POLAR O+ FLUX)
C
C     DATA PPOLAR/+1.E+8/
C     DATA PPOLAR/+2.5E+7/
      DATA PPOLAR/0./
C
C
C     ****     RADIANS TO DEGREES
C
      RTD = 180./pi
C
C     ****     CALCULATE O+ FLUX AT UPPER BOUNDARY
C
C     ****     CALCULATE SOLAR ZENITH ANGLE, CHI
C
      RLAT=-.5*pi+(FLOAT(J-1)+.5)*dphi
      COSLAT=COS(RLAT)
      SINLAT=SIN(RLAT)
      DO 1 I = 1,LEN1
        T2(I)=FLOAT(I-3)
    1 CONTINUE
      DO 2 I=1,LEN1
C
C     ****     T2 = LOCAL TIME
C
        T2(I)=AMOD(SECS/3600.+(T2(I)*dlamda+pi)*12./pi,24.)
C
C     ****     T2 = CHI
C
        T2(I)=ACOS(sin_sundec*SINLAT+cos_sundec*COSLAT*
     |    COS(pi*(T2(I)-12.)
     1  /12.))
C
C     ****
C     ****     T3 = A = .5*(1.+SIN(PI*(ABS(RLATM)-PI/6.)/(PI/3.)))
C     ****               FOR  ABS(RLATM).LT.PI/3.
C     ****          A = 1.    FOR  ABS(RLATM).GE.PI/3
C     ****
C       if (ABS(RLATM(I,J))-pi/12.>=0.) then
	if (ABS(RLATM(I,J))-pi/24.>=0.) then
          t3(i) = 1.
        else
C         T3(I)=.5*(1.+SIN(pi*(ABS(RLATM(I,J))-pi/
C    |      24.)/(pi/12.)))
          T3(I)=.5*(1.+SIN(pi*(ABS(RLATM(I,J))-pi/
     |      48.)/(pi/24.)))
C         T3(I) = AMAX1(T3(I),0.1)
	  T3(I) = AMAX1(T3(I),0.05)
        endif
C
C     ****     T4 = FED       T5 = FEN
C
        T4(I) = PHID*T3(I)
        T5(I) = PHIN*T3(I)
C
C     ****     T1 = FE
C
        if (T2(I)-.5*pi>=0.) then
          T7(I) = T5(I)
        else
          T7(I) = T4(I)
        endif
        if ((T2(I)*RTD-80.)*(T2(I)*RTD-100.) < 0.)
     |    T7(I) = .5*(T4(I)+T5(I))+.5*(T4(I)-T5(I))*COS(pi*
     |      (T2(I)*RTD-80.)/20.)
C
C       ****     ADD PPOLAR IF MAGNETIC LATITUDE .GE. 60. DEG
C
        if (ABS(RLATM(I,J))-pi/3.>=0.) T7(I) = T7(I)+PPOLAR
    2 CONTINUE
      RETURN
      END
C