#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"
#include "grid.h"
!
! Local:
      real :: phid,phin,ppolar,rtd,rlat,coslat,sinlat
      integer :: i
      real,dimension(zimxp,zkmxp) :: 
     |  rlatm_ik,                     ! for debug (redundant in vertical)
     |  t2_ik,t3_ik,t4_ik,t5_ik,t7_ik ! redundant in vertical and latitude
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/
      DATA PPOLAR/+2.5E+7/
C     DATA PPOLAR/0./

! Save rlatm_ik (redundant in vertical):
      do i=1,zimxp
        rlatm_ik(i,:) = rlatm(i,j)
      enddo
!     call addfsech('RLATM',' ',' ',rlatm_ik,zimxp,zkmxp,zkmx,j)
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     ****
	if (ABS(RLATM(I,J))-pi/24.>=0.) then
          t3(i) = 1.
        else
          T3(I)=.5*(1.+SIN(pi*(ABS(RLATM(I,J))-pi/
     |      48.)/(pi/24.)))
	  T3(I) = AMAX1(T3(I),0.05)

!         write(6,"('opflux: j=',i3,' lat=',f10.3,' i=',i3,' t3(i)=',
!    |      e12.4,' rlatm(i,j)=',e12.4)") j,glat(j),i,t3(i),rlatm(i,j)

        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

!     write(6,"('opflux: j=',i3,' lat=',f10.3,' t3=',/,(6e12.4))") 
!    |  j,glat(j),t3

! Save t2-5,t7 (redundant in vertical and latitude)
      do i=1,zimxp
        t2_ik(i,:) = t2(i)
        t3_ik(i,:) = t3(i)
        t4_ik(i,:) = t4(i)
        t5_ik(i,:) = t5(i)
        t7_ik(i,:) = t7(i)
      enddo
!     call addfsech('T2',' ',' ',t2_ik,zimxp,zkmxp,zkmx,j)
!     call addfsech('T3',' ',' ',t3_ik,zimxp,zkmxp,zkmx,j)
!     call addfsech('T4',' ',' ',t4_ik,zimxp,zkmxp,zkmx,j)
!     call addfsech('T5',' ',' ',t5_ik,zimxp,zkmxp,zkmx,j)
!     call addfsech('T7',' ',' ',t7_ik,zimxp,zkmxp,zkmx,j)

      RETURN
      END
C