      SUBROUTINE TRAHALO_RUN(MDOC
     * ,NAMEF,NAMED,NAMEO,NAMEOA,NAMEA,NAMEHA,NAMECPH
     * ,RESMINR,RESMAXR,STR_PARAM,RPARAM_TRAH
     * ,NP,POOL,MEMORY,ISYM,IPRSYM
     * ,RESULT,RESULTP,RESULTB,RESULTA,RESULTC,N_FINAL
     * ,IPRMAX,IPRROW,IPRNPR,IPRMAXD,NCOMM,IERR)
C --------------------------------------------------------------------
      INTEGER*2 ISYM(5,3,IPRSYM)
      REAL      POOL(MEMORY)
CMS$LARGE: POOL
C ----
C     SLIM == SMOOTHRR
C     CLIM == DIST_R
C 
C     !!!!!  IPRROWW must be iqual IPRROW
C ---
      COMMON/COMDIF/ NUMB_CFF,SLIM_PATT,BADD_COEF,DEL_LIM,DELA_LIM
C ---
      COMMON/RESDER/ NO_DER,RES_DER
      PARAMETER (IPRMAXDD = 71)
      PARAMETER (IPRROWW  = 21)
      PARAMETER (MAXDER   =  4)
      REAL      RES_DER(IPRMAXDD,IPRROWW,MAXDER)

      REAL      RESULTD(IPRMAXDD,IPRROWW)
C ---  
C      main_trahalo.f
c      PARAMETER  ( IPRMAXD = 71      )
c      PARAMETER  ( IPRMAX  = 21      )
c      PARAMETER  ( IPRROW  = 21      )
c      PARAMETER  ( IPRNPR  = 21      )
C      SUBR   = 'N'
C      NO_DER = 0
C      NPST   = 0
C      IERR   = 0
C ---
C -------------
C      SIR.F
C       level of messages
C       IERR = IPLEVEL
C -------------
C  NAMEA   - intput file of ABCD 
C  NAMECPH - intput file atoms of ABCD 
C  NAMEHA  - intput file with known heavy atoms
C ---------------------------------------------------------
      REAL      RESULT (IPRMAX ,IPRROW)
      REAL      RESULTB(IPRMAXD,IPRROW)
      REAL      RESULTC(IPRMAXD,IPRROW)
      REAL      RESULTA(IPRMAXD,IPRROW,IPRNPR)

      REAL      RESULTP(21,21)
      REAL      RPARAM_TRAH(10)      
C ---
      PARAMETER ( NCRDMAX = 100 )
      PARAMETER ( IPRORG=24)
      REAL      SOR(3,IPRORG)
C ------------------------------------------------
      PARAMETER ( IPRPST=1)
      COMMON/COMPST/ NPST,SPST,RESEFF,BRES,BEFF
      REAL      SPST(3,IPRPST)
C ------------------------------------------------
      PARAMETER ( IPRNEQ=24)
      INTEGER   INVERS1(IPRNEQ),IORIGN1(IPRNEQ),ISYMOP1(IPRNEQ)
      REAL      DISTN1(IPRNEQ)
      INTEGER   INVERS2(IPRNEQ),IORIGN2(IPRNEQ),ISYMOP2(IPRNEQ)
      REAL      DISTN2(IPRNEQ)
C ---
      CHARACTER STR_PARAM*(*)


      CHARACTER NAMEFS*80, NAMEC2*80,NAMEFS2*80,NAMEPH2*80, NAMEF2*80
      CHARACTER NAMEF*(*), NAMED*(*),NAMEO*(*), NAMEA*(*),  NAMET*80
      CHARACTER NAMES*80,  NAMEP*80, NAMEC*80,  NAME_SLF*80,NAMEPH*80
      CHARACTER NAMEHA*(*),NAMECPH*(*),NAMEOA*(*)
      CHARACTER NAMECA*80, NAME_SLFA*80,NAME_SLFN*80,NAMEAS*80

      CHARACTER LINE*80,TTT*80,ITYPE*4,MODE*1,RFAC*1,GAUSS*1,ANOM*1
      CHARACTER MOD2*1,FUNC*1,ACCUM*1,REMOV*1,INVER*1,NAMES2*80
      CHARACTER STOP_T*1,MSGA*1,ASYM*1,LAP*1,FOUR*1,WAT*1,P1*1
      CHARACTER MSCL*1,LSORT*1,TRANS*1,FAST*1,SUBR*1,FASTR*1,MISS*1
      CHARACTER ALT*1,SETT*1,REST*1,HATM*1,NUM*1,NOCC*1,INVERT*1
      CHARACTER PACK*1,SPEC*1,SPEC_P*1,COMB*1,MSGF*1,FOCC*1

      CHARACTER NAME(10)*80,FILE1*80,FILE2*80,NAMEC3*80,OPT1*1
      CHARACTER PDB*1,STR_TITLE*40,STR_DATE*9,PATT*1,SOLV*1,ANISO*1
      CHARACTER MODE_F*4,MODE_CC*1,AANOM*1

      REAL      AM(3,3),SH2(3,10),SH1(3),PERC2(10)
C -------------------------------------------------------------
      PI  = 4.0*ATAN(1.0)
C     level of messages
      IPLEVEL = IERR
      IERR    = 0
      N_FINAL = 0
C ---
      MD  = -ABS(MDOC)-1
      M   = 99
      MDD = MD
      IF(IPLEVEL.GT.0) MDD = 99
C ---
      TDEL_LIM  = RPARAM_TRAH(1) 
      SLIM_PATT = RPARAM_TRAH(2)
      BOFF      = RPARAM_TRAH(3)
      BADDR     = RPARAM_TRAH(4)
      RADR      = RPARAM_TRAH(5)
      TIMES_RMS = RPARAM_TRAH(6)
      SLIM_SAD  = RPARAM_TRAH(7)
      BADD_TRP  = RPARAM_TRAH(8)
      BADD_FOUR = RPARAM_TRAH(9)

      BADD_COEF = BADDR
      SMOOTH    = BADD_TRP
      RAD       = RADR

C ----
      INVERT = STR_PARAM(1:1)
C mscl doesn't used
      MSCL   = STR_PARAM(2:2)
      SPEC   = STR_PARAM(3:3)
      PACK   = STR_PARAM(4:4)
      FASTR  = STR_PARAM(5:5)
      TRANS  = STR_PARAM(6:6)
      GAUSS  = STR_PARAM(7:7)
      LSORT  = STR_PARAM(8:8)
      REMOV  = STR_PARAM(9:9)
C     comb = N/A/Y
      COMB   = STR_PARAM(10:10)
      ANISO  = STR_PARAM(11:11)
C     remove HA of input abcd == N
      OPT1   = STR_PARAM(12:12)
      SUBR   = STR_PARAM(13:13)
C ??????
      IF(SUBR.EQ.'N') MDD = MDOC
C ??????
      STOP_T   = 'N'
      FOCC     = 'N'
C ---
      BADD_REF  = 0.0
      BADD      = 0.0
C ///
      DIST_R    = 0.0
      DIST_LIM  = 0.0
      SLIM      = 0.0
C ///
 
      DO I=1,IPRMAX
      DO J=1,IPRROW
        RESULT(I,J) = 0.0
      ENDDO
      ENDDO

      DO I=1,IPRMAXDD
      DO J=1,IPRROWW
        RESULTD (I,J) = 0.0
      ENDDO
      ENDDO

      DO I=1,IPRMAXD
      DO J=1,IPRROW
        RESULTB (I,J) = 0.0
        RESULTC (I,J) = 0.0
        DO K=1,IPRNPR
          RESULTA(I,J,K) = 0.0
        ENDDO
      ENDDO
      ENDDO

      NCOMM = 0

C  -----------------------------------------------------------
C     check Fobs file
      IERR  = 0
      MTN   = 10
      ITYPE = 'FF  '
      CALL LENSTR_BL(NAMEF,LEN)
      IF(LEN.GT.0.AND.NAMEF(1:1).NE.','.AND.NAMEF(1:1).NE.' ') THEN
        CALL ORFILE(M,MTN,NAMEF,ITYPE,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) THEN
          CALL MSGERR(MDOC,' ERR: OPEN INPUT-FILE-FOBS')
          RETURN
        ENDIF
        CALL GET_TITLE(M,A,B,C,AL,BE,GA
     *  ,F1,F2,RN1,RX1,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
        CLOSE(MTN)
        CALL C_RD_CR      
        ANOM = 'N'
        IF(SUBR.EQ.'N') THEN
          IF(LEN.GT.60) LEN=60
          LINE='Input file Fobs:'//NAMEF(1:LEN)
          CALL MSGDOC(MD,LINE)
        ENDIF
      ELSE
        MTN       = 0
        ANOM      = 'Y'
        NAMEF     = ' '
        SLIM_PATT = SLIM_SAD
      ENDIF
C -------------------------------------------------------------
C     check Fder file
      IERR  = 0
      MTD   = 11
      ITYPE = 'FF  '
      CALL ORFILE(M,MTD,NAMED,ITYPE,NSYM,ISYM,IPRSYM,IERR)
      IF(IERR.NE.0) THEN
        CALL MSGERR(MDOC,' ERR: OPEN INPUT-FILE-FDER')
        RETURN
      ENDIF

      IF(SUBR.EQ.'N') THEN
        CALL LENSTR_BL(NAMED,LEN)
        IF(LEN.GT.60) LEN=60
        LINE='Input file Fder:'//NAMED(1:LEN)
        CALL MSGDOC(MD,LINE)
      ENDIF

      CALL GET_TITLE(M,A,B,C,AL,BE,GA
     *  ,F1,F2,RN,RX,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
      CLOSE(MTD)
      IF(MTN.EQ.0) CALL C_RD_CR      
C --------------------------------------------------------------
C     check input file ABCD
      CALL LENSTR_BL(NAMEA,LEN)
      IF(LEN.GT.0.AND.NAMEA(1:1).NE.','.AND.NAMEA(1:1).NE.' ') THEN
        IF(SUBR.EQ.'N') THEN
          CALL LENSTR_BL(NAMEA,LEN)
          IF(LEN.GT.60) LEN=60
          LINE='Input file ABCD:'//NAMEA(1:LEN)
          CALL MSGDOC(MD,LINE)
        ENDIF
        FOUR = 'Y'
      ELSE
        FOUR = 'N'
      ENDIF      
C --------------------------------------------------------------
      IF(FASTR.NE.'F'.AND.FASTR.NE.'M'.AND.FASTR.NE.'S'.AND.
     *   FASTR.NE.'Y'.AND.FASTR.NE.'N') FASTR = 'Y'

      FAST = FASTR

      IF(FAST.EQ.'F') FAST = 'Y'
      IF(FAST.EQ.'M') FAST = 'N'
      IF(FAST.EQ.'S') FAST = 'N'

      IF(FAST.NE.'N') FAST = 'Y'

      IF(FAST.EQ.'Y') THEN
        TRANS = 'N'
        IF(NP.EQ.0) NP = 6
      ENDIF

      RAD_DEF = 2.5
C ----------------------------------------------------
C     resolution
C 
      RES_DEF = 3.0     
      IF(FAST.EQ.'Y') RES_DEF = 5.0

C     resolution as subroutine's parameters
      RESMIN = RESMINR
      RESMAX = RESMAXR
   
C     resolution from Fder
      RESMIN_D = RN
      RESMAX_D = RX

      IF(MTN.GT.0) THEN
C       resolution from Fobs if there is this file
        IF(RN1.LT.RESMIN_D)   RESMIN_D = RN1
        IF(RX1.GT.RESMAX_D)   RESMAX_D = RX1
      ENDIF    

      IF((RESMIN_D.LT.RESMIN).OR.RESMIN.LE.0.0000001)
     *   RESMIN = RESMIN_D
      IF(RESMAX_D.GT.RESMAX)   RESMAX = RESMAX_D


      IF(RESMAX.LT.RES_DEF .AND.RESMAXR.LE.0.0000001) 
     *   RESMAX = RES_DEF

      IF(RESMIN.GT.20.0.AND.RESMINR.LE.0.0) RESMIN = 20.0
C ----------------------------------------------------

      IF(BOFF.EQ.0.0) BOFF = 400.0
      IF(BOFF.LT.0.0) BOFF = 0.0

      IF(RAD.LE.0.0)  RAD = RAD_DEF

      IF(NP.EQ.0) NP = 20
      IF(NP.LT.0) NP = 0
      IF(NP.GE.IPRNPR) NP = IPRNPR-1

      NAMES     = ' '
      NAMES2    = ' '
      NAMEC2    = 'scratch'
      NAMEC     = 'trahalo_cff'
      NAMECA    = 'trahalo_cffa'
      NAMEC3    = 'trahalo_fft'
      NAMEP     = 'trahalo_dns'
      NAMEFS    = 'trahalo_fsc'
      NAME_SLF  = 'trahalo_slf'
      NAME_SLFA = 'trahalo_slfa'
      NAME_SLFN = 'trahalo_slfn'
      NAMEPH    = 'trahalo_ph'
      NAMET     = 'trahalo_trp'
      NAMEF2    = 'trahalo_f2'
      NAMEFS2   = 'trahalo_fs2'
      NAMEPH2   = 'trahalo_ph2'
      NAMEAS    = NAMEA
C ----
      CALL CALC_ORIG_PS(MDOC,NORIG,SOR,IPRORG
     *   ,NSYM,ISYM,IPRSYM,IERR)

C     IF(NPST.GT.0) THEN
C       CALL CHECK_LIST_ORIG
C     ENDIF
C ---

      IF(DIST_LIM.LE.0.0) THEN
C       2.05.97
        DIST_LIM = RESMAX/2.0
C       DIST_LIM=RESMAX
C       till 2.05.97
C       DIST_LIM=3.0
      ENDIF

      RESMX = RESMAX/2.0

      IF(FAST.EQ.'Y') RESMX = RESMAX*0.75
      CALL LIMHKL(RESMX,IHMAX,IKMAX,ILMAX,NXMIN,NYMIN,NZMIN)

      CALL CHECK_REF(NSYM,ISYM,IFX,IFY,IFZ)             
C     IFX = 0  - not fixed
C ----
      CALL MSGDOC(MD,'grid spacing for calculation SF:') 
      WRITE(LINE,1000) NXMIN,NYMIN,NZMIN
 1000 FORMAT(' NX=',I3,', NY=',I3,', NZ=',I3) 
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,*) ' Resolution   :',RESMIN,RESMAX
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,*) ' GAUSS,RAD    : ',GAUSS,' ',RAD
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,*) ' Boff         :',BOFF
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,*) ' Badd_C,SMOOTH:',BADD_COEF,SMOOTH
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,*) ' SLIM         :',SLIM_PATT
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,*) ' CLIM         :',TDEL_LIM
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,*) ' DIST_LIM     :',DIST_LIM
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,*) ' FAST_mode    : ',FASTR
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,*) 
     * ' SORT,PACK,SPEC,TRANS,REMOV,COMB: ',LSORT,' ',PACK,' ',SPEC,' '
     *   ,TRANS,' ',REMOV,' ',COMB
      CALL MSGDOC(MD,LINE)
C ---
      CALL LENSTR_BL(NAMEO,LEN)
      IF(LEN.GT.0.AND.NAMEO(1:1).NE.','.AND.NAMEO(1:1).NE.' ') THEN      

      ELSE
        NAMEO  ='trahalo'
      ENDIF
C ====================================================================
C --------------------------------------------
      NAMEFS = NAMED
      IF(MTN.NE.0) THEN
        IF(BADDR.LT.0) THEN
          B = BRES-BEFF
          IF(B.LT.0.0) B = 0.0
          BADD = B   
          WRITE(LINE,'(''B_res,B_eff,Badd_new:'',3F10.3)')
     *                 BRES,BEFF,BADD
          CALL MSGDOC(MDOC,LINE)
          IF(SMOOTHR.LT.0.0) THEN
            CALL MSGDOC(MDOC,' SMOOTH_new = Badd_new')
            SMOOTH = BADD
          ENDIF
        ENDIF    
        IF(DIST_R.LT.0) THEN
          DIST_LIM = RESEFF
          CALL MSGDOC(MDOC,' Now DIST_LIM = Effective Resolution')
        ENDIF
      ENDIF
C --------------------------------------------------
C --------------------------------------------------
      ISIM   = 0
      SCALEF = 1.0
      LAP    = 'N'
      FOMMIN = 0.0

      IF(MTN.NE.0) THEN
c        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
C        CALL MSGDOC(MDOC,'--- coef for Diff Patterson ---')
c        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        MSGA  = 'N'
        FILE1 = NAMEF
        FILE2 = NAMEFS
      ELSE
c        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
C        CALL MSGDOC(MDOC,'--- coef for Anom. Diff Patterson ---')
c        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        MSGA  = 'Y'
        FILE1 = NAMED
        FILE2 = ' '
      ENDIF
      NAMES  = ' '
C ---

      BADDD    = BADD_COEF
      SLIM     = SLIM_PATT
      DRMS     = 0.0
      DRMS_A   = 0.0
C
      IF(COMB.NE.'N') MSGA = 'Y'

      DEL_LIM  = 0.0
      DELA_LIM = 0.0
      BADDD    = 0.0
      SLIM     = 0.0
      DRMS     = 0.0
      DRMS_A   = 0.0

      CALL COEF_FFT(M,FILE1,FILE2,NAMES,NAMES,NAMES
     * ,FOMMIN,RESMIN,RESMAX,BOFF,SLIM
     * ,DEL_LIM,DELA_LIM,DRMS,DRMS_A,LAP
     * ,SCALEF,BADDD,MSGA,ISIM,ISYM,IPRSYM,RESULT,IERR)

      FP_MAX_C    = RESULT( 1,9)
      FD_MAX_C    = RESULT( 2,9)
      DEL_MAX_C   = RESULT( 3,9)
      DEL_AVER_C  = RESULT( 4,9)
      DEL_RMS_C   = RESULT( 5,9)
      NTOT_DEL_C  = RESULT( 6,9)
      DELA_MAX_C  = RESULT( 7,9)
      DELA_AVER_C = RESULT( 8,9)
      DELA_RMS_C  = RESULT( 9,9)
      NTOT_DELA_C = RESULT(10,9)
      NUMB_C      = RESULT(11,9)
      NNOUT_C     = RESULT(12,9)
      NBSLIM_C    = RESULT(13,9)
      NBDLIM_C    = RESULT(14,9) 
      NBDLIMA_C   = RESULT(14,9) 
      NBFOM_C     = RESULT(15,9) 

      DEL_LIM = 0.0
      IF(TDEL_LIM.GT.0.0001) THEN
        DEL_LIM = TDEL_LIM*DEL_RMS_C + DEL_AVER_C 
      ENDIF  
      WRITE(LINE,'(''DEL_LIM      :'',F10.2)') 
     * DEL_LIM
      CALL MSGDOC(MD,LINE)
      DELA_LIM = 0.0
      IF(TDEL_LIM.GT.0.0001) THEN
        DELA_LIM = TDEL_LIM*DELA_RMS_C + DELA_AVER_C 
      ENDIF  
      WRITE(LINE,'(''DEL_LIM_anom :'',F10.2)') 
     * DELA_LIM
      CALL MSGDOC(MD,LINE)
C ------------------------------------------------------

C     input MSGA = 'Y': output MSGA = 'N' if no anom.signal

      IF(COMB.NE.'N'.AND.MSGA.NE.'N') ANOM = 'Y'

C ==================================================
C     calc coefficients for Patterson
C
      ISIM   = 0
      SCALEF = 1.0
      LAP    = 'N'
      FOMMIN = 0.0

C     input MSGA = 'Y': output MSGA = 'N' if no anom.signal

      IF(COMB.NE.'N'.AND.MSGA.NE.'N') ANOM = 'Y'
 
C      write(*,*) ' -->msga,comb,anom:',msga,comb,anom
     
      IF(MTN.NE.0) THEN
        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDOC,'--- coef for Diff Patterson ---')
        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')         
        MSGA  = 'N'
        IF(ANOM.EQ.'Y') THEN
          MSGA = 'Y'
          IF(COMB.EQ.'Y') MSGA = 'C'
        ENDIF
        FILE1 = NAMEF
        FILE2 = NAMEFS
      ELSE
        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDOC,'--- coef for Anom. Diff Patterson ---')
        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        MSGA  = 'Y'
        IF(COMB.EQ.'Y') MSGA = 'C'
        FILE1 = NAMED
        FILE2 = ' '
      ENDIF
C ---
      BADDD   = BADD_COEF
      SLIM    = SLIM_PATT
      DRMS    = 0.0
      DRMS_A  = 0.0
      NAMES   = ' '
C--   NAMEC   = 'trahalo_cff'

C      write(*,*) ' msga,anom:',msga,anom

      CALL COEF_FFT(M,FILE1,FILE2,NAMES,NAMEC,NAMECA
     *         ,FOMMIN,RESMIN,RESMAX,BOFF,SLIM
     *         ,DEL_LIM,DELA_LIM,DRMS,DRMS_A,LAP
     *         ,SCALEF,BADDD,MSGA,ISIM,ISYM,IPRSYM,RESULT,IERR)
      IF(IERR.NE.0) GO TO 900

      FP_MAX_C    = RESULT( 1,9)
      FD_MAX_C    = RESULT( 2,9)
      DEL_MAX_C   = RESULT( 3,9)
      DEL_AVER_C  = RESULT( 4,9)
      DEL_RMS_C   = RESULT( 5,9)
      NTOT_DEL_C  = RESULT( 6,9)
      DELA_MAX_C  = RESULT( 7,9)
      DELA_AVER_C = RESULT( 8,9)
      DELA_RMS_C  = RESULT( 9,9)
      NTOT_DELA_C = RESULT(10,9)
      NUMB_C      = RESULT(11,9)
      NNOUT_C     = RESULT(12,9)
      NBSLIM_C    = RESULT(13,9)
      NBDLIM_C    = RESULT(14,9) 
      NBDLIMA_C   = RESULT(15,9) 
      NBFOM_C     = RESULT(16,9) 

      WRITE(LINE,'(''N_coef,N_res_out,N_bad_sig     :'',3I6)') 
     * NUMB_C,NNOUT_C,NBSLIM_C
      CALL MSGDOC(MD,LINE)

      WRITE(LINE,'(''Fnat_max,Fder_max (or F+ F-)   :'',2F10.2)') 
     * FP_MAX_C,FD_MAX_C
      CALL MSGDOC(MD,LINE)

      WRITE(LINE,'(''DEL     : N,N_big,max,aver,rms :'',2I6,3F10.2)') 
     * NTOT_DEL_C,NBDLIM_C,DEL_MAX_C,DEL_AVER_C,DEL_RMS_C
      CALL MSGDOC(MD,LINE)

      WRITE(LINE,'(''DEL_anom: N,N_big,max,aver,rms :'',2I6,3F10.2)') 
     * NTOT_DELA_C,NBDLIMA_C,DELA_MAX_C,DELA_AVER_C,DELA_RMS_C
      CALL MSGDOC(MD,LINE)

C     RESULT(I,1) = RESOL
C     RESULT(I,2) = DEL = !FOBS -FDER!
C     RESULT(I,3) = FOBS 
C     RESULT(I,4) = RFAC
C     RESULT(I,5) = N
C     RESULT(I,6) = DEL_ano 
C     RESULT(I,7) = Nano
C     RESULT(I,8) = F/S or D/S
C     RESULT(1,9) = FP_MAX
C     RESULT(2,9) = FD_MAX
C     RESULT(3,9) = DEL_MAX

      CALL MSGDOC(MD,' ')
      WRITE(LINE,1700)  (RESULT(I,1),I=1,9)
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,1701)  (RESULT(I,2),I=1,9)
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,1703)  (RESULT(I,6),I=1,9)
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,1702)  (RESULT(I,4),I=1,9)
      CALL MSGDOC(MD,LINE)
      CALL MSGDOC(MD,' ')
      WRITE(LINE,1700)  (RESULT(I,1),I=10,18)
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,1701)  (RESULT(I,2),I=10,18)
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,1703)  (RESULT(I,6),I=10,18)
      CALL MSGDOC(MD,LINE)
      WRITE(LINE,1702)  (RESULT(I,4),I=10,18)
      CALL MSGDOC(MD,LINE)
      CALL MSGDOC(MD,' ')
 1700 FORMAT('Res :',9F8.2)
 1701 FORMAT('DEL :',3X,9(1X,F7.3))
 1703 FORMAT('Dano:',3X,9(1X,F7.3))
 1702 FORMAT('Rfac:',3X,9(1X,F7.3))


C      CALL COEF(M,FILE1,FILE2,NAMES,NAMEC
C     *         ,FOMMIN,RESMIN,RESMAX,BOFF,SLIM,LAP
C     *         ,SCALEF,BADDD,MSGA,ISIM,ISYM,IPRSYM,IERR)
C      IF(IERR.NE.0) GO TO 900

C     RESULT(I,1) = RESOL
C     RESULT(I,2) = DEL = !FOBS -FDER!
C     RESULT(I,3) = FOBS 
C     RESULT(I,4) = RFAC
C     RESULT(I,5) = N
C     RESULT(I,6) = DEL_ano 
C     RESULT(I,7) = Nano

      IF(COMB.EQ.'Y') THEN
        IF(MSGA.NE.'N') THEN
          CALL MSGDOC(MDOC,'    Anom.data of derivative was used')
        ENDIF
      ENDIF
      IF(ANOM.EQ.'Y'.AND.MSGA.EQ.'N') THEN
        CALL MSGERR(MDOC,
     *  ' ERROR: No anom.signal in the derivative')
        IERR=1
        RETURN
      ENDIF

C =======================================================
C ===============================================================

      NPEAKS_PATT = 0
      NKNOWN      = 0

      IF(FOUR.EQ.'Y'.AND.SUBR.EQ.'Y'.AND.NO_DER.GT.0) THEN
C
C --- restore self-list
C
        DO I=1,IPRMAXDD
        DO J=1,IPRROWW
          RESULTD(I,J) = RES_DER(I,J,NO_DER)  
        ENDDO
        ENDDO
        NPEAKS      = RESULTD(1,12)
        NPEAKS_PATT = RESULTD(1,12)
        NX          = NXMIN
        NY          = NYMIN
        NZ          = NZMIN
C ========================
C
C --- check known HA file      
C      
        CALL LENSTR_BL(NAMEHA,LEN)
        IF(LEN.GT.0.AND.NAMEHA(1:1).NE.','.AND.
     *    NAMEHA(1:1).NE.' ') THEN
          IF(SUBR.EQ.'N') THEN
            CALL LENSTR_BL(NAMEHA,LEN)
            IF(LEN.GT.56) LEN=56
            LINE='Input file known HA:'//NAMEHA(1:LEN)
            CALL MSGDOC(MD,LINE)
          ENDIF
C
C         add known atoms in the top of self-list
C         maximum 9
C
          NATOMR = 0
          NAMES  = ' '
          PDB    = 'C'
          CALL CHECK_CRD(M,NAMEHA,ISYM,IPRSYM
     *    ,STR_TITLE,STR_DATE,NAMES,NATOMR
     *    ,RESULTB,IPRMAXD,IPRROW,PDB,IERR)
          NHATOM = RESULTB(11,12)
          IF(NHATOM.GT.0) THEN
            IF(NHATOM.GT.9) NHATOM=9
            DO I=1,NHATOM
              IND           = I+10
              RESULTB(I,1 ) = RESULTB(IND,9 )*NX
              RESULTB(I,2 ) = RESULTB(IND,10)*NY
              RESULTB(I,3 ) = RESULTB(IND,11)*NZ
              RESULTB(I,4 ) = RESULTB(IND,4 )
              RESULTB(I,5 ) = RESULTB(IND,5 )
              RESULTB(I,6 ) = RESULTB(IND,6 )
              RESULTB(I,7 ) = RESULTB(IND,7 )
              RESULTB(I,8 ) = RESULTB(IND,8 )
              RESULTB(I,9 ) = RESULTB(IND,9 )
              RESULTB(I,10) = RESULTB(IND,10)
              RESULTB(I,11) = RESULTB(IND,11)
              RESULTB(I,13) = 0.0
              RESULTB(I,15) = I
              RESULTB(I,16) = RESULTB(IND,4 )

              RESULTB(1,12) = NHATOM

            ENDDO
      CALL MSGDOC(MDOC,'====================')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(MDOC,'---   List of known heavy atoms   ---')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
            NPEAK  = RESULTB(1,12)
            IPRINT = 1
            CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *      ,RESULTB,IPRMAXD,IPRROW,IERR)
      CALL MSGDOC(MDOC,'====================')
      
            DMAX = 0.0
            DO I=1,NPEAKS
              IF(RESULTD(I,4).GT.DMAX) DMAX = RESULTD(I,4)
            ENDDO

            NP = NPEAKS+NHATOM
            IF(NP.GE.IPRMAXDD) NP = IPRMAXDD-1

            DO I=NP,1,-1
              IF(I.LE.NHATOM) THEN
                DO J=1,IPRROW
                  RESULTD(I,J) = RESULTB(I,J)  
                ENDDO
                RESULTD(I,4) = RESULTB(I,4)*DMAX
              ELSE
                II = I-NHATOM
                DO J=1,IPRROW
                  RESULTD(I,J) = RESULTD(II,J)  
                ENDDO
              ENDIF
            ENDDO
            NPEAKS        = NP
            RESULTD(1,12) = NP
            NPEAKS_PATT   = RESULTD(1,12)

      CALL MSGDOC(MDD,'====================')
      CALL MSGDOC(MDD ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(MDD,'---   Now Self_patt List ---')
            NPEAK  = RESULTD(1,12)
            IPRINT = 1
      CALL PR_RES_PS(MDD,NPEAK,IPRINT
     *      ,RESULTD,IPRMAXDD,IPRROW,IERR)
      CALL MSGDOC(MDD,'====================')

          ENDIF
C
C         end to add known atoms 
C
        ENDIF
C ========================

        IF(NPEAKS.LE.0) THEN
          IERR=20
          CALL MSGERR(MDOC,' ERROR: input number of self_peaks = 0')
          GO TO 900
        ENDIF
C
C ---   go to calc diff fourie
C
        GO TO 700
      ENDIF
C
      IF(NO_DER.EQ.-2) THEN
        NX = NXMIN
        NY = NYMIN
        NZ = NZMIN
C ---   go to calc diff fourie
        GO TO 700
      ENDIF

      IF(NO_DER.LT.0) NO_DER = 0
     
C -------------------------
      IF(MTN.NE.0) THEN
        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDOC,'--- trpack: search in Diff. Patterson ---')
        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
      ELSE
        CALL MSGDOC(MD  
     *  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
        CALL MSGDOC(MDOC
     *  ,'--- trpack: search in Anom.Diff. Patterson ---')
        CALL MSGDOC(MD  
     *  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      ENDIF
      IOP1 = 0
      IOP2 = 0
      MOD2 = 'N'
      MODE = 'M'
      FUNC = 'Y'

      IF(TRANS.EQ.'N') MODE='T'

      IF(PACK.EQ.'N') FUNC='N'
       
      ACCUM = 'M'

      INVER    = ' '
      RAD1     = RAD
      RAD2     = RAD
      PERC2(1) = 0.0
      NX       = NXMIN
      NY       = NYMIN      
      NZ       = NZMIN

      IZMIN    = 0
C     IZMAX    = NZMIN-1
      IZMAX    = NZ/2+1
      IF(IFZ.EQ.0) IZMAX = 0 
      NHATOM   = 0
C -----------------

      MDDD = MD

      CALL LENSTR_BL(NAMEHA,LEN)
      IF(LEN.GT.0.AND.NAMEHA(1:1).NE.','.AND.NAMEHA(1:1).NE.' ') THEN

        IF(SUBR.EQ.'N') THEN
          CALL LENSTR_BL(NAMEHA,LEN)
          IF(LEN.GT.56) LEN=60
          LINE='Input file known HA:'//NAMEHA(1:LEN)
          CALL MSGDOC(MD,LINE)
        ENDIF
C
C       fixed model is all known atoms
C       maximum 9
C
        NATOMR = 0
        NAMES  = ' '
        PDB    = 'C'
        CALL CHECK_CRD(M,NAMEHA,ISYM,IPRSYM
     *  ,STR_TITLE,STR_DATE,NAMES,NATOMR
     *  ,RESULTB,IPRMAXD,IPRROW,PDB,IERR)
        IF(IERR.NE.0) GO TO 900
        NHATOM = RESULTB(11,12)
        IF(NHATOM.GT.0) THEN
          IF(NHATOM.GT.9) NHATOM = 9

          PERC2(1) = 50.0

          IF(NHATOM.GT.1) THEN
            WRITE(MOD2,'(I1)') NHATOM
          ELSE
            MOD2 = 'Y'
          ENDIF
          DO I=1,NHATOM
            IND=I+10
            RESULTB(I,1 ) = RESULTB(IND,9 )*NX
            RESULTB(I,2 ) = RESULTB(IND,10)*NY
            RESULTB(I,3 ) = RESULTB(IND,11)*NZ
            RESULTB(I,4 ) = RESULTB(IND,4 )
            RESULTB(I,5 ) = RESULTB(IND,5 )
            RESULTB(I,6 ) = RESULTB(IND,6 )
            RESULTB(I,7 ) = RESULTB(IND,7 )
            RESULTB(I,8 ) = RESULTB(IND,8 )
            RESULTB(I,9 ) = RESULTB(IND,9 )
            RESULTB(I,10) = RESULTB(IND,10)
            RESULTB(I,11) = RESULTB(IND,11)
            RESULTB(I,13) = 0.0
            RESULTB(I,15) = I
            RESULTB(I,16) = RESULTB(IND,4 )

            RESULTB(1,12) = NHATOM

          ENDDO
      CALL MSGDOC(MDOC,'====================')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(MDOC,'---   List of known heavy atoms   ---')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
          NPEAK  = RESULTB(1,12)
          IPRINT = 1
          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *    ,RESULTB,IPRMAXD,IPRROW,IERR)
      CALL MSGDOC(MDOC,'====================')

          SH1(1) = RESULTB(1,9 )
          SH1(2) = RESULTB(1,10)
          SH1(3) = RESULTB(1,11)
          IF(NHATOM.GT.1) THEN
            DO I=2,NHATOM
              II        = I-1
              SH2(1,II) = RESULTB(I,9 )
              SH2(2,II) = RESULTB(I,10)
              SH2(3,II) = RESULTB(I,11)
            ENDDO
          ENDIF
          IZMIN = 0
          IZMAX = NZMIN-1

        ELSE

          SH1(1) = 0.0
          SH1(2) = 0.0
          SH1(3) = 0.0

          SH2(1,1) = 0.0
          SH2(2,1) = 0.0
          SH2(3,1) = 0.0
C
C         for P1 fixed model is one atom (0,0,0)
C
          IF(NSYM.EQ.1) THEN
            MOD2  = 'Y'
            IZMIN = 0
            IZMAX = NZMIN-1
          ENDIF

        ENDIF
        MDDD = MD

      ELSE
C
C       without fixed model
C
        SH1(1)   = 0.0
        SH1(2)   = 0.0
        SH1(3)   = 0.0

        SH2(1,1) = 0.0
        SH2(2,1) = 0.0
        SH2(3,1) = 0.0
C
C       for P1 fixed model is one atom (0,0,0)
C 
        IF(NSYM.EQ.1) THEN
          MOD2  = 'Y'
          IZMIN = 0
          IZMAX = NZMIN-1
        ENDIF
      ENDIF


C ???
c      MODE_CC = 'N'
c      IFIRST  = 1
c      NAMES   = ' '
c      NA_REF  = 1
c        CALL TEST_REF(MDOC,MODE_CC,IFIRST   
c     *          ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
c     *          ,RESMIN,RESMAX,BOFF,BADD
c     *          ,RESULTC,IPRMAXD,RESULT,IPRMAX,IPRROW
c     *          ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
c        IF(IERR.NE.0) GO TO 900
C ???




C -------------------------------------

      NAME(1)  = ' '
      NAME(2)  = ' '
      NAME(3)  = ' '
      NAME(4)  = ' '
C                NAMEC - file of coeff of Diff. Patt.
C                NAMEC = 'trahalo_cff'
      NAME(5)  = NAMEC
      NAME(6)  = ' '
      NAME(7)  = NAMET
      NAME(7)  = 'patt_pack'
      NAME(8)  = 'trahalo.scr'
      NAME(9)  = NAME_SLF
      NAME(10) = ' '

      BADDD = SMOOTH
C ---
      CALL TRPACK(M,NAME,MODE,MOD2,FUNC,REMOV,INVER
     * ,ACCUM,NX,NY,NZ,IZMIN,IZMAX,IOP1,IOP2,BADDD,BOFF
     * ,RAD1,RAD2,SH2,SH1,GAUSS,PERC2,RESMIN,RESMAX
     * ,POOL,MEMORY,ISYM,IPRSYM,IERR)
      IF(IERR.EQ.10) THEN
        CALL MSGERR(MDOC,
     *  ' ERROR: bad Packing Function')
      ENDIF
      IF(IERR.NE.0) GO TO 900

C ---
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(MDOC,'---  peaksrch for pattsrch function  ---')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      ALEVEL = 0.0

      NPMAX = 1
      CALL PEAKSRCH(M,NAME_SLF,NAMEO,ALEVEL,DIST_LIM,SPEC,NPMAX
     *     ,POOL,MEMORY,ISYM,IPRSYM,RESULTD,IPRMAXDD,IPRROW,IERR)
      IF(IERR.NE.0) GO TO 900
C                                
C  RESULT(i,1) IX 
C  RESULT(i,2) IY
C  RESULT(i,3) IZ
C  RESULT(i,4) Dens
C  RESULT(i,5) sigma
C  RESULT(i,6) alpha      Xort
C  RESULT(i,7) beta   or  Yort
C  RESULT(i,8) gamma      Zort
C  RESULT(i,9)  Xfrac
C  RESULT(i,10) Yfrac
C  RESULT(i,11) Zfrac
C
C  RESULT(1,12)  number of peaks
C
      NPEAKS = RESULTD(1,12)

      IF(NPEAKS.GE.IPRMAXDD) NPEAKS = NPEAKS-1
      RESULTD(1,12)        = NPEAKS
      RESULTD(NPEAKS+1,4)  = 0.0
      RESULTD(NPEAKS+1,13) = 0.0
      RESULTD(NPEAKS+1,16) = 0.0

      IF(NPEAKS.LE.0) THEN
        CALL MSGERR(MDOC,' ERR: number of peaks = 0')
        IERR = 1
        GO TO 900
      ENDIF      

      IF(NSYM.EQ.1.AND.NHATOM.LE.0) THEN
C
C      put atom (0,0,0) to self-list for P1
C  
       IF(NPEAKS+1.GE.IPRMAXDD) NPEAKS = NPEAKS-1

        DO IP=NPEAKS-1,1,-1
          DO I=1,IPRROW
            RESULTD(IP+1,I) = RESULTD(IP,I)
          ENDDO
        ENDDO

        DO J=1,IPRROW
          RESULTD(1,J) = 0.0  
        ENDDO
C
C       set dens. for 0,0,0 as second peak
C
        RESULTD(1,4)  = RESULTD(2,4)
        RESULTD(1,16) = RESULTD(2,16)

        NPEAKS        = NPEAKS+1
        RESULTD(1,12) = NPEAKS

        NP = NP + 1
        IF(NP.GE.IPRMAXDD) NP = IPRMAXDD-1
C
C       now number of peaks to check NP = NP_input + 1
C
      ENDIF
      
      IF(NHATOM.GT.0) THEN
C
C       put known atoms to self-list
C       with dens(i) = OCC(i) * DENS_max 
C 
        DMAX = 0.0
        DO I=1,NPEAKS
          IF(RESULTD(I,4).GT.DMAX) DMAX = RESULTD(I,4)
        ENDDO

        NPP = NPEAKS+NHATOM
        IF(NPP.GE.IPRMAXDD) NPP = IPRMAXDD-1

        DO I=NPP,1,-1
          IF(I.LE.NHATOM) THEN
            DO J=1,IPRROW
              RESULTD(I,J) = RESULTB(I,J)  
            ENDDO
            RESULTD(I,4)  = RESULTB(I,4)*DMAX
            RESULTD(I,16) = RESULTB(I,16)
          ELSE
            II = I-NHATOM
            DO J=1,IPRROW
              RESULTD(I,J) = RESULTD(II,J)  
            ENDDO
          ENDIF
        ENDDO

        NPEAKS        = NPP
        RESULTD(1,12) = NPP

        MDDD          = MD
       
        NP = NP+NHATOM
        IF(NP.GE.IPRMAXDD) NP = IPRMAXDD-1
C
C       now number of peaks to check NP = NP_input + N_known
C
      ENDIF

      CALL MSGDOC(M ,'====================')
      CALL MSGDOC(M ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(M ,'--- Self_patt List before refinement  ---')
      NPEAK  = RESULTD(1,12)
      IPRINT = 1


      CALL PR_RES_PS(M,NPEAK,IPRINT
     *    ,RESULTD,IPRMAXDD,IPRROW,IERR)
      CALL MSGDOC(MDD,'====================')

      II  = 1
      IST = 2      
      IF(NHATOM.GT.0) THEN
        IST = NHATOM+1
        II  = NHATOM
      ELSE IF(NSYM.EQ.1) THEN
        IST = 3
        II  = 2
      ENDIF
C
C     remove equal peaks
C
      DO I=IST,NPEAKS

        IX1 = RESULTD(I,1)
        IY1 = RESULTD(I,2)
        IZ1 = RESULTD(I,3)

        DO IC2=1,I-1

          IX2 = RESULTD(IC2,1)
          IY2 = RESULTD(IC2,2)
          IZ2 = RESULTD(IC2,3)

          CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
     *    ,INV,IOR,ISOP,DIST
     *    ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *    ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)

          IF(DIST.LT.DIST_LIM) THEN

                  WRITE(LINE,
     *       '('' Reject peak :'',I5,'' because = :'',I5)') 
     *            I,IC2
                  CALL MSGDOC(M,LINE)

            GO TO 440
          ENDIF
        ENDDO
        II = II+1
        DO J=1,IPRROW
          RESULTD(II,J) = RESULTD(I,J)  
        ENDDO
 440    CONTINUE
      ENDDO
            
      RESULTD(1,12) = II
      NPEAKS        = RESULTD(1,12)

      IF(NPEAKS.LE.0) THEN
        CALL MSGERR(MDOC,' ERROR: number of peaks = 0')
        IERR=1
        GO TO 900
      ENDIF

C
C      CALL MSGDOC(MDD,'--- refinement ---')
C     
      II     = 0
      NKNOWN = 0
      IF(NHATOM.GT.0) THEN
        NKNOWN = NHATOM
      ELSE IF(NSYM.EQ.1) THEN
        NKNOWN = 1
      ENDIF

      IF(NKNOWN.GT.0) THEN
        DO IPN=1,NKNOWN
          II = II+1
          DO J=1,IPRROW
            RESULTC(II,J) = RESULTD(IPN,J)
          ENDDO
C    ---  set the same weight --
          RESULTC(II,4) = 1.0
        ENDDO

        IFIRST        = 1
        NA_REF        = II
        RESULTC(1,12) = NA_REF
        NAMES         = ' '
        MODE_CC = 'N'
        CALL CHECK_AT_REF(M,MODE_CC,IFIRST   
     *          ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *          ,RESMIN,RESMAX,BOFF,BADD
     *          ,RESULTC,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *          ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) GO TO 900

        DO I=1,NKNOWN
          RESULTD(I,13) = RESULTC(1,13)
        ENDDO

        DO I=1,NKNOWN
          RESULTD(I,16) = RESULTC(I,16)
        ENDDO

      ENDIF
      
      NSTART = NKNOWN+1
      IFIRST = 0

C            RESULTD(1,9 ) = .0791
C            RESULTD(1,10) = .63438
C            RESULTD(1,11) = .17299
C            RESULTD(1,1 ) = RESULTD(1,9 )*NX
C            RESULTD(1,2 ) = RESULTD(1,10)*NY
C            RESULTD(1,3 ) = RESULTD(1,11)*NZ
C            RESULTD(1,6 ) = RESULTD(1,9 )*52.35
C            RESULTD(1,7 ) = RESULTD(1,10)*57.73
C            RESULTD(1,8 ) = RESULTD(1,11)*82.57
C            RESULTD(1,13) = 0.0
C            RESULTD(1,15) = 1
C            RESULTD(1,16) = RESULTD(1,4)
C
C            RESULTD(2,9 ) = .63767
C            RESULTD(2,10) = .36383
C            RESULTD(2,11) = .16412
C            RESULTD(2,1 ) = RESULTD(2,9 )*NX
C            RESULTD(2,2 ) = RESULTD(2,10)*NY
C            RESULTD(2,3 ) = RESULTD(2,11)*NZ
C            RESULTD(2,6 ) = RESULTD(2,9 )*52.35
C            RESULTD(2,7 ) = RESULTD(2,10)*57.73
C            RESULTD(2,8 ) = RESULTD(2,11)*82.57
C            RESULTD(2,13) = 0.0
C            RESULTD(2,15) = 1
C            RESULTD(2,16) = RESULTD(2,4)


      DO IPK=NSTART,NPEAKS
        II = 0
        IF(NKNOWN.GT.0) THEN
          DO IPN=1,NKNOWN
            II = II+1
            DO J=1,IPRROW
              RESULTC(II,J) = RESULTD(IPN,J)
            ENDDO
C      ---  set the same weight --
            RESULTC(II,4) = 1.0
          ENDDO
        ENDIF

        NA_REF = II+1
        II     = II+1
        DO J=1,IPRROW
          RESULTC(II,J) = RESULTD(IPK,J)
        ENDDO
C   --- set the same weight --
        RESULTC(II,4) = 1.0
        RESULTC(1,12) = NA_REF
        NAMES         = ' '
        MODE_CC = 'N'
        IFIRST  = IFIRST + 1
        CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *          ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *          ,RESMIN,RESMAX,BOFF,BADD
     *          ,RESULTC,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *          ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) GO TO 900

        RESULTD(IPK, 6) = RESULTC(1, 6)
        RESULTD(IPK, 7) = RESULTC(1, 7)
        RESULTD(IPK, 8) = RESULTC(1, 8)
        RESULTD(IPK, 9) = RESULTC(1, 9)
        RESULTD(IPK,10) = RESULTC(1,10)
        RESULTD(IPK,11) = RESULTC(1,11)
        RESULTD(IPK, 1) = RESULTD(IPK, 9)*NX
        RESULTD(IPK, 2) = RESULTD(IPK,10)*NY
        RESULTD(IPK, 3) = RESULTD(IPK,11)*NZ
        RESULTD(IPK,13) = RESULTC(1,13)
        RESULTD(IPK,16) = RESULTC(1,16)

      ENDDO


C         WRITE(LINE,
C     *       '('' P1,P2 :'',2F16.5)') 
C     *   RESULTD(1,13), RESULTD(2,13)          
C                  CALL MSGDOC(MDOC,LINE)
C         RESULTD(1,13) = 4.1
C         RESULTD(2,13) = 4.0

C --- SORT RESULTD
      IST = 1+NKNOWN
      IF(LSORT.EQ.'P') THEN
        CALL SORT_PEAK(MDOC,IST,RESULTD,IPRMAXDD,IPRROW,IERR)
      ELSE IF(LSORT.EQ.'M') THEN
        CALL SORT_PEAKM(MDOC,IST,RESULTD,IPRMAXDD,IPRROW,IERR)
      ENDIF
C ----------
      NPEAKS = RESULTD(1,12)
      IF(NPEAKS.GT.0) THEN
        DO I=1,NPEAKS
          RESULTD(I,15) = I
        ENDDO
      ENDIF
C ---
      CALL MSGDOC(MD,'====================')
      CALL MSGDOC(MD ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(MD,'---   List of expected heavy atoms   ---')
      CALL MSGDOC(MD ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(MD,'---   Self_patt List ---')
      IF(NKNOWN.LE.0) THEN
        CALL MSGDOC(MD,
     *   ' Power for current peak after unphased refinement ')
      ELSE
        CALL MSGDOC(MD,
     *  ' Power for current peak plus known atoms ')
        CALL MSGDOC(MD,
     *  ' after unphased refinement with the same weight')
      ENDIF
      NPEAK  = RESULTD(1,12)
      IPRINT = 1
      CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *    ,RESULTD,IPRMAXDD,IPRROW,IERR)
      CALL MSGDOC(MDD,'====================')

      NPEAKS      = RESULTD(1,12)
      NPEAKS_PATT = RESULTD(1,12)

      IF(NPEAKS.LE.0) THEN
        IERR=13
        RETURN
      ENDIF

      IF(LSORT.EQ.'M') LSORT='P' 
C     IF(NP.LE.0)     GO TO 900
C
      IF(FOUR.EQ.'Y') GO TO 700
C
C ========== if four ---> ==================
C
      IF(NP.LE.0.AND.NKNOWN.LE.0) THEN
C ---   case when input NP < 0
        IC     = NP   
        NPEAKS = 1  
        RESULTD(1,12) = 1
        GO TO 310
C ---   go to the final analysis
      ENDIF
C ----
      IF(NO_DER.GT.0) THEN
C
C       store self-list
C
        DO I=1,IPRMAXDD
        DO J=1,IPRROWW
          RES_DER(I,J,NO_DER) = RESULTD(I,J)  
        ENDDO
        ENDDO
      ENDIF
C ----
      IF(NKNOWN.GT.0) THEN
        NATOMS = NKNOWN
        TEST   = RESULTD(NKNOWN,13)*1.01
        IF(RESULTD(NKNOWN+1,13).GE.TEST) THEN 
          NATOMS = NATOMS+1
        ENDIF

        RESULTD(1,12) = NATOMS
        NA_REF        = NATOMS

        DO IPN=1,NATOMS

C    ---  set the same weight --
c          RESULTD(IPN,4)=1.0

          RESULTD(IPN,4)=RESULTD(IPN,16)
        ENDDO

        MODE_CC = 'N'
        NAMES   = ' '
        IFIRST  = 1
        CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *            ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *            ,RESMIN,RESMAX,BOFF,BADD
     *            ,RESULTD,IPRMAXDD,RESULT,RESULTP,IPRMAX,IPRROW
     *            ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) RETURN

C ---
        IF(NATOMS.GT.0) THEN
          CALL MSGDOC(MDOC,'====================')
          CALL MSGDOC(MDOC,'---  List of expected heavy atoms  ---')
          NPEAK  = NATOMS   
          NPMAX  = NATOMS
          IPRINT = 1         
          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *      ,RESULTD,IPRMAXDD,IPRROW,IERR)
          CALL MSGDOC(MDOC,'====================')

          DO IPN=1,NATOMS      
            RESULTD(IPN,4) = RESULTD(IPN,16)
          ENDDO

          CALL WR_RES_PS(M,NAMEO,NPMAX,RESULTD,IPRMAXDD,IPRROW
     *    ,NSYM,ISYM,IPRSYM,IERR)
          N_FINAL = NATOMS 

        ENDIF

        RETURN

      ENDIF
C =================================================================      
C  --- search with fixed model_2 ----


 710  CONTINUE
      NPEAKS = RESULTD(1,12)
C
C --------- patt. search --- 
C
      IF(NPEAKS.LE.1) GO TO 310
C --- go to the final analysis

        IF(NPEAKS.LT.NP) NP = NPEAKS
        IC  = 0
        IF(NPEAKS.GE.IPRMAXDD) NPEAKS = NPEAKS-1

C  --- cycle of search with fixed model ----
        CALL MSGDOC(MD
     * ,'*-*-*-*-**-*-*-*-*-*--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDOC,' ')
        CALL MSGDOC(MDOC
     * ,' --- Search with one peak as fixed model ---')
        CALL MSGDOC(MDOC,' ')
        CALL MSGDOC(MD
     * ,'*-*-*-*-**-*-*-*-*-*--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')

 200    CONTINUE

        IC = IC + 1
        CALL MSGDOC(MD
     * ,'*-*-*-*-*-*--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        WRITE(LINE,'('' --- peak '',I2,'' ('',3F8.3,
     *  '') is fixed model_2 ---'')') 
     *  IC,RESULTD(IC, 9),RESULTD(IC,10),RESULTD(IC,11)
        CALL MSGDOC(MDOC,LINE)

        IF(MTN.NE.0) THEN
          CALL MSGDOC(MDOC,'--- trpack: search in Diff. Patterson ---')
        ELSE
          CALL MSGDOC(MDOC
     * ,'--- trpack: search in Anom. Diff. Patterson ---')
        ENDIF
        CALL MSGDOC(MD
     * ,'*-*-*-*-*-*--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        IOP1 = 0
        IOP2 = 0
        MOD2 = 'Y'
        MODE = 'M'
        FUNC = 'Y'

        IF(TRANS.EQ.'N') MODE = 'T'

        IF(PACK.EQ.'N') FUNC = 'N'
C       IF(PACK.EQ.'P') MODE = 'P'

        ACCUM    = 'M'
        INVER    = ' '
        RAD1     = RAD
        RAD2     = RAD
        PERC2(1) = 50.0
        SH2(1,1) = 0.0
        SH2(2,1) = 0.0
        SH2(3,1) = 0.0
        SH1(1)   = RESULTD(IC, 9)
        SH1(2)   = RESULTD(IC,10)
        SH1(3)   = RESULTD(IC,11)
        NX       = NXMIN
        NY       = NYMIN
        NZ       = NZMIN
        IZMIN    = 0
        IZMAX    = NZMIN-1
        NAME(1)  = ' '
        NAME(2)  = ' '
        NAME(3)  = ' '
        NAME(4)  = ' '
        NAME(5)  = NAMEC
        NAME(6)  = ' '
        NAME(7)  = NAMET
        NAME(8)  = 'trahalo.scr'
        NAME(9)  = NAMEP
        NAME(10) = ' '

        BADDD    = SMOOTH
C ---
        CALL TRPACK(M,NAME,MODE,MOD2,FUNC,REMOV,INVER
     *   ,ACCUM,NX,NY,NZ,IZMIN,IZMAX,IOP1,IOP2,BADDD,BOFF
     *   ,RAD1,RAD2,SH2,SH1,GAUSS,PERC2,RESMIN,RESMAX
     *   ,POOL,MEMORY,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) THEN 
          IERR = 0
          GO TO 300
        ENDIF
C ---
        CALL MSGDOC(MDD,  '*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDD,'--- peaksrch for pattsrch function ---')
        CALL MSGDOC(MDD,  '*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        ALEVEL = 0.0
        NPMAX  = 0
        SPEC_P = SPEC
        IF(SPEC.EQ.'N') SPEC_P = 'P'
        CALL PEAKSRCH(M,NAMEP,NAMES,ALEVEL,DIST_LIM,SPEC_P,NPMAX
     *  ,POOL,MEMORY,ISYM,IPRSYM,RESULTA(1,1,IC),IPRMAXD,IPRROW,IERR)
        IF(IERR.NE.0) THEN
          IERR = 0
          GO TO 300
        ENDIF
        NPEAKSA = RESULTA(1,12,IC)

        IF(NPEAKSA.GT.30) THEN          
          NPEAKSA = 30
          RESULTA(1,12,IC) = NPEAKSA
        ENDIF

        IF(NPEAKSA.GT.0) THEN          
C ----    set number of peak in array RESULTD ---

          CALL MSGDOC(MDD,'--- analysis ---')

          DO IPK=1,NPEAKSA

            IX1 = RESULTA(IPK,1,IC)
            IY1 = RESULTA(IPK,2,IC)
            IZ1 = RESULTA(IPK,3,IC)

            IF(IFX.EQ.0) IX1 = 0
            IF(IFY.EQ.0) IY1 = 0
            IF(IFZ.EQ.0) IZ1 = 0

            DO ICD=1,NPEAKS

              IX2 = RESULTD(ICD,1)
              IY2 = RESULTD(ICD,2)
              IZ2 = RESULTD(ICD,3)

              IF(IFX.EQ.0) IX2 = 0
              IF(IFY.EQ.0) IY2 = 0
              IF(IFZ.EQ.0) IZ2 = 0

              CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
     *        ,INV,IOR,ISOP,DIST
     *        ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *        ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)
 
              IF(DIST.LT.DIST_LIM) THEN
                RESULTA(IPK,15,IC) = ICD
                GO TO 320
              ENDIF
            ENDDO  
            RESULTA(IPK,15,IC) = 0
 320        CONTINUE
          ENDDO
C -----
          CALL MSGDOC(MDD,'--- refinement ---')

          IFIRST = 0

          DO IPK=1,NPEAKSA

            DO J=1,IPRROW
              RESULTC(1,J) = RESULTD(IC,J)
              RESULTC(2,J) = RESULTA(IPK,J,IC)
            ENDDO
            RESULTC(1,12) = 2

C ---       set the same weight --
            RESULTC(2,4) = RESULTC(1,4)

            NA_REF  = 2
            NAMES   = ' '
            MODE_CC = 'N'
            IFIRST  = IFIRST + 1
            CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *          ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *          ,RESMIN,RESMAX,BOFF,BADD
     *          ,RESULTC,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *          ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
            IF(IERR.NE.0) RETURN

            RESULTA(IPK, 6,IC) = RESULTC(2, 6)
            RESULTA(IPK, 7,IC) = RESULTC(2, 7)
            RESULTA(IPK, 8,IC) = RESULTC(2, 8)
            RESULTA(IPK, 9,IC) = RESULTC(2, 9)
            RESULTA(IPK,10,IC) = RESULTC(2,10)
            RESULTA(IPK,11,IC) = RESULTC(2,11)
            RESULTA(IPK, 1,IC) = RESULTA(IPK, 9,IC)*NX
            RESULTA(IPK, 2,IC) = RESULTA(IPK,10,IC)*NY
            RESULTA(IPK, 3,IC) = RESULTA(IPK,11,IC)*NZ

            RESULTA(IPK,13,IC) = RESULTC(1,13)
            RESULTA(IPK,16,IC) = RESULTC(2,16)

          ENDDO
          IF(LSORT.EQ.'P') THEN
C ---       sort RESULTA(1,1,IC)
            IST=1
            CALL SORT_PEAK(MDOC,IST,RESULTA(1,1,IC)
     *      ,IPRMAXD,IPRROW,IERR)
          ENDIF
C --
          CALL MSGDOC(MDD,'====================')
          CALL MSGDOC(MDD,
     *'Power for pair current peak and fixed model peak after unphased r
     *efinement')
          CALL MSGDOC(MDD,
     *'with the same weight. No - peak number in Self_patt list')
          NPEAK  = RESULTA(1,12,IC)
          IPRINT = 1
          CALL PR_RES_PS(MDD,NPEAK,IPRINT
     *    ,RESULTA(1,1,IC),IPRMAXD,IPRROW,IERR)
          CALL MSGDOC(MDD,'====================')

        ENDIF
C ---
 300  CONTINUE
      IF(IC.LT.NP) GO TO 200

C -------------------------------------------------------------------
C ------ additional peak -----------

      P_MAX  = 0.0
      IP_MAX = 0      
C --- choice of additional peak ---
      DO IC=1,NP
        NPEAKSA = RESULTA(1,12,IC)

        D_MAX    = 0.0
        POWER_D  = 0.0
        ID_MAX   = 0
        IC_MAXI  = 0
        IPK_MAXI = 0
        PINIT    = RESULTD(IC,13)
        
        IF(NPEAKSA.GT.0) THEN  
          DO IPK=1,NPEAKSA

            IX = RESULTA(IPK,1,IC)
            IY = RESULTA(IPK,2,IC)
            IZ = RESULTA(IPK,3,IC)

            CALL COMP_SELF_PNT(MDOC,IX,IY,IZ,DIST_LIM,ISPEC
     *      ,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *      ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)
C
C      ISPEC=1 SPECIAL POSITION
C
            ID = RESULTA(IPK,15,IC)
            D  = RESULTA(IPK, 4,IC)
            P  = RESULTA(IPK,13,IC)
            IF(ID.GT.0.AND.ID.NE.IC.AND.ID.GT.NP) THEN
              IF(P.GT.POWER_D.AND.P.GT.PINIT) THEN
                D_MAX    = D
                POWER_D  = P
                ID_MAX   = ID
                IC_MAXI  = IC
                IPK_MAXI = IPK
              ENDIF
            ELSE IF(ID.NE.IC.AND.ID.EQ.0) THEN
              IF(P.GT.POWER_D.AND.P.GT.PINIT) THEN
                D_MAX    = D
                POWER_D  = P
                ID_MAX   = -1
                IC_MAXI  = IC
                IPK_MAXI = IPK
              ENDIF
            ENDIF
          ENDDO
        ENDIF

        IF(POWER_D.GT.P_MAX) THEN
          P_MAX   = POWER_D
          IP_MAX  = ID_MAX
          IPK_MAX = IPK_MAXI
          IC_MAX  = IC_MAXI
        ENDIF

      ENDDO
C ------
C --- checking of additional peak ---
C ///

      IF(IP_MAX.NE.0) THEN           
C ///
        IF(IP_MAX.GT.0) THEN           

          NP    = NP + 1
          NPNEW = IP_MAX
          DO J=1,IPRROW

            RESULTC(1,J)      = RESULTD(NP,J)
            RESULTD(NP,J)     = RESULTD(IP_MAX,J)
            RESULTD(IP_MAX,J) = RESULTC(1,J)

          ENDDO

          IF(NP.GT.NPEAKS) THEN
            NPEAKS        = NP
            RESULTD(1,12) = NPEAKS
          ENDIF

          WRITE(LINE,'('' Additional peak : ('',3F8.3,'')'')') 
     *    RESULTD(NP, 9),RESULTD(NP,10),RESULTD(NP,11)
     *    
          CALL MSGDOC(MDOC,LINE)
          WRITE(LINE,'('' old number of Additional peak :'',I4,
     *    '' new number :'',I4)') 
     *    NPNEW,NP
          CALL MSGDOC(MDOC,LINE)
          WRITE(LINE,'('' and change number: old:'',I4,'' --> new:''
     *    ,I4)')
     *    NP,NPNEW
          CALL MSGDOC(MDOC,LINE)
          
          IF(NO_DER.GT.0) THEN
            RES_DER(1,12,NO_DER) = NPEAKS
            DO J=1,IPRROWW
              RES_DER(NP,J,NO_DER) = RESULTD(NP,J)  
            ENDDO
          ENDIF

          RESULTD(NP   ,15) = NPNEW
          RESULTD(NPNEW,15) = NP

          DO IC=1,NP-1
            NPEAKSA = RESULTA(1,12,IC)
            IF(NPEAKSA.GT.0) THEN  
              DO IPK=1,NPEAKSA
                IP  = RESULTA(IPK,15,IC)
                IF(IP.EQ.NP) THEN
                  RESULTA(IPK,15,IC) = NPNEW
                ELSE IF(IP.EQ.NPNEW) THEN
                  RESULTA(IPK,15,IC) = NP
                ENDIF
              ENDDO  
            ENDIF
          ENDDO

        ELSE IF(IP_MAX.LT.0) THEN

          NP = NP + 1
          IF(NPNEW.GT.IPRMAXDD) NPNEW = IPRMAXDD
          NPNEW = NPEAKS+1
          IF(NPNEW.GT.IPRMAXDD) NPNEW = IPRMAXDD
          DO J=1,IPRROW
            RESULTC(1,J)     = RESULTD(NP,J)
            RESULTD(NPNEW,J) = RESULTC(1,J)
            RESULTD(NP,J)    = RESULTA(IPK_MAX,J,IC_MAX)
          ENDDO
          RESULTD(NP,4) = RESULTD(IC_MAX,4)*0.75
          NPEAKS        = NPNEW
          RESULTD(1,12) = NPEAKS

          WRITE(LINE,'('' Additional peak : ('',3F8.3,'')'')') 
     *    RESULTD(NP, 9),RESULTD(NP,10),RESULTD(NP,11) 
          CALL MSGDOC(MDOC,LINE)
          WRITE(LINE,'('' old number of Additional peak :'',I4,''(''
     *    ,I4,'') new number :'',I4)') IC_MAX,IPK_MAX,NP
          CALL MSGDOC(MDOC,LINE)
          WRITE(LINE,'('' and change number : old:'',I4,'' --> new:''
     *    ,I4)') NP,NPNEW
          CALL MSGDOC(MDOC,LINE)

          RESULTD(NP   ,15) = NPNEW
          RESULTD(NPNEW,15) = NP

          IF(NO_DER.GT.0) THEN
            RES_DER(1,12,NO_DER) = NPEAKS
            DO J=1,IPRROWW
              RES_DER(NP,J,NO_DER) = RESULTD(NP,J)  
            ENDDO
          ENDIF

          DO IC=1,NP-1
            NPEAKSA = RESULTA(1,12,IC)
            IF(NPEAKSA.GT.0) THEN  
              DO IPK=1,NPEAKSA

                IX1 = RESULTA(IPK,1,IC)
                IY1 = RESULTA(IPK,2,IC)
                IZ1 = RESULTA(IPK,3,IC)
                IP  = RESULTA(IPK,15,IC)

                IF(IFX.EQ.0) IX1 = 0
                IF(IFY.EQ.0) IY1 = 0
                IF(IFZ.EQ.0) IZ1 = 0

                IX2 = RESULTD(NP,1)
                IY2 = RESULTD(NP,2)
                IZ2 = RESULTD(NP,3)

                IF(IFX.EQ.0) IX2 = 0
                IF(IFY.EQ.0) IY2 = 0
                IF(IFZ.EQ.0) IZ2 = 0

                IF(IP.LE.0) THEN

                  CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
     *            ,INV,IOR,ISOP,DIST
     *            ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *            ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)

                  IF(DIST.LT.DIST_LIM) THEN
                    RESULTA(IPK,15,IC) = NP
                  ENDIF
 
                ELSE IF(IP.EQ.NP) THEN

                  RESULTA(IPK,15,IC) = NPNEW

                ENDIF

              ENDDO  
            ENDIF
          ENDDO
C          IP_MAX=NP
        ENDIF

        IC               = NP
        RESULTA(1,12,IC) = 0
        CALL MSGDOC(MD
     * ,'*-*-*-*-*-*--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        WRITE(LINE,'('' --- peak '',I2,'' is fixed model_2 ---'')') 
     *  IC
        CALL MSGDOC(MDOC,LINE)

        IF(MTN.NE.0) THEN
          CALL MSGDOC(MDOC,'--- trpack: search in Diff. Patterson ---')
        ELSE
          CALL MSGDOC(MDOC
     * ,'--- trpack: search in Anom. Diff. Patterson ---')
        ENDIF
        CALL MSGDOC(MD
     * ,'*-*-*-*-*-*--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        IOP1 = 0
        IOP2 = 0
        MOD2 = 'Y'
        MODE = 'M'
        FUNC = 'Y'

        IF(TRANS.EQ.'N') MODE = 'T'

        IF(PACK.EQ.'N') FUNC = 'N'
        ACCUM    = 'M'

        INVER    = ' '
        RAD1     = RAD
        RAD2     = RAD
        PERC2(1) = 50.0
        SH2(1,1) = 0.0
        SH2(2,1) = 0.0
        SH2(3,1) = 0.0
        SH1(1)   = RESULTD(IC, 9)
        SH1(2)   = RESULTD(IC,10)
        SH1(3)   = RESULTD(IC,11)
        NX       = NXMIN
        NY       = NYMIN
        NZ       = NZMIN
        IZMIN    = 0
        IZMAX    = NZMIN-1
        NAME(1)  = ' '
        NAME(2)  = ' '
        NAME(3)  = ' '
        NAME(4)  = ' '
        NAME(5)  = NAMEC
        NAME(6)  = ' '
        NAME(7)  = NAMET
        NAME(8)  = 'trahalo.scr'
        NAME(9)  = NAMEP
        NAME(10) = ' '

        BADDD    = SMOOTH
C ---
        CALL TRPACK(M,NAME,MODE,MOD2,FUNC,REMOV,INVER
     *   ,ACCUM,NX,NY,NZ,IZMIN,IZMAX,IOP1,IOP2,BADDD,BOFF
     *   ,RAD1,RAD2,SH2,SH1,GAUSS,PERC2,RESMIN,RESMAX
     *   ,POOL,MEMORY,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) THEN 
          IERR = 0
C ---     go to the final analysis
          GO TO 310
        ENDIF
C ---
        CALL MSGDOC(MDD,  '*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDD,'--- peaksrch for pattsrch function ---')
        CALL MSGDOC(MDD,  '*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        ALEVEL = 0.0
        NPMAX  = 0
        SPEC_P = SPEC
        IF(SPEC.EQ.'N') SPEC_P='P'
        CALL PEAKSRCH(M,NAMEP,NAMES,ALEVEL,DIST_LIM,SPEC_P,NPMAX
     *  ,POOL,MEMORY,ISYM,IPRSYM,RESULTA(1,1,IC),IPRMAXD,IPRROW,IERR)
        IF(IERR.NE.0) THEN
          IERR = 0
C ---     go to the final analysis
          GO TO 310
        ENDIF
        NPEAKSA = RESULTA(1,12,IC)

        IF(NPEAKSA.GT.30) THEN          
          NPEAKSA = 30
          RESULTA(1,12,IC) = NPEAKSA
        ENDIF

        IF(NPEAKSA.GT.0) THEN

C ---- set number of peak in array RESULTD ---

          DO IPK=1,NPEAKSA

            IX1 = RESULTA(IPK,1,IC)
            IY1 = RESULTA(IPK,2,IC)
            IZ1 = RESULTA(IPK,3,IC)


            IF(IFX.EQ.0) IX1 = 0
            IF(IFY.EQ.0) IY1 = 0
            IF(IFZ.EQ.0) IZ1 = 0

            DO ICD=1,NPEAKS

              IX2 = RESULTD(ICD,1)
              IY2 = RESULTD(ICD,2)
              IZ2 = RESULTD(ICD,3)

              IF(IFX.EQ.0) IX2 = 0
              IF(IFY.EQ.0) IY2 = 0
              IF(IFZ.EQ.0) IZ2 = 0


              CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
     *        ,INV,IOR,ISOP,DIST
     *        ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *        ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)
 
              IF(DIST.LT.DIST_LIM) THEN
                RESULTA(IPK,15,IC) = ICD
                GO TO 330
              ENDIF
            ENDDO  
            RESULTA(IPK,15,IC) = 0
 330        CONTINUE
          ENDDO
C -----
          IFIRST = 0
          DO IPK=1,NPEAKSA

            DO J=1,IPRROW
              RESULTC(1,J) = RESULTD(IC,J)
              RESULTC(2,J) = RESULTA(IPK,J,IC)
            ENDDO
            RESULTC(1,12) = 2

C ---       set the same weight --
            RESULTC(2,4) = RESULTC(1,4)

            NA_REF  = 2
            NAMES   = ' '
            MODE_CC = 'N'
            IFIRST  = IFIRST + 1
            CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *          ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *          ,RESMIN,RESMAX,BOFF,BADD
     *          ,RESULTC,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *          ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
            IF(IERR.NE.0) RETURN

            RESULTA(IPK, 6,IC) = RESULTC(2, 6)
            RESULTA(IPK, 7,IC) = RESULTC(2, 7)
            RESULTA(IPK, 8,IC) = RESULTC(2, 8)
            RESULTA(IPK, 9,IC) = RESULTC(2, 9)
            RESULTA(IPK,10,IC) = RESULTC(2,10)
            RESULTA(IPK,11,IC) = RESULTC(2,11)
            RESULTA(IPK, 1,IC) = RESULTA(IPK, 9,IC)*NX
            RESULTA(IPK, 2,IC) = RESULTA(IPK,10,IC)*NY
            RESULTA(IPK, 3,IC) = RESULTA(IPK,11,IC)*NZ

            RESULTA(IPK,13,IC) = RESULTC(1,13)
            RESULTA(IPK,16,IC) = RESULTC(2,16)

          ENDDO
          IF(LSORT.EQ.'P') THEN
C ---       sort RESULTA(1,1,IC)
            IST = 1
            CALL SORT_PEAK(MDOC,IST,RESULTA(1,1,IC)
     *      ,IPRMAXD,IPRROW,IERR)
          ENDIF
C --
          CALL MSGDOC(MDD,'====================')
          NPEAK  = RESULTA(1,12,IC)
          IPRINT = 1
          CALL PR_RES_PS(MDD,NPEAK,IPRINT
     *    ,RESULTA(1,1,IC),IPRMAXD,IPRROW,IERR)
          CALL MSGDOC(MDD,'====================')

        ENDIF
      ENDIF

C ---- end add. peaks --
C --------------------------------------------------------------------
C ---- search for best single, pair, triple

 310  CONTINUE
      NPEAKS  = RESULTD(1,12)
      IBEST2  =    0
      IBEST3  =    0
      POWER2  = -1.0
      POWER3  = -1.0

      IF(NPEAKS.GT.1) THEN

        CALL ANALYS_PK(MDOC,NP,RESULTD,RESULTA,RESULTB,RESULTC
     *          ,IPRMAXDD,IPRROWW,IPRMAXD,IPRROW,IPRNPR
     *          ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX
     *          ,INVERS1,IORIGN1,ISYMOP1,DISTN1
     *          ,INVERS2,IORIGN2,ISYMOP2,DISTN2,IPRNEQ,IPRORG
     *          ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)

C --------------------------------------
C      SUBROUTINE ANALYS_PK(MDOC,NP,RESULT,RESULTA,RESULTB,RESULTC
C     *          ,IPRMAXDD,IPRROWW,IPRMAX,IPRROW,IPRNPR
C     *          ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX
C     *          ,INVERS1,IORIGN1,ISYMOP1,DISTN1
C     *          ,INVERS2,IORIGN2,ISYMOP2,DISTN2,IPRNEQ,IPRORG
C     *          ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)
C --------------------------------------
C      input NPEAKS = RESULT(1,12)
C        IF(NPEAKS.LT.NP) NP = NPEAKS
C        RESULT (IPRMAXDD,IPRROWW)          IX1 = RESULT(IC1,1)
C        RESULTA(IPRMAX,IPRROW,IPRNPR)      JX1  =RESULTA(I,1,IC1)
C     
C      output
C
C        RESULTB(IPRMAX,IPRROW)  NPAIRS*2   = RESULTB(1,12)
C        RESULTC(IPRMAX,IPRROW)  NTRIPLES*3 = RESULTC(1,12) 
C
C      PARAMETER ( IPRPAIR = 60      )
C      PARAMETER ( IPRROW1 = 21      )
C                  ' ERR: IPRROW1 < IPRROW'
C                    IF(NATOMS.GT.IPRMAX) NATOMS = (IPRMAX/2)*2 
C
C --------------------------------------
        IF(IERR.NE.0) RETURN
C///\\\
C ---
        CALL MSGDOC(MDOC,'====================')
        CALL MSGDOC(MDOC
     *    ,'--- best single heavy atom  ---')
        NPEAK  = 1
        IPRINT = 1        
        POWER1 = RESULTD(1,13) 
        CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *  ,RESULTD,IPRMAXDD,IPRROW,IERR)

        CALL MSGDOC(MDOC,'====================')
C ---
        NPAIR = RESULTB(1,12)
        IF(NPAIR.GT.0) THEN
          IBEST2  =    0
          POWER2  = -1.0
          IFIRST  = 0
          DO IT=1,NPAIR,2

            DO J=1,IPRROW
              RESULTA(1,J,1) = RESULTB(IT  ,J)
              RESULTA(2,J,1) = RESULTB(IT+1,J)
            ENDDO

C ---       set the same weight --
            RESULTA(1,4,1) = 1.0
            RESULTA(2,4,1) = 1.0

            RESULTA(1,12,1)=2

            NA_REF  = 2
            NAMES   = ' '
            MODE_CC = 'N'
            IFIRST  = IFIRST + 1
            CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *          ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *          ,RESMIN,RESMAX,BOFF,BADD
     *          ,RESULTA,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *          ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
            IF(IERR.NE.0) RETURN

            RESULTB(IT  ,13) = RESULTA(1,13,1)
            RESULTB(IT+1,13) = RESULTA(1,13,1)
            RESULTB(IT  ,16) = RESULTA(1,16,1)
            RESULTB(IT+1,16) = RESULTA(2,16,1)

            IF(RESULTB(IT  ,13).GT.POWER2) THEN
              IBEST2 = IT
              POWER2 = RESULTB(IT  ,13)
            ENDIF
          ENDDO
          CALL MSGDOC(MDOC,'====================')
          CALL MSGDOC(MDOC
     *    ,'---  all pairs of heavy atoms  ---')
          CALL MSGDOC(MDOC,
     *    'Power for pair after unphased refinement')
          NPEAK  = NPAIR       
          IPRINT = 2        
          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *    ,RESULTB,IPRMAXD,IPRROW,IERR)
          CALL MSGDOC(MDOC,'====================')
C
C --- version : best pair ---
C

C --- end version : best pair ---

        ENDIF
C ----
C VVVV
        NTRIPLE = RESULTC(1,12)
        IF(NTRIPLE.GT.0) THEN
          IBEST3  =    0
          POWER3  = -1.0
          IFIRST  = 0
          DO IT=1,NTRIPLE,3

            DO J=1,IPRROW
              RESULTA(1,J,1) = RESULTC(IT  ,J)
              RESULTA(2,J,1) = RESULTC(IT+1,J)
              RESULTA(3,J,1) = RESULTC(IT+2,J)
            ENDDO
C ---       set the same weight --
            RESULTA(1,4,1) = 1.0
            RESULTA(2,4,1) = 1.0
            RESULTA(3,4,1) = 1.0

            RESULTA(1,12,1) = 3

            NA_REF  = 3
            NAMES   = ' '
            MODE_CC = 'N'
            IFIRST  = IFIRST + 1
            CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *          ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *          ,RESMIN,RESMAX,BOFF,BADD
     *          ,RESULTA,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *          ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
            IF(IERR.NE.0) RETURN
            RESULTC(IT  ,13) = RESULTA(1,13,1)
            RESULTC(IT+1,13) = RESULTA(1,13,1)
            RESULTC(IT+2,13) = RESULTA(1,13,1)
            IF(RESULTC(IT  ,13).GT.POWER3) THEN
              IBEST3 = IT
              POWER3 = RESULTC(IT  ,13)
            ENDIF
          ENDDO
          CALL MSGDOC(MDOC,'====================')
          CALL MSGDOC(MDOC
     *    ,'---  all triples of heavy atoms  ---')
          CALL MSGDOC(MDOC,
     *    'Power for triple after unphased refinement')
          NPEAK  = NTRIPLE       
          IPRINT = 3        
          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *    ,RESULTC,IPRMAXD,IPRROW,IERR)
          CALL MSGDOC(MDOC,'====================')
        ENDIF
C ^^^^^^
C ---------
 311    CONTINUE

        TEST1 = POWER1*0.999*0.999
        TEST2 = POWER2*0.999

        IF(TEST1.GT.TEST2.AND.TEST1.GT.POWER3) THEN
          NATOMS = 1
          DO K=1,IPRROW
            RESULTB(1,K) = RESULTD(1,K)
          ENDDO
          RESULTB(1,12) = 1

          POWER_PRV = POWER1
          NPEAKS_NK = 1

        ELSE IF(TEST2.GT.POWER3) THEN
          NATOMS = 2
          DO K=1,IPRROW
            RESULTB(1,K) = RESULTB(IBEST2  ,K)
            RESULTB(2,K) = RESULTB(IBEST2+1,K)
          ENDDO
          RESULTB(1,12) = 2

          POWER_PRV = POWER2
          NPEAKS_NK = 2

        ELSE
          NATOMS = 3
          DO K=1,IPRROW
            RESULTB(1,K) = RESULTC(IBEST3  ,K)
            RESULTB(2,K) = RESULTC(IBEST3+1,K)
            RESULTB(3,K) = RESULTC(IBEST3+2,K)
          ENDDO
          RESULTB(1,12) = 3

          POWER_PRV = POWER3
          NPEAKS_NK = 3

        ENDIF

        NPEAKS_SLF = RESULTD(1,12)


        IF(NPEAKS_NK.GT.1) THEN

 313      CONTINUE

          RESULTB(1,12) = NPEAKS_NK
          NPMAX         = NPEAKS_NK
          NPEAK         = NPEAKS_NK


          DO IPN=1,NPEAKS_NK
            RESULTB(IPN,4) = 1.0
          ENDDO

          NA_REF  = NPEAKS_NK
          NAMES   = ' '
          MODE_CC = 'N'
          IFIRST  = 1
          CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *        ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *        ,RESMIN,RESMAX,BOFF,BADD
     *        ,RESULTB,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *        ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
          IF(IERR.NE.0) THEN
            RETURN
          ENDIF

          CALL MSGDOC(MDOC,'====================')
          CALL MSGDOC(MDOC,'---  expected heavy atom.  ---')
          IPRINT=1        
          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *    ,RESULTB,IPRMAXD,IPRROW,IERR)
          CALL MSGDOC(MDOC,'====================')

          DO IPN=1,NPEAKS_NK
            RESULTB(IPN,4) = RESULTB(IPN,16)
          ENDDO

          CALL WR_RES_PS(M,NAMEO,NPMAX,RESULTB,IPRMAXD,IPRROW
     *    ,NSYM,ISYM,IPRSYM,IERR)
          N_FINAL = NPEAKS_NK   


          NPEAKS_ADD_MAX = 9
          IF(FAST.EQ.'Y') NPEAKS_ADD_MAX = 5 
          IF(NPEAKS_NK.GE.NPEAKS_ADD_MAX) THEN
            RESULTD(1,12) = NPEAKS_NK
            GO TO 9000
          ENDIF 
          POWER_PRV = RESULTB(1,13)
C ----

          IOP1 = 0
          IOP2 = 0
          WRITE(MOD2,'(I1)') NPEAKS_NK
           
          MODE = 'M'
          FUNC = 'Y'

          IF(TRANS.EQ.'N') MODE = 'T'

          IF(PACK.EQ.'N') FUNC = 'N'
          ACCUM   = 'M'

          INVER   = ' '
          RAD1    = RAD
          RAD2    = RAD
          PERC2(1)= 50.0

          SH1(1)  = RESULTB(1, 9)
          SH1(2)  = RESULTB(1,10)
          SH1(3)  = RESULTB(1,11)
          DO I=2,NPEAKS_NK
            SH2(1,I-1) = RESULTB(I, 9)
            SH2(2,I-1) = RESULTB(I,10)
            SH2(3,I-1) = RESULTB(I,11)
          ENDDO 

          NX       = NXMIN
          NY       = NYMIN
          NZ       = NZMIN
          IZMIN    = 0
          IZMAX    = NZMIN-1
          NAME(1)  = ' '
          NAME(2)  = ' '
          NAME(3)  = ' '
          NAME(4)  = ' '
          NAME(5)  = NAMEC
          NAME(6)  = ' '
          NAME(7)  = NAMET
          NAME(8)  = 'trahalo.scr'
          NAME(9)  = NAMEP
          NAME(10) = ' '

          BADDD    = SMOOTH

          IF(MTN.NE.0) THEN
            CALL MSGDOC(MDOC
     *     ,'--- trpack: search in Diff. Patterson         ---')
          ELSE
            CALL MSGDOC(MDOC
     *     ,'--- trpack: search in Anom. Diff. Patterson   ---')
          ENDIF
          CALL MSGDOC(MDOC
     *    ,'---          for expected atoms as fixed model ---')
          CALL MSGDOC(MD
     *    ,'*-*-*-*-*-*--*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
C ---
          CALL TRPACK(M,NAME,MODE,MOD2,FUNC,REMOV,INVER
     *    ,ACCUM,NX,NY,NZ,IZMIN,IZMAX,IOP1,IOP2,BADDD,BOFF
     *    ,RAD1,RAD2,SH2,SH1,GAUSS,PERC2,RESMIN,RESMAX
     *    ,POOL,MEMORY,ISYM,IPRSYM,IERR)
          IF(IERR.NE.0) THEN 
            IERR = 0
            GO TO 9000
          ENDIF
C ---
        CALL MSGDOC(MDD,  '*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDD,'--- peaksrch for pattsrch function ---')
        CALL MSGDOC(MDD,  '*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
          ALEVEL = 0.0
          NPMAX  = 0
          SPEC_P = SPEC
          IF(SPEC.EQ.'N') SPEC_P = 'P'
        CALL PEAKSRCH(M,NAMEP,NAMES,ALEVEL,DIST_LIM,SPEC_P,NPMAX
     *  ,POOL,MEMORY,ISYM,IPRSYM,RESULTA(1,1,1),IPRMAXD,IPRROW,IERR)
          IF(IERR.NE.0) THEN
            IERR = 0
            GO TO 9000
          ENDIF
          NPEAKSA = RESULTA(1,12,1)

          IF(NPEAKSA.GT.30) THEN          
            NPEAKSA = 30
            RESULTA(1,12,1) = NPEAKSA
          ENDIF

          IF(NPEAKSA.GT.0) THEN

C ----      set number of peak in array RESULTD ---

            CALL MSGDOC(MDD,'--- analysis ---')
            DO IPK=1,NPEAKSA

              IX1 = RESULTA(IPK,1,1)
              IY1 = RESULTA(IPK,2,1)
              IZ1 = RESULTA(IPK,3,1)

              IF(IFX.EQ.0) IX1 = 0
              IF(IFY.EQ.0) IY1 = 0
              IF(IFZ.EQ.0) IZ1 = 0

              DO ICD=1,NPEAKS

                IX2 = RESULTD(ICD,1)
                IY2 = RESULTD(ICD,2)
                IZ2 = RESULTD(ICD,3)

                IF(IFX.EQ.0) IX2 = 0
                IF(IFY.EQ.0) IY2 = 0
                IF(IFZ.EQ.0) IZ2 = 0

                CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
     *          ,INV,IOR,ISOP,DIST
     *          ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *          ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)
 
                IF(DIST.LT.DIST_LIM) THEN
                  RESULTA(IPK,15,1) = ICD
                  GO TO 321
                ENDIF
              ENDDO  
              RESULTA(IPK,15,1) = 0
 321          CONTINUE

            ENDDO
C -----
            CALL MSGDOC(MDD,'--- refinement ---')

            IFIRST = 0
            DO IPK=1,NPEAKSA

              DO K=1,NPEAKS_NK 
                DO J=1,IPRROW
                  RESULTC(K,J) = RESULTB(K,J  )
                ENDDO
C ---           set the same weight --
                RESULTC(K,4) = 1.0
              ENDDO

              DO J=1,IPRROW
                RESULTC(NPEAKS_NK+1,J) = RESULTA(IPK,J,1)
              ENDDO
C ---         set the same weight --
              RESULTC(NPEAKS_NK+1,4) = 1.0

              RESULTC(1,12) = NPEAKS_NK+1

              NA_REF  = NPEAKS_NK+1
              NAMES   = ' '
              MODE_CC = 'N'
              IFIRST  = IFIRST + 1
              CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *          ,NAMEF,NAMEFS,NAMECA,NAMES,NAMES,NA_REF,ANOM
     *          ,RESMIN,RESMAX,BOFF,BADD
     *          ,RESULTC,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *          ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
              IF(IERR.NE.0) RETURN
              RESULTA(IPK,13,1) = RESULTC(1,13)

C             OCC1             = RESULTC(1,16)
C             OCC2             = RESULTC(2,16)
              RESULTA(IPK,16,1) = RESULTC(3,16)

            ENDDO
            IF(LSORT.EQ.'P') THEN
C ---         sort RESULTA(1,1,1)
              IST = 1
              CALL SORT_PEAK(MDOC,IST,RESULTA(1,1,1)
     *        ,IPRMAXD,IPRROW,IERR)
            ENDIF
C --
          ELSE
            GO TO 9000
          ENDIF
C --
          CALL MSGDOC(MDD,'====================')
          CALL MSGDOC(MDD,
     *'Power current peak and exp.atoms after unphased refinement')
          CALL MSGDOC(MDD,
     *'with the same weight. No - peak number in Self_patt list')
          NPEAK  = RESULTA(1,12,1)
          IPRINT = 1
          CALL PR_RES_PS(MDD,NPEAK,IPRINT
     *    ,RESULTA(1,1,1),IPRMAXD,IPRROW,IERR)
          CALL MSGDOC(MDD,'====================')

          I_ADD = 0
          DO I=1,NPEAKSA

            POWER = RESULTA(I,13,1)
            NO    = RESULTA(I,15,1)
            TEST  = POWER

            IF(TEST.LE.POWER_PRV) GO TO 9000

            DO J=1,NPEAKS_NK

              NO_NK = RESULTB(J,15)

              IF(NO.EQ.NO_NK) THEN
                GO TO 322
              ELSE IF(NO.LE.0) THEN

                IX1 = RESULTA(I,1,1)
                IY1 = RESULTA(I,2,1)
                IZ1 = RESULTA(I,3,1)

                IX2 = RESULTB(J,1)
                IY2 = RESULTB(J,2)
                IZ2 = RESULTB(J,3)

                CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
     *          ,INV,IOR,ISOP,DIST
     *          ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *          ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)
 
                IF(DIST.LT.DIST_LIM) THEN
                  RESULTA(I,15,1) = NO_NK
                  GO TO 322
                ENDIF
              ENDIF

            ENDDO
            I_ADD = I
            IF(NO.LE.0) THEN
              NPEAKS_SLF      = NPEAKS_SLF+1
              RESULTA(I,15,1) = NPEAKS_SLF
              NO              = NPEAKS_SLF  
            ENDIF
            GO TO 323
 322        CONTINUE
          ENDDO

          GO TO 9000

 323      CONTINUE

          POWER_PRV = POWER
          NPEAKS_NK = NPEAKS_NK+1

          DO J=1,IPRROW
            RESULTB(NPEAKS_NK,J) = RESULTA(I_ADD,J,1)
          ENDDO
          RESULTB(1,12) = NPEAKS_NK
C
C    add to self_list if peak is new
C
          N = RESULTD(1,12)
          DO I=1,N              
            NO_SLF = RESULTD(I,15)
            IF(NO_SLF.EQ.NO) GO TO 313
          ENDDO
          DO I=1,NPEAKS_NK-1              
            NO_NK = RESULTB(I,15)
            IF(NO_NK.EQ.NO) GO TO 313
          ENDDO

          IF(N.LT.IPRMAXD) THEN

            NPNEW = N + 1
            DO J=1,IPRROW
              RESULTD(NPNEW,J) = RESULTB(NPEAKS_NK,J)
            ENDDO
            RESULTD(1,12) = NPNEW

            IF(NO_DER.GT.0) THEN
              RES_DER(1,12,NO_DER) = NPNEW
              DO J=1,IPRROWW
                RES_DER(NPNEW,J,NO_DER) = RESULTD(NPNEW,J)  
              ENDDO
            ENDIF
            WRITE(LINE,'('' New peak : ('',3F8.3,
     *      '') was added to self_list'')') 
     *      RESULTD(NPNEW, 9),RESULTD(NPNEW,10),RESULTD(NPNEW,11) 
            CALL MSGDOC(MDOC,LINE)
            WRITE(LINE,'('' with number:'',I4)') NPNEW
            CALL MSGDOC(MDOC,LINE)

          ELSE   
               
            DO I=N,1,-1              
              NO_SLF = RESULTD(I,15)
              DO J=1,NPEAKS_NK-1
                NO_NK = RESULTB(J,15)
                IF(NO_NK.EQ.NO_SLF) GO TO 324 
              ENDDO
              NPNEW = I
              GO TO 325
 324          CONTINUE
            ENDDO
            GO TO 313

 325        CONTINUE
            DO J=1,IPRROW
              RESULTD(NPNEW,J) = RESULTB(NPEAKS_NK,J)
            ENDDO
            RESULTB(NPEAKS_NK,15) = NPNEW
            RESULTD(NPNEW,15)     = NPNEW

            IF(NO_DER.GT.0) THEN
              DO J=1,IPRROW
                RES_DER(NPNEW,J,NO_DER) = RESULTD(NPNEW,J)  
              ENDDO
            ENDIF

            WRITE(LINE,'(
     *      '' Program replaced peak with number:'',I4)') NPNEW
            CALL MSGDOC(MDOC,LINE)
            WRITE(LINE,'(''to new peak : ('',3F8.3,
     *      '') in the self_list'')') 
     *      RESULTD(NPNEW, 9),RESULTD(NPNEW,10),RESULTD(NPNEW,11) 
            CALL MSGDOC(MDOC,LINE)

          ENDIF

          GO TO 313

        ENDIF

        NATOMS        = RESULTB(1,12)
        RESULTD(1,12) = NATOMS
C ---
        IF(NATOMS.GT.0) THEN
          CALL MSGDOC(MDOC,'====================')
          CALL MSGDOC(MDOC,'---  List of expected heavy atoms  ---')
          NPEAK  = NATOMS   
          NPMAX  = NATOMS
          IPRINT = 1         
          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *      ,RESULTB,IPRMAXD,IPRROW,IERR)
          CALL MSGDOC(MDOC,'====================')

          DO IPN=1,NATOMS      
            RESULTB(IPN,4) = RESULTB(IPN,16)
          ENDDO

          CALL WR_RES_PS(M,NAMEO,NPMAX,RESULTB,IPRMAXD,IPRROW
     *    ,NSYM,ISYM,IPRSYM,IERR)
          N_FINAL = NATOMS
        ENDIF
C --
      ELSE

        NATOMS        = 1
        RESULTD(1,12) = NATOMS
        NPMAX         = NATOMS
        CALL MSGDOC(MDOC,'====================')
        CALL MSGDOC(MDOC,'---  expected heavy atom  ---')
        NPEAK  = NATOMS   
        NPMAX  = NATOMS
        IPRINT = 1        
        CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *  ,RESULTD,IPRMAXD,IPRROW,IERR)
        CALL MSGDOC(MDOC,'====================')

c        RESULTD(IPN,4)=RESULTD(IPN,16)

        CALL WR_RES_PS(MDOC,NAMEO,NPMAX,RESULTD,IPRMAXDD,IPRROW
     *  ,NSYM,ISYM,IPRSYM,IERR)
        N_FINAL = NATOMS

      ENDIF

 9000 CONTINUE

      WRITE(LINE,'(I5,'' atoms were found'')') N_FINAL
      CALL MSGDOC(MDOC,LINE)

      IF(N_FINAL.GT.0) THEN
        NATOMR = 0
        NAMES  = ' '
        PDB    = 'C'
        CALL CHECK_CRD(M,NAMEO,ISYM,IPRSYM
     *    ,STR_TITLE,STR_DATE,NAMES,NATOMR
     *    ,RESULTB,IPRMAXD,IPRROW,PDB,IERR)
        IF(IERR.NE.0) GO TO 900
        NATOMR = RESULTB(1,2)
        IF(NATOMR.GT.0) THEN
          DO I=1,NATOMR
            IND           = I+10
            RESULTB(I,1 ) = RESULTB(IND,9 )*NX
            RESULTB(I,2 ) = RESULTB(IND,10)*NY
            RESULTB(I,3 ) = RESULTB(IND,11)*NZ
C           occ
            RESULTB(I,4 ) = RESULTB(IND,4 )
C ---         set the same weight --
            RESULTB(I,4)  = 1.0
C           
            RESULTB(I,5 ) = RESULTB(IND,5 )
C           ort
            RESULTB(I,6 ) = RESULTB(IND,6 )
            RESULTB(I,7 ) = RESULTB(IND,7 )
            RESULTB(I,8 ) = RESULTB(IND,8 )
C           fract
            RESULTB(I,9 ) = RESULTB(IND,9 )
            RESULTB(I,10) = RESULTB(IND,10)
            RESULTB(I,11) = RESULTB(IND,11)
            RESULTB(I,13) = 0.0
            RESULTB(I,15) = I
            RESULTB(I,16) = RESULTB(IND,4 )

            RESULTB(1,12) = NATOMR

          ENDDO

          NPMAX_PST = NATOMR

          IF(NPST.GT.0) THEN
            DO IPR=1,NATOMR
              DO J=1,NPST
                NPMAX_PST = NPMAX_PST +1
                DO K=1,IPRROW
                  RESULTB(NPMAX_PST,K) = RESULTB(IPR,K) 
                ENDDO
                RESULTB(NPMAX_PST, 9) = RESULTB(IPR, 9) + SPST(1,J)
                RESULTB(NPMAX_PST,10) = RESULTB(IPR,10) + SPST(2,J)
                RESULTB(NPMAX_PST,11) = RESULTB(IPR,11) + SPST(3,J)
              ENDDO
            ENDDO
            RESULTB(1,12) = NPMAX_PST
          ENDIF

          N_FINAL = NPMAX_PST

c      CALL MSGDOC(MDOC,'====================')
c      CALL MSGDOC(MD  ,'--*-*-*-*-*-*-*-*-*-*-*-*-*-*')
c      CALL MSGDOC(MDOC,'----  List of heavy atoms ---')
c      CALL MSGDOC(MD  ,'--*-*-*-*-*-*-*-*-*-*-*-*-*-*')
c          NPEAK  = RESULTB(1,12)
c          IPRINT = 1
c          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
c     *      ,RESULTB,IPRMAXD,IPRROW,IERR)
c      CALL MSGDOC(MDOC,'====================')

        ENDIF

        NA_REF  = N_FINAL
        MODE_CC = 'N'
        NAMEAS  = NAMEOA
        IFIRST  = 1
        CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *      ,NAMEF,NAMEFS,NAMECA,NAMEAS,NAMEO,NA_REF,ANOM
     *      ,RESMIN,RESMAX,BOFF,BADD
     *      ,RESULTB,IPRMAXD,RESULT,RESULTP,IPRMAX,IPRROW
     *      ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) RETURN


        CALL MSGDOC(MDOC,'====================')
        CALL MSGDOC(MD  ,'--*-*-*-*-*-*-*-*-*-*-*-*-*-*')
        CALL MSGDOC(MDOC,'----  List of heavy atoms ---')
        CALL MSGDOC(MD  ,'--*-*-*-*-*-*-*-*-*-*-*-*-*-*')
        NPEAK  = RESULTB(1,12)
        IPRINT = 1
        CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *      ,RESULTB,IPRMAXD,IPRROW,IERR)
        CALL MSGDOC(MDOC,'====================')

C        NPMAX = NPEAK
C        CALL WR_RES_PS(MDOC,NAMEO,NPMAX,RESULTB,IPRMAXD,IPRROW
C     *  ,NSYM,ISYM,IPRSYM,IERR)
C

        GO TO 700

      ENDIF

      RETURN

C =================================================================      
C --- four. search ---


 700  CONTINUE

      CALL LENSTR_BL(NAMEAS,LEN)
      IF(LEN.GT.0.AND.NAMEAS(1:1).NE.','.AND.NAMEAS(1:1).NE.' ') THEN      
C ---
        CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDOC,'--- ABCD --> PH  ---')
        CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-')
        IPRINT = 0
        SCCOMB = 0.0
        NAMES2 = ' '
        STOP_T = 'N'
        CALL ABCDPH(M,NAMEAS,NAMES,NAMEPH,NAMES2,SCCOMB
     *  ,RESMIN,RESMAX,IPRINT,STOP_T
     *  ,ISYM,IPRSYM,RESULT,IPRROW,IERR)
C       RESULT(1) - number of refls /in output file/
C       RESULT(2) - <FOM>
        IF(IERR.NE.0) GO TO 900
        WRITE(LINE,'(''figure-of-merit of input ABCD:'',F10.3)') 
     *  RESULT(2,1)  
        CALL MSGDOC(MDOC,LINE)
C ---

        ISIM   = 0
        SCALEF = 1.0
        LAP    = 'N'
        FOMMIN = 0.0

        IF(MTN.NE.0) THEN
          CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
          CALL MSGDOC(MDOC,'---    coef for Diff Fourie    ---')
          CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
          MSGA  = 'N'
c          IF(ANOM.EQ.'Y') THEN
c            MSGA = 'Y'
c            IF(COMB.EQ.'Y') MSGA = 'C'
c          ENDIF
          FILE1 = NAMEF
          FILE2 = NAMEFS
        ELSE
          CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
          CALL MSGDOC(MDOC,'--- coef for Anom. Diff Fourie ---')
          CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-')
          MSGA  = 'Y'
c          IF(COMB.EQ.'Y') MSGA = 'C'
          FILE1 = NAMED
          FILE2 = ' '
        ENDIF
C ---
        SLIM   = 0.0
        SLIM   = SLIM_PATT
C              = BADD_FOUR ??? but slim and del_lim
        BADDD  = BADD_COEF
        NAMES  = ' ' 
        CALL COEF_FFT(M,FILE1,FILE2,NAMEPH,NAMEC3,NAMES
     *         ,FOMMIN,RESMIN,RESMAX,BOFF,SLIM
     *         ,DEL_LIM,DELA_LIM,DRMS,DRMS_A,LAP
     *         ,SCALEF,BADDD,MSGA,ISIM,ISYM,IPRSYM,RESULT,IERR)
        IF(IERR.NE.0) GO TO 900
        IF(COMB.EQ.'Y') THEN
          IF(MSGA.NE.'N') THEN
            CALL MSGDOC(MDOC,'    Anom.data of derivative was used')
          ENDIF
        ENDIF
        IF(ANOM.EQ.'Y'.AND.MSGA.EQ.'N'.AND.MTN.EQ.0) THEN
          CALL MSGERR(MDOC,
     *    ' ERROR: No anom.signal in the derivative')
          IERR=1
          RETURN
        ENDIF
C ---
        CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDOC,'---     cfft     ---')
        CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-')
        F000  = 0.0
        IZMIN = 0
        IZMAX =-1
        CALL FFT(M,POOL,NX,NY,NZ,IZMIN,IZMAX,F000,NAMEC3
     *  ,NAMEP,MEMORY,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) GO TO 900
C ---
        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
        CALL MSGDOC(MDOC,'--- peak search in  diff. fourie ---')
        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
        ALEVEL = 0.0
        NPMAX  = 0
        CALL PEAKSRCH(M,NAMEP,NAMES,ALEVEL,DIST_LIM,SPEC,NPMAX
     *  ,POOL,MEMORY,ISYM,IPRSYM,RESULTA,IPRMAXD,IPRROW,IERR)
        IF(IERR.NE.0) GO TO 900
C                                
C  RESULT(i,1) IX 
C  RESULT(i,2) IY
C  RESULT(i,3) IZ
C  RESULT(i,4) Dens
C  RESULT(i,5) sigma
C  RESULT(i,6) alpha      Xort
C  RESULT(i,7) beta   or  Yort
C  RESULT(i,8) gamma      Zort
C  RESULT(i,9)  Xfrac
C  RESULT(i,10) Yfrac
C  RESULT(i,11) Zfrac
C
C  RESULT(1,12)  number of peaks
C
        PK_LIM   = 0.0
        NPEAKS_F = RESULTA(1,12,1)


        IF(NPEAKS_F.GT.30) THEN          
          NPEAKS_F = 30
          RESULTA(1,12,1) = NPEAKS_F
        ENDIF


        IF(NPEAKS_F.LE.0) THEN
          IERR = 0
          CALL MSGERR(MDOC,
     *    'WARNING: can not found any peaks in diff.fourie map.')
          RESULTD(1,12) = 0
          IERR = 14
          GO TO 900
        ELSE

          DPK1   = RESULTA(1,4,1)
          DRMS1  = RESULTA(1,5,1)
          IF(DRMS1.LT.0.0000001) DRMS1 = 1.0
          DRMS   = DPK1/DRMS1

          PRMS = 0.0
          DO IC1=1,NPEAKS_F
            PK1               = RESULTA(IC1,4,1)
            PRMS              = PRMS+PK1*PK1
            RESULTA(IC1,15,1) = IC1
            RESULTA(IC1,13,1) = 1.0
            RESULTA(IC1,16,1) = 1.0
          ENDDO
          T    = PRMS/NPEAKS_F
          PRMS = SQRT(ABS(T))

      CALL MSGDOC(MDOC,'====================')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(MDOC,'--- List of peaks of diff. fourie map  ---')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
            NPEAK  = RESULTA(1,12,1)
            IPRINT = 1
            CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *      ,RESULTA,IPRMAXD,IPRROW,IERR)
      CALL MSGDOC(MDOC,'====================')

c          WRITE(LINE,'('' D1,DRMS1,RMS:'',3G16.4)') 
c     *    DPK1,DRMS1,DRMS
c          CALL MSGDOC(MD,LINE)

C          WRITE(LINE,'('' Dens_rms,Npeaks:'',G16.4,I6)') 
C     *    PRMS,NPEAKS_F
C          CALL MSGDOC(MD,LINE)

C          PK_LIM = PRMS*0.9

          PK_LIM = DRMS * TIMES_RMS

          IF(NPEAKS_F.LE.2) PK_LIM = PK_LIM*0.0
          WRITE(LINE,'('' Peak_limit :'',G16.4)') PK_LIM
          CALL MSGDOC(MDOC,LINE)
            
        ENDIF
C -----
C -----------------
        NCATOM     = 0
        CALL LENSTR_BL(NAMECPH,LEN)

        IF(LEN.GT.0.AND.NAMECPH(1:1).NE.','.AND.NAMECPH(1:1).NE.' ') 
     *  THEN

C        IF(LEN.GT.0.AND.NAMECPH(1:1).NE.','.AND.NAMECPH(1:1).NE.' '.AND.
C     *     OPT1.NE.'N') THEN


          IF(SUBR.EQ.'N') THEN
            CALL LENSTR_BL(NAMECPH,LEN)
            IF(LEN.GT.54) LEN=54
            LINE='Input file HA of phases:'//NAMECPH(1:LEN)
            CALL MSGDOC(MD,LINE)
          ENDIF
          NATOMR = 0
          NAMES  = ' '
          PDB    = 'C'
          CALL CHECK_CRD(M,NAMECPH,ISYM,IPRSYM
     *    ,STR_TITLE,STR_DATE,NAMES,NATOMR
     *    ,RESULTB,IPRMAXD,IPRROW,PDB,IERR)
          IF(IERR.NE.0) GO TO 900
          NCATOM = RESULTB(11,12)
          IF(NCATOM.GT.0) THEN
            DO I=1,NCATOM
              IND           = I+10
              RESULTB(I,1 ) = RESULTB(IND,9 )*NX
              RESULTB(I,2 ) = RESULTB(IND,10)*NY
              RESULTB(I,3 ) = RESULTB(IND,11)*NZ
              RESULTB(I,4 ) = RESULTB(IND,4 )
              RESULTB(I,5 ) = RESULTB(IND,5 )
              RESULTB(I,6 ) = RESULTB(IND,6 )
              RESULTB(I,7 ) = RESULTB(IND,7 )
              RESULTB(I,8 ) = RESULTB(IND,8 )
              RESULTB(I,9 ) = RESULTB(IND,9 )
              RESULTB(I,10) = RESULTB(IND,10)
              RESULTB(I,11) = RESULTB(IND,11)
              RESULTB(I,13) = 0.0
              RESULTB(I,15) = I
              RESULTB(I,16) = RESULTB(IND,4 )

              RESULTB(1,12) = NCATOM

            ENDDO
      CALL MSGDOC(MDOC,'====================')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
      CALL MSGDOC(MDOC,'---  List of heavy atoms of file ABCD  ---')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
            NPEAK  = RESULTB(1,12)
            IPRINT = 1
            CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *      ,RESULTB,IPRMAXD,IPRROW,IERR)
      CALL MSGDOC(MDOC,'====================')

            II    = 0
            IPRMS = 0
            PRMS  = 0.0
            
            DO IC1=1,NPEAKS_F

              IX1 = RESULTA(IC1,1,1)
              IY1 = RESULTA(IC1,2,1)
              IZ1 = RESULTA(IC1,3,1)
              PK1 = RESULTA(IC1,4,1)

              IFLAG_RMS = 0

              DO IC2=1,NCATOM
                IX2 = RESULTB(IC2,1)
                IY2 = RESULTB(IC2,2)
                IZ2 = RESULTB(IC2,3)
                PK2 = RESULTB(IC2,4)

                CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
     *          ,INV,IOR,ISOP,DIST
     *          ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *          ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)

                IF(DIST.LT.DIST_LIM) THEN

                  IF(OPT1.NE.'N') THEN

                    RESULTA(IC1, 4,1) = 0.0
                    RESULTA(IC1,13,1) = 0.0
                    RESULTA(IC1,15,1) = 0.0
           

                    WRITE(LINE,
     *         '('' Reject peak :'',I5,'' because = HA-peak :'',I5)') 
     *              IC1,IC2
                    CALL MSGDOC(MD,LINE)

                    GO TO 410
                  ELSE
                    IFLAG_RMS = 1
                  ENDIF

                ENDIF
              ENDDO
              II = II + 1

c              DO J=1,IPRROW
c                RESULTA(II,J,1)=RESULTA(IC1,J,1)
c              ENDDO

              PK1  = RESULTA(IC1,4,1)
              IF(IFLAG_RMS.EQ.0) THEN
                PRMS  = PRMS+PK1*PK1
                IPRMS = IPRMS + 1
              ENDIF

 410          CONTINUE
            ENDDO            

            IF(II.GT.0) THEN
              PK_LIM = PK_LIM * 0.8
              LINE = 
     *' Peak_limit will be changed because there are common atoms'
              CALL MSGDOC(MDOC,LINE)
              WRITE(LINE,
     *'('' in Self_list and in Input file HA of phases N_com ='',I4)') 
     *        II
              CALL MSGDOC(MDOC,LINE)
              WRITE(LINE,'('' and now Peak_limit :'',G16.4)') PK_LIM
              CALL MSGDOC(MDOC,LINE)
            ENDIF

C //////////////////
C            IF(II.LE.0) THEN
C              IERR = 0
C              CALL MSGERR(MDOC,
C     *        'WARNING: can not found any peaks in diff.fourie map.')
C              RESULTD(1,12) = 0
C              IERR = 14
C              GO TO 900
C            ENDIF
C
C            IF(IPRMS.LE.0) IPRMS = 1
C            T    = PRMS/IPRMS
C            PRMS = SQRT(ABS(T))
C            WRITE(LINE,'('' Dens_rms,Npeaks:'',G16.4,I6)') 
C     *      PRMS,II
C --
C            IF(NPEAKS_F.GT.1) THEN
C              II = 1
C              DO IC1=2,NPEAKS_F
C
C                IX1 = RESULTA(IC1,1,1)
C                IY1 = RESULTA(IC1,2,1)
C                IZ1 = RESULTA(IC1,3,1)
C                PK1 = RESULTA(IC1,4,1)
C                IF(PK1.LE.0.0) GO TO 430
C                IF(PK1.LT.PK_LIM) THEN
C                  IF(PK1.GT.0.0) THEN
C                    WRITE(LINE,
C     *         '('' Reject peak :'',I5,'' because < rms '')') 
C     *              IC1
C                    CALL MSGDOC(M,LINE)
C                  ENDIF
C                  RESULTA(IC1,4,1)  = 0.0
C                  RESULTA(IC1,13,1) = 0.0
C                  RESULTA(IC1,15,1) = 0.0
C                  GO TO 430
C                ENDIF
C                DO IC2=1,II
C                  IX2 = RESULTA(IC2,1,1)
C                  IY2 = RESULTA(IC2,2,1)
C                  IZ2 = RESULTA(IC2,3,1)
C                  PK2 = RESULTA(IC2,4,1)
C                  IF(PK2.LE.0.0) GO TO 430
C
C                  CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
C     *            ,INV,IOR,ISOP,DIST
C     *            ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
C     *            ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)
C                  IF(DIST.LT.DIST_LIM) THEN
C                    WRITE(LINE,
C     *       '('' Reject peak :'',I5,'' because = peak :'',I5)') 
C     *              IC1,IC2
C                    CALL MSGDOC(MD,LINE)
C                    RESULTA(IC1,4,1)  = 0.0
C                    RESULTA(IC1,13,1) = 0.0
C                    RESULTA(IC1,15,1) = 0.0
C                    GO TO 430
C                  ENDIF
C 431              CONTINUE
C                ENDDO
C                II = II + 1
C
C 430            CONTINUE
C              ENDDO            
C
C              IF(II.LE.0) THEN
C                IERR=0
C                CALL MSGERR(MDOC,
C     *     'WARNING: can not found any peaks in diff.fourie map.')
C                RESULTD(1,12) = 0
C                IERR          = 14
C                GO TO 900
C              ENDIF
C
C            ENDIF
C ///////////////
C --
          ENDIF

        ENDIF
C ----------
C ----------
        NCOMM = 0
        DO IC1=1,NPEAKS_F

          IX1 = RESULTA(IC1,1,1)
          IY1 = RESULTA(IC1,2,1)
          IZ1 = RESULTA(IC1,3,1)
          PK1 = RESULTA(IC1,4,1)
          POW = RESULTA(IC1,13,1)
          NOA = RESULTA(IC1,15,1)
          IF(PK1.LE.0.0) GO TO 400

          IF(PK1.LT.PK_LIM) THEN
            IF(PK1.GT.0.0) THEN
              WRITE(LINE,
     *        '('' Reject peak :'',I5,'' because < rms '')') 
     *        IC1
              CALL MSGDOC(M,LINE)
            ENDIF 
            RESULTA(IC1,15,1) = 0.0
            RESULTA(IC1,13,1) = 0.0
            GO TO 400
          ENDIF

          IF(IC1.GT.1) THEN
            DO IC2=1,IC1-1
              IX2 = RESULTA(IC2,1,1)
              IY2 = RESULTA(IC2,2,1)
              IZ2 = RESULTA(IC2,3,1)
              PK2 = RESULTA(IC2,4,1)
              IF(PK2.GT.0.0) THEN

                CALL COMP1_PNT(MDOC,IX1,IY1,IZ1,IX2,IY2,IZ2
     *          ,INV,IOR,ISOP,DIST
     *          ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *          ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)


C ???????????????????
C                IF(DIST.LT.DIST_LIM) THEN
C                  RESULTA(IC1, 4,1) = 0.0
C                  RESULTA(IC1,13,1) = 0.0
C                  RESULTA(IC1,15,1) = 0.0
C                    WRITE(LINE,
C     *       '('' Reject peak :'',I5,'' because = peak :'',I5)') 
C     *              IC1,IC2
C                  CALL MSGDOC(MD,LINE)
C                  RESULTA(IC1,15,1) = 0.0
C                  RESULTA(IC1,13,1) = 0.0
C                  GO TO 400
C                ENDIF
C ???????????????????
              ENDIF
            ENDDO
          ENDIF

          IF(NPEAKS_PATT.GT.0) THEN
            DO IC2=1,NPEAKS_PATT

              IX2  = RESULTD(IC2,1)
              IY2  = RESULTD(IC2,2)
              IZ2  = RESULTD(IC2,3)
              PK2  = RESULTD(IC2,4)

              IXX1 = IX1
              IYY1 = IY1
              IZZ1 = IZ1

              IF(NSYM.GT.1) THEN
                IF(IFX.EQ.0) IXX1 = 0
                IF(IFY.EQ.0) IYY1 = 0
                IF(IFZ.EQ.0) IZZ1 = 0

                IF(IFX.EQ.0) IX2 = 0
                IF(IFY.EQ.0) IY2 = 0
                IF(IFZ.EQ.0) IZ2 = 0
              ENDIF

              IF(PK2.LE.0.0) GO TO 401
              POWD = RESULTD(IC2,13)
              PK12 = PK1*PK2

              CALL COMP1_PNT(MDOC,IXX1,IYY1,IZZ1,IX2,IY2,IZ2
     *        ,INV,IOR,ISOP,DIST
     *        ,DIST_LIM,NX,NY,NZ,IZMIN,IZMAX,IPRORG
     *        ,SOR,NORIG,NSYM,ISYM,IPRSYM,IERR)
 
              IF(DIST.LT.DIST_LIM) THEN
                NCOMM             = NCOMM+1
                RESULTA(IC1,15,1) = RESULTD(IC2,15)
                WRITE(LINE,
     *'(I4,'' Common peaks :'',I5,'' /four/ ,'',I5,'' /self/'')') 
     *          NCOMM,IC1,IC2
                CALL MSGDOC(MD,LINE)
                GO TO 400
              ENDIF
 401          CONTINUE
            ENDDO  

            RESULTA(IC1,15,1) = 0.0
            RESULTA(IC1,13,1) = 0.0

          ENDIF

 400      CONTINUE
        ENDDO

        IF(NPEAKS_PATT.GT.0) THEN
          WRITE(LINE,
     *    '('' Npeaks_common in self_patt and fourie maps.:'',I5)') 
     *    NCOMM
          CALL MSGDOC(MDOC,LINE)
        ENDIF

        NPF_MAX = 10
        N       = 0
        DO K=1,NPEAKS_F
          IX = RESULTA(K,1,1)
          IY = RESULTA(K,2,1)
          IZ = RESULTA(K,3,1)
          DI = RESULTA(K,13,1)
          IF(DI.LE.0.0) GO TO 420
          IF(NCOMM.EQ.0) DI = RESULTA(K,4,1)
          IF(N.GT.0.AND.NCOMM.GT.0) THEN
            DO L=1,N
              IF(RESULTA(K,15,1).EQ.RESULTA(L,15,1)) GO TO 420             
            ENDDO
          ENDIF
          IF(DI.GT.0.0) THEN
            N = N+1
            DO L=1,IPRROW
              RESULTA(N,L,1) = RESULTA(K,L,1)
            ENDDO
            IF(N.GE.NPF_MAX) GO TO 460
          ENDIF
 420      CONTINUE
        ENDDO
 460    CONTINUE

        NATOMS          = N
        NPEAKS_F        = N
        RESULTA(1,12,1) = NATOMS

        IF(N.GT.0) THEN

C          CALL MSGDOC(MDOC,'====================')
C          CALL MSGDOC(MDOC,
C     * '---  List of expected heavy atoms from diff. fourier map ---')
C          NPEAK  = NATOMS        
C          IPRINT = 1        
C          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
C     *    ,RESULTA,IPRMAXD,IPRROW,IERR)
C          CALL MSGDOC(MDOC,'====================')

          DO I=1,NATOMS
            DO L=1,IPRROW
              RESULTD(I,L) = RESULTA(I,L,1)
            ENDDO
          ENDDO
          RESULTD(1,12) = NATOMS
          NPMAX         = NATOMS

          DO IPN=1,NATOMS      
C            RESULTD(IPN,4)=RESULTD(IPN,16)
C ---           set the same weight --
            RESULTD(IPN,4) = 1.0
          ENDDO

          NA_REF  = NATOMS
          MODE_CC = 'N'
          NAMEAS  = NAMEOA
          IFIRST  = 1

          CALL CHECK_AT_REF(M,MODE_CC,IFIRST
     *        ,NAMEF,NAMEFS,NAMECA,NAMEAS,NAMEO,NA_REF,ANOM
     *        ,RESMIN,RESMAX,BOFF,BADD
     *        ,RESULTD,IPRMAXDD,RESULT,RESULTP,IPRMAX,IPRROW
     *        ,POOL,MEMORY,NSYM,ISYM,IPRSYM,IERR)
          IF(IERR.NE.0) RETURN

          CALL MSGDOC(MDOC,'====================')
          CALL MSGDOC(MDOC,
     * '---  List of expected heavy atoms from diff. fourier map ----')
          NPEAK  = NATOMS        
          IPRINT = 1        
          CALL PR_RES_PS(MDOC,NPEAK,IPRINT
     *    ,RESULTD,IPRMAXDD,IPRROW,IERR)
          CALL MSGDOC(MDOC,'====================')

C          CALL WR_RES_PS(MDOC,NAMEO,NPMAX,RESULTD,IPRMAXDD,IPRROW
C     *    ,NSYM,ISYM,IPRSYM,IERR)

          N_FINAL = NATOMS
          IF(NCOMM.LE.0) THEN
            CALL MSGDOC(MDOC,
     *      ' no common peaks in self_patt and fourie maps .')
            IERR=13
          ENDIF

        ELSE

          CALL MSGERR(MDOC,
     *    ' no common peaks in self_patt and fourie maps.')
          RESULTD(1,12) = 0
          N_FINAL       = 0
          IERR          = 12

        ENDIF

      ELSE

        CALL MSGERR(MDOC,
     *  'ERROR: wrong name of ABCD_file.')
         RESULTD(1,12) = 0
         N_FINAL       = 0
         IERR          = 1

      ENDIF
C --------------------------------------------
  900 CONTINUE
      RETURN
      END
C --- TRAHALO.FTN ---
C --------------------------------------------------------------------
      PARAMETER ( MEMORY = 4000000   )
      REAL      POOL(MEMORY)
C --- IPRSYM - maximal number of symmetry operators
      PARAMETER ( IPRSYM = 48       )
      INTEGER*2 ISYM(5,3,IPRSYM)
C ----
C     PARAMETER ( IPRMAXD= 41      )
      PARAMETER  ( IPRMAXD = 71      )
      PARAMETER  ( IPRMAX  = 21      )
      PARAMETER  ( IPRROW  = 21      )
      PARAMETER  ( IPRNPR  = 21      )
      REAL      RESULT (IPRMAX,IPRROW)
      REAL      RESULTB(IPRMAXD,IPRROW)
      REAL      RESULTC(IPRMAXD,IPRROW)
      REAL      RESULTA(IPRMAXD,IPRROW,IPRNPR)
C
      REAL      SPST(3,1)
C ---
      CHARACTER NAMEF*80,NAMED*80,NAMEO*80,NAMEA*80,SORT*1,SUBR*1
      CHARACTER NAMEHA*80,NAMECPH*80,COMB*1,ANISO*1
      CHARACTER LINE*80,TTT*80,ITYPE*4,GAUSS*1,INVERT*1,MSCL*1
      CHARACTER PACK*1,SPEC*1,FAST*1,TRANS*1,STR_PARAM*16,REMOV*1
      CHARACTER PROG*80,MSG*1
C ---
      INCLUDE 'trahalo_version.fh'
C --------------------------------------------------------------------
C --------------------------------------------------------------------
      M    = 0
      PROG = 'trahalo'
      CALL START(0,M,PROG)
C ---------------------------
      CALL MSGDOC(M,' ')
      CALL MSGDOC(M,VERSION)
      CALL MSGDOC(M,' ')
C -------------------------------------------------------------
      CALL MSGDOC(M,
     *' Do you want to have FILE-DOCUMENT /trahalo.doc/ ? /<N>/Y/A :')
      CALL MSGDOC(M,'   N - means without DOC-file')
      CALL MSGDOC(M,'   Y - with new contents')
      CALL MSGDOC(M,
     *'   A - means to keep old contents and add new information')
      CALL MSGDOC(M,
     *'       with DOC-file program creates batch file: trahalo.bat')
      CALL DISPL_BL('_DOC:')
      READ(*,'(A)') LINE

      IBATCH = 0
      IF(LINE(1:4).EQ.'_DOC') THEN
        CALL LENSTR_BL(LINE,LEN)
        IF(LEN.GE.6) THEN
          DO I=6,LEN
            IF(LINE(I:I).NE.' ') THEN
              MSG=LINE(I:I)
              GO TO 10
            ENDIF
          ENDDO
        ENDIF
  10    CONTINUE
        IBATCH = 1
      ELSE
        MSG = LINE(1:1)
      ENDIF

      IF(MSG.EQ.'y') MSG = 'Y'
      IF(MSG.EQ.'a') MSG = 'A'
      IF(MSG.NE.'Y'.AND.MSG.NE.'A'.AND.MSG.NE.'#') MSG='N'

      IF(MSG.NE.'N'.OR.MSG.EQ.'#') THEN
        IF(MSG.EQ.'Y') MDOC = 999
        IF(MSG.EQ.'A') MDOC = 998
        IF(MSG.EQ.'#') MDOC = 997
        CALL START(0,MDOC,'trahalo')
      ELSE
        MDOC=0
      ENDIF
      IF(IBATCH.EQ.1) CALL SET_BATCH(IBATCH)
      CALL MSGDOC(M,' ')
C =============================================================
      MD   = 0
      M    = 99
 400  IERR = 0
      CALL MSGDOC(MD,'input file F_derivative')
      MTN   = 10
      ITYPE = 'FF  '
      LINE  = 'FILE_D'
      CALL ASK(LINE)
      NAMED=LINE
      CALL CORR_NAME_BL(NAMED)
      CALL ORFILE(MD,MTN,NAMED,ITYPE,NSYM,ISYM,IPRSYM,IERR)
      IF(IERR.NE.0) THEN
        CALL MSGERR(MD,' ERR: OPEN INPUT-FILE-Fder')
        GO TO 400
      ENDIF
      CALL GET_TITLE(MD,A,B,C,AL,BE,GA
     *  ,F1,F2,RESMIN,RESMAX,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
      CLOSE(MTN)
C --------------
      CALL MSGDOC(MD,' ')
 100  IERR = 0
      CALL MSGDOC(MD
     * ,'input file Fobs /for Anom. Diff. Patteson search press "CR"/')
      MTN   = 10
      ITYPE = 'FF  '
      LINE  = 'FILE_F'
      CALL ASK(LINE)
      NAMEF = LINE
      CALL LENSTR_BL(NAMEF,LEN)
      IF(LEN.GT.0.AND.NAMEF(1:1).NE.','.AND.NAMEF(1:1).NE.' ') THEN
        CALL CORR_NAME_BL(NAMEF)
        CALL ORFILE(MD,MTN,NAMEF,ITYPE,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) THEN
          CALL MSGERR(MD,' ERR: OPEN INPUT-FILE-FOBS')
          GO TO 100
        ENDIF
        CALL GET_TITLE(MD,A,B,C,AL,BE,GA
     *    ,F1,F2,RESMIN,RESMAX,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
        CLOSE(MTN)
      ENDIF
C -------------------------------------------------------------
 300  CONTINUE
      IERR   = 0
      NAMECPH= ' '
      NAMEHA = ' '
      RESN   = 0.0
      RESX   = 0.0
      BOFF   = 0.0
      BADD   = 0.0
      RAD    = 0.0
      NP     = 0
      CLIM   = 0.0
      SLIM   = 0.0
      NPST   = 0
      SPST(1,1) = 0.0
      SPST(2,1) = 0.0
      SPST(3,1) = 0.0
      CALL SET_DEFAULT(' ')
      CALL SET_DEFAULT('     Keywords:')
      CALL SET_DEFAULT(' ')
      CALL SET_DEFAULT(
     *'RESOL: MIN,MAX    - resolution <default from FILE_Fobs>')
      CALL SET_DEFAULT(
     *'NP:      <20>     - number of peaks to check, MAX =20, =6 for FAS
     *T')
      CALL SET_DEFAULT(
     *'GAUSS:   <Y>,N    - N means a sphere, Y - gaussian.')
      CALL SET_DEFAULT(
     *'RAD:     <2.5>    - radius of sphere of model of heavy atom or wi
     *dth')
      CALL SET_DEFAULT(
     *'                    of gaussian /in angstroms/.')
      CALL SET_DEFAULT(
     *'BOFF:    <400>    - 400 corresponds to RESmin=10A, BOFF=4*RESmin^
     *2')
      CALL SET_DEFAULT(
     *'BADD:     <0>       BOFF and BADD mean :')
      CALL SET_DEFAULT(
     *'                    !F!new = !F!input *EXP(-BADD*RSQ)*(1-EXP(-BOF
     *F*RSQ)')
      CALL SET_DEFAULT(
     *'                    BADD = -1 means automatical choice')
C      CALL SET_DEFAULT(
C     *'DIST:     <0>     - minimal distance between two peaks (in angstr
C     *om)') 
c      CALL SET_DEFAULT(
c     *'                    default is RESMAX/2, -1 means DIST = Effectiv
c     *e Resolution')
      CALL SET_DEFAULT(
     *'SORT:  P/D/<M>    - sort the list of atom positions according to:
     *')
      CALL SET_DEFAULT(
     *'                    P - phasing power, D - density , M - mixed.
     *')
      CALL SET_DEFAULT(
     *'SPEC:    <Y>/N    - Y means to use peaks in special  position.')
      CALL SET_DEFAULT(
     *'PACK:    <Y>/N    - N means Translation function without Packing 
     *function.') 
      CALL SET_DEFAULT(
     *'PST:     <N>/Y    - Y use pseudo-translation vector VPST')
      CALL SET_DEFAULT(
     *'VPST:  <0,0,0>    - vector of pseudo-translation /in fract./')
      CALL SET_DEFAULT(
     *'SLIM:    <1.0>    - limit for DEL=Fder-Fnat. !DEL!>sd(DEL)*SLIM w
     *ill used')
      CALL SET_DEFAULT(
     *'CLIM:    <3.0>    - limit for DEL=Fder-Fnat. !DEL!< DEL_MAX will
     *used')
      CALL SET_DEFAULT(
     *'                    DEL_MAX = <DEL> + DEL_rms * CLIM, 0 means no 
     *this limit')
      CALL SET_DEFAULT(
     *'TRANS:   <Y>/N    - Y means advanced translation function, N stan
     *dard')
      CALL SET_DEFAULT(
     *'FAST:    <F>/S    - F fast mode, automatical choice parameter for
     *FAST run:')
      CALL SET_DEFAULT(
     *'                    NP=6, TRANS="N", res_max=5')
      CALL SET_DEFAULT(
     *'SCALE:   <Y>/N    - N means without scaling Fnat and Fder.')
      CALL SET_DEFAULT(
     *'ANISO:   <N>/Y    - Y means anisoscaling Fnat and Fder.')
      CALL SET_DEFAULT(
     *'REMOV:   <N>/Y    - remove origin peak')
      CALL SET_DEFAULT(
     *'COMB:  N/A/<Y>    - Y use combined diff. patterson (D_iso + Kemp  
     ** D_ano)')
      CALL SET_DEFAULT(
     *'                  - A use Anomalous diff. patterson') 
      CALL SET_DEFAULT(
     *'FILE_A:   < >     - intput file of ABCD ')
C      CALL SET_DEFAULT(
C     *'FILE_T:   < >     - intput file atoms of ABCD ')
C      CALL SET_DEFAULT(
C     *'FILE_H:   < >     - intput file with known heavy atoms')
      CALL SET_DEFAULT(
     *'FILE_O:   name    - output file of heavy atoms, default: trahalo.
     *crd')
      CALL SET_DEFAULT(' ')
      CALL SET_DEFAULT(
     *'           trahalo_abc.dat - output file of ABCD_coefficients.')
      CALL SET_DEFAULT(
     *'           Use program ABCDPH to compute phases from ABCD or to c
     *ombine ABCD.')
C ---
      CALL SET_DEFAULT(' ')
      CALL SET_DEFAULT('?')
C -------------------------------------------
      LINE='RESOL'
      CALL ASK(LINE)
      READ(LINE,*) RESN,RESX
C ----
      LINE='BOFF'
      CALL ASK(LINE)
      READ(LINE,*) BOFF
C --
      LINE='BADD'
      CALL ASK(LINE)
      READ(LINE,*) BADD
C --
      LINE='NP'
      CALL ASK(LINE)
      READ(LINE,*) NP
C --
      LINE='RAD'
      CALL ASK(LINE)
      READ(LINE,*) RAD
C ----
      LINE='PST'
      CALL ASK(LINE)
      NPST = 0
      IF(LINE(1:1).EQ.'Y'.OR.LINE(1:1).EQ.'y') NPST = 1
C --
      LINE='VPST'
      CALL ASK(LINE)
      READ(LINE,*) SPST(1,1),SPST(2,1),SPST(3,1)
C --
      LINE='RAD'
      CALL ASK(LINE)
      READ(LINE,*) RAD
C ----
      LINE='SLIM'
      CALL ASK(LINE)
      READ(LINE,*) SLIM
      IF(LINE(1:1).EQ.',') SLIM=1.0
C --
      LINE='CLIM'
      CALL ASK(LINE)
      READ(LINE,*) CLIM
      IF(LINE(1:1).EQ.',') CLIM=3.0
C ----
      LINE='FILE_O'
      CALL ASK(LINE)
      NAMEO=LINE
      CALL CORR_NAME_BL(NAMEO)
C ----
C      LINE='FILE_T'
C      CALL ASK(LINE)
C      NAMECPH=LINE
C      CALL CORR_NAME_BL(NAMECPH)
C ----
      LINE='GAUSS'
      CALL ASK(LINE)
      GAUSS=LINE(1:1)
C ----
C      LINE='INVERT'
C      CALL ASK(LINE)
C      INVERT=LINE(1:1)
C ----
      LINE='FILE_A'
      CALL ASK(LINE)
      NAMEA=LINE
      CALL CORR_NAME_BL(NAMEA)
C ----
C      LINE='FILE_H'
C      CALL ASK(LINE)
C      NAMEHA=LINE
C      CALL CORR_NAME_BL(NAMEHA)
C ---
      LINE='SCALE'
      CALL ASK(LINE)
      MSCL=LINE(1:1)
C ---
      LINE='SORT'
      CALL ASK(LINE)
      SORT=LINE(1:1)
C ---
      LINE='SPEC'
      CALL ASK(LINE)
      SPEC=LINE(1:1)
C ---
      LINE='PACK'
      CALL ASK(LINE)
      PACK=LINE(1:1)
C ---
      LINE='REMOV'
      CALL ASK(LINE)
      REMOV=LINE(1:1)
C ---
      LINE='COMB'
      CALL ASK(LINE)
      COMB=LINE(1:1)
C     COMB='N'
C ---
      LINE='FAST'
      CALL ASK(LINE)
      FAST=LINE(1:1)
C --
      LINE='TRANS'
      CALL ASK(LINE)
      TRANS=LINE(1:1)
C -----
      LINE='ANISO'
      CALL ASK(LINE)
      ANISO=LINE(1:1)
C -----
C      GAUSS  ='N'
      INVERT ='N' 
C --
      STR_PARAM(1:1)   = INVERT
      STR_PARAM(2:2)   = MSCL
      STR_PARAM(3:3)   = SPEC
      STR_PARAM(4:4)   = PACK
      STR_PARAM(5:5)   = FAST
      STR_PARAM(6:6)   = TRANS
      STR_PARAM(7:7)   = GAUSS
      STR_PARAM(8:8)   = SORT
      STR_PARAM(9:9)   = REMOV
      STR_PARAM(10:10) = COMB
      STR_PARAM(11:11) = ANISO
      STR_PARAM(12:12) = ' '

      SUBR = 'N'

      STR_PARAM(13:13) = SUBR

C ------------------------------------------
C      CALL CLOSE_BATCH(M)
C -------------------------------------------
      CALL TRAHALO(MDOC
     * ,NAMEF,NAMED,NAMEO,NAMEA,NAMEHA,NAMECPH
     * ,RESN,RESX,CLIM,STR_PARAM,SLIM
     * ,NP,BOFF,BADD,RAD,POOL,MEMORY,ISYM,IPRSYM
     * ,NPST,SPST
     * ,RESULT,RESULTB,RESULTA,RESULTC
     * ,IPRMAX,IPRROW,IPRNPR,IPRMAXD,IERR)
C -------------
  900 CALL FINISH
      END

      SUBROUTINE TRAHALO(MDOC
     * ,NAMEF,NAMED,NAMEO,NAMEA,NAMEHA,NAMECPH
     * ,RESMINR,RESMAXR,CLIM_R,STR_PARAM,SLIMR
     * ,NP,BOFF,BADDR,RADR,POOL,MEMORY,ISYM,IPRSYM
     * ,NPSTR,VSPST
     * ,RESULT,RESULTB,RESULTA,RESULTC
     * ,IPRMAX,IPRROW,IPRNPR,IPRMAXD,IERR)

C --------------------------------------------------------------------
      INTEGER*2 ISYM(5,3,IPRSYM)
      REAL      POOL(MEMORY)
CMS$LARGE: POOL
C ----
C ---
C      COMMON/RESDER/ NO_DER,RES_DER
C      PARAMETER (IPRMAXDD = 71)
C      PARAMETER (IPRROWW  = 21)
C      PARAMETER (MAXDER   =  4)
C      REAL      RES_DER(IPRMAXDD,IPRROWW,MAXDER)
C ---  
C 
      REAL      RESULT (IPRMAX,IPRROW)
      REAL      RESULTP(21,21)
C 
      REAL      RESULTB(IPRMAXD,IPRROW)
      REAL      RESULTC(IPRMAXD,IPRROW)
      REAL      RESULTA(IPRMAXD,IPRROW,IPRNPR)
C ------------------------------------------------
      PARAMETER ( IPRPST=1)
      COMMON/COMPST/ NPST,SPST,RESEFF,BRES,BEFF
      REAL      SPST(3,IPRPST)
C -------------------------------
      REAL      VSPST(3,1)
C ------------------------------------------------
C ---
      PARAMETER ( NCRDMAX = 100 )
      PARAMETER ( IPRORG=24)
      REAL      SOR(3,IPRORG)
C ---
      PARAMETER ( IPRNEQ=24)
      INTEGER   INVERS1(IPRNEQ),IORIGN1(IPRNEQ),ISYMOP1(IPRNEQ)
      REAL      DISTN1(IPRNEQ)
      INTEGER   INVERS2(IPRNEQ),IORIGN2(IPRNEQ),ISYMOP2(IPRNEQ)
      REAL      DISTN2(IPRNEQ)
C ---
      REAL      RPARAM_TRAH(10),RHR(11)
C ---
      CHARACTER NAMEF*(*),NAMED*(*),NAMEO*(*),NAMEA*(*),NAMEHA*(*)
      CHARACTER NAMECPH*(*),NAME*80,SUBR*1,MSCL*1,ANISO*1
      CHARACTER LINE*80,STR_PARAM*(*),NAMEFS*80,NAMEDD*80,NAMEOA*80
      CHARACTER INVERT*1,SPEC*1,PACK*1,FAST*1,TRANS*1,GAUSS*1,SORT*1
     *          ,REMOV*1,COMB*1
C ---
      INCLUDE 'trahalo_version.fh'
C -------------------------------------------------------------
C  NAMEA   - intput file of ABCD 
C  NAMECPH - intput file atoms of ABCD 
C  NAMEHA  - intput file with known heavy atoms
C -------------------------------------------------------------
      IF(ABS(MDOC).GE.997) THEN
        MDOC_C=1
        IF(MDOC.LT.0) THEN
          MDOC=0
        ELSE      
C   
C         MDOC=999
C         IF(MSG.EQ.'Y') MDOC=999 do not keep old contents
C         IF(MSG.EQ.'A') MDOC=998 keep old contents
C         IF(MSG.EQ.'#') MDOC=997 batch mode, "_DOC " will be read
C          
C         for subroutine case:
C
C         MDOC=999 do not keep old contents
C         MDOC=998 keep old contents
C
C         MDOC=-999 without DOC-file
C
        ENDIF
        CALL START(0,MDOC,'trahalo')
      ELSE
        MDOC_C=0
      ENDIF
C =============================================================
      IF(MDOC.GT.0) THEN
        M=-1
        CALL MSGDOC(M,' ')
        CALL MSGDOC(M,VERSION)
        CALL MSGDOC(M,' ')
      ENDIF
      CALL CLOSE_BATCH(M)
C -------------------------------------------------------------
      IERR   = 0
      MD     = -ABS(MDOC)-1
      NO_DER = 0
C -----------------------------
      NPST      = 0
      SPST(1,1) = 0.0
      SPST(2,1) = 0.0
      SPST(3,1) = 0.0
C se
C      SPST(1,1) = 0.5
C      SPST(2,1) = 0.5
C      SPST(3,1) = 0.0
C      
      IF(NPSTR.GT.0) THEN
        CALL MSGDOC(MDOC,' ')
        CALL MSGDOC(MDOC,
     *  ' --- Pseudo-Translation vector will be used ---')
        CALL MSGDOC(MDOC,' ')
        NPST      = NPSTR
        SPST(1,1) = VSPST(1,1)
        SPST(2,1) = VSPST(2,1)
        SPST(3,1) = VSPST(3,1)
        WRITE(LINE,'(''  Vector:'',3F8.3)')
     *  SPST(1,1),SPST(2,1),SPST(3,1)
        CALL MSGDOC(MDOC,LINE)
      ENDIF
C ---------------------------------------------------
      INVERT = STR_PARAM(1:1)
      MSCL   = STR_PARAM(2:2)
      SPEC   = STR_PARAM(3:3)
      PACK   = STR_PARAM(4:4)
      FAST   = STR_PARAM(5:5)
      TRANS  = STR_PARAM(6:6)
      GAUSS  = STR_PARAM(7:7)
      SORT   = STR_PARAM(8:8)
      REMOV  = STR_PARAM(9:9)
      COMB   = STR_PARAM(10:10)
      ANISO  = STR_PARAM(11:11)
      SUBR   = STR_PARAM(13:13)

      IF(MSCL .EQ.'n') MSCL  = 'N'
      IF(MSCL .NE.'N') MSCL  = 'Y'
      IF(ANISO.EQ.'y') ANISO = 'Y'
      IF(ANISO.NE.'Y') ANISO = 'N'
      IF(FAST .EQ.'f') FAST  = 'F'
      IF(FAST .EQ.'m') FAST  = 'M'
      IF(FAST .EQ.'s') FAST  = 'S'
      IF(FAST .EQ.'y') FAST  = 'Y'
      IF(FAST .EQ.'n') FAST  = 'N'

      IF(FAST.NE.'F'.AND.FAST.NE.'M'.AND.FAST.NE.'S'.AND.
     *   FAST.NE.'Y'.AND.FAST.NE.'N') FAST = 'Y'

      IF(GAUSS.EQ.'n') GAUSS = 'N'
      IF(GAUSS.NE.'N') GAUSS = 'Y'

      IF(TRANS .EQ.'n') TRANS  = 'N'
      IF(TRANS .EQ.'y') TRANS  = 'Y'
      IF(TRANS .NE.'Y') TRANS  = 'N'
      IF(INVERT.EQ.'y') INVERT = 'Y'
      IF(INVERT.NE.'Y') INVERT = 'N'

      IF(REMOV.EQ.'y') REMOV = 'Y'
      IF(REMOV.NE.'Y') REMOV = 'N'

      IF(COMB.EQ.'n') COMB = 'N'
      IF(COMB.EQ.'a') COMB = 'A'
      IF(COMB.NE.'N'.AND.COMB.NE.'A') COMB = 'Y'

      IF(SORT.EQ.'d') SORT = 'D'
      IF(SORT.EQ.'p') SORT = 'P'
      IF(SORT.NE.'D'.AND.SORT.NE.'P') SORT = 'M'

      IF(PACK.EQ.'n') PACK = 'N'
      IF(PACK.NE.'N') PACK = 'Y'

      IF(SPEC.EQ.'n') SPEC = 'N'
      IF(SPEC.NE.'N') SPEC = 'Y'

      NAMEFS = 'trahalo_fsc'

      CALL FSCALE_TRH(MDOC,NAMEF,NAMED,NAMEFS,SUBR,ANISO,MSCL
     *  ,RESULT,IPRMAX,IPRROW,ISYM,IPRSYM,IERR)      
      IF(IERR.NE.0) GO TO 900                

      STR_PARAM(1:1)   = INVERT
      STR_PARAM(2:2)   = 'N'
      STR_PARAM(3:3)   = SPEC
      STR_PARAM(4:4)   = PACK
      STR_PARAM(5:5)   = FAST
      STR_PARAM(6:6)   = TRANS
      STR_PARAM(7:7)   = GAUSS
      STR_PARAM(8:8)   = SORT
      STR_PARAM(9:9)   = REMOV
      STR_PARAM(10:10) = COMB
      STR_PARAM(11:11) = ANISO
      STR_PARAM(12:12) = 'N'
      STR_PARAM(13:13) = SUBR

      SLIM_PATT      = SLIMR
      TIMES_RMS      = 6.0

      RPARAM_TRAH(1) = CLIM_R
      RPARAM_TRAH(2) = SLIMR 
      RPARAM_TRAH(3) = BOFF
      RPARAM_TRAH(4) = BADDR 
      RPARAM_TRAH(5) = RADR
      RPARAM_TRAH(6) = TIMES_RMS

      SLIM_SAD       = SLIM_PATT
      BADD_TRP       = 0.0
      BADD_FOUR      = 0.0

      RPARAM_TRAH( 7) = SLIM_SAD  
      RPARAM_TRAH( 8) = BADD_TRP 
      RPARAM_TRAH( 9) = BADD_FOUR 
      RPARAM_TRAH(10) = 0.0 

      NAMEDD =  NAMEFS
      NAMEOA =  'trahalo_abc'
C ---------------------------------------------------
      CALL TRAHALO_RUN(MDOC
     * ,NAMEF,NAMEDD,NAMEO,NAMEOA,NAMEA,NAMEHA,NAMECPH
     * ,RESMINR,RESMAXR,STR_PARAM,RPARAM_TRAH
     * ,NP,POOL,MEMORY,ISYM,IPRSYM
     * ,RESULT,RESULTP,RESULTB,RESULTA,RESULTC,N_FINAL
     * ,IPRMAX,IPRROW,IPRNPR,IPRMAXD,NCOMM,IERR)

      IF(IERR.EQ.0) THEN
        CALL MSGDOC(MD,' ')
        CALL MSGDOC(MD,
     *  'Lack of closure of error: E = rms(!!Fph(exp)!-!Fph(calc)!!)')
        CALL MSGDOC(MD,
     *  'rms heavy atoms         : rmsHA= rms(!Fh!)')
        CALL MSGDOC(MD,
     *  'Phasing power           : P = rmsHA/E')
        CALL MSGDOC(MDOC,' ')
        WRITE(LINE,'(''LOC            :'',F12.4)') RESULTP(5,1)
        CALL MSGDOC(MDOC,LINE)
        WRITE(LINE,'(''rmsHA          :'',F12.4)') RESULTP(6,1)
        CALL MSGDOC(MDOC,LINE)
        P = 0.0
        IF(RESULT(5,1).GT.0.0) P=RESULTP(6,1)/RESULTP(5,1)
        WRITE(LINE,'(''Power          :'',F10.3)') P
        CALL MSGDOC(MDOC,LINE)
        WRITE(LINE,'(''figure-of-merit:'',F10.3)') RESULTP(2,1)
        CALL MSGDOC(MDOC,LINE)

        DRANG  = RESULTP(9,1)
        RHR(1) = RESULTP(8,1)
        DO     I=2,11
          T      = 1./RHR(I-1)+DRANG
          RHR(I) = 1./T
        ENDDO
        CALL MSGDOC(MDOC,' ')
        WRITE(LINE,'(''R:'',F6.1,10(1X,F6.1))') RHR
        CALL MSGDOC(MDOC,LINE)
        WRITE(LINE,'(''Power'',10F7.2)') (RESULTP(10+K,1),K=1,10)
        CALL MSGDOC(MDOC,LINE)
        CALL MSGDOC(MDOC,' ')

      ENDIF
C -------------------------------------------------------------
 900   CONTINUE

       MT=10
       NAME='trahalo_cff.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=701)
       CLOSE(MT,STATUS='DELETE',ERR=701)
 701   CONTINUE

       NAME='trahalo.scr'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=702)
       CLOSE(MT,STATUS='DELETE',ERR=702)
 702   CONTINUE

       NAME='trahalo_dns.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=703)
       CLOSE(MT,STATUS='DELETE',ERR=703)
 703   CONTINUE

       NAME='trahalo_fft.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=704)
       CLOSE(MT,STATUS='DELETE',ERR=704)
 704   CONTINUE

       NAME='trahalo_fsc.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=705)
       CLOSE(MT,STATUS='DELETE',ERR=705)
 705   CONTINUE

       NAME='trahalo_ph.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=706)
       CLOSE(MT,STATUS='DELETE',ERR=706)
 706   CONTINUE

       NAME='trahalo_slf.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=707)
       CLOSE(MT,STATUS='DELETE',ERR=707)
 707   CONTINUE

       NAME='trahalo_trp.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=708)
       CLOSE(MT,STATUS='DELETE',ERR=708)
 708   CONTINUE

       NAME='scratch.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=608)
       CLOSE(MT,STATUS='DELETE',ERR=608)
 608   CONTINUE

       NAME='scratch.crd'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=609)
       CLOSE(MT,STATUS='DELETE',ERR=609)
 609   CONTINUE

       NAME='refine.crd'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=610)
       CLOSE(MT,STATUS='DELETE',ERR=610)
 610   CONTINUE

       NAME='phase_abc.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=611)
       CLOSE(MT,STATUS='DELETE',ERR=611)
 611   CONTINUE

       NAME='patt_pack.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=612)
       CLOSE(MT,STATUS='DELETE',ERR=612)
 612   CONTINUE

       MT=10
       NAME='trahalo_cffa.dat'
       OPEN(UNIT=MT,FILE=NAME,ACCESS='SEQUENTIAL',ERR=613)
       CLOSE(MT,STATUS='DELETE',ERR=613)
 613   CONTINUE
C -------------------------------------------------------------
      IF(MDOC_C.EQ.1) CALL FINISH
      RETURN      
      END

      SUBROUTINE SIRR(MDOC
     * ,NAMEF,NAMED,NAMEC,NAMEA,NAMEHA,NAMEO,NAMECPH
     * ,RESMINR,RESMAXR,PARAM_SIR,SUBR
     * ,NP,RADR,BOFF,BADD,POOL,MEMORY,ISYM,IPRSYM
     * ,RESULT,RESULTB,RESULTC,RESULTA,IPRMAX,IPRMAXD,IPRROW
     * ,IPRNPR,NCOMM,IPLEVEL,IERR)
C --------------------------------------------------------------------
C     NAMEA   - intput file of ABCD /for search and refinement/
C     NAMECPH - intput file atoms of ABCD 
C     NAMEHA  - intput file with known heavy atoms
C     NAMEC   - intput file of heavy atoms /to refine and calc phases only/
C     NAMEO   - output file of HL_coefficients.<sir_abc> 
C     sir.crd - output file of refined coordinates        
C --------------------------------------------------------------------
      INTEGER*2 ISYM(5,3,IPRSYM)
      REAL      POOL(MEMORY)
      CHARACTER NAMEF*(*),NAMED*(*),NAMEC*(*),NAMEO*(*),NAMEA*(*)
      CHARACTER NAMECPH*(*),NAMEHA*(*),SUBR*1,PARAM_SIR*(*)
      CHARACTER NAMEOA*80
C ----
C ---  
C      main_trahalo.f
c      PARAMETER  ( IPRMAXD = 71      )
c      PARAMETER  ( IPRMAX  = 21      )
c      PARAMETER  ( IPRROW  = 21      )
c      PARAMETER  ( IPRNPR  = 21      )
      REAL      RESULTP(21,21)

      REAL      RESULT (IPRMAX,IPRROW)
      REAL      RESULTB(IPRMAXD,IPRROW)
      REAL      RESULTC(IPRMAXD,IPRROW)
      REAL      RESULTA(IPRMAXD,IPRROW,IPRNPR)
      REAL      RHR(11)
C ------------------------------------------------
      PARAMETER ( IPRPST=1)
      COMMON/COMPST/ NPST,SPST,RESEFF,BRES,BEFF
      REAL      SPST(3,IPRPST)
      REAL      RPARAM_TRAH(10)      
C ------------------------------------------------
C ---
      COMMON/COMDIF/ NUMB_CFF,SLIM_PATT,BADD_COEF,DEL_LIM,DELA_LIM

C ---
C     !!!!!  IPRROWW must be iqual IPRROW
C ---
C      COMMON/RESDER/ NO_DER,RES_DER
C      PARAMETER (IPRMAXDD = 71)
C      PARAMETER (IPRROWW  = 21)
C      PARAMETER (MAXDER = 4)
C      REAL      RES_DER(IPRMAXDD,IPRROWW,MAXDER)
C ---
CMS$LARGE: POOL
      CHARACTER NAMEFS*80,NAMES2*80,FASTR*1
      CHARACTER NAMES*80,NAMER*80,NAMET*80,NAMEPH*80
      CHARACTER LINE*80,TTT*80,ITYPEF*4,ITYPED*4,ITYPE*4,COMB*1
      CHARACTER MODE*1,MISS*1,REST*1,RFAC*1,PATT*1,INVERT*1
      CHARACTER REFS*1,REFA*1,REFX*1,REFB*1,REFO*1,REFC*1
      CHARACTER STOP_T*1,GAUSS*1,ANOM*1,CENT*1,REF*1,LSORT*1
      CHARACTER SPEC*1,PACK*1,FAST*1,TRANS*1,RMOD*1
      CHARACTER MOD2*1,FUNC*1,ACCUM*1,REMOV*1,INVER*1,MSCL*1
      CHARACTER STR_PARAM*16,APATT*1,SOLV*1,ANISO*1,OPT1*1,SPEC_R*1
C -------------------------------------------------------------
      MD =-ABS(MDOC)-1
      M  = 99
C -------------------------
      CDEL_DEF = 3.0
      SLIM_DEF = 1.0
C     CDEL_DEF = 4.0
C     SLIM_DEF = 0.5
C -------------------------
      IERR   = 0
      MTN    = 0
      ANOM   = 'N'
      CALL LENSTR_BL(NAMEF,LEN)
      IF(LEN.GT.0.AND.NAMEF(1:1).NE.','.AND.NAMEF(1:1).NE.' ') THEN
        MTN    = 10
        ITYPEF = 'FF  '
        CALL ORFILE(M,MTN,NAMEF,ITYPEF,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) THEN
          CALL MSGERR(MDOC,' ERR: OPEN INPUT-FILE-FOBS')
          RETURN
        ENDIF
        CALL GET_TITLE(M,A,B,C,AL,BE,GA
     *    ,F1,F2,RN1,RX1,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
        CLOSE(MTN)
        CALL C_RD_CR     
      ELSE
        ANOM = 'Y'
      ENDIF 
C --------------------------------------------------------------
      IERR   = 0
      MTD    = 11
      ITYPED = 'FF  '
      CALL ORFILE(M,MTD,NAMED,ITYPED,NSYM,ISYM,IPRSYM,IERR)
      IF(IERR.NE.0) THEN
        CALL MSGERR(MDOC,' ERR: OPEN INPUT-FILE-FDER')
        RETURN
      ENDIF
      CALL GET_TITLE(M,A,B,C,AL,BE,GA
     *  ,F1,F2,RN,RX,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
      CLOSE(MTD)
      IF(MTN.EQ.0) THEN
        CALL C_RD_CR
        RN1 = RN
        RX1 = RX     
      ENDIF
C ---------------------------------------------------------------
      SLIM_PATT = 0.0
      BADD_COEF = 0.0
      DEL_LIM   = 0.0
      DELA_LIM  = 0.0
      RAD_DEF   = 2.5
      RES_DEF   = 3.0     
C ---------------------
      TRANS  = PARAM_SIR(1:1)
      PATT   = PARAM_SIR(2:2)
C     COMB   = N/A/Y
      COMB   = PARAM_SIR(3:3)
      CENT   = PARAM_SIR(4:4)
      REF    = PARAM_SIR(5:5)
      FASTR  = PARAM_SIR(6:6)
      RMOD   = PARAM_SIR(7:7)
C     mscl doesn't used
      MSCL   = PARAM_SIR(8:8)
      REMOV  = PARAM_SIR(9:9)
      ANISO  = PARAM_SIR(10:10)
C     == N
      OPT1   = PARAM_SIR(11:11)
      SPEC_R = PARAM_SIR(12:12)

      IF(COMB.NE.'N') ANOM = 'Y'

      RAD    = RADR

      RESMIN = RESMINR
      RESMAX = RESMAXR

      FAST = FASTR

C      IF(FAST.EQ.'M') FAST='N'
C      IF(FAST.EQ.'S') FAST='N'

      IF(FAST.NE.'F'.AND.FAST.NE.'Y') FAST = 'S'

      IF(FAST.EQ.'Y'.OR.FAST.EQ.'F') RES_DEF = 5.0

      RESMIN_D = RN1
      RESMAX_D = RX1

      IF(RN .LT.RESMIN_D)   RESMIN_D = RN
      IF(RX .GT.RESMAX_D)   RESMAX_D = RX
        
      IF((RESMIN_D.LT.RESMIN).OR.RESMIN.LE.0.0001) RESMIN = RESMIN_D
      IF(RESMAX_D.GT.RESMAX)                    RESMAX = RESMAX_D

      IF(RESMAX.LT.RES_DEF.AND.RESMAXR.LE.0.00001)  RESMAX = RES_DEF
      IF(RESMIN.GT.20.0   .AND.RESMINR.LE.0.00001)  RESMIN = 20.0

C ---
      IF(BOFF.EQ.0.0) BOFF = 400.0
      IF(BOFF.LT.0.0) BOFF = 0.0
      IF(BADD.LE.0.0) BADD = 0.0

      IF(FAST.EQ.'Y'.OR.FAST.EQ.'F') THEN
        IF(RESMAX.LT.RES_DEF.AND.RESMAXR.LE.0.0000001) RESMAX=RES_DEF
      ENDIF

      IF((FAST.EQ.'Y'.OR.FAST.EQ.'F').AND.NP.LE.0) THEN
        IF(NP.LE.0) NP = 6
      ELSE 
        IF(NP.LE.0) NP = 20
      ENDIF

      IF(NP.GT.IPRNPR) NP = IPRNPR

      IF(RAD.LE.0.0) RAD = RAD_DEF

C ---
      NAMES  = ' '
      NAMES2 = ' '
      NAMER  = 'sir'
      NAMEFS = 'sir_fsc'

      DO I=1,IPRMAXD
      DO J=1,IPRROW
        RESULTB(I,J) = 0.0
        RESULTC(I,J) = 0.0
        DO K=1,IPRNPR
          RESULTA(I,J,K) = 0.0
        ENDDO
      ENDDO
      ENDDO

      DO I=1,IPRMAX
      DO J=1,IPRROW
        RESULT (I,J) = 0.0
      ENDDO
      ENDDO

      DO I=1,21
      DO J=1,21
        RESULTP(I,J) = 0.0
      ENDDO
      ENDDO
      CALL LENSTR_BL(NAMEO,LEN)
      IF(LEN.GT.0.AND.NAMEO(1:1).NE.','.AND.NAMEO(1:1).NE.' ') THEN      

      ELSE
        NAMEO  ='sir_abc'
      ENDIF
C --------------------------------------------
      CALL LENSTR_BL(NAMEC,LEN)
      IF(LEN.GT.0.AND.NAMEC(1:1).NE.','.AND.NAMEC(1:1).NE.' ') THEN      
        NAMET = NAMEC
        MTC   = 12
        CALL ORCRD(M,MTC,NAMET,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) THEN
          CALL MSGERR(MDOC,' ERR: OPEN INPUT COORD_FILE')
          RETURN
        ENDIF
        CLOSE(MTC)
      ELSE

c        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-')
c        CALL MSGDOC(MDOC,'---  patt search  ---')
c        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-')
        GAUSS  = 'Y'
        INVERT = 'N'
C /////
C       LSORT  = 'P'
        LSORT  = 'M'
C ///
c        SPEC   = 'Y'
c        SPEC   = 'N'

        TIMES_RMS = 6.0


        SPEC = SPEC_R

        PACK   = 'Y'
        IF(FASTR.EQ.'M') TRANS='N'
        NAMET  = 'trahalo'     
        NAMEOA = 'trahalo_abc'     
        IF(PATT.EQ.'Y') THEN
          NAMET  = 'sir' 
          NAMEOA = NAMEO     
        ENDIF
C ---
        STR_PARAM(1:1)   = INVERT
C       mscl doesn't used
        STR_PARAM(2:2)   = MSCL
        STR_PARAM(3:3)   = SPEC
        STR_PARAM(4:4)   = PACK
        STR_PARAM(5:5)   = FAST
        STR_PARAM(6:6)   = TRANS
        STR_PARAM(7:7)   = GAUSS
        STR_PARAM(8:8)   = LSORT
        STR_PARAM(9:9)   = REMOV
C       comb = N/A/Y
        STR_PARAM(10:10) = COMB
        STR_PARAM(11:11) = ANISO
C        == N , remove HA of input abcd.
C        STR_PARAM(12:12) = OPT1
        STR_PARAM(12:12) = 'N'
        STR_PARAM(13:13) = SUBR

C /////
C CDEL_DEF times of drms : del_lim = <del> + t * drms
C /////

        TIMES_RMS = 6.0

        RPARAM_TRAH(1)  = CDEL_DEF
        RPARAM_TRAH(2)  = SLIM_DEF
        RPARAM_TRAH(3)  = BOFF
        RPARAM_TRAH(4)  = BADD
        RPARAM_TRAH(5)  = RAD
        RPARAM_TRAH(6)  = TIMES_RMS

        SLIM_SAD        = SLIM_DEF
        BADD_TRP        = 0.0
        BADD_FOUR       = 0.0

        RPARAM_TRAH( 7) = SLIM_SAD  
        RPARAM_TRAH( 8) = BADD_TRP 
        RPARAM_TRAH( 9) = BADD_FOUR 
        RPARAM_TRAH(10) = 0.0 
C -------------
C       level of messages
        IERR = IPLEVEL
C --
        CALL TRAHALO_RUN(MDOC
     *  ,NAMEF,NAMED,NAMET,NAMEOA,NAMEA,NAMEHA,NAMECPH
     *  ,RESMIN,RESMAX,STR_PARAM,RPARAM_TRAH
     *  ,NP,POOL,MEMORY,ISYM,IPRSYM
     *  ,RESULT,RESULTP,RESULTB,RESULTA,RESULTC,NPEAKS
     *  ,IPRMAX,IPRROW,IPRNPR,IPRMAXD,NCOMM,IERR)
C
C ierr= 13 - no self peaks
C ierr= 14 - no diff peaks
C ierr= 12 - no common peaks
C

c        NPEAKS = RESULTD(1,12)
        WRITE(LINE,'(''Number of heavy atoms :'',I5)') NPEAKS
        CALL MSGDOC(MDOC,LINE)
        IF(NPEAKS.LE.0) IERR = 1

        IF(PATT.EQ.'Y') THEN

          IF(NPEAKS.GT.0) THEN
      CALL MSGDOC(MD,' ')
      CALL MSGDOC(MD,
     *'Lack of closure of error: E = rms(!!Fph(exp)!-!Fph(calc)!!)')
      CALL MSGDOC(MD,
     *'rms heavy atoms         : rmsHA= rms(!Fh!)')
      CALL MSGDOC(MD,
     *'Phasing power           : P = rmsHA/E')
      CALL MSGDOC(MDOC,' ')
      WRITE(LINE,'(''LOC            :'',F12.4)') RESULTP(5,1)
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,'(''rmsHA          :'',F12.4)') RESULTP(6,1)
      CALL MSGDOC(MDOC,LINE)
            P = 0.0
            IF(RESULT(5,1).GT.0.0) P=RESULTP(6,1)/RESULTP(5,1)
            WRITE(LINE,'(''Power          :'',F10.3)') P
            CALL MSGDOC(MDOC,LINE)
            WRITE(LINE,'(''figure-of-merit:'',F10.3)') RESULTP(2,1)
            CALL MSGDOC(MDOC,LINE)

            DRANG  = RESULTP(9,1)
            RHR(1) = RESULTP(8,1)
            DO     I=2,11
              T      = 1./RHR(I-1)+DRANG
              RHR(I) = 1./T
            ENDDO
            CALL MSGDOC(MDOC,' ')
            WRITE(LINE,'(''R:'',F6.1,10(1X,F6.1))') RHR
            CALL MSGDOC(MDOC,LINE)
            WRITE(LINE,'(''Power'',10F7.2)') (RESULTP(10+K,1),K=1,10)
            CALL MSGDOC(MDOC,LINE)
            CALL MSGDOC(MDOC,' ')

            DO I=1,21
            DO J=1,21
              RESULT(I,J) = RESULTP(I,J)
            ENDDO
            ENDDO
          ENDIF
          GO TO 900
        ENDIF
        IF(IERR.NE.13.AND.IERR.NE.0) GO TO 900
        IERR = 0
      ENDIF
C --------------------------------------------------------------
      MTOUT = 13
      ITYPE = 'ABCD'
      CALL C_CR_WR      
      CALL SET_TITLE(M,A,B,C,AL,BE,GA
     *  ,F1,F2,RESMIN,RESMAX,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
      CALL OWFILE(M,MTOUT,NAMEO,ITYPE,NSYM,ISYM,IPRSYM,IERR)
      IF(IERR.NE.0) THEN
        CALL MSGERR(MDOC,' ERR: OPEN OUTPUT-FILE-ABCD')
        RETURN
      ENDIF
      CLOSE(MTOUT)

      CALL LENSTR_BL(NAMEF,LEN)
      IF(LEN.LE.0.OR.NAMEF(1:1).EQ.','.OR.NAMEF(1:1).EQ.' ') THEN
C        MSCL = 'N'
        ANOM = 'Y'
      ENDIF
C --------------------------------------------
      NAMEFS = NAMED
C -----
      IF(RMOD.EQ.'P') THEN
C ---
        CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDOC,'--- ABCD --> PH  ---')
        CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-')
        IPRINT = 0
        SCCOMB = 0.0
        STOP_T = 'N'
        NAMES  = ' '
        NAMES2 = ' '
        NAMEPH = 'sir_ph'
        CALL ABCDPH(M,NAMEA,NAMES,NAMEPH,NAMES2,SCCOMB
     *  ,RESMIN,RESMAX,IPRINT,STOP_T
     *  ,ISYM,IPRSYM,RESULT,IPRROW,IERR)
C       RESULT(1) - number of refls /in output file/
C       RESULT(2) - <FOM>
        IF(IERR.NE.0) GO TO 900
        WRITE(LINE,'(''figure-of-merit:'',F10.3)') RESULT(2,1)  
        CALL MSGDOC(MDOC,LINE)
C ---
      ELSE
        NAMEPH = ' '
      ENDIF
C --------------------------------------------
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*')
      LINE = '---  refine : mode "'//RMOD//'" ---'
      CALL MSGDOC(MDOC,LINE)
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*')
      IF(RMOD.EQ.'P') THEN
        IF(REF.EQ.'S') THEN
          REFS = 'Y'
          REFC = 'S'
          REFA = 'N'
          REFX = 'y'
          REFB = 'N'
          REFO = 'Y'
        ELSE IF(REF.EQ.'Y') THEN
          REFS = 'Y'
          REFC = 'N'
          REFA = 'N'
          REFX = 'y'
          REFB = 'N'
          REFO = 'Y'
        ELSE
          REFS = 'Y'
          REFC = 'N'
          REFA = 'N'
          REFX = 'y'
          REFB = 'Y'
          REFO = 'N'
        ENDIF
      ELSE
        IF(REF.EQ.'S') THEN
          REFS = 'N'
          REFC = 'S'
          REFA = 'N'
          REFX = 'y'
          REFB = 'N'
          REFO = 'Y'
        ELSE IF(REF.EQ.'Y') THEN
          REFS = 'N'
          REFC = 'N'
          REFA = 'N'
          REFX = 'y'
          REFB = 'N'
          REFO = 'Y'
        ELSE
          REFS = 'N'
          REFC = 'N'
          REFA = 'N'
          REFX = 'y'
          REFB = 'Y'
          REFO = 'N'
        ENDIF
      ENDIF
      STOP_T = 'N'
      NCYCL  = 5

C      SIGMIN = 0.0
C      DMAX   = 0.0
C      DAMAX  = 0.0
C      BBADD  = 0.0

      SIGMIN   = SLIM_PATT
      DMAX     = DEL_LIM
      DAMAX    = DELA_LIM
      BBADD    = BADD_COEF

      DMIN     = 0.0
      DAMIN    = 0.0

      BBOFF    = BOFF
      CALL REFINE(M,NAMET,NAMEF,NAMEFS,NAMER,NAMEPH
     * ,BBOFF,BBADD,ANOM,DMIN,DMAX,DAMIN,DAMAX
     * ,SIGMIN,RESMIN,RESMAX,NCYCL,CENT
     * ,REFS,REFC,REFA,REFX,REFB,REFO,STOP_T
     * ,RESULT,IPRROW,IERR)
      IF(IERR.NE.0) GO TO 900
C -------------------
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-')
      CALL MSGDOC(MDOC,'---  phase  ---')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-')

C      SIGMIN = 0.0
C      DMAX   = 0.0
C      DAMAX  = 0.0
C      BBADD  = 0.0

      SIGMIN   = SLIM_PATT
      DMAX     = DEL_LIM
      DAMAX    = DELA_LIM
      BBADD    = BADD_COEF

      DMIN     = 0.0
      DAMIN    = 0.0

      BBOFF    = BOFF
      CALL PHASE(M,NAMER,NAMEF,NAMEFS,NAMEO,NAMEPH
     * ,ANOM,BBOFF,BBADD,DMIN,DMAX,DAMIN,DAMAX,CENT
     * ,SIGMIN,RESMIN,RESMAX,STOP_T,RESULT,IPRROW,IERR)
      IF(IERR.NE.0) GO TO 900
C     RESULT(1) - number of refls /in output file/
C     RESULT(2) - <FOM>
C     RESULT(3) - R-factor
C     RESULT(4) - R_ST = SUM(!(Fph(exp)-Fph(calc)!)/SUM(!Fph(exp)-Fp!)
C     RESULT(5) - Lack of closure of error: E = rms(!!Fph(exp)!-!Fph(calc)!!)
C     RESULT(6) - rms heavy atom = rms(!Fh!)
C     RESULT(7)= NA number of atoms
c      RHMIN = 1.0/RESULT(8)
c      DRANG = RESULT(9)
c      RHR(1)= RESULT(8)
c      DO     I=2,11
c        T=1./RHR(I-1)+DRANG
c        RHR(I)=1./T
c      ENDDO
C      E     - (WRZNZ(K),K=1,10)
C      rmsHA - (FSZ(K)  ,K=1,10)
C      P=FSZ(I)/WRZNZ(I),I=1,10
C        RESULT(8)  = RESMIN
C        RESULT(9)  = DRANG
C        RESULT(10) = CORRCOEF
C        RESULT(10+I)=P,I=1,10

      CALL MSGDOC(MD,' ')
      CALL MSGDOC(MD,
     *'Lack of closure of error: E = rms(!!Fph(exp)!-!Fph(calc)!!)')
      CALL MSGDOC(MD,
     *'rms heavy atoms         : rmsHA= rms(!Fh!)')
      CALL MSGDOC(MD,
     *'Phasing power           : P = rmsHA/E')
      CALL MSGDOC(MDOC,' ')
      WRITE(LINE,'(''LOC            :'',F12.4)') RESULT(5,1)
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,'(''rmsHA          :'',F12.4)') RESULT(6,1)
      CALL MSGDOC(MDOC,LINE)
      P = 0.0
      IF(RESULT(5,1).GT.0.0) P=RESULT(6,1)/RESULT(5,1)
      WRITE(LINE,'(''Power          :'',F10.3)') P
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,'(''figure-of-merit:'',F10.3)') RESULT(2,1)
      CALL MSGDOC(MDOC,LINE)

      DRANG  = RESULT(9,1)
      RHR(1) = RESULT(8,1)
      DO     I=2,11
        T      = 1./RHR(I-1)+DRANG
        RHR(I) = 1./T
      ENDDO
      CALL MSGDOC(MDOC,' ')
      WRITE(LINE,'(''R:'',F6.1,10(1X,F6.1))') RHR
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,'(''Power'',10F7.2)') (RESULT(10+K,1),K=1,10)
      CALL MSGDOC(MDOC,LINE)
      CALL MSGDOC(MDOC,' ')
C
C --------------------------------------------
  900 RETURN
      END

      SUBROUTINE SIRR(MDOC
     * ,NAMEF,NAMED,NAMEC,NAMEA,NAMEHA,NAMEO,NAMECPH
     * ,RESMINR,RESMAXR,PARAM_SIR,SUBR
     * ,NP,RADR,BOFF,BADD,POOL,MEMORY,ISYM,IPRSYM
     * ,RESULT,RESULTB,RESULTC,RESULTA,IPRMAX,IPRMAXD,IPRROW
     * ,IPRNPR,NCOMM,IPLEVEL,IERR)
C --------------------------------------------------------------------
C     NAMEA   - intput file of ABCD /for search and refinement/
C     NAMECPH - intput file atoms of ABCD 
C     NAMEHA  - intput file with known heavy atoms
C     NAMEC   - intput file of heavy atoms /to refine and calc phases only/
C     NAMEO   - output file of HL_coefficients.<sir_abc> 
C     sir.crd - output file of refined coordinates        
C --------------------------------------------------------------------
      INTEGER*2 ISYM(5,3,IPRSYM)
      REAL      POOL(MEMORY)
      CHARACTER NAMEF*(*),NAMED*(*),NAMEC*(*),NAMEO*(*),NAMEA*(*)
      CHARACTER NAMECPH*(*),NAMEHA*(*),SUBR*1,PARAM_SIR*(*)
      CHARACTER NAMEOA*80
C ----
C ---  
C      main_trahalo.f
c      PARAMETER  ( IPRMAXD = 71      )
c      PARAMETER  ( IPRMAX  = 21      )
c      PARAMETER  ( IPRROW  = 21      )
c      PARAMETER  ( IPRNPR  = 21      )
      REAL      RESULTP(21,21)

      REAL      RESULT (IPRMAX,IPRROW)
      REAL      RESULTB(IPRMAXD,IPRROW)
      REAL      RESULTC(IPRMAXD,IPRROW)
      REAL      RESULTA(IPRMAXD,IPRROW,IPRNPR)
      REAL      RHR(11)
C ------------------------------------------------
      PARAMETER ( IPRPST=1)
      COMMON/COMPST/ NPST,SPST,RESEFF,BRES,BEFF
      REAL      SPST(3,IPRPST)
      REAL      RPARAM_TRAH(10)      
C ------------------------------------------------
C ---
      COMMON/COMDIF/ NUMB_CFF,SLIM_PATT,BADD_COEF,DEL_LIM,DELA_LIM

C ---
C     !!!!!  IPRROWW must be iqual IPRROW
C ---
C      COMMON/RESDER/ NO_DER,RES_DER
C      PARAMETER (IPRMAXDD = 71)
C      PARAMETER (IPRROWW  = 21)
C      PARAMETER (MAXDER = 4)
C      REAL      RES_DER(IPRMAXDD,IPRROWW,MAXDER)
C ---
CMS$LARGE: POOL
      CHARACTER NAMEFS*80,NAMES2*80,FASTR*1
      CHARACTER NAMES*80,NAMER*80,NAMET*80,NAMEPH*80
      CHARACTER LINE*80,TTT*80,ITYPEF*4,ITYPED*4,ITYPE*4,COMB*1
      CHARACTER MODE*1,MISS*1,REST*1,RFAC*1,PATT*1,INVERT*1
      CHARACTER REFS*1,REFA*1,REFX*1,REFB*1,REFO*1,REFC*1
      CHARACTER STOP_T*1,GAUSS*1,ANOM*1,CENT*1,REF*1,LSORT*1
      CHARACTER SPEC*1,PACK*1,FAST*1,TRANS*1,RMOD*1
      CHARACTER MOD2*1,FUNC*1,ACCUM*1,REMOV*1,INVER*1,MSCL*1
      CHARACTER STR_PARAM*16,APATT*1,SOLV*1,ANISO*1,OPT1*1,SPEC_R*1
C -------------------------------------------------------------
      MD =-ABS(MDOC)-1
      M  = 99
C -------------------------
      CDEL_DEF = 3.0
      SLIM_DEF = 1.0
C     CDEL_DEF = 4.0
C     SLIM_DEF = 0.5
C -------------------------
      IERR   = 0
      MTN    = 0
      ANOM   = 'N'
      CALL LENSTR_BL(NAMEF,LEN)
      IF(LEN.GT.0.AND.NAMEF(1:1).NE.','.AND.NAMEF(1:1).NE.' ') THEN
        MTN    = 10
        ITYPEF = 'FF  '
        CALL ORFILE(M,MTN,NAMEF,ITYPEF,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) THEN
          CALL MSGERR(MDOC,' ERR: OPEN INPUT-FILE-FOBS')
          RETURN
        ENDIF
        CALL GET_TITLE(M,A,B,C,AL,BE,GA
     *    ,F1,F2,RN1,RX1,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
        CLOSE(MTN)
        CALL C_RD_CR     
      ELSE
        ANOM = 'Y'
      ENDIF 
C --------------------------------------------------------------
      IERR   = 0
      MTD    = 11
      ITYPED = 'FF  '
      CALL ORFILE(M,MTD,NAMED,ITYPED,NSYM,ISYM,IPRSYM,IERR)
      IF(IERR.NE.0) THEN
        CALL MSGERR(MDOC,' ERR: OPEN INPUT-FILE-FDER')
        RETURN
      ENDIF
      CALL GET_TITLE(M,A,B,C,AL,BE,GA
     *  ,F1,F2,RN,RX,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
      CLOSE(MTD)
      IF(MTN.EQ.0) THEN
        CALL C_RD_CR
        RN1 = RN
        RX1 = RX     
      ENDIF
C ---------------------------------------------------------------
      SLIM_PATT = 0.0
      BADD_COEF = 0.0
      DEL_LIM   = 0.0
      DELA_LIM  = 0.0
      RAD_DEF   = 2.5
      RES_DEF   = 3.0     
C ---------------------
      TRANS  = PARAM_SIR(1:1)
      PATT   = PARAM_SIR(2:2)
C     COMB   = N/A/Y
      COMB   = PARAM_SIR(3:3)
      CENT   = PARAM_SIR(4:4)
      REF    = PARAM_SIR(5:5)
      FASTR  = PARAM_SIR(6:6)
      RMOD   = PARAM_SIR(7:7)
C     mscl doesn't used
      MSCL   = PARAM_SIR(8:8)
      REMOV  = PARAM_SIR(9:9)
      ANISO  = PARAM_SIR(10:10)
C     == N
      OPT1   = PARAM_SIR(11:11)
      SPEC_R = PARAM_SIR(12:12)

      IF(COMB.NE.'N') ANOM = 'Y'

      RAD    = RADR

      RESMIN = RESMINR
      RESMAX = RESMAXR

      FAST = FASTR

C      IF(FAST.EQ.'M') FAST='N'
C      IF(FAST.EQ.'S') FAST='N'

      IF(FAST.NE.'F'.AND.FAST.NE.'Y') FAST = 'S'

      IF(FAST.EQ.'Y'.OR.FAST.EQ.'F') RES_DEF = 5.0

      RESMIN_D = RN1
      RESMAX_D = RX1

      IF(RN .LT.RESMIN_D)   RESMIN_D = RN
      IF(RX .GT.RESMAX_D)   RESMAX_D = RX
        
      IF((RESMIN_D.LT.RESMIN).OR.RESMIN.LE.0.0001) RESMIN = RESMIN_D
      IF(RESMAX_D.GT.RESMAX)                    RESMAX = RESMAX_D

      IF(RESMAX.LT.RES_DEF.AND.RESMAXR.LE.0.00001)  RESMAX = RES_DEF
      IF(RESMIN.GT.20.0   .AND.RESMINR.LE.0.00001)  RESMIN = 20.0

C ---
      IF(BOFF.EQ.0.0) BOFF = 400.0
      IF(BOFF.LT.0.0) BOFF = 0.0
      IF(BADD.LE.0.0) BADD = 0.0

      IF(FAST.EQ.'Y'.OR.FAST.EQ.'F') THEN
        IF(RESMAX.LT.RES_DEF.AND.RESMAXR.LE.0.0000001) RESMAX=RES_DEF
      ENDIF

      IF((FAST.EQ.'Y'.OR.FAST.EQ.'F').AND.NP.LE.0) THEN
        IF(NP.LE.0) NP = 6
      ELSE 
        IF(NP.LE.0) NP = 20
      ENDIF

      IF(NP.GT.IPRNPR) NP = IPRNPR

      IF(RAD.LE.0.0) RAD = RAD_DEF

C ---
      NAMES  = ' '
      NAMES2 = ' '
      NAMER  = 'sir'
      NAMEFS = 'sir_fsc'

      DO I=1,IPRMAXD
      DO J=1,IPRROW
        RESULTB(I,J) = 0.0
        RESULTC(I,J) = 0.0
        DO K=1,IPRNPR
          RESULTA(I,J,K) = 0.0
        ENDDO
      ENDDO
      ENDDO

      DO I=1,IPRMAX
      DO J=1,IPRROW
        RESULT (I,J) = 0.0
      ENDDO
      ENDDO

      DO I=1,21
      DO J=1,21
        RESULTP(I,J) = 0.0
      ENDDO
      ENDDO
      CALL LENSTR_BL(NAMEO,LEN)
      IF(LEN.GT.0.AND.NAMEO(1:1).NE.','.AND.NAMEO(1:1).NE.' ') THEN      

      ELSE
        NAMEO  ='sir_abc'
      ENDIF
C --------------------------------------------
      CALL LENSTR_BL(NAMEC,LEN)
      IF(LEN.GT.0.AND.NAMEC(1:1).NE.','.AND.NAMEC(1:1).NE.' ') THEN      
        NAMET = NAMEC
        MTC   = 12
        CALL ORCRD(M,MTC,NAMET,NSYM,ISYM,IPRSYM,IERR)
        IF(IERR.NE.0) THEN
          CALL MSGERR(MDOC,' ERR: OPEN INPUT COORD_FILE')
          RETURN
        ENDIF
        CLOSE(MTC)
      ELSE

c        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-')
c        CALL MSGDOC(MDOC,'---  patt search  ---')
c        CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-')
        GAUSS  = 'Y'
        INVERT = 'N'
C /////
C       LSORT  = 'P'
        LSORT  = 'M'
C ///
c        SPEC   = 'Y'
c        SPEC   = 'N'

        TIMES_RMS = 6.0


        SPEC = SPEC_R

        PACK   = 'Y'
        IF(FASTR.EQ.'M') TRANS='N'
        NAMET  = 'trahalo'     
        NAMEOA = 'trahalo_abc'     
        IF(PATT.EQ.'Y') THEN
          NAMET  = 'sir' 
          NAMEOA = NAMEO     
        ENDIF
C ---
        STR_PARAM(1:1)   = INVERT
C       mscl doesn't used
        STR_PARAM(2:2)   = MSCL
        STR_PARAM(3:3)   = SPEC
        STR_PARAM(4:4)   = PACK
        STR_PARAM(5:5)   = FAST
        STR_PARAM(6:6)   = TRANS
        STR_PARAM(7:7)   = GAUSS
        STR_PARAM(8:8)   = LSORT
        STR_PARAM(9:9)   = REMOV
C       comb = N/A/Y
        STR_PARAM(10:10) = COMB
        STR_PARAM(11:11) = ANISO
C        == N , remove HA of input abcd.
C        STR_PARAM(12:12) = OPT1
        STR_PARAM(12:12) = 'N'
        STR_PARAM(13:13) = SUBR

C /////
C CDEL_DEF times of drms : del_lim = <del> + t * drms
C /////

        TIMES_RMS = 6.0

        RPARAM_TRAH(1)  = CDEL_DEF
        RPARAM_TRAH(2)  = SLIM_DEF
        RPARAM_TRAH(3)  = BOFF
        RPARAM_TRAH(4)  = BADD
        RPARAM_TRAH(5)  = RAD
        RPARAM_TRAH(6)  = TIMES_RMS

        SLIM_SAD        = SLIM_DEF
        BADD_TRP        = 0.0
        BADD_FOUR       = 0.0

        RPARAM_TRAH( 7) = SLIM_SAD  
        RPARAM_TRAH( 8) = BADD_TRP 
        RPARAM_TRAH( 9) = BADD_FOUR 
        RPARAM_TRAH(10) = 0.0 
C -------------
C       level of messages
        IERR = IPLEVEL
C --
        CALL TRAHALO_RUN(MDOC
     *  ,NAMEF,NAMED,NAMET,NAMEOA,NAMEA,NAMEHA,NAMECPH
     *  ,RESMIN,RESMAX,STR_PARAM,RPARAM_TRAH
     *  ,NP,POOL,MEMORY,ISYM,IPRSYM
     *  ,RESULT,RESULTP,RESULTB,RESULTA,RESULTC,NPEAKS
     *  ,IPRMAX,IPRROW,IPRNPR,IPRMAXD,NCOMM,IERR)
C
C ierr= 13 - no self peaks
C ierr= 14 - no diff peaks
C ierr= 12 - no common peaks
C

c        NPEAKS = RESULTD(1,12)
        WRITE(LINE,'(''Number of heavy atoms :'',I5)') NPEAKS
        CALL MSGDOC(MDOC,LINE)
        IF(NPEAKS.LE.0) IERR = 1

        IF(PATT.EQ.'Y') THEN

          IF(NPEAKS.GT.0) THEN
      CALL MSGDOC(MD,' ')
      CALL MSGDOC(MD,
     *'Lack of closure of error: E = rms(!!Fph(exp)!-!Fph(calc)!!)')
      CALL MSGDOC(MD,
     *'rms heavy atoms         : rmsHA= rms(!Fh!)')
      CALL MSGDOC(MD,
     *'Phasing power           : P = rmsHA/E')
      CALL MSGDOC(MDOC,' ')
      WRITE(LINE,'(''LOC            :'',F12.4)') RESULTP(5,1)
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,'(''rmsHA          :'',F12.4)') RESULTP(6,1)
      CALL MSGDOC(MDOC,LINE)
            P = 0.0
            IF(RESULT(5,1).GT.0.0) P=RESULTP(6,1)/RESULTP(5,1)
            WRITE(LINE,'(''Power          :'',F10.3)') P
            CALL MSGDOC(MDOC,LINE)
            WRITE(LINE,'(''figure-of-merit:'',F10.3)') RESULTP(2,1)
            CALL MSGDOC(MDOC,LINE)

            DRANG  = RESULTP(9,1)
            RHR(1) = RESULTP(8,1)
            DO     I=2,11
              T      = 1./RHR(I-1)+DRANG
              RHR(I) = 1./T
            ENDDO
            CALL MSGDOC(MDOC,' ')
            WRITE(LINE,'(''R:'',F6.1,10(1X,F6.1))') RHR
            CALL MSGDOC(MDOC,LINE)
            WRITE(LINE,'(''Power'',10F7.2)') (RESULTP(10+K,1),K=1,10)
            CALL MSGDOC(MDOC,LINE)
            CALL MSGDOC(MDOC,' ')

            DO I=1,21
            DO J=1,21
              RESULT(I,J) = RESULTP(I,J)
            ENDDO
            ENDDO
          ENDIF
          GO TO 900
        ENDIF
        IF(IERR.NE.13.AND.IERR.NE.0) GO TO 900
        IERR = 0
      ENDIF
C --------------------------------------------------------------
      MTOUT = 13
      ITYPE = 'ABCD'
      CALL C_CR_WR      
      CALL SET_TITLE(M,A,B,C,AL,BE,GA
     *  ,F1,F2,RESMIN,RESMAX,IDM,IDM,IDM,IDM,IDM,TTT,IERR)
      CALL OWFILE(M,MTOUT,NAMEO,ITYPE,NSYM,ISYM,IPRSYM,IERR)
      IF(IERR.NE.0) THEN
        CALL MSGERR(MDOC,' ERR: OPEN OUTPUT-FILE-ABCD')
        RETURN
      ENDIF
      CLOSE(MTOUT)

      CALL LENSTR_BL(NAMEF,LEN)
      IF(LEN.LE.0.OR.NAMEF(1:1).EQ.','.OR.NAMEF(1:1).EQ.' ') THEN
C        MSCL = 'N'
        ANOM = 'Y'
      ENDIF
C --------------------------------------------
      NAMEFS = NAMED
C -----
      IF(RMOD.EQ.'P') THEN
C ---
        CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-')
        CALL MSGDOC(MDOC,'--- ABCD --> PH  ---')
        CALL MSGDOC(MD  ,'*-*-*-*-*-*-*-*-*-*-')
        IPRINT = 0
        SCCOMB = 0.0
        STOP_T = 'N'
        NAMES  = ' '
        NAMES2 = ' '
        NAMEPH = 'sir_ph'
        CALL ABCDPH(M,NAMEA,NAMES,NAMEPH,NAMES2,SCCOMB
     *  ,RESMIN,RESMAX,IPRINT,STOP_T
     *  ,ISYM,IPRSYM,RESULT,IPRROW,IERR)
C       RESULT(1) - number of refls /in output file/
C       RESULT(2) - <FOM>
        IF(IERR.NE.0) GO TO 900
        WRITE(LINE,'(''figure-of-merit:'',F10.3)') RESULT(2,1)  
        CALL MSGDOC(MDOC,LINE)
C ---
      ELSE
        NAMEPH = ' '
      ENDIF
C --------------------------------------------
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*')
      LINE = '---  refine : mode "'//RMOD//'" ---'
      CALL MSGDOC(MDOC,LINE)
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-*-*-*-*-*-*')
      IF(RMOD.EQ.'P') THEN
        IF(REF.EQ.'S') THEN
          REFS = 'Y'
          REFC = 'S'
          REFA = 'N'
          REFX = 'y'
          REFB = 'N'
          REFO = 'Y'
        ELSE IF(REF.EQ.'Y') THEN
          REFS = 'Y'
          REFC = 'N'
          REFA = 'N'
          REFX = 'y'
          REFB = 'N'
          REFO = 'Y'
        ELSE
          REFS = 'Y'
          REFC = 'N'
          REFA = 'N'
          REFX = 'y'
          REFB = 'Y'
          REFO = 'N'
        ENDIF
      ELSE
        IF(REF.EQ.'S') THEN
          REFS = 'N'
          REFC = 'S'
          REFA = 'N'
          REFX = 'y'
          REFB = 'N'
          REFO = 'Y'
        ELSE IF(REF.EQ.'Y') THEN
          REFS = 'N'
          REFC = 'N'
          REFA = 'N'
          REFX = 'y'
          REFB = 'N'
          REFO = 'Y'
        ELSE
          REFS = 'N'
          REFC = 'N'
          REFA = 'N'
          REFX = 'y'
          REFB = 'Y'
          REFO = 'N'
        ENDIF
      ENDIF
      STOP_T = 'N'
      NCYCL  = 5

C      SIGMIN = 0.0
C      DMAX   = 0.0
C      DAMAX  = 0.0
C      BBADD  = 0.0

      SIGMIN   = SLIM_PATT
      DMAX     = DEL_LIM
      DAMAX    = DELA_LIM
      BBADD    = BADD_COEF

      DMIN     = 0.0
      DAMIN    = 0.0

      BBOFF    = BOFF
      CALL REFINE(M,NAMET,NAMEF,NAMEFS,NAMER,NAMEPH
     * ,BBOFF,BBADD,ANOM,DMIN,DMAX,DAMIN,DAMAX
     * ,SIGMIN,RESMIN,RESMAX,NCYCL,CENT
     * ,REFS,REFC,REFA,REFX,REFB,REFO,STOP_T
     * ,RESULT,IPRROW,IERR)
      IF(IERR.NE.0) GO TO 900
C -------------------
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-')
      CALL MSGDOC(MDOC,'---  phase  ---')
      CALL MSGDOC(MD  ,'-*-*-*-*-*-*-*-')

C      SIGMIN = 0.0
C      DMAX   = 0.0
C      DAMAX  = 0.0
C      BBADD  = 0.0

      SIGMIN   = SLIM_PATT
      DMAX     = DEL_LIM
      DAMAX    = DELA_LIM
      BBADD    = BADD_COEF

      DMIN     = 0.0
      DAMIN    = 0.0

      BBOFF    = BOFF
      CALL PHASE(M,NAMER,NAMEF,NAMEFS,NAMEO,NAMEPH
     * ,ANOM,BBOFF,BBADD,DMIN,DMAX,DAMIN,DAMAX,CENT
     * ,SIGMIN,RESMIN,RESMAX,STOP_T,RESULT,IPRROW,IERR)
      IF(IERR.NE.0) GO TO 900
C     RESULT(1) - number of refls /in output file/
C     RESULT(2) - <FOM>
C     RESULT(3) - R-factor
C     RESULT(4) - R_ST = SUM(!(Fph(exp)-Fph(calc)!)/SUM(!Fph(exp)-Fp!)
C     RESULT(5) - Lack of closure of error: E = rms(!!Fph(exp)!-!Fph(calc)!!)
C     RESULT(6) - rms heavy atom = rms(!Fh!)
C     RESULT(7)= NA number of atoms
c      RHMIN = 1.0/RESULT(8)
c      DRANG = RESULT(9)
c      RHR(1)= RESULT(8)
c      DO     I=2,11
c        T=1./RHR(I-1)+DRANG
c        RHR(I)=1./T
c      ENDDO
C      E     - (WRZNZ(K),K=1,10)
C      rmsHA - (FSZ(K)  ,K=1,10)
C      P=FSZ(I)/WRZNZ(I),I=1,10
C        RESULT(8)  = RESMIN
C        RESULT(9)  = DRANG
C        RESULT(10) = CORRCOEF
C        RESULT(10+I)=P,I=1,10

      CALL MSGDOC(MD,' ')
      CALL MSGDOC(MD,
     *'Lack of closure of error: E = rms(!!Fph(exp)!-!Fph(calc)!!)')
      CALL MSGDOC(MD,
     *'rms heavy atoms         : rmsHA= rms(!Fh!)')
      CALL MSGDOC(MD,
     *'Phasing power           : P = rmsHA/E')
      CALL MSGDOC(MDOC,' ')
      WRITE(LINE,'(''LOC            :'',F12.4)') RESULT(5,1)
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,'(''rmsHA          :'',F12.4)') RESULT(6,1)
      CALL MSGDOC(MDOC,LINE)
      P = 0.0
      IF(RESULT(5,1).GT.0.0) P=RESULT(6,1)/RESULT(5,1)
      WRITE(LINE,'(''Power          :'',F10.3)') P
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,'(''figure-of-merit:'',F10.3)') RESULT(2,1)
      CALL MSGDOC(MDOC,LINE)

      DRANG  = RESULT(9,1)
      RHR(1) = RESULT(8,1)
      DO     I=2,11
        T      = 1./RHR(I-1)+DRANG
        RHR(I) = 1./T
      ENDDO
      CALL MSGDOC(MDOC,' ')
      WRITE(LINE,'(''R:'',F6.1,10(1X,F6.1))') RHR
      CALL MSGDOC(MDOC,LINE)
      WRITE(LINE,'(''Power'',10F7.2)') (RESULT(10+K,1),K=1,10)
      CALL MSGDOC(MDOC,LINE)
      CALL MSGDOC(MDOC,' ')
C
C --------------------------------------------
  900 RETURN
      END
