#include "dims.h"
      SUBROUTINE COMP
      use input_module,only: difhor
      use cons_module,only: len1,len3,imax,imaxh,imaxp2,imaxp4,kmax,
     |  kmaxp1,difk,expz,t0,kut,rmass,rmassinv,shapiro,dtsmooth,
     |  dtsmooth_div2,expzmid_inv,expzmid,dtx2inv
      implicit none
C     ****     ADVANCE PSI IN TIME
#include "params.h"
#include "fcom.h"
#include "ccomp.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
#include "cmpbnd.h"
#include "crates_tdep.h"
!
! Local:
      real :: AK0(2,2),PHI(2,3),DELTA(2,2),tau,t00,small
      integer :: i,l,ll,nps1k,km,kp,npsl,m,ntk,k,kb,nwk,nrjk,npsk,
     |  npjp1k,npsnmk,npjm1k,npjm2k,npsdhk,kk,npsm,npsnpk,npsmnk,len,k1,
     |  k2,n,nps1mk,nps2mk,nps2k,npjp2k
      real,dimension(zimxp,zkmxp) :: o2smooth,o1smooth
      integer :: ii
      PHI(:,1)=(/0.,0.673/)
      PHI(:,2)=(/1.35,0./)
      PHI(:,3)=(/1.11,0.769/)
      TAU=1.86e+3
      DELTA(:,1)=(/1.,0./)
      DELTA(:,2)=(/0.,1./)
      T00=273.
!     SMALL = 1.e-20
      SMALL = 1.e-10 ! modsrc.eq
!     SMALL = 1.e-6  ! orig tgcm15
!
      DO 34 I=1,LEN1
      T1(I)=1.
   34 CONTINUE
      IF (difhor.EQ.1) CALL DFACT(T1,J,IMAXP4)

!     write(6,"('comp: lat=',i2,' dfactor=',/,(6e12.4))")
!    |  j,t1

C     ****     CALCULATE EMBAR
      NPS1K=NJ+NPS
      NPS2K=NPS1K+KMAXP1
      DO 1 I=1,LEN3
        EMBAR(I,1)=1./(F(I,NPS1K)*rmassinv(1)+F(I,NPS2K)*rmassinv(2)+
     |    (1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3))
    1 CONTINUE

!     call addfsech('EMBAR',' ',' ',embar,zimxp,zkmxp,zkmxp,j)

C     CALCULATE PS0 AND EMBAR0
      DO 30 I=1,LEN1
        PS0(I,1) = B(I,1,1)*F(I,NPS1K)+B(I,1,2)*F(I,NPS2K)+FB(I,1)
        PS0(I,2) = B(I,2,1)*F(I,NPS1K)+B(I,2,2)*F(I,NPS2K)+FB(I,2)
        EMBAR0(I)=1./(PS0(I,1)*rmassinv(1)+PS0(I,2)*rmassinv(2)+
     |    (1.-PS0(I,1)-PS0(I,2))*rmassinv(3))
C     ****     WKS4=.5*(DMBAR/DZ)/MBAR
        WKS4(I)=(EMBAR(I,1)-EMBAR0(I))/(dz*(EMBAR0(I)+EMBAR(I,1)))
   30 CONTINUE

!     write(6,"(/,'comp: lat=',i2,' ps0(:,1)=',/,(6e12.4))") 
!    |  j,ps0(:,1)
!     write(6,"('ps0(:,2)=',/,(6e12.4))") ps0(:,2)
!     write(6,"('embar0(:)=',/,(6e12.4))") embar0(:)

C     ****     CALCULTE EP AND AK AT LEVEL 1/2, SET Z(0) TO ZERO
      KM=1
      KP=2
C     ****     EP(1/2)
      DO 2 L=1,2
      DO 2 I=1,LEN1
      EP(I,L,KP)=1.-(2./(EMBAR0(I)+EMBAR(I,1)))*(RMASS(L)+(EMBAR(I,1)-
     AEMBAR0(I))/dz)
      Z(I,1,L)=0.
    2 CONTINUE

!     write(6,"('ep(:,1,kp)=',/,(6e12.4))") ep(:,1,kp)
!     write(6,"('ep(:,2,kp)=',/,(6e12.4))") ep(:,2,kp)

      DO 3 L=1,2
      LL=3-L
      NPSL=NJ+NPS+(L-1)*KMAXP1
      DO 3 M=1,2
      DO 3 I=1,LEN1
      AK(I,L,M,KP)=-DELTA(L,M)*(PHI(LL,3)+(PHI(LL,L)-PHI(LL,3))*.5*
     A(PS0(I,L)+F(I,NPSL)))-(1.-DELTA(L,M))*(PHI(L,M)-PHI(L,3))*
     B.5*(PS0(I,L)+F(I,NPSL))
    3 CONTINUE

!     do m=1,2
!       write(6,"('lat=',i2,' m=',i2,' ak(io2)=',/,(6e12.4))") 
!    |    j,m,ak(:,1,m,kp) 
!       write(6,"('lat=',i2,' m=',i2,' ak(io1)=',/,(6e12.4))")
!    |    j,m,ak(:,2,m,kp)
!     enddo

C     ****     WKS1=MBAR/M3*(T00/(T0+T))*0.25/(TAU*DET) AT LEVEL 1/2
      NTK=NJ+NT+KMAX
      DO 4 I=1,LEN1
        WKS1(I)=0.5*(EMBAR0(I)+EMBAR(I,1))*rmassinv(3)*(T00/(T0(1)+
     |    F(I,NTK)))**0.25/(TAU*(AK(I,1,1,KP)*AK(I,2,2,KP)-AK(I,1,2,KP)*
     |    AK(I,2,1,KP)))
    4 CONTINUE

!     write(6,"('comp: lat=',i3,' wks1=',/,(6e12.4))") j,wks1

C     ****     COMPLETE CALCULATION OF AK(1/2)
      DO 5 L=1,2
      DO 5 M=1,2
      DO 5 I=1,LEN1
      AK(I,L,M,KP)=AK(I,L,M,KP)*WKS1(I)
C     ****     SET GAMA(0) TO 0.
      GAMA(I,1,L,M)=0.
    5 CONTINUE

!     do m=1,2
!       write(6,"('lat=',i2,' m=',i2,' ak(io2)=',/,(6e12.4))") 
!    |    j,m,ak(:,1,m,kp) 
!       write(6,"('lat=',i2,' m=',i2,' ak(io1)=',/,(6e12.4))")
!    |    j,m,ak(:,2,m,kp)
!     enddo

C     ****     MAIN HEIGHT DO LOOP
      DO 6 K=1,KMAX
C     ****     CYCLE K LEVELS
      KB=KM
      KM=KP
      KP=KB
C     **** EP(K+1/2)
      DO 7 L=1,2
      DO 7 I=1,LEN1
      EP(I,L,KP)=1.-(2./(EMBAR(I,K)+EMBAR(I,K+1)))*(RMASS(L)+
     A(EMBAR(I,K+1)-EMBAR(I,K))/dz)
    7 CONTINUE

!     write(6,"('lat=',i3,' k=',i3,' ep(io2)=',/,(6e12.4))") j,k,
!    |  ep(:,1,kp)
!     write(6,"('ep(io1)=',/,(6e12.4))") ep(:,2,kp)

C     ****     AK(K+1/2)
      DO 8 L=1,2
      LL=3-L
      NPSL=NJ+NPS+(L-1)*KMAXP1+K
      DO 8 M=1,2
      DO 8 I=1,LEN1
      AK(I,L,M,KP)=-DELTA(L,M)*(PHI(LL,3)+(PHI(LL,L)-PHI(LL,3))*.5*
     A(F(I,NPSL-1)+F(I,NPSL)))-(1.-DELTA(L,M))*(PHI(L,M)-PHI(L,3))*
     B.5*(F(I,NPSL-1)+F(I,NPSL))
    8 CONTINUE
C     ****     WKS1=MBAR/M3*(T00/(T0+T))**0.25/(TAU*DET(ALFA))
      NTK=NJ+NT+K
      DO 9 I=1,LEN1
      WKS1(I)=0.5*(EMBAR(I,K)+EMBAR(I,K+1))*rmassinv(3)*(T00/(T0(K+1)+
     |  .5*(F(I,NTK-1)+F(I,NTK))))**0.25/(TAU*(AK(I,1,1,KP)*
     |  AK(I,2,2,KP)-AK(I,1,2,KP)*AK(I,2,1,KP)))
C     ****     EDDY DIFFUSION TERMS IN WKS3 AND WKS4
      WKS3(I)=WKS4(I)
      WKS4(I)=(EMBAR(I,K+1)-EMBAR(I,K))/(dz*(EMBAR(I,K)+EMBAR(I,K+1)))
    9 CONTINUE
C     ****     FINISH CALCULATING AK(K+1/2) AND GENERATE PK, QK, RK
      NWK=NJ+NW+K
      NRJK=NRJ+K
      DO 1000 M=1,2
      DO 1001 L=1,2
      DO 1002 I=1,LEN1

      AK(I,L,M,KP)=AK(I,L,M,KP)*WKS1(I)

      PK(I,L,M)=(AK(I,L,M,KM)*(1./dz+EP(I,M,KM)/2.)-expz(K)
     A*(expzmid_inv*DIFK(K)*T1(I)*(1./dz-WKS3(I))+0.25*(F(I,NWK-1)+
     BF(I,NWK)))*DELTA(L,M))/dz

      RK(I,L,M)=(AK(I,L,M,KP)*(1./dz-EP(I,M,KP)/2.)-expz(K)
     A*(expzmid*DIFK(K+1)*T1(I)*(1./dz+WKS4(I))-0.25*(F(I,NWK-1)+
     BF(I,NWK)))*DELTA(L,M))/dz

      QK(I,L,M)=-(AK(I,L,M,KM)*(1./dz-EP(I,M,KM)/2.)+AK(I,L,M,KP)*
     A(1./dz+EP(I,M,KP)/2.))/dz+expz(K)*(((expzmid*DIFK(K+1)*
     B(1./dz-WKS4(I))+expzmid_inv*DIFK(K)*(1./dz+WKS3(I)))*T1(I)/dz+
     Bdtx2inv)*DELTA(L,M)-FS(I,K,L,M))

 1002 CONTINUE
!      write(6,"(/,'comp: lat=',i3,' m=',i2,' isp=',i2,' k=',i3)") 
!    |   j,m,l,k
!      write(6,"('ak=',/,(6e12.4))") ak(:,l,m,km) 
!      write(6,"('ak=',/,(6e12.4))") ak(:,l,m,kp) 

!      write(6,"('ep=',/,(6e12.4))") ep(:,m,km) 
!      write(6,"('ep=',/,(6e12.4))") ep(:,m,kp) 

!      write(6,"('pk=',/,(6e12.4))") pk(:,l,m)
!      write(6,"('rk=',/,(6e12.4))") rk(:,l,m)
!      write(6,"('qk=',/,(6e12.4))") qk(:,l,m)

 1001 CONTINUE
 1000 CONTINUE
C     ****     CALCULATE HORIZONTAL ADVECTION IN FK ARRAY
      DO 11 L=1,2
      NPSL=NPS+(L-1)*KMAXP1+K-1
      ii = 0 
      if (l==1) ii = 1
      CALL ADVECL(NPSL,FK(1,L),K,ii)
   11 CONTINUE

!     write(6,"(/,'comp: lat=',i3,' k=',i2)") j,k
!     write(6,"('fk(:,io2)=',/,(6e12.4))") fk(3:imaxp2,1)
!     write(6,"('fk(:,io1)=',/,(6e12.4))") fk(3:imaxp2,2)

C     ****
C     ****     ADD EXPLICIT SOURCE TERMS  TO FK
C     ****
      DO 211 I = 1,LEN1
        FK(I,1) = FK(I,1)-FS(I,K,1,0)
        FK(I,2) = FK(I,2)-FS(I,K,2,0)
  211 CONTINUE

!     write(6,"(/,'comp: lat=',i3,' k=',i2)") j,k
!     write(6,"('fk(:,io2)=',/,(6e12.4))") fk(3:imaxp2,1)
!     write(6,"('fk(:,io1)=',/,(6e12.4))") fk(3:imaxp2,2)

C     ****     SHAPIRO SMOOTHER
      DO 111 L=1,2
      NPSK=NPSNM+(L-1)*KMAXP1+K-1
      NPJP2K=NJP2+NPSK
      NPJP1K=NJP1+NPSK
      NPSNMK=NJ+NPSK
      NPJM1K=NJM1+NPSK
      NPJM2K=NJM2+NPSK
      DO 111 I=1,LEN1
      WKV1(I,L)=F(I,NPSNMK)-shapiro*(F(I,NPJP2K)+F(I,NPJM2K)-4.*
     A(F(I,NPJP1K)+F(I,NPJM1K))+6.*F(I,NPSNMK))
  111 CONTINUE
      DO 112 L=1,2
      DO 112 I=3,IMAXP2
      WKV2(I,L)=WKV1(I,L)-shapiro*(WKV1(I+2,L)+WKV1(I-2,L)-4.*
     A(WKV1(I+1,L)+WKV1(I-1,L))+6.*WKV1(I,L))

        o2smooth(i,k) = wkv2(i,1)
        o1smooth(i,k) = wkv2(i,2)
  112 CONTINUE
C     ****     COMPLETE CALCULATION OF RHS IN FK
      DO 12 L=1,2
!     NPSDHK=NPSDH+(L-1)*KMAX+K-1
      NPSDHK=NPSDH+(L-1)*KMAXP1+K-1
      DO 12 I=3,IMAXP2
      FK(I,L)=expz(K)*(WKV2(I,L)*dtx2inv-FK(I,L)+F(I,NPSDHK))
   12 CONTINUE

!     write(6,"(/,'comp: lat=',i3,' k=',i2)") j,k
!     write(6,"('fk(:,io2)=',/,(6e12.4))") fk(3:imaxp2,1)
!     write(6,"('fk(:,io1)=',/,(6e12.4))") fk(3:imaxp2,2)

C     ****     INSERT PERIODIC POINTS
      DO 121 L=1,2
      DO 121 I=1,2
      FK(I,L)=FK(I+IMAX,L)
      FK(I+IMAXP2,L)=FK(I+2,L)
  121 CONTINUE
      IF(K.EQ.1) GO TO 13
      IF(K.EQ.KMAX) GO TO 14
      GO TO 15
   13 CONTINUE
C     ****     LOWER BNDRY
      DO 16 L=1,2
      DO 16 M=1,2
      DO 16 KK=1,2
      DO 16 I=1,LEN1

      QK(I,L,M) = QK(I,L,M)+PK(I,L,KK)*B(I,KK,M)

!     if (l==1) then
!       write(6,"('after: lat=',i2,' i=',i2,' kk=',i2,' m=',i2,
!    |    ' qk(lev0)=',e12.4)") j,i,kk,m,qk(i,l,m)
!     endif

   16 CONTINUE
      DO 33 L=1,2
      DO 33 M=1,2
      DO 33 I=1,LEN1
      FK(I,L) = FK(I,L)-PK(I,L,M)*FB(I,M)
      PK(I,L,M)=0.
   33 CONTINUE

!     do m=1,2 
!       write(6,"('comp lbc: m=',i2,' lat=',i2)") m,j
!       write(6,"('qk(io2)=',/,(6e12.4))") qk(:,1,m)
!       write(6,"('qk(io1)=',/,(6e12.4))") qk(:,2,m)
!     enddo ! m=1,2
!     write(6,"('fk(io2)=',/,(6e12.4))") fk(:,1)
!     write(6,"('fk(io1)=',/,(6e12.4))") fk(:,2)

      GO TO 15
   14 CONTINUE
C     ****     UPPER BNDRY
      DO 17 L=1,2
      DO 17 M=1,2
      DO 17 I=1,LEN1

!     if (l==1) 
!    |  write(6,"('comp: i=',i2,' m=',i2,' j=',i2,' ep=',e12.4,' rk=',
!    |    e12.4,' qk=',e12.4)") i,m,j,ep(i,m,kp),
!    |    rk(i,l,m),qk(i,l,m)

      QK(I,L,M)=QK(I,L,M)+(1.+.5*EP(I,M,KP)*dz)/(1.-.5*EP(I,M,KP)*
     Adz)*RK(I,L,M)

      RK(I,L,M)=0.
   17 CONTINUE

!     do m=1,2
!       write(6,"('comp ubc: m=',i2,' lat=',i2)") m,j
!       write(6,"('qk(io2)=',/,(6e12.4))") qk(:,1,m)
!       write(6,"('qk(io1)=',/,(6e12.4))") qk(:,2,m)
!     enddo ! m=1,2

   15 CONTINUE

C     ****     QK=ALFAK=QK-PK*GAMA(K-1)
      DO 1800 L=1,2
      DO 1801 M=1,2
      DO 1802 KK=1,2
      DO 1803 I=1,LEN1

!     if (l==1) then
!       write(6,"('comp: i=',i2,' kk=',i2,' m=',i2,' k=',i2,
!    |    ' lat=',i2)") i,kk,m,k,j
!       write(6,"('qk=',e12.4,' pk=',e12.4,' gama=',e12.4))") 
!    |    qk(i,l,m),pk(i,l,kk),gama(i,k,kk,m)
!     endif

      QK(I,L,M)=QK(I,L,M)-PK(I,L,KK)*GAMA(I,K,KK,M)
 1803 CONTINUE
 1802 CONTINUE
 1801 continue
 1800 continue

!     do m=1,2
!       write(6,"('comp: m=',i2,' k=',i2,' lat=',i2)") m,k,j
!       write(6,"('qk(io2)=',/,(6e12.4))") qk(:,1,m)
!       write(6,"('qk(io1)=',/,(6e12.4))") qk(:,2,m)
!     enddo


!     do m=1,2
!       write(6,"('comp: m=',i2,' k=',i2,' lat=',i2,
!    |    ' qk(:,k,io2,m,lat)=',/,(6e12.4))") m,k,j,qk(:,1,m)
!     enddo

C     ****     WKS1=DET(ALFA)
      DO 19 I=1,LEN1
      WKS1(I)=QK(I,1,1)*QK(I,2,2)-QK(I,1,2)*QK(I,2,1)
   19 CONTINUE
C     ****     WKM1=ALFAI
      DO 20 L=1,2
      LL=3-L
      DO 20 M=1,2
      DO 20 I=1,LEN1
      WKM1(I,L,M)=(DELTA(L,M)*QK(I,LL,LL)-(1.-DELTA(L,M))*QK(I,L,M))/
     AWKS1(I)
   20 CONTINUE
C     ****     GAMA(K+1)=ALFAI*RK
C     ****     WKV1=FK-PK*Z(K)
      DO 21 L=1,2
      DO 22 I=1,LEN1
      WKV1(I,L)=FK(I,L)
   22 CONTINUE
      DO 21 M=1,2
      DO 23 I=1,LEN1
      GAMA(I,K+1,L,M)=0.
      WKV1(I,L)=WKV1(I,L)-PK(I,L,M)*Z(I,K,M)
   23 CONTINUE
      DO 21 KK=1,2
      DO 21 I=1,LEN1
      GAMA(I,K+1,L,M)=GAMA(I,K+1,L,M)+WKM1(I,L,KK)*RK(I,KK,M)
   21 CONTINUE
C     ****     Z(K+1)=WKM1*WKV1
      DO 231 L=1,2
      DO 232 I=1,LEN1
      Z(I,K+1,L)=0.
  232 CONTINUE
      DO 231 M=1,2
      DO 231 I=1,LEN1
      Z(I,K+1,L)=Z(I,K+1,L)+WKM1(I,L,M)*WKV1(I,M)
  231 CONTINUE

!     do m=1,2 
!       write(6,"('comp: m=',i2,' k=',i2,' lat=',i2)") m,k,j 
!       write(6,"('gama(k+1,io2)=',/,(6e12.4))") gama(:,k+1,1,m) 
!       write(6,"('gama(k+1,io1)=',/,(6e12.4))") gama(:,k+1,2,m)   
!     enddo 
!     write(6,"('zz(k+1,io2)=',/,(6e12.4))") z(:,k+1,1)
!     write(6,"('zz(k+1,io1)=',/,(6e12.4))") z(:,k+1,2)

    6 CONTINUE

!     call addfsech('O2SMOOTH',' ',' ',o2smooth,zimxp,zkmxp,zkmxp,j)
!     call addfsech('O1SMOOTH',' ',' ',o1smooth,zimxp,zkmxp,zkmxp,j)
!     call addfsech('ZZ_O2'   ,' ',' ',z(:,:,1),zimxp,zkmxp,zkmxp,j)
!     call addfsech('ZZ_O1'   ,' ',' ',z(:,:,2),zimxp,zkmxp,zkmxp,j)
!     call addfsech('GAMAO2M1',' ',' ',gama(:,:,1,1),zimxp,zkmxp,zkmxp,
!    |  j)
!     call addfsech('GAMAO2M2',' ',' ',gama(:,:,1,2),zimxp,zkmxp,zkmxp,
!    |  j)
!     call addfsech('GAMAO1M1',' ',' ',gama(:,:,2,1),zimxp,zkmxp,zkmxp,
!    |  j)
!     call addfsech('GAMAO1M2',' ',' ',gama(:,:,2,2),zimxp,zkmxp,zkmxp,
!    |  j)

C     ****     SET PSNP(KMAXP1) TO ZERO
      DO 24 L=1,2
      NPSL=NJNP+NPS+(L-1)*KMAXP1+KMAX
      DO 24 I=1,LEN1
      F(I,NPSL)=0.
   24 CONTINUE
C     ****     DOWNWARD SWEEP
      DO 25 KK=1,KMAX
      K=KMAX+1-KK
      DO 25 L=1,2
      NPSL=NJNP+NPS+(L-1)*KMAXP1+(K-1)
      DO 26 I=1,LEN1
      F(I,NPSL)=Z(I,K+1,L)
   26 CONTINUE
      DO 25 M=1,2
      NPSM=NJNP+NPS+(M-1)*KMAXP1+K
      DO 25 I=1,LEN1
      F(I,NPSL)=F(I,NPSL)-GAMA(I,K+1,L,M)*F(I,NPSM)
   25 CONTINUE
C     ****     INSERT VALUE OF PS(KMAXP1) USING BNDRY CONDITION
      DO 27 L=1,2
      NPSL=NJNP+NPS+(L-1)*KMAXP1+KMAX
      DO 27 I=1,LEN1
      F(I,NPSL)=(1.+.5*EP(I,L,KP)*dz)/(1.-.5*EP(I,L,KP)*dz)*
     AF(I,NPSL-1)
   27 CONTINUE

!     call addfsech('O2_SOLV',' ',' ',f(1,njnp+nps),
!    |  zimxp,zkmxp,zkmxp,j)
!     call addfsech('O1_SOLV',' ',' ',f(1,njnp+nps2),
!    |  zimxp,zkmxp,zkmxp,j)

C     ****     FILTER PSI
      IF(KUT(J).GE.IMAXH) GO TO 28
      NPSK=NJNP+NPS
      CALL FILTER(NPSK,2*KMAXP1,KUT(J))
   28 CONTINUE

!     call addfsech('O2_FILT',' ',' ',f(1,njnp+nps),
!    |  zimxp,zkmxp,zkmxp,j)
!     call addfsech('O1_FILT',' ',' ',f(1,njnp+nps2),
!    |  zimxp,zkmxp,zkmxp,j)

C     ****     TIME SMOOTHING OF PSI
      NPSNMK=NJ+NPSNM
      NPSK=NJ+NPS
      NPSNPK=NJNP+NPS
      NPSMNK=NJNP+NPSNM
      LEN=2.*LEN3
      DO 29 I=1,LEN
      F(I,NPSMNK)=dtsmooth*F(I,NPSK)+dtsmooth_div2*(F(I,NPSNMK)+
     |  F(I,NPSNPK))
   29 CONTINUE

!     call addfsech('O2NM_OUT',' ',' ',f(1,njnp+npsnm),
!    |  zimxp,zkmxp,zkmxp,j)
!     call addfsech('O1NM_OUT',' ',' ',f(1,njnp+nps2nm),
!    |  zimxp,zkmxp,zkmxp,j)

C     ****     INSERT PERIODIC POINTS
      K1=NJNP+NPS
      K2=K1+2*KMAXP1-1
      DO 31 N=1,2
      DO 32 I=1,2
      DO 32 K=K1,K2
      F(I,K)=F(I+IMAX,K)
      F(I+IMAXP2,K)=F(I+2,K)
   32 CONTINUE
      K1=NJNP+NPSNM
      K2=K1+2*KMAXP1-1
   31 CONTINUE
C     ****     INSURE NON-NEGATIVE PSI
      NPS1K=NJNP+NPS
      NPS2K=NJNP+NPS2
      NPS1MK=NJNP+NPSNM
      NPS2MK=NJNP+NPS2NM
      DO 41 I=1,LEN3
        if (f(i,nps1k) < small) f(i,nps1k) = small
        if (f(i,nps2k) < small) f(i,nps2k) = small
        if (f(i,nps1mk) < small) f(i,nps1mk) = small
        if (f(i,nps2mk) < small) f(i,nps2mk) = small
        if (1.-SMALL-F(I,NPS1K)-F(I,NPS2K) < 0.) then
          f(i,nps1k) = f(i,nps1k)*((1.-SMALL)/(F(I,NPS1K)+F(I,NPS2K)))
          f(i,nps2k) = f(i,nps2k)*((1.-SMALL)/(F(I,NPS1K)+F(I,NPS2K)))
        endif
        if (1.-SMALL-F(I,NPS1MK)-F(I,NPS2MK) < 0.) then
          f(i,nps1mk) = f(i,nps1mk)*((1.-SMALL)/(F(I,NPS1MK)+
     |                  F(I,NPS2MK)))
          f(i,nps2mk) = f(i,nps2mk)*((1.-SMALL)/(F(I,NPS1MK)+
     |                  F(I,NPS2MK)))
        endif
   41 CONTINUE

!     call addfsech('O2_FINAL',' ',' ',f(1,njnp+nps),
!    |  zimxp,zkmxp,zkmxp,j)
!     call addfsech('O1_FINAL',' ',' ',f(1,njnp+nps2),
!    |  zimxp,zkmxp,zkmxp,j)
!     call addfsech('O2MFINAL',' ',' ',f(1,njnp+npsnm),
!    |  zimxp,zkmxp,zkmxp,j)
!     call addfsech('O1MFINAL',' ',' ',f(1,njnp+nps2nm),
!    |  zimxp,zkmxp,zkmxp,j)

      RETURN
      END
C