SUBROUTINE OPFLUX
      implicit none
C     ****
C     ****     CALCULATES UPWARD O+ NUMBER FLUX IN T7
C     ****
      include "params.h"
      include "vscr.h"
      include "cons.h"
      include "phys.h"
      include "strt.h"
      real :: rlatm,rlonm,dipmag,decmag,sndec,csdec,sn2dec,sncsdc,
     |  dumdum
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),
     1  DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),
     2  CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX),
     3  DUMDUM(ZJMX,3)
!
! Local:
      real :: phid,phin,ppolar,rtd,rlat,coslat,sinlat
      integer :: i
C     ****
C     ****     SET PHID AND PHIN
C     ****
C     DATA PHID,PHIN/1.5E+8,-1.5E+8/
C     DATA PHID,PHIN/0.,-3.0E+8/
C     DATA PHID,PHIN/0.,-1.5E+8/
      DATA PHID,PHIN/-1.5E+8,-3.0E+8/
C     DATA PHID,PHIN/-2.5E+8,-3.0E+8/
C     DATA PHID,PHIN/-5.0E+8,-6.0E+8/
C     DATA PHID,PHIN/-7.0E+8,-3.0E+8/
C     DATA PHID,PHIN/-9.0E+8,-3.0E+8/
C
C     ****     SET PPOLAR (POLAR O+ FLUX)
C
C     DATA PPOLAR/+1.E+8/
C     DATA PPOLAR/+5.E+7/
      DATA PPOLAR/0./
C
C
C     ****     RADIANS TO DEGREES
C
      RTD = 180./C(110)
C
C     ****     CALCULATE O+ FLUX AT UPPER BOUNDARY
C
C     ****     CALCULATE SOLAR ZENITH ANGLE, CHI
C
      RLAT=-.5*C(110)*FLOAT(1-NHEMI)+(FLOAT(J-1)+.5)*C(2)
      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)*C(1)+C(110))*12./C(110),24.)
C
C     ****     T2 = CHI
C
        T2(I)=ACOS(C(95)*SINLAT+C(96)*COSLAT*COS(C(110)*(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     ****
        T3(I)=merge(1.,.5*(1.+SIN(C(110)*(ABS(RLATM(I,J))-C(110)/
     1   24.)/(C(110)/12.))),ABS(RLATM(I,J))-C(110)/12.>=0.)
C
C     ****     T4 = FED       T5 = FEN
C
        T4(I) = PHID*T3(I)
        T5(I) = PHIN*T3(I)
C
C     ****     T1 = FE
C
        T7(I) = merge(T5(I),T4(I),T2(I)-.5*C(110)>=0.)
        T7(I) = merge(T7(I),.5*(T4(I)+T5(I))+.5*(T4(I)-T5(I))
     1  *COS(C(110)*(T2(I)*RTD-80.)/20.),(T2(I)*RTD-80.)
     2  *(T2(I)*RTD-100.)>=0.)
C
C       ****     ADD PPOLAR IF MAGNETIC LATITUDE .GE. 60. DEG
C
        T7(I) = merge(T7(I)+PPOLAR,T7(I),ABS(RLATM(I,J))-C(110)/3.>=0.)
    2 CONTINUE
      RETURN
      END
C