      SUBROUTINE OPLUS
C     ****
C     ****     ADVANCES N(O+) BY ONE TIME STEP
C     ****
C     ****     ON ENTRY, T7 = UPWARD NUMBER FLUX OF O+ AT TOP
C     ****
      include "params.h"
      include "blnk.h"
      include "vscr.h"
      include "cons.h"
      include "index.h"
      include "buff.h"
      include "phys.h"
      include "crates.h"
C     ****     COMMON BLOCK DESCRIBING MAGNETIC FIELD
      COMMON/MAGFLD/BXM(ZIMXP,2),BX(ZIMXP,ZJMX),BXP(ZIMXP,4),
     1  BY(ZIMXP,ZJMX),BYP(ZIMXP,4),BZ(ZIMXP,ZJMX),BZP(ZIMXP,4),
     2  BMOD(ZIMXP,ZJMX),BMODP(ZIMXP,2)
      DIMENSION KUTT(ZJMX)
      DATA KUTT/1,2,3,5,6,7,9,10,11,13,14,15,17,10*17,17,15,14,13,11,10,
     19,7,6,5,3,2,1/
C     ****     MOLECULAR WEIGHTS OF O+ AND N2D
      COMMON/MASS/RMN4S,RMN2D,RMNO,BRN2D,CEE
      DATA RMOP/16./
      DATA EXPLIC/1./
C     ****
C     ****      T1 = DIV(B)
C     ****
      CALL DIVB(T1,J)
C     ****
C     ****      S4 = TP(J-1), S5 = TP(J), S6 = TP(J+1)
C     ****
      NTJM1K = NJM1+NT
      NTK = NJ+NT
      NTJP1K = NJP1+NT
      NTEM1K = NJM1+NTE
      NTEK = NJ+NTE
      NTEP1K = NJP1+NTE
      NTIM1K = NJM1+NTI
      NTIK = NJ+NTI
      NTIP1K = NJP1+NTI
      DO 40 I = 1,LEN2
        S4(I,1) = .5*(F(I,NTEM1K)+F(I,NTIM1K))
        S5(I,1) = .5*(F(I,NTEK)+F(I,NTIK))
        S6(I,1) = .5*(F(I,NTEP1K)+F(I,NTIP1K))
   40 CONTINUE
C     ****
C     ****      S13 = D(J-1),  S14 = D(J),  S15 = D(J+1)
C     ****
      NMSM1K = NJM1+NMS
      NMSK = NJ+NMS
      NMSP1K = NJP1+NMS
      NPS1MK = NJM1+NPS
      NPS1K = NJ+NPS
      NPS1PK = NJP1+NPS
      NPS2MK = NJM1+NPS2
      NPS2K = NJ+NPS2
      NPS2PK = NJP1+NPS2
      CALL RRK(F(1,NTJM1K),F(1,NMSM1K),F(1,NPS1MK),F(1,NPS2MK),S4,S13)
      CALL RRK(F(1,NTK),F(1,NMSK),F(1,NPS1K),F(1,NPS2K),S5,S14)
      CALL RRK(F(1,NTJP1K),F(1,NMSP1K),F(1,NPS1PK),F(1,NPS2PK),S6,S15)
C     ****
C     ****      S4 = 2.*TP(J-1), S5 = 2.*TP(J), S6 = 2.*TP(J+1)
C     ****
      DO 41 I = 1,LEN2
        S4(I,1) = 2.*S4(I,1)
        S5(I,1) = 2.*S5(I,1)
        S6(I,1) = 2.*S6(I,1)
   41 CONTINUE
C     ****
C     ****      S10 = H(J-1), S11 = H(J), S12 = H(J+1)
C     ****
      DO 42 I = 1,LEN2
        S10(I,1) = C(57)*F(I,NTJM1K)/(.5*(F(I,NMSM1K)+F(I,NMSM1K+1))*
     1  C(54))
        S11(I,1) = C(57)*F(I,NTK)/(.5*(F(I,NMSK)+F(I,NMSK+1))*
     1  C(54))
        S12(I,1) = C(57)*F(I,NTJP1K)/(.5*(F(I,NMSP1K)+F(I,NMSP1K+1))*
     1  C(54))
   42 CONTINUE
C     ****
C     ****      S7 = (D/(H*DZ)*2.*TP+M*G/R)*N(O+)   (J-1)
C     ****      S8 = (D/(H*DZ)*2.*TP+M*G/R)*N(O+)     (J)
C     ****      S9 = (D/(H*DZ)*2.*TP+M*G/R)*N(O+)   (J+1)
C     ****
      NOPM1K = NJM1+NOP
      NOPK = NJ+NOP
      NOPP1K = NJP1+NOP
      CALL DIFFUS(S4,F(1,NOPM1K),S10,S7)
      CALL DIFFUS(S5,F(1,NOPK),S11,S8)
      CALL DIFFUS(S6,F(1,NOPP1K),S12,S9)
C     ****
C     ****      S9 = (B(H).DEL(H))*(D/(H*DZ)*TP+M*G/R)*N(O+)   (J)
C     ****
      CALL BDOTDH(S7,S8,S9,J,S9)
C     ****
C     ****      S9 = D*BZ*S9         (J)
C     ****
      DO 43 K = 1,KMAX
        DO 44 I = 1,LEN1
          S9(I,K) = S14(I,K)*BZ(I,J)*S9(I,K)
   44   CONTINUE
   43 CONTINUE
C     ****
C     ****      S3 = 2.*TP*N(O+)    (J-2)
C     ****      S4 = 2.*TP*N(O+)    (J-1)
C     ****      S5 = 2.*TP*N(O+)    (J)
C     ****      S6 = 2.*TP*N(O+)    (J+1)
C     ****      S7 = 2.*TP*N(O+)    (J+2)
C     ****
      NTEM2K = NJM2+NTE
      NTEP2K = NJP2+NTE
      NOPM2K = NJM2+NOP
      NOPP2K = NJP2+NOP
      NTIM2K = NJM2+NTI
      NTIP2K = NJP2+NTI
      DO 45 I = 1,LEN2
        S3(I,1) = F(I,NOPM2K)*(F(I,NTEM2K)+F(I,NTIM2K))
        S4(I,1) = S4(I,1)*F(I,NOPM1K)
        S5(I,1) = S5(I,1)*F(I,NOPK)
        S6(I,1) = S6(I,1)*F(I,NOPP1K)
        S7(I,1) = F(I,NOPP2K)*(F(I,NTEP2K)+F(I,NTIP2K))
   45 CONTINUE
C     ****
C     ****      S1 = (B(H).DEL(H))*2.*TP*N(O+)    (J-1)
C     ****      S2 = (B(H).DEL(H))*2.*TP*N(O+)      (J)
C     ****      S3 = (B(H).DEL(H))*2.*TP*N(O+)    (J+1)
C     ****
      CALL BDOTDH(S3,S4,S5,J-1,S1)
      CALL BDOTDH(S4,S5,S6,J,S2)
      CALL BDOTDH(S5,S6,S7,J+1,S3)
C     ****
C     ****      S1 = S1*D(J-1), S2 = S2*D(J), S3 = S3*D(J+1)
C     ****
      DO 46 I = 1,LEN2
        S1(I,1) = S1(I,1)*S13(I,1)
        S2(I,1) = S2(I,1)*S14(I,1)
        S3(I,1) = S3(I,1)*S15(I,1)
   46 CONTINUE
C     ****
C     ****      S8 = (B(H).DEL(H))*D*(B(H).DEL(H))*2.*TP*N(O+)   (J)
C     ****
      CALL BDOTDH(S1,S2,S3,J,S8)
C     ****
C     ****      S7 = (BZ*D/(H*DZ)+DIV(*B))*S2
C     ****
      CALL BDZDVB(S2,T1,S11,J,S7)
C     ****
C     ****      COLLECT NEW EXPLICIT TERMS IN S4
C     ****        ( RHS OF TRIDIAGONALSYSTEM )
C     ****
      DO 47 I = 1,LEN2
        S4(I,1) = -EXPLIC*(S7(I,1)+S8(I,1)+S9(I,1))
   47 CONTINUE
C     ****
C     ****      S9 = H(J-1), S10 = H(J), S11 = H(J+1)
C     ****
      DO 48 I=1,LEN2
        S9(I,1) = S10(I,1)
        S10(I,1) = S11(I,1)
        S11(I,1) = S12(I,1)
   48 CONTINUE
C     ****
C     ****      S6 = (B.U)*N(O+)    (J-1)
C     ****      S7 = (B.U)*N(O+)      (J)
C     ****      S8 = (B.U)*N(O+)    (J+1)
C     ****
      NUK = NJM1+NU-1
      NVK = NJM1+NV-1
      NWK = NJM1+NW-1
      NOPK = NJM1+NOP-1
      DO 49 K = 1,KMAX
        NUK = NUK+1
        NVK = NVK+1
        NWK = NWK+1
        NOPK = NOPK+1
        DO 50 I = 1,LEN1
          S6(I,K)=(BX(I,J-1)*F(I,NUK)+BY(I,J-1)*F(I,NVK)+S9(I,K)*
     1    BZ(I,J-1)*.5*(F(I,NWK)+F(I,NWK+1)))*F(I,NOPK)
   50   CONTINUE
   49 CONTINUE
      NUK = NJ+NU-1
      NVK = NJ+NV-1
      NWK = NJ+NW-1
      NOPK = NJ+NOP-1
      DO 51 K = 1,KMAX
        NUK = NUK+1
        NVK = NVK+1
        NWK = NWK+1
        NOPK = NOPK+1
        DO 52 I = 1,LEN1
          S7(I,K)=(BX(I,J)*F(I,NUK)+BY(I,J)*F(I,NVK)+S10(I,K)*
     1    BZ(I,J)*.5*(F(I,NWK)+F(I,NWK+1)))*F(I,NOPK)
   52   CONTINUE
   51 CONTINUE
      NUK = NJP1+NU-1
      NVK = NJP1+NV-1
      NWK = NJP1+NW-1
      NOPK = NJP1+NOP-1
      DO 53 K = 1,KMAX
        NUK = NUK+1
        NVK = NVK+1
        NWK = NWK+1
        NOPK = NOPK+1
        DO 54 I = 1,LEN1
          S8(I,K)=(BX(I,J+1)*F(I,NUK)+BY(I,J+1)*F(I,NVK)+S11(I,K)*
     1    BZ(I,J+1)*.5*(F(I,NWK)+F(I,NWK+1)))*F(I,NOPK)
   54   CONTINUE
   53 CONTINUE
C     ****
C     ****     NOW ADD TERMS TO F IN S4
C     ****
      NOPMK=NJM1+NOP-1
      NOPK=NJ+NOP-1
      NOPPK=NJP1+NOP-1
      NUIK=NDJ+NUI-1
      NVIK=NDJ+NVI-1
      DO 24 K=1,KMAX
      NOPMK=NOPMK+1
      NOPK=NOPK+1
      NOPPK=NOPPK+1
      NUIK=NUIK+1
      NVIK=NVIK+1
      DO 24 I = 3,IMAXP2
      S4(I,K)=S4(I,K)+1./(2.*C(51))*(1./(CS(J)*C(1))*
     1(BX(I,J)*(S7(I+1,K)-S7(I-1,K))+.5*(F(I,NUIK)+F(I,NUIK+1))*
     2BMOD(I,J)**2*(F(I+1,NOPK)/BMOD(I+1,J)**2-F(I-1,NOPK)/BMOD(I-1,J)
     3**2))+1./C(2)*(BY(I,J)*(S8(I,K)-S6(I,K))+.5*(F(I,NVIK)+F(I,NVIK+1)
     4)*BMOD(I,J)**2*(F(I,NOPPK)/BMOD(I,J+1)**2-F(I,NOPMK)/BMOD(I,J-1)**
     52)))
   24 CONTINUE
C     ****     PERIODIC POINTS
      DO 33 I=1,2
      DO 33 K=1,KMAX
      S4(I,K)=S4(I+IMAX,K)
      S4(I+IMAXP2,K)=S4(I+2,K)
   33 CONTINUE
C     ****
C     ****     T1 = DIV(B)/BZ
C     ****
      DO 1 I=1,LEN1
      T1(I)=T1(I)/BZ(I,J)
    1 CONTINUE
C     ****
C     ****      S7=D(J-1), S8=D(J), S9=D(J+1)
C     ****
      DO 2 I=1,LEN2
        S7(I,1)=S13(I,1)
        S8(I,1)=S14(I,1)
        S9(I,1)=S15(I,1)
    2 CONTINUE
C     ****
C     ****     S15 = 1./(H*DZ),      S14 = TP = .5*(TI+TE)
C     ****
      NTIK=NJ+NTI
      NTEK=NJ+NTE
      DO 3 I=1,LEN2
      S15(I,1)=1./(S10(I,1)*C(3))
      S14(I,1)=.5*(F(I,NTIK)+F(I,NTEK))
    3 CONTINUE
C     ****
C     ****     S13 = 2.*TP/(H(-)*DZ)+GMR     (K = 1,KMAXP1)
C     ****     S12 = 2.*TP/(H(+)*DZ)-GMR     (K = 0,KMAX)
C     ****
C     ****     GMR = G*M(0+)/(2.*R)
      GMR=C(54)*RMOP/(2.*C(57))
      DO 4 I=1,LEN2-LEN1
      S13(I,2)=2.*S14(I,2)*(.5*(S15(I,1)+S15(I,2)))+GMR
      S12(I,2)=2.*S14(I,1)*(.5*(S15(I,1)+S15(I,2)))-GMR
    4 CONTINUE
C     ****     UPPER AND LOWER BOUNDARIES
      DO 5 I=1,LEN1
      S13(I,1)=2.*S14(I,1)*(1.5*S15(I,1)-0.5*S15(I,2))+GMR
      S13(I,KMAXP1)=2.*(2.*S14(I,KMAX)-S14(I,KMAXM1))*(1.5*S15(I,KMAX)
     1-0.5*S15(I,KMAXM1))+GMR
      S12(I,1)=2.*(2.*S14(I,1)-S14(I,2))*(1.5*S15(I,1)-0.5*S15(I,2))-GMR
      S12(I,KMAXP1)=2.*S14(I,KMAX)*(1.5*S15(I,KMAX)-0.5*S15(I,KMAXM1))
     1-GMR
    5 CONTINUE
C     ****
C     ****     S11 = D(K+1/2)
C     ****
      DO 8 I=1,LEN2-LEN1
      S11 (I,2)=.5*(S8(I,1)+S8(I,2))
    8 CONTINUE
C     ****     UPPER AND LOWER BOUNDARIES
      DO 9 I=1,LEN1
      S11(I,1)=(1.5*S8(I,1)-0.5*S8(I,2))
      S11(I,KMAXP1)=(1.5*S8(I,KMAX)-.5*S8(I,KMAXM1))
    9 CONTINUE
C     ****
C     ****      S7 = (DIV(B)+(DH*D*BZ)/(D*BZ)
C     ****
C     ****      PERIODIC POINTS
C     ****
      DO 91 I = 1,2
        DO 92 K = 1,3*KMAXP1
          S9(I,K) = S9(I+IMAX,K)
          S9(I+IMAXP2,K) = S9(I+2,K)
   92   CONTINUE
   91 CONTINUE
      DO 93 K = 1,KMAX
        DO 94 I = 3,IMAXP2
          S7(I,K) = T1(I)+1./(C(51)*S8(I,K)*BZ(I,J)**2)*(BX(I,J)/CS(J)
     A              *(S8(I+1,K)*BZ(I+1,J)-S8(I-1,K)*BZ(I-1,J))/(2.*C(1))
     B              +BY(I,J)*(S9(I,K)*BZ(I,J+1)-S7(I,K)*BZ(I,J-1))
     C              /(2.*C(2)))
   94   CONTINUE
   93 CONTINUE
C     ****
C     ****      PERIODIC POINTS
C     ****
      DO 96 I = 1,2
        DO 97 K = 1,KMAXP1
          S7(I,K) = S7(I+IMAX,K)
          S7(I+IMAXP2,K) = S7(I+2,K)
   97   CONTINUE
   96 CONTINUE
C     ****
C     ****     S10 = (1./(H*DZ)-(DIV(B)+DH*D*BZ/(D*BZ))/(2*BZ))*BZ**2
C     ****     S9 = (1./(H*DZ)+(DIV(B)+DH*D*BZ/(D*BZ))/(2*BZ))*BZ**2
C     ****
      DO 10 K=1,KMAX
      DO 10 I=1,LEN1
      S10(I,K)=(S15(I,K)-.5*S7(I,K))*BZ(I,J)**2
      S9(I,K)=(S15(I,K)+.5*S7(I,K))*BZ(I,J)**2
   10 CONTINUE
C     ****
C     ****     SHAPIRO SMOOTHER
C     ****     SET S4 TO N(O+)/(2*DT)  (N-1)
C     ****
      NOJP2K=NJP2+NOPNM
      NOJP1K=NJP1+NOPNM
      NOJNMK=NJ+NOPNM
      NOJM1K=NJM1+NOPNM
      NOJM2K=NJM2+NOPNM
      DO 36 I=1,LEN2
      S1(I,1)=F(I,NOJNMK)-C(26)*(F(I,NOJP2K)+F(I,NOJM2K)-4.*(F(I,NOJP1K)
     1+F(I,NOJM1K))+6.*F(I,NOJNMK))
   36 CONTINUE
      DO 37 I=3,LEN2-2
      S4(I,1)=S4(I,1)-(S1(I,1)-C(26)*(S1(I+2,1)+S1(I-2,1)-4.*(S1(I+1,1)+
     1S1(I-1,1))+6.*S1(I,1)))*C(7)
   37 CONTINUE
C     ****
C     ****     BEGIN COEFFS P,Q,R IN S1,S2,S3
C     ****
      DO 11 I=1,LEN2
      S1(I,1)=S10(I,1)*S11(I,1)*S12(I,1)
      S2(I,1)=-(S9(I,1)*S11(I,2)*S12(I,2)+S10(I,1)*S11(I,1)*S13(I,1))
      S3(I,1)=S9(I,1)*S11(I,2)*S13(I,2)
   11 CONTINUE
C     ****
C     ****     S7 = (B.U)
C     ****
      DO 12 I=1,LEN2
C     ****     S6 = H
      S6(I,1)=1./(S15(I,1)*C(3))
   12 CONTINUE
      NUK=NJ+NU-1
      NVK=NJ+NV-1
      NWK=NJ+NW-1
      DO 13 K=1,KMAX
      NUK=NUK+1
      NVK=NVK+1
      NWK=NWK+1
      DO 13 I=1,LEN1
C     ****     S7 = B.U
      S7(I,K)=BX(I,J)*F(I,NUK)+BY(I,J)*F(I,NVK)+S6(I,K)*BZ(I,J)*.5*
     1(F(I,NWK)+F(I,NWK+1))
   13 CONTINUE
C     ****
C     ****     FINISH P AND R IN S1 AND S3
C     ****
      NWIK=NDJ+NWI-1
      DO 14 K=1,KMAXM1
      NWIK=NWIK+1
      DO 14 I=1,LEN1
      S1(I,K+1)=S1(I,K+1)+(BZ(I,J)*S7(I,K)+.5*(F(I,NWIK+1)+F(I,NWIK+2)))
     1*.5*S15(I,K+1)
      S2(I,K)=S2(I,K)-.5*(F(I,NWIK)+F(I,NWIK+1))*6./C(51)
      S3(I,K)=S3(I,K)-(BZ(I,J)*S7(I,K+1)+.5*(F(I,NWIK)+F(I,NWIK+1)))*.5*
     1S15(I,K)
   14 CONTINUE
C     ****     UPPER AND LOWER BOUNDARIES
      NWIK=NDJ+NWI
      DO 15 I=1,LEN1
      S1(I,1)=S1(I,1)+(BZ(I,J)*(2.*S7(I,1)-S7(I,2))+.5*(F(I,NWIK)+
     1F(I,NWIK+1)))*.5*S15(I,1)
      S2(I,KMAX)=S2(I,KMAX)-.5*(F(I,NWIK+KMAX)+F(I,NWIK+KMAXM1))*6./
     AC(51)
      S3(I,KMAX)=S3(I,KMAX)-(BZ(I,J)*(2.*S7(I,KMAX)-S7(I,KMAXM1))+.5*
     1(F(I,NWIK+KMAX)+F(I,NWIK+KMAXM1)))*.5*S15(I,KMAX)
   15 CONTINUE
C     ****
C     ****     ADDITIONS TO Q IN S2
C     ****
      DO 16 K=1,KMAX
      DO 16 I=1,LEN1
      S2(I,K)=S2(I,K)-S7(I,K)*T1(I)*BZ(I,J)-C(7)
   16 CONTINUE
C     ****
C     ****     UPPER BOUNDARY CONDITION
C     ****
C     ****     T2 = A, T3= B
      NWIK=NDJ+NWI+KMAX
      DO 17 I=1,LEN1
C     ****
C     ****     MODS FOR OPLUS -- UPPER BOUNDARY CONDITION
C     ****
C     T2(I) = .5*(2.*S7(I,KMAX)-S7(I,KMAXM1))*BZ(I,J)+.5*F(I,NWIK)
C     T2(I) = .5*(2.*S7(I,KMAX)-S7(I,KMAXM1))*BZ(I,J)
C     T2(I)=0. IS DIFFUSIVE EQUILIBRIUM
      T2(I) = 0.
C     T2(I) = 0.5*F(I,NWIK) IS OPFLUX AND WIND COMPONENT
C     T2(I) = 0.5*F(I,NWIK)
      T3(I) = -BZ(I,J)**2*S11(I,KMAXP1)*S12(I,KMAXP1)-T2(I)
      T2(I) = -BZ(I,J)**2*S11(I,KMAXP1)*S13(I,KMAXP1)+T2(I)
C     ****     S2 = Q = Q+B/A*R
      S2(I,KMAX)=S2(I,KMAX)+T3(I)/T2(I)*S3(I,KMAX)
C     ****     S4 = F = F -R/A*PHI
      S4(I,KMAX)=S4(I,KMAX)-T7(I)*S3(I,KMAX)/T2(I)
      S3(I,KMAX)=0.
   17 CONTINUE
C     ****
C     ****
C     ****     SOURCES AND SINKS
C     ****
      NTK=NJ+NT-1
      NMSK=NJ+NMS-1
      NPS1K=NJ+NPS-1
      NPS2K=NJ+NPS2-1
      NPN2DK=NJ+NPN2D-1
      NQOP2PK=NQOP2P-1
      NQOP2DK=NQOP2D-1
      NEK=NJ+NE-1
      NOPK=NJ+NOP-1
      DO 26 K=1,KMAX
      NTK=NTK+1
      NMSK=NMSK+1
      NPS1K=NPS1K+1
      NPS2K=NPS2K+1
      NPN2DK=NPN2DK+1
      NQOP2PK=NQOP2PK+1
      NQOP2DK=NQOP2DK+1
      NEK=NEK+1
      NOPK=NOPK+1
      DO 26 I=1,LEN1
C     ****     CALCULATE LOSS IN S13, SUBSTRACT FROM S2
      XIOP2P(I,K) = .5*(F(I,NQOP2PK)+F(I,NQOP2PK+1))/((C(81)*EXPS(K)
     1  *.5*(F(I,NMSK)+F(I,NMSK+1))/(C(84)*F(I,NTK))*((RK16(I,K)
     2  +RK17(I,K))*(1.-F(I,NPS1K)-F(I,NPS2K))/RMASS(3)+RK18(I,K)
     3  *F(I,NPS2K)/RMASS(2)))+(RK19(I,K)+RK20(I,K))*F(I,NEK)+
     4  RK21(I,K)+RK22(I,K))
      XIOP2D(I,K) = (.5*(F(I,NQOP2DK)+F(I,NQOP2DK+1))+(RK20(I,K)*
     1  F(I,NEK)+RK22(I,K))*XIOP2P(I,K))/((C(81)*EXPS(K)*.5*(F(I,NMSK)
     2  +F(I,NMSK+1))/(C(84)*F(I,NTK))*(RK23(I,K)*(1.-F(I,NPS1K)
     3  -F(I,NPS2K))/RMASS(3)+RK24(I,K)*F(I,NPS2K)/RMASS(2)+RK26(I,K)
     4  *F(I,NPS1K)/RMASS(1)))+RK25(I,K)*F(I,NEK)+RK27(I,K))
      S13(I,K)=C(81)*EXPS(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(C(84)*
     1F(I,NTK))*(RK1(I,K,1)*F(I,NPS1K)/RMASS(1)+RK2(I,K)*(1.-F(I,NPS1K)
     2-F(I,NPS2K))/RMASS(3)+RK10(I,K)*F(I,NPN2DK)/RMN2D)
      S2(I,K)=S2(I,K)-S13(I,K)
   26 CONTINUE
C     ****     ADD SOURCE TERM TO RHS (F) IN S4
      NQOPK=NQOP-1
      NEK=NJ+NE-1
      NPS2K=NJ+NPS2-1
      NMSK=NJ+NMS-1
      NTK=NJ+NT-1
      DO 27 K=1,KMAX
      NQOPK=NQOPK+1
      NEK=NEK+1
      NPS2K=NPS2K+1
      NMSK=NMSK+1
      NTK=NTK+1
      DO 27 I=1,LEN1
      S4(I,K)=S4(I,K)-.5*(F(I,NQOPK)+F(I,NQOPK+1))-(RK19(I,K)*
     1  F(I,NEK)+RK21(I,K))*XIOP2P(I,K)-(RK25(I,K)*F(I,NEK)+
     2  RK27(I,K))*XIOP2D(I,K)-(RK18(I,K)*XIOP2P(I,K)+RK24(I,K)*
     3  XIOP2D(I,K))*
     4  F(I,NPS2K)/RMASS(2)*
     5  C(81)*EXPS(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(C(84)*
     6  F(I,NTK))
   27 CONTINUE
C     ****     LOWER BOUNDARY CONDITION
C     ****
C     ****     N(O+) = Q/L
      NQOPK=NQOP
      DO 28 I=1,LEN1
      S2(I,1)=S2(I,1)-S1(I,1)
      S4(I,1)=S4(I,1)-2.*S1(I,1)*F(I,NQOPK)/(1.5*S13(I,1)-.5*S13(I,2))
      S1(I,1)=0.
   28 CONTINUE
C     ****
C     ****     SOLVE TRIDIAGONAL SYSTEM
C     ****
      NOPK=NJNP+NOP
      CALL TRSOLV(S1,S2,S3,S4,S7,F(1,NOPK),LEN1,3,LEN1-2,KMAXP1,1,KMAX,
     11)
C     ****
C     ****     FILTER N(O+)
C     ****
      IF(KUTT(J).LT.IMAXH)THEN
          NOPK=NJNP+NOP
          CALL FILTER2(NOPK,KMAX,KUTT(J))
      ENDIF
C     ****     TIME SMOOTHING
      NOPNMK=NJ+NOPNM
      NOPK=NJ+NOP
      NOPNPK=NJNP+NOP
      NOPMNK=NJNP+NOPNM
      DO 29 I=1,LEN2
      F(I,NOPMNK)=C(30)*F(I,NOPK)+C(31)*(F(I,NOPNMK)+F(I,NOPNPK))
   29 CONTINUE
C     ****     PERIODIC POINTS
      K1=NJNP+NOP
      K2=K1+KMAXM1
      DO 30 N=1,2
      DO 31 I=1,2
      DO 31 K=K1,K2
      F(I,K)=F(I+IMAX,K)
      F(I+IMAXP2,K)=F(I+2,K)
   31 CONTINUE
      K1=NJNP+NOPNM
      K2=K1+KMAXM1
   30 CONTINUE
C     ****     INSURE NON-NEGATIVE N(O+)
      NOPK=NJNP+NOP
      NOPNMK=NJNP+NOPNM
      DO 32 I=1,LEN2
      F(I,NOPK)=merge(F(I,NOPK),1.E-5,F(I,NOPK)-1.E-5>=0.)
      F(I,NOPNMK)=merge(F(I,NOPNMK),1.E-5,F(I,NOPNMK)-1.E-5>=0.)
   32 CONTINUE
      RETURN
      END
C
