C++*********************************************************************
C
C  REPCG.F         USED PROJT FOR BCKCQ CALL      FEB 2000 ARDEAN LEITH
C                  OPFILEC                        FEB   03 ARDEAN LEITH
C                  VERBOSE                        FEB   06 ARDEAN LEITH
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
C=*                                                                    *
C=* This program is free software; you can redistribute it and/or      *
C=* modify it under the terms of the GNU General Public License as     *
C=* published by the Free Software Foundation; either version 2 of the *
C=* License, or (at your option) any later version.                    *
C=*                                                                    *
C=* This program is distributed in the hope that it will be useful,    *
C=* but WITHOUT ANY WARRANTY; without even the implied warranty of     *
C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C
C  REPROJECTIONS 3D, CONJUGATE GRADIENTS METHOD, 
C  RECONSTRUCTION KEPT IN THE SQUARE TO INTRODUCE OTHER CONSTRAINTS.
C  AVERAGE OUTSIDE THE WINDOW IS SUBTRACTED
C  SYMMETRIES IMPOSED ... - NOT YET!  INTERPOLATION CHANGED !
C  CORRECTIONS APPLIED ON THE VOLUME SIDE.      11/09/98
C  COMPRESSION OF ANGLES - 08/14/96
C
C  REPCG(MAXMEM)
C  REDPRCG
C  RPRQ(N,B,CUBE,IPCUBE,NN,PHI,THETA,PSI,DM,RI)
C  ASTASQ(X,N,RI,ABA,KLP)
C  PREPCUB_S(N,NN,IPCUBE,RI)
C  BCKCQ(CUBE,LTC,DM,B,N,IPCUBE,NN)
C  PRJCQ(CUBE,LTC,DM,B,N,IPCUBE,NN)
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C--*********************************************************************

        SUBROUTINE REPCG

        INCLUDE 'CMBLOCK.INC'
        INCLUDE 'CMLIMIT.INC'
        INCLUDE 'F90ALLOC.INC'

        REAL, DIMENSION(:,:), POINTER ::        PANG,PSYM
        REAL, ALLOCATABLE, DIMENSION(:,:,:) ::  CB
        REAL, ALLOCATABLE, DIMENSION(:,:) ::    ANG,DM
        REAL, ALLOCATABLE, DIMENSION(:) ::      BCKN

        INTEGER, ALLOCATABLE, DIMENSION(:,:) :: QM
        INTEGER, ALLOCATABLE, DIMENSION(:) ::   ILIST,LB

        LOGICAL           MD 
        COMMON  /F_SPEC/  FINPAT,NLET,FINPIC
        CHARACTER*80      FINPIC,FINFO,FINPAT
        CHARACTER(LEN=MAXNAM)   ::  ANGDOC

        CHARACTER*1       NULL,ANS
        COMMON /PAR/      LDP,NM,LDPNM
        DATA  INPIC/99/

#ifdef USE_MPI
        include 'mpif.h'
        INTEGER  MYPID, COMM, MPIERR, J
C
        COMM = MPI_COMM_WORLD
        CALL MPI_COMM_RANK(COMM, MYPID, MPIERR)
#else
        MYPID = -1
#endif
        NULL = CHAR(0)
        NILMAX = NIMAX
     
C       N - LINEAR DIMENSION OF PROJECTIONS AND RESTORED CUBE
C       NANG - NUMBER OF ANGLES (PROJECTIONS)

        ALLOCATE (ILIST(NILMAX), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'REPCG, ILIST',IER)
           RETURN
        ENDIF

        CALL FILELIST(.TRUE.,INPIC,FINPAT,NLET,ILIST,NILMAX,NANG,
     &                 'ENTER TEMPLATE FOR 2-D IMAGES',IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           CALL ERRT(27,'REPCG ',NE)
           DEALLOCATE(ILIST)          
           RETURN
        ENDIF  
        
        MAXNUM = MAXVAL(ILIST(1:NANG))

        IF (MYPID .LE. 0) WRITE(NOUT,2001) NANG
2001    FORMAT('  NUMBER OF IMAGES =',I5)
C       NANG - TOTAL NUMBER OF IMAGES

        CALL  RDPRM(RI,NOT_USED,'RADIUS OF RECONSTRUCTED OBJECT')

C       RETRIEVE ARRAY WITH ANGLES DATA IN IT
        MAXXT = 4
        MAXYT = MAXNUM
        CALL GETDOCDAT('ANGLES DOC',.TRUE.,ANGDOC,77,.FALSE.,MAXXT,
     &                       MAXYT,PANG,IRTFLG)
 
        IF (IRTFLG .NE. 0) THEN
           CALL ERRT(4,'REPCG ',NE)
           DEALLOCATE (ILIST)
           RETURN
        ENDIF       

C       OPEN FIRST IMAGE FILE TO DETERMINE NSAM, NROW, NSL
        CALL FILGET(FINPAT,FINPIC,NLET,ILIST(1),IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           CALL ERRT(4,'REPCG ',NE)
           DEALLOCATE (ILIST)
           RETURN
        ENDIF

        MAXIM = 0
        CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM,NSAM,NROW,NSL,
     &               MAXIM,'DUMMY',.FALSE.,IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           DEALLOCATE (ILIST)
           RETURN
        ENDIF

        N=NSAM
        M=N

        NM=(N-M)/2
        LDP=M/2+1
        LDPNM=LDP+NM


        CALL RDPRMC (ANS, NLETI, .TRUE.,
     &     'DOES YOUR VOLUME HAVE SYMMETRIES? (Y/N)',NULL,IRT)
        IF (ANS.EQ.'Y')  THEN
           WRITE(NOUT,*)  '*** OPTION NOT IMPLEMENTED YET'
           DEALLOCATE(PANG)
           DEALLOCATE(ILIST)
           CLOSE(INPIC) 
           RETURN
        ENDIF

       NSYM = 0
       IF (ANS .EQ. 'Y')  THEN
           MAXXT = 0
           MAXYT = 0
           CALL GETDOCDAT('SYMMETRIES ANGLES DOC',.TRUE.,FINPIC,INPIC,
     &         .TRUE.,MAXXT,MAXYT,PSYM,IRTFLG)
           IF (IRTFLG.NE.0) THEN
              CALL ERRT(4,'REPCG ',NE)
              GOTO 9120
           ENDIF
           NSYM = MAXYT
           
           WRITE(NOUT,2021) NSYM
2021       FORMAT(/,' NUMBER OF SYMMETRIES:',I5,/)
        ELSE 
C          PSYM IS NOW A DUMMY POINTER
           PSYM => PANG
        ENDIF

C       DUM IS A DUMMY VARIABLE, VALUE OF NN IS DETERMINED
        
        MD = .FALSE.
        CALL  PREPCUB_S(N,NN,DUM,RI,MD)

C       USE NN TO ALLOCATE (QM) 
        ALLOCATE (QM(5,NN), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'REPS, QM',IER)
           GOTO 9118
        ENDIF

        MD = .TRUE.
        CALL  PREPCUB_S(N,NN,QM,RI,MD)
c       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        NMAT=N*N*N
        LSD=N+2-MOD(N,2)
C       IN THIS VERSION TOTAL MEMORY IS THREE VOLUMES BCKN, BCKE AND CUBE
C       PLUS 2 2D PROJECTIONS.
C
C       CUBE - KEEPS BACK-PROJECTED ORIGINAL PROJECTIONS, READ FROM THE DISK
C       BCKE - WORKING VOLUME
C       BCKN - CURRENT RECONSTRUCTION

C       FIND NUMBER OF OMP THREADS
        CALL GETTHREADS(NUMTH)

        MAXIM = 0
        IFORM = 3
        CALL OPFILEC(0,.TRUE.,FINPIC,INPIC,'U',IFORM,M,M,M,
     &             MAXIM,'RECONSTRUCTED 3-D',.FALSE.,IRTFLG)
        IF (IRTFLG .NE. 0) GOTO 9116

        ALLOCATE (CB(N,N,N), ANG(3,NANG),DM(9,NANG),STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'REPCG, CB, ANG & DM',IER)
           GOTO 9116
        ENDIF

        CALL REDPRCG(N,NANG,CB,ANG,ILIST,QM,NN,DM,RI,PANG)

C       COMPRESS ANGLES - IT CHANGES NANG !!

        ALLOCATE (LB(NANG), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'REPCG, LB',IER)
           GOTO 9110
        ENDIF

        CALL  HIANG(ANG,NANG,DM,LB,LO)
        NANG=LO
        IF (MYPID .LE. 0) WRITE(NOUT,2027) NANG
2027    FORMAT('  EFFECTIVE NUMBER OF ANGLES: ',I6)

        NSYMT=MAX0(1,NSYM)

        ALLOCATE (BCKN(N*N*N), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'REPCG, BCKN',IER)
           GOTO 9108
        ENDIF

        CALL  REPRCG3_Q(CB,BCKN,LSD,N,QM,NN,DM,LB,NANG,RI,
     &                  PSYM,NSYMT,NSYM,NUMTH)

        DO  K=1,M*M
           MI=(K-1)*M
           CALL WRTLIN(INPIC,BCKN(1+MI),M,K)
        ENDDO

        DEALLOCATE(BCKN)
9108    DEALLOCATE(LB)
9110    DEALLOCATE(DM)
9112    DEALLOCATE(ANG)
9114    DEALLOCATE(CB)
9116    DEALLOCATE(QM)
9118    IF (ANS.EQ.'Y') DEALLOCATE(PSYM) 
9120    DEALLOCATE(PANG)
        DEALLOCATE(ILIST)                
        CLOSE(INPIC)   
        END

C++*********************************************************************
C--*********************************************************************

        SUBROUTINE  REPRCG3_Q(CB,BCKN,LSD,N,IPCUBE,NN,DM,LB,NANG,RI,
     &          SYM,NSYMT,NSYM,NUMTH)

C       NUMTH - NUMTHREDS FOR MP; =NUMTHREADS() FOR MP, OTHERWISE=1.
        INCLUDE 'CMBLOCK.INC'

        DIMENSION   CB(N,N,N),BCKN(N,N,N)
        DIMENSION   SYM(3,NSYMT),DM(9,NANG),PROJ(N*N,NUMTH)
        DIMENSION   projt(4,n*n)
        DIMENSION   IPCUBE(5,NN),LB(NANG)
        REAL, ALLOCATABLE, DIMENSION(:,:,:) :: BCKE,BCKP

        COMMON /NORMB/ BNORM,CHI2
        COMMON /PAR/  LDP,NM,LDPNM

#ifdef USE_MPI
        include 'mpif.h'
        INTEGER MYPID, COMM, NPROCS, MPIERR, MPISTAT(MPI_STATUS_SIZE)
        INTEGER NANGLOC, IP, DMTAG, LBTAG, MASTER, N2, IZ, JZ
        REAL     ALA, SQ, DELSQ
        INTEGER, ALLOCATABLE, DIMENSION(:)    :: PSIZE, NBASE
        INTEGER, ALLOCATABLE, DIMENSION(:)    :: LBLOC
        REAL,    ALLOCATABLE, DIMENSION(:,:)  :: DMLOC
        REAL,    ALLOCATABLE, DIMENSION(:,:,:):: BCKE_SUM
        REAL,    ALLOCATABLE, DIMENSION(:,:,:):: RESID
C
        COMM = MPI_COMM_WORLD
        CALL MPI_COMM_RANK(COMM, MYPID , MPIERR)
        CALL MPI_COMM_SIZE(COMM, NPROCS, MPIERR)
C
        ALLOCATE(PSIZE(NPROCS), STAT=IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           WRITE(0,*) 'REPRCG3_Q: FAILED TO ALLOCATE PSIZE.'
           STOP
        END IF
        ALLOCATE(NBASE(NPROCS), STAT=IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           WRITE(0,*) 'REPRCG3_Q: FAILED TO ALLOCATE NBASE.'
           STOP
        END IF
        CALL SETPART(NANG, PSIZE, NBASE)
        NANGLOC = PSIZE(MYPID+1)
#endif
        
        CALL  RDPRM2(ERRM,CHIM,NOT_USED,'ERROR LIMIT, CHI^2 LIMIT')
        CALL  RDPRMI(MAXIT,MODE,NOT_USED,'ITERATION LIMIT, MODE')
        CALL  RDPRM(ALA,NOT_USED,'LAMBDA')
        DI=N/2
        IF (MODE.EQ.1.AND.RI+1.0.GT.DI
     &          .OR. MODE.EQ.2.AND.RI+2.0.GT.DI
     &          .OR. MODE.EQ.3.AND.RI+3.0.GT.DI)  THEN
           CALL ERRT(101,'RADIUS TOO LARGE',NE)
           RETURN
        ENDIF   

        ALLOCATE (BCKE(N,N,N), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'REPRCG3_Q, BCKE',IER)
           RETURN
        ENDIF

        ALLOCATE (BCKP(N,N,N), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'REPRCG3_Q, BCKP',IER)
           DEALLOCATE (BCKE)
           RETURN
        ENDIF

#ifdef USE_MPI
        ALLOCATE (BCKE_SUM(N,N,N), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN
           CALL ERRT(46,'REPRCG3_Q, BCKE_SUM',IER)
           RETURN
        ENDIF
#endif

        NMAT=N*N*N
        LTB=N*N
c$omp   parallel do private(k,j,i)
        DO K=1,N
           DO J=1,N
              DO I=1,N
C                X=0
                 BCKN(I,J,K)=0.0
C                CDEL=0
                 BCKP(I,J,K)=0.0
#ifdef USE_MPI
                 BCKE_sum(i,j,k) = 0.0
#endif
              ENDDO
           ENDDO
        ENDDO

        ERR=1.0
        ITER=0
        Q=0.0

#ifdef USE_MPI

C       DISTRIBUTE LB AND DM

        ALLOCATE(LBLOC(NANGLOC), STAT = IRTFLAG)
        IF (IRTFLG.NE.0) THEN
           CALL ERRT(46,'REPRCG3_Q, LBLOC',IER)
           RETURN
        ENDIF
        ALLOCATE(DMLOC(9, NANGLOC), STAT = IRTFLAG)
        IF (IRTFLG.NE.0) THEN
           CALL ERRT(46,'REPRCG3_Q, DMLOC',IER)
           RETURN
        ENDIF

        DMTAG  = 1
        LBTAG  = 2
        MASTER = 0
        IF (MYPID .EQ. 0) THEN

C          PROCESS 0 DISTRIBUTES DM() AND LB() TO OTHER PROCESSORS

           DO IP = 2, NPROCS
              IBEGIN = NBASE(IP) + 1
              CALL MPI_SEND(DM(1,IBEGIN), 9*PSIZE(IP), MPI_REAL,
     &                      IP-1        , DMTAG      , COMM    ,
     &                      IERR)
              CALL MPI_SEND(LB(IBEGIN)  , PSIZE(IP)  , MPI_INTEGER,
     &                      IP-1        , LBTAG      , COMM       ,
     &                      IERR)
           END DO
           CALL SCOPY(9*NANGLOC, DM, 1, DMLOC, 1)
           CALL SCOPY(NANGLOC  , LB, 1, LBLOC, 1)
        ELSE

C          === SLAVES RECEIVE DATA FROM THE MASTER ===

           CALL MPI_RECV(DMLOC, 9*NANGLOC, MPI_REAL, MASTER,
     &                   DMTAG, COMM     , MPISTAT , IERR)
           CALL MPI_RECV(LBLOC, NANGLOC  , MPI_INTEGER, MASTER,
     &                   LBTAG, COMM     , MPISTAT , IERR)
        END IF
#ifdef MPI_DEBUG
        WRITE(6,*) 'REPRCG3_Q: DATA DISTRIBUTION COMPLETED, MYPID = ',
     &              MYPID
#endif

        DO  ITER=1,MAXIT
           DELSQ=0.0
c$omp      parallel do private(kn,j,k,i),reduction(+:DELSQ)
           DO KN=1,NN
              J=IPCUBE(4,KN)
              K=IPCUBE(5,KN)
              DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 DELSQ=DELSQ+CB(I,J,K)*CB(I,J,K)
              ENDDO
           ENDDO
           Q=Q*DELSQ

c$omp      parallel do private(kn,j,k,i)
           DO    KN=1,NN
              J=IPCUBE(4,KN)
              K=IPCUBE(5,KN)
              DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 BCKP(I,J,K)=CB(I,J,K)-Q*BCKP(I,J,K)
              ENDDO
           ENDDO
           Q=-1.0/DELSQ

c$omp      parallel do private(k,j,i)
           DO    K=1,N
              DO    J=1,N
                 DO    I=1,N
C                   TEMP=0
                    BCKE(I,J,K)=0.0
                    bcke_sum(I,J,K) = 0.0
                 ENDDO
              ENDDO
           ENDDO

C          BCKP -> PROJ -> BCKE

C          LOOP OVER PROJECTIONS

           DO K=1,nangloc
              L_EN=MIN0(NANGloc,K+NUMTH-1)
              L_NUM=MIN0(NUMTH,NANGloc-K+1)
c$omp         parallel do private(i,j)
              DO J=1,L_NUM
                 DO I=1,LTB
                    PROJ(I,J)=0.0
                 ENDDO
              ENDDO

              CALL PRJCQ(BCKP,NMAT,dmloc(1,K),
     &                   PROJ,N,IPCUBE,NN)

C             HERE BCKCQ ITSELF IS MP  
              DO  L_TH=1,L_NUM
C                MULTIPLY PROJECTIONS BY THEIR WEIGHTS
                 IF (lbloc(K+L_TH-1).GT.1)  THEN
c$omp               parallel do private(i)
                    DO  I=1,N*N
                       PROJ(I,L_TH)=PROJ(I,L_TH)*lbloc(K+L_TH-1)
                    ENDDO
                 ENDIF

c$omp            parallel do private(i,pt)
                 do i = 1,n*n - n - 1
                    pt         = proj(i, L_TH)
                    projt(1,i) = pt
                    projt(2,i) = proj(i+n,  L_TH) - pt
                    projt(3,i) = proj(i+1,  L_TH) - pt 
                    projt(4,i) = proj(i+n+1,L_TH) - proj(i+1,L_TH) - 
     &                           projt(2,i)
                 enddo

C                BACKPROJECT
                 CALL  BCKCQ(BCKE,NMAT,DMloc(1,K+L_TH-1),
     &                       projt,N,IPCUBE,NN)
              ENDDO
           ENDDO
C
           CALL MPI_ALLREDUCE(BCKE   , BCKE_SUM,
     &                        NMAT   , MPI_REAL,
     &                        MPI_SUM, COMM    ,
     &                        IERR)
           IF (IERR .NE. 0) THEN
              WRITE(0,*) 'REPRCG3_Q: FAILED AT ALLREDUCE BCKE_SUM'
              STOP
           END IF

C          === IGNORE MODE > 0 FOR THE TIME BEING ===

           IF (MODE .EQ. 1)  THEN
              CALL  FIXEDGE1(BCKP,NMAT,BCKP,N,IPCUBE,NN)      
c              CALL  BFIRSTS(BCKE,BCKP,N,ALA,IPCUBE,NN,RI) sep 01 al
              CALL  BFIRSTS(BCKE_sum,BCKP,N,ALA,IPCUBE,NN)

           ELSEIF(MODE.EQ.2)  THEN
              CALL  FIXEDGE2(BCKP,NMAT,BCKP,N,IPCUBE,NN)      
              CALL  BSECOND(BCKE_sum,BCKP,N,ALA,IPCUBE,NN)

           ELSEIF(MODE.EQ.3)  THEN
              CALL  FIXEDGE3(BCKP,NMAT,BCKP,N,IPCUBE,NN)      
              CALL  BTHIRD(BCKE_sum,BCKP,N,ALA,IPCUBE,NN)
           ENDIF

           AKDEN=0.0

c$omp      parallel do private(kn,j,k,i),reduction(+:AKDEN)
           DO KN=1,NN
              J=IPCUBE(4,KN)
              K=IPCUBE(5,KN)
              DO I=IPCUBE(3,KN),
     &              IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 AKDEN=AKDEN+BCKP(I,J,K)*BCKE_sum(I,J,K)
              ENDDO
           ENDDO

           P=DELSQ/AKDEN

c$omp      parallel do private(kn,j,k,i)
           DO KN=1,NN
              J=IPCUBE(4,KN)
              K=IPCUBE(5,KN)
              DO I=IPCUBE(3,KN),
     &           IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 BCKN(I,J,K) = BCKN(I,J,K)+P*BCKP(I,J,K)
                 CB(I,J,K)   = CB(I,J,K)-P*BCKE_sum(I,J,K)
              ENDDO
           ENDDO
 
           ERR  = DELSQ/BNORM
           CHI2 = CHI2-P*DELSQ
           IF (MYPID .EQ. 0) WRITE(NOUT,2041) ITER,ERR,CHI2
2041       FORMAT('  ITERATION:', I3,'  ERROR:',1PE12.4,
     &            '  CHISQ:',1PE12.4)

C          STOPPING CRITERIA
           IF (ABS(ERR).LE.ERRM .OR. ABS(CHI2).LE.CHIM) THEN  
              DEALLOCATE (BCKE)
              DEALLOCATE (BCKP)
              DEALLOCATE (BCKE_SUM)
              RETURN
           ENDIF
        ENDDO
#else
        DO  ITER=1,MAXIT
           DELSQ=0.0
c$omp      parallel do private(kn,j,k,i),reduction(+:DELSQ)
           DO KN=1,NN
              J=IPCUBE(4,KN)
              K=IPCUBE(5,KN)
              DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 DELSQ=DELSQ+CB(I,J,K)*CB(I,J,K)
              ENDDO
           ENDDO
           Q=Q*DELSQ

c$omp      parallel do private(kn,j,k,i)
           DO    KN=1,NN
              J=IPCUBE(4,KN)
              K=IPCUBE(5,KN)
              DO I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 BCKP(I,J,K)=CB(I,J,K)-Q*BCKP(I,J,K)
              ENDDO
           ENDDO
           Q=-1.0/DELSQ

c$omp      parallel do private(k,j,i)
           DO K=1,N
              DO J=1,N
                 DO I=1,N
                    BCKE(I,J,K)=0.0
                 ENDDO
              ENDDO
           ENDDO

C          BCKP -> PROJ -> BCKE

C          LOOP OVER PROJECTIONS
           DO K=1,NANG,NUMTH
              L_EN=MIN0(NANG,K+NUMTH-1)
              L_NUM=MIN0(NUMTH,NANG-K+1)

c$omp         parallel do private(i,j)
              DO J=1,L_NUM
                 DO I=1,LTB
                    PROJ(I,J)=0.0
                 ENDDO
              ENDDO

c$omp       parallel do private(l_th),shared(nmat,n,nn),schedule(static)
              DO  L_TH=1,L_NUM
                 CALL  PRJCQ(BCKP,NMAT,DM(1,K+L_TH-1),
     &                       PROJ(1,L_TH),N,IPCUBE,NN)
              ENDDO

C             HERE BCKCQ ITSELF IS MP  
              DO  L_TH=1,L_NUM
C                MULTIPLY PROJECTIONS BY THEIR WEIGHTS
                 IF (LB(K+L_TH-1).GT.1)  THEN
c$omp               parallel do private(i)
                    DO  I=1,N*N
                       PROJ(I,L_TH)=PROJ(I,L_TH)*LB(K+L_TH-1)
                    ENDDO
                 ENDIF

c$omp            parallel do private(i,pt)
                 DO I = 1,N*N - N - 1
                    PT         = PROJ(I, L_TH)
                    PROJT(1,I) = PT
                    PROJT(2,I) = PROJ(I+N,  L_TH) - PT
                    PROJT(3,I) = PROJ(I+1,  L_TH) - PT 
                    PROJT(4,I) = PROJ(I+N+1,L_TH) - PROJ(I+1,L_TH) - 
     &                           PROJT(2,I)
                 ENDDO

C                BACKPROJECT
                 CALL  BCKCQ(BCKE,NMAT,DM(1,K+L_TH-1),
     &                       PROJT,N,IPCUBE,NN)
              ENDDO
           ENDDO

           IF (MODE .EQ. 1)  THEN
              CALL  FIXEDGE1(BCKP,NMAT,BCKP,N,IPCUBE,NN)      
c             CALL  BFIRSTS(BCKE,BCKP,N,ALA,IPCUBE,NN,RI) sep 01 al
              CALL  BFIRSTS(BCKE,BCKP,N,ALA,IPCUBE,NN)

           ELSEIF(MODE.EQ.2)  THEN
              CALL  FIXEDGE2(BCKP,NMAT,BCKP,N,IPCUBE,NN)      
              CALL  BSECOND(BCKE,BCKP,N,ALA,IPCUBE,NN)

           ELSEIF(MODE.EQ.3)  THEN
              CALL  FIXEDGE3(BCKP,NMAT,BCKP,N,IPCUBE,NN)      
              CALL  BTHIRD(BCKE,BCKP,N,ALA,IPCUBE,NN)
           ENDIF

           AKDEN=0.0

c$omp      parallel do private(kn,j,k,i),reduction(+:AKDEN)
           DO KN=1,NN
              J=IPCUBE(4,KN)
              K=IPCUBE(5,KN)
              DO I=IPCUBE(3,KN),
     &              IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 AKDEN=AKDEN+BCKP(I,J,K)*BCKE(I,J,K)
              ENDDO
           ENDDO
           P=DELSQ/AKDEN

c$omp      parallel do private(kn,j,k,i)
           DO KN=1,NN
              J=IPCUBE(4,KN)
              K=IPCUBE(5,KN)
              DO I=IPCUBE(3,KN),
     &           IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 BCKN(I,J,K) = BCKN(I,J,K)+P*BCKP(I,J,K)
                 CB(I,J,K)   = CB(I,J,K)-P*BCKE(I,J,K)
              ENDDO
           ENDDO

           ERR=DELSQ/BNORM
           CHI2=CHI2-P*DELSQ

           IF (VERBOSE) THEN
C            SKIP THIS IF NOT VERBOSE
             WRITE(NOUT,2041) ITER,ERR,CHI2
2041         FORMAT('  ITERATION: 'I3,' ERROR:',1PE12.4,'
     &                 CHISQ: ',1PE12.4)
           ENDIF

C          STOPPING CRITERIA
           IF (ABS(ERR).LE.ERRM .OR. ABS(CHI2).LE.CHIM) THEN  
              DEALLOCATE (BCKE)
              DEALLOCATE (BCKP)            
              IF (.NOT. VERBOSE) THEN
C                WRITE FINAL ERROR 
                 WRITE(NOUT,2041) ITER,ERR,CHI2
              ENDIF
              RETURN
           ENDIF
       ENDDO
#endif
C      NUMBER OF ITERATIONS EXCEEDED
 
       DEALLOCATE (BCKE)
       DEALLOCATE (BCKP)
#ifdef USE_MPI
       IF (ALLOCATED(BCKE_SUM)) DEALLOCATE(BCKE_SUM)
#endif

       END

C++*********************************************************************
C--*********************************************************************

        SUBROUTINE REDPRCG(N,NANG,CB,ANG,ILIST,IPCUBE,NN,DM,RI,ANGBUF)
 
        INCLUDE 'CMBLOCK.INC'

        DIMENSION     CB(N,N,N),PROJ(N,N),ANG(3,NANG)
        DIMENSION     ILIST(NANG),IPCUBE(5,NN),DM(9,NANG)
        DIMENSION     ANGBUF(4,NANG)
        CHARACTER*80  FINPIC
        CHARACTER*80  FINPAT

        COMMON  /F_SPEC/  FINPAT,NLET,FINPIC
        COMMON /NORMB/    BNORM,CHI2
        DOUBLE PRECISION  ABA,SUS,SSQ

        DATA  INPROJ/98/

#ifdef USE_MPI
        include 'mpif.h'
        INTEGER MYPID, COMM, MPIERR, NPROCS
        INTEGER NANGLOC, NREM, IP, KGLB, NLOC, IPROC
        REAL, ALLOCATABLE, DIMENSION(:,:,:) :: CB_LOC, PRJBUF, PRJLOC
        REAL, ALLOCATABLE, DIMENSION(:,:)   :: ANG_LOC
        INTEGER, ALLOCATABLE, DIMENSION(:)  :: PSIZE
        INTEGER, ALLOCATABLE, DIMENSION(:)  :: NBASE
        DOUBLE PRECISION      ABA_LOC, SUS_LOC, SSQ_LOC
        INTEGER               KLP_LOC, KLS_LOC
C
        COMM = MPI_COMM_WORLD
        CALL MPI_COMM_RANK(COMM, MYPID , MPIERR)
        CALL MPI_COMM_SIZE(COMM, NPROCS, MPIERR)
        ALLOCATE(PSIZE(NPROCS), STAT=IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           WRITE(0,*) ' REPRCG3_Q: FAILED TO ALLOCATE PSIZE.'
           STOP
        ENDIF
        ALLOCATE(NBASE(NPROCS), STAT=IRTFLG)
        IF (IRTFLG .NE. 0) THEN
           WRITE(0,*) ' REPRCG3_Q: FAILED TO ALLOCATE NBASE.'
           STOP
        ENDIF
#endif
        ABA=0.0D0
        KLP=0
        SUS=0.0D0
        SSQ=0.0D0
        KLS=0
        
        ABA_LOC = 0.0D0 
        KLP_LOC = 0
        KLS_LOC = 0
        SUS_LOC = 0.0D0
        SSQ_LOC = 0.0D0

c$omp   parallel do private(i,j,k)
        DO    K=1,N
           DO    J=1,N
              DO    I=1,N
                 CB(I,J,K)=0.0
              ENDDO
           ENDDO
        ENDDO

#ifdef USE_MPI
        ALLOCATE (CB_LOC(N,N,N), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN
           CALL ERRT(46,'REPRCG3_Q, CB_LOC',IER)
           RETURN
        ENDIF
        CB_LOC = 0.0
        CALL SETPART(NANG, PSIZE, NBASE)
        NANGLOC = PSIZE(MYPID+1)
C
        ALLOCATE (ANG_LOC(3,NANGLOC), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN
           CALL ERRT(46,'REPRCG3_Q, ANG_LOC',IER)
           RETURN
        ENDIF
        ANG_LOC = 0.0
C
        ALLOCATE (PRJBUF(N,N,PSIZE(1)),PRJLOC(N,N,NANGLOC),
     &            STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN
           CALL ERRT(46,'REPRCG3_Q, PRJBUF, PRJLOC',IER)
           RETURN
        ENDIF
C
        DO IPROC = 1, NPROCS
           NLOC = PSIZE(IPROC)

C          READ A SUBSET OF IMAGES (ONLY ONE PROCESSOR READS)

           DO K=1,NLOC
              KGLB = K + NBASE(IPROC)
              CALL  FILGET(FINPAT,FINPIC,NLET,ILIST(KGLB),IRTFLG)
              IF (IRTFLG .NE. 0) RETURN
              MAXIM = 0
              CALL OPFILEC(0,.FALSE.,FINPIC,INPROJ,'O',IFORM,
     &                     LSAM,LROW,NSL,
     &                     MAXIM,'DUMMY',.FALSE.,IRTFLG)
              IF (IRTFLG .NE. 0) RETURN
C
              DO K2=1,N
                 CALL  REDLIN1P(INPROJ,PRJBUF(1,K2,K),N,K2)
              ENDDO
              IF (MYPID .EQ. 0) CLOSE(INPROJ)
           ENDDO

C          DISTRIBUTE IMAGES

           IF (IPROC .GT. 1) THEN
              IF (MYPID .EQ. 0) THEN
                 CALL MPI_SEND(PRJBUF , N*N*NLOC, MPI_REAL,
     &                         IPROC-1, IPROC-1 , COMM    ,
     &                         MPIERR)
                 IF (MPIERR .NE. 0) THEN
                     WRITE(6,*) 'REDPRS: SEND FAILED'
                     STOP
                 ENDIF
              ELSE IF (MYPID .EQ. IPROC-1) THEN
                 CALL MPI_RECV(PRJLOC, N*N*NLOC   , MPI_REAL,
     &                         0     , MPI_ANY_TAG, COMM    ,
     &                         ISTAT , MPIERR)
                 IF (MPIERR .NE. 0) THEN
                    WRITE(6,*) 'REDPRS: RECV FAILED'
                    STOP
                 ENDIF
              ENDIF
           ELSE IF (MYPID .EQ. 0) THEN
              CALL  SCOPY(N*N*NLOC,PRJBUF,1,PRJLOC,1)
           ENDIF 
        ENDDO

        DO K = 1, NANGLOC
           KGLB = K + NBASE(MYPID+1)

C          ORDER IN DOCUMENT FILE IS PSI, THETA, PHI AND ANGLES ARE IN
C          DEGREES!
C          IN ANG ARRAY IT IS OTHER WAY AROUND

           ITMP     = ILIST(KGLB)
           ICOUNT   = ANGBUF(1,ITMP)
           ANG_LOC(3,K) = ANGBUF(2,ITMP)
           ANG_LOC(2,K) = ANGBUF(3,ITMP)
           ANG_LOC(1,K) = ANGBUF(4,ITMP)

           CALL  ASTASQ(PRJLOC(1,1,K), N, RI, ABA_LOC, KLP_LOC, 
     &                  SUS_LOC, SSQ_LOC, KLS_LOC)
 
           CALL  RPRQ(N,PRJLOC(1,1,K),CB_LOC,IPCUBE,NN,
     &     ANG_LOC(1,K),ANG_LOC(2,K),ANG_LOC(3,K),DM(1,KGLB),RI)
        ENDDO
        IF (ALLOCATED(PRJBUF)) DEALLOCATE(PRJBUF) 
        IF (ALLOCATED(PRJLOC)) DEALLOCATE(PRJLOC) 

C       GATHER ANG_LOC INTO ANG

        DO J = 1, NPROCS
           NBASE(J) = 3*NBASE(J)
           PSIZE(J) = 3*PSIZE(J)
        END DO
        CALL MPI_ALLGATHERV(ANG_LOC,  3*NANGLOC, MPI_REAL,
     &                      ANG    ,  PSIZE    , NBASE   ,
     &                      MPI_REAL, COMM     , IERR)    

        DO K = 1, NANG
           IF (VERBOSE) THEN
              IF (MYPID .EQ. 0) THEN
                 WRITE(NOUT,3331)  K,(ANG(J,K),J=3,1,-1)
              ENDIF
3331          FORMAT(' PROJECTION #',I6,
     &            '; PSI=',F6.1,' THETA=',F6.1,' PHI=',F6.1)
           ENDIF
        END DO

        N3 = N*N*N
        IERR = 0
        CALL MPI_ALLREDUCE(CB_LOC , CB  , N3,  MPI_REAL,
     &                     MPI_SUM, COMM, IERR)
        IF (IERR .NE. 0) THEN
           WRITE(0,*) 'REPRCG3_Q: FAILED AT ALLREDUCE CB'
           STOP
        END IF
        CALL MPI_ALLREDUCE(ABA_LOC, ABA , 1,   MPI_DOUBLE_PRECISION,
     &                     MPI_SUM, COMM, IERR)
        IF (IERR .NE. 0) THEN
           WRITE(0,*) 'REPRCG3_Q: FAILED AT ALLREDUCE ABA'
           STOP
        END IF
        CALL MPI_ALLREDUCE(KLP_LOC, KLP , 1,   MPI_INTEGER,
     &                     MPI_SUM, COMM, IERR)
        IF (IERR .NE. 0) THEN
           WRITE(0,*) 'REPRCG3_Q: FAILED AT ALLREDUCE KLP'
           STOP
        END IF
        CALL MPI_ALLREDUCE(SSQ_LOC, SSQ , 1,   MPI_DOUBLE_PRECISION,
     &                     MPI_SUM, COMM, IERR)
        IF (IERR .NE. 0) THEN
           WRITE(0,*) 'REPRCG3_Q: FAILED AT ALLREDUCE ABA'
           STOP
        END IF
        CALL MPI_ALLREDUCE(SUS_LOC, SUS , 1,   MPI_DOUBLE_PRECISION,
     &                     MPI_SUM, COMM, IERR)
        IF (IERR .NE. 0) THEN
           WRITE(0,*) 'REPRCG3_Q: FAILED AT ALLREDUCE ABA'
           STOP
        END IF
        CALL MPI_ALLREDUCE(KLS_LOC, KLS , 1,   MPI_INTEGER,
     &                     MPI_SUM, COMM, IERR)
        IF (IERR .NE. 0) THEN
           WRITE(0,*) 'REPRCG3_Q: FAILED AT ALLREDUCE KLP'
           STOP
        END IF
        IF (ALLOCATED(CB_LOC))  DEALLOCATE(CB_LOC)
        IF (ALLOCATED(ANG_LOC)) DEALLOCATE(ANG_LOC)
        IF (ALLOCATED(PSIZE))   DEALLOCATE(PSIZE)
        IF (ALLOCATED(NBASE))   DEALLOCATE(NBASE)
#else
        DO    K=1,NANG
        CALL  FILGET(FINPAT,FINPIC,NLET,ILIST(K),IRTFLG)
           IF (IRTFLG .NE. 0) RETURN

           MAXIM = 0
           CALL OPFILEC(0,.FALSE.,FINPIC,INPROJ,'O',IFORM,
     &                  LSAM,LROW,NSL,
     &                  MAXIM,'DUMMY',.FALSE.,IRTFLG)
           IF (IRTFLG .NE. 0) RETURN

           DO K2=1,N
              CALL  REDLIN(INPROJ,PROJ(1,K2),N,K2)
           ENDDO
           CLOSE(INPROJ)

C          ORDER IN DOCUMENT FILE IS PSI, THETA, PHI AND ANGLES ARE IN 
C          DEGREES! IN ANG ARRAY IT IS OTHER WAY AROUND

           ITMP     = ILIST(K)
           ICOUNT   = ANGBUF(1,ITMP)
           IF (ICOUNT .LE. 0) THEN
C             MISSING KEY
              CALL ERRT(102,'MISSING ANGLE FOR IMAGE',ITMP)
              RETURN
           ENDIF

           ANG(3,K) = ANGBUF(2,ITMP)
           ANG(2,K) = ANGBUF(3,ITMP)
           ANG(1,K) = ANGBUF(4,ITMP)
           IF (VERBOSE) THEN
              WRITE(NOUT,3331)  K,(ANG(J,K),J=3,1,-1)
3331          FORMAT(' PROJECTION #',I6,
     &            '; PSI=',F6.1,' THETA=',F6.1,' PHI=',F6.1)
           ENDIF

           CALL  ASTASQ(PROJ,N,RI,ABA,KLP,SUS,SSQ,KLS)

           CALL  RPRQ(N,PROJ,CB,IPCUBE,NN,
     &     ANG(1,K),ANG(2,K),ANG(3,K),DM(1,K),RI)
        ENDDO
#endif

C       CLOSE DOCUMENT FILE
        IF (MYPID .LE. 0) CLOSE(77)
        ABA=ABA/KLP
C       PRINT STATISTICS
        IF (MYPID .LE. 0) WRITE(NOUT,2044)   KLP,ABA
2044    FORMAT
     &      ('  TOTAL NUMBER OF POINTS IN PROJECTIONS =',I10,
     &     /,'  AVERAGE OUTSIDE THE WINDOW            =',1PE10.3,/)

C       SUBTRACT THE AVERAGE AND WRITE TO THE FILE

        BNORM=0.0
        QT=ABA*NANG

c$omp   parallel do private(kn,j,k,i),reduction(+:BNORM)
        DO KN=1,NN
           J=IPCUBE(4,KN)
           K=IPCUBE(5,KN)
           DO    I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
              CB(I,J,K)=CB(I,J,K)-QT
              BNORM=BNORM+CB(I,J,K)*CB(I,J,K)
           ENDDO
        ENDDO

C       CALCULATE SUM OF SQUARES OF PROJECTIONS WITHIN CIRCLES AFTER
C       SUBTRACTION OF THE AVERAGE
        SSQ=SSQ-ABA*(2*SUS-ABA*KLS)
        CHI2=SSQ
        IF (MYPID .LE. 0) WRITE(NOUT,2045)   BNORM,CHI2
2045    FORMAT ('  SQUARED BP          =',1PE10.3,
     &        /,'  CHI2 =',1PE10.3,/)

        END

C       --------------------- ASTASQ -------------------------------

        SUBROUTINE  ASTASQ(X,N,RI,ABA,KLP,SUS,SSQ,KLS)

        DIMENSION   X(N,N)
        DOUBLE PRECISION  ABA,SUS,SSQ

C       ESTIMATE AVERAGE OUTSIDE THE CIRCLE
        R=RI*RI
        NC=N/2+1
        DO   J=1,N
           T=J-NC
           XX=T*T
           DO   I=1,N
              T=I-NC
              IF(XX+T*T.GT.R)    THEN
                 ABA=ABA+DBLE(X(I,J))
                 KLP=KLP+1
              ELSE
                 SSQ=SSQ+X(I,J)*DBLE(X(I,J))
                 SUS=SUS+X(I,J)
                 KLS=KLS+1
              ENDIF
           ENDDO
        ENDDO
        END

C++*********************************************************************
C--*********************************************************************

        SUBROUTINE  FIXEDGE1(BCKP,NMAT,BCK3,N,IPCUBE,NN)

        DIMENSION    BCKP(NMAT),BCK3(N,N,N),IPCUBE(5,NN)

C-------------------------------------------------------
C   PUT ZEROS OUTSIDE
C-------------------------------------------------------
        NT=1
        DO    I=1,NMAT
           IF(NT.GT.NN)  THEN
              BCKP(I)=0.0
           ELSEIF(I.LT.IPCUBE(1,NT))  THEN
              BCKP(I)=0.0
           ELSEIF(I.EQ.IPCUBE(2,NT))  THEN
              NT=NT+1
           ENDIF
        ENDDO
C  ADD PIXELS ON THE EDGE
C FIX THE EDGES IN BCKP
        DO  K=1,N
           DO J=1,N
              DO  I=1,N-1
                 IF(BCK3(I+1,J,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I+1,J,K)
                    GOTO  101
                 ENDIF
              ENDDO
101     CONTINUE
              DO  I=N,2,-1
                 IF(BCK3(I-1,J,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I-1,J,K)
                    GOTO  102
                 ENDIF
              ENDDO
102     CONTINUE
           ENDDO
        ENDDO
        DO  K=1,N
           DO I=1,N
              DO  J=1,N-1
                 IF(BCK3(I,J+1,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J+1,K)
                    GOTO  201
                 ENDIF
              ENDDO
201     CONTINUE
              DO  J=N,2,-1
                 IF(BCK3(I,J-1,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J-1,K)
                    GOTO  202
                 ENDIF
              ENDDO
202     CONTINUE
           ENDDO
        ENDDO
        DO  J=1,N
           DO I=1,N
              DO  K=1,N-1
                 IF(BCK3(I,J,K+1).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J,K+1)
                    GOTO  301
                 ENDIF
              ENDDO
301     CONTINUE
              DO  K=N,2,-1
                 IF(BCK3(I,J,K-1).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J,K-1)
                    GOTO  302
                 ENDIF
              ENDDO
302     CONTINUE
           ENDDO
        ENDDO
        END

C++*********************************************************************
C--*********************************************************************

        SUBROUTINE  FIXEDGE2(BCKP,NMAT,BCK3,N,IPCUBE,NN)
        DIMENSION    BCKP(NMAT),BCK3(N,N,N),IPCUBE(5,NN)

C-------------------------------------------------------
C   PUT ZEROS OUTSIDE
C-------------------------------------------------------
        NT=1
        DO    I=1,NMAT
           IF(NT.GT.NN)  THEN
              BCKP(I)=0.0
           ELSEIF(I.LT.IPCUBE(1,NT))  THEN
              BCKP(I)=0.0
           ELSEIF(I.EQ.IPCUBE(2,NT))  THEN
              NT=NT+1
           ENDIF
        ENDDO
C  ADD PIXELS ON THE EDGE
C FIX THE EDGES IN BCKP
        DO  K=1,N
           DO J=1,N
              DO  I=2,N-1
                 IF(BCK3(I+1,J,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I+1,J,K)
                    BCK3(I-1,J,K)=BCK3(I+1,J,K)
                    GOTO  101
                 ENDIF
              ENDDO
101     CONTINUE
              DO  I=N-1,2,-1
                 IF(BCK3(I-1,J,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I-1,J,K)
                    BCK3(I+1,J,K)=BCK3(I-1,J,K)
                    GOTO  102
                 ENDIF
             ENDDO
102     CONTINUE
           ENDDO
        ENDDO
        DO  K=1,N
           DO I=1,N
              DO  J=2,N-1
                 IF(BCK3(I,J+1,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J+1,K)
                    BCK3(I,J-1,K)=BCK3(I,J+1,K)
                    GOTO  201
                 ENDIF
              ENDDO
201           CONTINUE
              DO  J=N-1,2,-1
                 IF(BCK3(I,J-1,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J-1,K)
                    BCK3(I,J+1,K)=BCK3(I,J-1,K)
                    GOTO  202
                 ENDIF
              ENDDO
202           CONTINUE
           ENDDO
        ENDDO
        DO  J=1,N
           DO I=1,N
              DO  K=2,N-1
                 IF(BCK3(I,J,K+1).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J,K+1)
                    BCK3(I,J,K-1)=BCK3(I,J,K+1)
                    GOTO  301
                 ENDIF
              ENDDO
301           CONTINUE
              DO  K=N-1,2,-1
                 IF(BCK3(I,J,K-1).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J,K-1)
                    BCK3(I,J,K+1)=BCK3(I,J,K-1)
                    GOTO  302
                 ENDIF
              ENDDO
302           CONTINUE
           ENDDO
        ENDDO
        END

C++*********************************************************************
C--*********************************************************************

        SUBROUTINE  FIXEDGE3(BCKP,NMAT,BCK3,N,IPCUBE,NN)

        DIMENSION    BCKP(NMAT),BCK3(N,N,N),IPCUBE(5,NN)

C-------------------------------------------------------
C   PUT ZEROS OUTSIDE
C-------------------------------------------------------
        NT=1
        DO    I=1,NMAT
           IF(NT.GT.NN)  THEN
              BCKP(I)=0.0
           ELSEIF(I.LT.IPCUBE(1,NT))  THEN
              BCKP(I)=0.0
           ELSEIF(I.EQ.IPCUBE(2,NT))  THEN
              NT=NT+1
           ENDIF
        ENDDO
C  ADD PIXELS ON THE EDGE
C FIX THE EDGES IN BCKP
        DO  K=1,N
           DO J=1,N
              DO  I=3,N-1
                 IF(BCK3(I+1,J,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I+1,J,K)
                    BCK3(I-1,J,K)=BCK3(I+1,J,K)
                    BCK3(I-2,J,K)=BCK3(I+1,J,K)
                    GOTO  101
                 ENDIF
              ENDDO
101           CONTINUE
              DO  I=N-2,2,-1
                 IF(BCK3(I-1,J,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I-1,J,K)
                    BCK3(I+1,J,K)=BCK3(I-1,J,K)
                    BCK3(I+2,J,K)=BCK3(I-1,J,K)
                    GOTO  102
                 ENDIF
              ENDDO
102           CONTINUE
           ENDDO
        ENDDO
        DO  K=1,N
           DO I=1,N
              DO  J=3,N-1
                 IF(BCK3(I,J+1,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J+1,K)
                    BCK3(I,J-1,K)=BCK3(I,J+1,K)
                    BCK3(I,J-2,K)=BCK3(I,J+1,K)
                    GOTO  201
                 ENDIF
              ENDDO
201           CONTINUE
              DO  J=N-2,2,-1
                 IF(BCK3(I,J-1,K).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J-1,K)
                    BCK3(I,J+1,K)=BCK3(I,J-1,K)
                    BCK3(I,J+2,K)=BCK3(I,J-1,K)
                    GOTO  202
                 ENDIF
              ENDDO
202           CONTINUE
           ENDDO
        ENDDO
        DO  J=1,N
           DO I=1,N
              DO  K=3,N-1
                 IF(BCK3(I,J,K+1).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J,K+1)
                    BCK3(I,J,K-1)=BCK3(I,J,K+1)
                    BCK3(I,J,K-2)=BCK3(I,J,K+1)
                    GOTO  301
                 ENDIF
              ENDDO
301           CONTINUE
              DO  K=N-2,2,-1
                 IF(BCK3(I,J,K-1).NE.0.0) THEN
                    BCK3(I,J,K)=BCK3(I,J,K-1)
                    BCK3(I,J,K+1)=BCK3(I,J,K-1)
                    BCK3(I,J,K+2)=BCK3(I,J,K-1)
                    GOTO  302
                 ENDIF
              ENDDO
302           CONTINUE
           ENDDO
        ENDDO
        END

C++*********************************************************************
C--*********************************************************************

        SUBROUTINE  BFIRSTS(BCKE,BCKP,N,ALA,IPCUBE,NN)

        DIMENSION  BCKE(N,N,N),BCKP(N,N,N),ipcube(5,nn)
c$omp   parallel do private(kn,j,k,i)
        DO    KN=1,NN
           J=IPCUBE(4,KN)
           K=IPCUBE(5,KN)
           DO    I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
              BCKE(I,J,K)=BCKE(I,J,K)+ALA*(6*BCKP(I,J,K)-(
     &        +BCKP(I+1,J,K)+BCKP(I,J+1,K)+BCKP(I,J,K+1)
     &        +BCKP(I-1,J,K)+BCKP(I,J-1,K)+BCKP(I,J,K-1)))
           ENDDO
        ENDDO
        END

C++*********************************************************************
C--*********************************************************************

        SUBROUTINE  BSECOND(BCKE,BCKP,N,ALA,IPCUBE,NN)

        DIMENSION  BCKE(N,N,N),BCKP(N,N,N),ipcube(5,nn)

c$omp parallel do private(kn,j,k,i)
        DO    KN=1,NN
           J=IPCUBE(4,KN)
           K=IPCUBE(5,KN)
           DO    I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 BCKE(I,J,K)=BCKE(I,J,K)+ALA*(18*BCKP(I,J,K)
     &          -4.0*BCKP(I+1,J,K)+BCKP(I+2,J,K)
     &          -4.0*BCKP(I-1,J,K)+BCKP(I-2,J,K)
     &          -4.0*BCKP(I,J+1,K)+BCKP(I,J+2,K)
     &          -4.0*BCKP(I,J-1,K)+BCKP(I,J-2,K)
     &          -4.0*BCKP(I,J,K+1)+BCKP(I,J,K+2)
     &          -4.0*BCKP(I,J,K-1)+BCKP(I,J,K-2)
     &          )
           ENDDO
        ENDDO
        END

C++*********************************************************************
C--*********************************************************************

        SUBROUTINE  BTHIRD(BCKE,BCKP,N,ALA,IPCUBE,NN)

        DIMENSION  BCKE(N,N,N),BCKP(N,N,N),ipcube(5,nn)

c$omp parallel do private(kn,j,k,i)
        DO    KN=1,NN
           J=IPCUBE(4,KN)
           K=IPCUBE(5,KN)
           DO    I=IPCUBE(3,KN),IPCUBE(3,KN)+IPCUBE(2,KN)-IPCUBE(1,KN)
                 BCKE(I,J,K)=BCKE(I,J,K)+ALA*(60*BCKP(I,J,K)
     &          -15.0*BCKP(I+1,J,K)+6.0*BCKP(I+2,J,K)-BCKP(I+3,J,K)
     &          -15.0*BCKP(I-1,J,K)+6.0*BCKP(I-2,J,K)-BCKP(I-3,J,K)
     &          -15.0*BCKP(I,J+1,K)+6.0*BCKP(I,J+2,K)-BCKP(I,J+3,K)
     &          -15.0*BCKP(I,J-1,K)+6.0*BCKP(I,J-2,K)-BCKP(I,J-3,K)
     &          -15.0*BCKP(I,J,K+1)+6.0*BCKP(I,J,K+2)-BCKP(I,J,K+3)
     &          -15.0*BCKP(I,J,K-1)+6.0*BCKP(I,J,K-2)-BCKP(I,J,K-3)
     &          )
           ENDDO
        ENDDO
        END