#include "dims.h"
      SUBROUTINE OPFLUX
      use cons_module,only: len1,pi,dlamda,dphi
      use init_module,only: secs,cos_sundec,sin_sundec
      implicit none
C     ****
C     ****     CALCULATES UPWARD O+ NUMBER FLUX IN T7
C     ****
#include "params.h"
#include "vscr.h"
#include "phys.h"
      real :: rlatm,rlonm,dipmag,decmag,sndec,csdec,sn2dec,sncsdc,
     |  rlatmp,rlonmp,dlons
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),
     1DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),
     2CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),
     3RLATMP(ZJMX),RLONMP(ZJMX),DLONS(ZJMX)
!
! Local:
      integer :: i
      real :: phid,phin,ppolar,rtd,rlat,coslat,sinlat
      real :: tmp(zimxp,zkmxp) ! to save T7
C     ****
C     ****     SET PHID AND PHIN
C     ****
! 11/30/98: Changed phid,phin as per modsrc.kibo:
      phid =  2.e8
      phin = -3.e8
!     phid = -9.0e+8
!     phin = -3.0e+8
C
C     ****     SET PPOLAR (POLAR O+ FLUX)
C
! 11/30/98: Changed ppolar as per modsrc.kibo:
!     ppolar = 1.e8
      ppolar = 5.0e+7
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.)/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     ****
        t3(i) = 1.
        if (abs(rlatm(i,j)) < pi/12.) 
     |    t3(i) = .5*(1.+SIN(pi*(ABS(RLATM(I,J))-pi/48.)/
     |            (pi/24.)))
        T3(I) = AMAX1(T3(I),0.05)
C
C     ****     T4 = FED       T5 = FEN
C
        T4(I) = PHID*T3(I)
        T5(I) = PHIN*T3(I)
C
C     ****     T1 = FE
C
        t7(i) = t5(i)
        if (t2(i) < .5*pi) t7(i) = t4(i)
        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       ****     ADD PPOLAR IF MAGNETIC LATITUDE .GE. 60. DEG
        if (abs(rlatm(i,j)) >= pi/3.)
     |    t7(i) = t7(i)+ppolar
    2 CONTINUE
      RETURN
      END