PROGRAM MINIMIZE_categ

! **********************************************************************
! *              .      .    .
! * PROGRAM:    MINIMIZE     
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2001-06-01
! *
! * ABSTRACT:   MINIMIZATION PACKAGE, WITH PRECONDITIONING
! *
! * 
! * OPTIONS:
! *
! *    (1) STEEPEST DESCENT,
! *
! *    (2) CONJUGATE GRADIENT (FLETCHER-REEVES, POLAK-RIBIERE)
! *
! *    (2) L-BFGS QUASI-NEWTON METHOD
! *        J. NOCEDAL (1980), WITH RESTART BY SHANNO (1985)
! *
! * PROGRAM LOG:
! *
! * 06/01/2001    ..... M. ZUPANSKI: Original version, variational data assimilation
! * 09/04/2003    ..... M. ZUPANSKI: Adopted for use with ensemble data assimilation
! * 09/04/2003    ..... M. ZUPANSKI: Added CG option
! *
! **********************************************************************
!   INPUT:      gradient (g)
!  OUTPUT:      descent direction (d)
!------------------------------------------------------------------------

      integer,parameter::isha=20       ! restart and iteration parameters
      integer,parameter::igrad=21      ! gradient file
      integer,parameter::idesc=51      ! descent direction file
  
!---define allocatable areas--------------------------------------------
      real,dimension(:),allocatable	:: g,d

      real 		g_inner,g_norm
      real 		d_inner,d_norm

      real 		SGITOT,GAMATOT,GAMANG
      integer 		ITERNO,ITEROLD

      integer 		iout,ioutm1,index
      integer 		NPAIRS,ioutmax
      integer 		NENS
      integer 		NENS_START
      integer 		NALL,ncateg
      character*5 	MINIM_ALG

      character*12 	gradname
      character*11 	descname

      NAMELIST /CURRENT_ITERATION/ iout,ioutm1,index
      NAMELIST /MINIMIZATION/ NPAIRS,ioutmax,MINIM_ALG
      NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS

!==============start calculation===================

       write(*,*) "start minimize"
       write(3,*) "start minimize"

!-- read current iteration 

       REWIND 13
       READ(13,CURRENT_ITERATION)

          write(*,*) "minimize: iout,ioutm1,index=",iout,ioutm1,index

!-- read minimization information

       REWIND 14
       READ(14,MINIMIZATION)

          write(*,*) "minimize: NPAIRS,ioutmax,MINIM_ALG=",NPAIRS,ioutmax,MINIM_ALG

!-- read ensemble size

       REWIND 15
       READ(15,ENSEMBLE_SIZE)

          write(*,*) "minimize: NENS_START, NENS=",NENS_START,NENS

!-- read number of categories of matrix A

      rewind 17
      read(17,*) ncateg

!- Total number of degrees of freedom
      NALL=NENS-NENS_START+1

      N_start=NENS_START
      N_end=NALL*ncateg
!-------------------------------------------------------------
        allocate(d(N_start:N_end))
        d=0.0
!---------------------------------

   IF(iout.eq.1) THEN

!-- first iteration
        ITERNO=1
        SGITOT=0.0

   ELSE

!--- read restart parameters

      close(isha)
      OPEN(UNIT=isha,FILE='shanno',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' ISHA OPEN UNIT ERROR IER=',IER
      read(isha) SGITOT,ITERNO

      ITERNO=ITERNO+1

   END IF

!--------------------------------------

        write(3,*) "SGITOT=",SGITOT
        write(3,*) "ITERNO=",ITERNO
        write(*,*) "SGITOT=",SGITOT
        write(*,*) "ITERNO=",ITERNO

 100  format('gradient_',i2.2)
 200  format('descent_',i2.2)
!-------------------------
! Read new gradient

        allocate(g(N_start:N_end))
        g=0.0

      write(gradname,100) iout 

      close(igrad)
      open(UNIT=igrad,FILE=gradname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*)' igrad OPEN UNIT ERROR IER=',IER
      read(igrad) (g(k),k=N_start,N_end)
      close(igrad,status='KEEP')

      ij=N_start-1
             do j=1,ncateg
            do i=NENS_START,NENS
      ij=ij+1
      if(j.eq.min(5,ncateg)) then
      write(*,*) " i=",i," j=",j," g(ij)=",g(ij)
      endif
            end do
             end do

!--- gradient norm and inner product ---------------------
                               call norm &
      (N_start,N_end,g,g_inner,g_norm)

        write(*,*) "CURRENT GRADIENT NORM sqrt(g*g): g_norm=",g_norm
        write(3,*) "CURRENT GRADIENT NORM sqrt(g*g): g_norm=",g_norm

               CALL  PREP_RESTRT &
      (SGITOT,g_inner,ITERNO,GAMATOT,GAMANG)

!--- minimization options 

      if(MINIM_ALG.eq.'STEEP') then
      write(*,*) "STEEPEST DESCENT MINIMIZATION ALGORITHM"
                               call STEEP &
      (N_start,N_end,g,d)
                               call norm &
      (N_start,N_end,d,d_inner,d_norm)

         else if(MINIM_ALG.eq.'CONJG') then
      write(*,*) "CONJUGATE GRADIENT MINIMIZATION ALGORITHM"

                               call  CONJUG &
      (N_start,N_end,iout,g,d)

                               call norm &
      (N_start,N_end,d,d_inner,d_norm)

                               call ANGLE &
      (N_start,N_end,g,d,g_norm,d_norm,ITERNO,iout,GAMATOT,GAMANG,SGITOT)

            else if(MINIM_ALG.eq.'QUASN') then
      write(*,*) "L-BFGS QUASIN-NEWTON MINIMIZATION ALGORITHM"

                               call  QUASIN &
      (N_start,N_end,iout,ITERNO,NPAIRS,ioutmax,g,d)

                               call norm &
      (N_start,N_end,d,d_inner,d_norm)

                               call ANGLE &
      (N_start,N_end,g,d,g_norm,d_norm,ITERNO,iout,GAMATOT,GAMANG,SGITOT)

      else
      write(*,*) "PROBLEM: MINIMIZATION ALGORITHM NOT FOUND !!!! ",MINIM_ALG
      stop
      end if

      deallocate(g)

!--- descent direction norm 

        write(*,*) "CURRENT DESCENT DIRECTION NORM sqrt(d*d): d_norm=",d_norm
        write(3,*) "CURRENT DESCENT DIRECTION NORM sqrt(d*d): d_norm=",d_norm

!---------------------------------
!--- write descent direction

      write(descname,200) iout 

      close(idesc)
      open(UNIT=idesc,FILE=descname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*)' idesc OPEN UNIT ERROR IER=',IER
      write(idesc) (d(k),k=N_start,N_end)
      close(idesc,status='KEEP')

      ij=N_start-1
             do j=1,ncateg
            do i=NENS_START,NENS
      ij=ij+1
      if(j.eq.min(5,ncateg)) then
      write(*,*) " i=",i," j=",j," d(ij)=",d(ij)
      endif
            end do
             end do

      deallocate(d)

!--- write restart parameters

      CLOSE(isha)
      OPEN(UNIT=isha,FILE='shanno',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' ISHA OPEN UNIT ERROR IER=',IER
      write(isha) SGITOT,ITERNO
      close(isha,status='KEEP')

!------------------------------------------
       write(*,*) "end minimize"
       write(3,*) "end minimize"
!------------------------------------------
 stop
 end

!=======================================================================
                               subroutine  STEEP &
      (N_start,N_end,g,d)

!---------------------------------------------------
!--- Preconditioned steepest descent minimization algorithm

     integer :: N_start,N_end
     real,dimension(N_start:N_end) :: g,d
     real,allocatable,dimension(:) :: g_prec
!---------------------------------------------------
!-- preconditioned g
      allocate(g_prec(N_start:N_end))
      call PRECOND (N_start,N_end,g,g_prec)

      d=-g_prec

      deallocate(g_prec)

 return
 end subroutine STEEP

!=======================================================================
                               subroutine  CONJUG &
      (N_start,N_end,iout,g,d)

!---------------------------------------------------
!--- Fletcher-Reeves or Polak-Ribiere conjugate-gradient minimization algorithms

     integer,parameter :: igrad_old=30
     integer,parameter :: idesc_old=31
     integer :: N_start,N_end
     integer :: iout
     real,dimension(N_start:N_end) :: g,d
     real,allocatable,dimension(:) :: g_old,d_old
     real,allocatable,dimension(:) :: g_prec,g_prec_old
     real :: sumu,sumd,BETA
     character*12 :: gold_name
     character*11 :: dold_name
!---------------------------------------------------
!-- preconditioned g
      allocate(g_prec(N_start:N_end))
      call PRECOND (N_start,N_end,g,g_prec)

                IF(iout.eq.1) THEN
      d=-g_prec
                ELSE
!-- read g_old
      allocate(g_old(N_start:N_end))
      iterm1=iout-1
      write(gold_name,100) iterm1
 100  format('gradient_',i2.2)

      close(igrad_old)
      open(UNIT=igrad_old,FILE=gold_name,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' Igrad_old OPEN UNIT ERROR IER=',IER
      read(igrad_old) (g_old(k),k=N_start,N_end)
      close(igrad_old,status='keep')

!-- preconditioned g_old
      allocate(g_prec_old(N_start:N_end))
      call PRECOND (N_start,N_end,g_old,g_prec_old)

!-- read d_old
   allocate(d_old(N_start:N_end))
      iterm1=iout-1
      write(dold_name,200) iterm1
 200  format('descent_',i2.2)

      close(idesc_old)
      open(UNIT=idesc_old,FILE=dold_name,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' idesc_old OPEN UNIT ERROR IER=',IER
      read(idesc_old) (d_old(k),k=N_start,N_end)
      close(idesc_old,status='keep')

      sumu=0.0
      sumd=0.0

           do i=N_start,N_end
            sumu=sumu+g(i)*g_prec(i)   !! Fletcher-Reeves
!!!!!       sumu=sumu+(g(i)-g_old(i))*g_prec(i)   !! Polak-Ribiere
            sumd=sumd+g_old(i)*g_prec_old(i)
           end do

       if(sumd.gt.0.) then
      BETA=sumu/sumd
       else
      write(*,*) "PROBLEM: NEGATIVE BETA in C-G: sumu,sumd=",sumu,sumd
      return
       end if

         do i=N_start,N_end
        d(i)=-g_prec(i)+BETA*d_old(i)
         end do

     deallocate(g_prec_old)
     deallocate(g_old) ; deallocate(d_old)

               END IF

     deallocate(g_prec)

 return
 end subroutine CONJUG

!=======================================================================
                               subroutine  QUASIN &
      (N_start,N_end,iout,ITERNO,NPAIRS,ioutmax,g,d)

!-------------------------------------------------------------------
!
!   THIS IS THE QUASI-NEWTON MINIMIZATION (L-BFGS)  
!                                                  
!   s(k)=alfa(k)*d(k)
!   y(k)=g(k+1)-g(k) (OLD) 
!                                                 
!   (J.NOCEDAL (1980), Math.Comp.,Vol.35,773-782)
!-------------------------------------------------------------------

     integer,parameter :: iHH=81
     integer :: N_start,N_end,N
     integer :: iout,ITERNO,NPAIRS,ioutmax
     integer :: nn_upper
     real,dimension(N_start:N_end) :: g,d
     real,dimension(:),allocatable	:: g_prec
     real,dimension(:),allocatable	:: Y,S
     real,dimension(:),allocatable	:: HH_upper
     real,dimension(:,:),allocatable		:: HH,HH_old
     real,dimension(:,:),allocatable		:: VK,SK
     integer 	IBOUND,IBNDMIN
     real SUM1,SUM2,SUMD
!====================================================================
        WRITE(3,*)"            "
        WRITE(3,*)"BEGIN QUASI-NEWTON "
        WRITE(3,*)"            "
!==================================================
!
!----------IBOUND = MAX NUMBER OF UPDATE FILES------------
!--- (correct if there was a restart during the course of minimization)
         IBOUND=ITERNO-1
         IBNDMIN=MAX(IBOUND-NPAIRS+1,iout-ITERNO+1)
!---------------------------------------------
          WRITE(3,*) "                         "
          WRITE(3,*) "NUMBER OF PAIRS (Y,S) IS ",NPAIRS
          WRITE(3,*) "IBNDMIN= ",IBNDMIN,"  IBOUND= ",IBOUND
          WRITE(3,*) "                         "
          WRITE(*,*) "NUMBER OF PAIRS (Y,S) IS ",NPAIRS
          WRITE(*,*) "IBNDMIN= ",IBNDMIN,"  IBOUND= ",IBOUND
          WRITE(*,*) "                         "
!--------------------------------------------------
      allocate(Y(N_start:N_end))
      allocate(S(N_start:N_end))

!---- start 1st main IB loop 
      allocate(HH_old(N_start:N_end,N_start:N_end))
      allocate(HH(N_start:N_end,N_start:N_end))
      allocate(VK(N_start:N_end,N_start:N_end))
      allocate(SK(N_start:N_end,N_start:N_end))

!-- initialize H_old=I -------------------------
      HH(:,:)=0.0
      do i=N_start,N_end
      HH(i,i)=1.0
      end do


      IF(iout.eq.1) THEN

      write(*,*) "START STEEPEST DESCENT: 1st iteration"
      allocate(g_prec(N_start:N_end))
      call PRECOND (N_start,N_end,g,g_prec)
      d=-g_prec
      deallocate(g_prec)

      ELSE

      write(*,*) "START QUASI-NEWTON: iout=",iout

            DO IB=IBNDMIN,IBOUND

      HH_old(:,:)=HH(:,:)

      CALL  READ_Ro_Y_S (IB,N_start,N_end,Y,S,Ro)

!---  calculate the VK-matrix = I-Ro*y*s**T
           do j=N_start,N_end
            do i=N_start,N_end
             VK(i,j)=-Ro*y(i)*s(j)
            end do
           end do
 
            do i=N_start,N_end
             VK(i,i)=1.0 +  VK(i,i)
            end do

!---  calculate the SK-matrix = Ro*s*s**T
           do j=N_start,N_end
            do i=N_start,N_end
             SK(i,j)= Ro*s(i)*s(j)
            end do
           end do
 
!-- (VK**T)*HH_old*VK+SK --------------------------

           do i=N_start,N_end
            do j=N_start,N_end
             sum1=0.0
             do k=N_start,N_end
             sum1=sum1+HH_old(i,k)*VK(k,j)
             end do
             HH(i,j)=sum1
             end do
            end do

           do i=N_start,N_end
            do j=N_start,N_end
             sum2=0.0
             do k=N_start,N_end
             sum2=sum2+VK(k,i)*HH(k,j)
             end do
             HH_old(i,j)=sum2
            end do
           end do

             HH(:,:)=HH_old(:,:)+SK(:,:)


                      ENDDO

       deallocate(VK) ; deallocate(SK)
       deallocate(HH_old)

       END IF

!-------- descent direction ----------------------
           do i=N_start,N_end
           sumd=0.0
            do j=N_start,N_end
             sumd=sumd+HH(i,j)*g(j)
            end do
           d(i)=-sumd
           end do

!--- HH matrix -------------
            do i=N_start,N_end
       write(*,110) i,(HH(i,j),j=1,10)
            end do
 110 format(I5,10E10.3)

!--- write upper-packed HH matrix -------------
      N=N_end-N_start+1
      nn_upper=N*(N+1)/2
      write(*,*) "iout=",iout," N=",N," nn_upper=",nn_upper
      allocate(HH_upper(1:nn_upper))
      HH_upper(:)=0.0

         Kindex=0
          do jj=N_start,N_end
          do ii=N_start,jj
         Kindex=Kindex+1
         HH_upper(Kindex)=HH(ii,jj)
          end do
          end do

      close(iHH)
      open(UNIT=iHH,FILE='HH_matrix',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*) iHH,' OPEN UNIT ERROR IER=',IER
      write(iHH) (HH_upper(i),i=1,nn_upper)
      close(iHH,status='keep')

       deallocate(HH_upper)
       deallocate(HH)
!=====================================================

        WRITE(3,*)"            "
        WRITE(3,*)"END QUASI-NEWTON"
        WRITE(3,*)"            "
!=========================================================

  return
  end subroutine QUASIN

!=======================================================================
                        subroutine PRECOND &
      (N_start,N_end,x,x_prec)

!----------------------------------------
! Currently, no Hessian preconditioning
!----------------------------------------

     integer :: N_start,N_end
     real,dimension(N_start:N_end) :: x,x_prec

      x_prec=x


 return
 end subroutine PRECOND

!==================================================================
                        subroutine ANGLE &
      (N_start,N_end,g,d,g_norm,d_norm,ITERNO,iout,GAMATOT,GAMANG,SGITOT)

!--------------------------------------------------------
!   THIS PROGRAM COMPUTES THE ANGLE BETWEEN g(k+1) AND d(k+1)     
!   AND DOES THE ANGLE RESTART TEST 
!--------------------------------------------------------
       real,parameter :: CSTR=180./3.1415
       real,parameter :: EPS_angle=1.E-05
       integer :: N_start,N_end
       integer :: ITERNO,iout
       real,dimension(N_start:N_end) :: g,d
       real :: g_norm,d_norm
       real :: g_inner
       real :: GAMATOT,GAMANG
       real :: TETATOT,TETANG,ANGND,GDTOT
!====================================================================
!           CALCULATE INNER PRODUCT (-(g,d))

       GDTOT =0.0 
         do i=N_start,N_end
       GDTOT=GDTOT - g(i)*d(i)
         end do

      IF(iout.gt.1) then
      write(*,*) "subroutine ANGLE: current minimization iteration=",iout
!-------------------------------------------------
           WRITE(3,2001) g_norm,d_norm
           WRITE(*,2001) g_norm,d_norm
!-------------------------------------------------
       IF((g_norm.GT.0.0).AND.(d_norm.GT.0.0)) THEN 
           WRITE(3,*) "GDTOT=",GDTOT
           WRITE(*,*) "GDTOT=",GDTOT
           WRITE(*,*) "g_norm,d_norm=",g_norm,d_norm
           gd_norm=g_norm*d_norm
           WRITE(*,*) "GDTOT,gd_norm=",GDTOT,gd_norm
             if(GDTOT.lt.EPS_angle) then
           write(*,*) "iout=",iout," old GDTOT=",GDTOT," old gd_norm=",gd_norm
           GDTOT=EPS_angle
           gd_norm=EPS_angle
           ANGND=1.0
             else
           write(*,*) "iout=",iout," GDTOT=",GDTOT
           ANGND=min(1.,GDTOT/gd_norm)
           write(*,*) "ANGND=",ANGND
             end if
       TETATOT=ANGND**2
       ELSE
           WRITE(*,*) "subroutine ANGLE: PROBLEM with negative norms !!!!!!"
           WRITE(*,*) "g_norm,d_norm=",g_norm,d_norm
       END IF
           WRITE(3,2103) TETATOT,GAMATOT
           WRITE(*,2103) TETATOT,GAMATOT
       ANGND =ACOS(ANGND )*CSTR
           WRITE(*,2003) ANGND
           WRITE(3,2003) ANGND
!-------------------------------------------------
 2001  FORMAT('TOTAL NORMS: GRAD,DESC ',2E15.5)
 2003  FORMAT('TOTAL ANGLES: GRAD-DESC ',F15.5)
 2103  FORMAT('TETA,GAMA TOTAL ',2E15.5)

       TETANG =ACOS(SQRT(TETATOT))*CSTR

          WRITE(3,7777) TETANG
          WRITE(*,7777) TETANG
 7777     FORMAT('TOTAL ANGLES: TETA (Q-N)  ',F12.5)
!----------------------------
!  MONITOR THE ANGLES (RESTART)
!  GAMANG = ANGLE BETWEEN G and D USING FLETCHER-REEVES CG
!  TETANG = ANGLE BETWEEN G and D USING QUASI-NEWTON MINIMIZATION
!  (SHOULD BE 0<TETANG<GAMANG FOR GOOD CONVERGENCE)

      WRITE(3,*)"COSTF ITERATION INDEX IS",iout

      im1=iout-1
      ICOST=301
      CLOSE(ICOST)
      OPEN(UNIT=ICOST,FILE='outangle',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' ICOST OPEN UNIT ERROR IER=',IER
      WRITE(ICOST,400) im1,GAMANG,TETANG
      CLOSE(ICOST,STATUS='KEEP')

 400  FORMAT(I2,2(1X,F10.2))
!----------------------------

!=================================================
!========THIS IS A RESTART TEST==================
!=================================================
       IF(TETATOT.LT.GAMATOT) THEN
!=================================================
       WRITE(3,*) "---------------------------"
       WRITE(3,1220) TETATOT,GAMATOT
       WRITE(3,*) "BEFORE RESTART: ITERNO=",ITERNO
 1220  FORMAT('RESTART PROCEDURE BEGINS: TETA,GAMA',2E15.5)
       WRITE(*,*) "---------------------------"
       WRITE(*,1220) TETATOT,GAMATOT
       WRITE(*,*) "BEFORE RESTART: ITERNO=",ITERNO

       ITERNO=1
       SGITOT=0.0

       WRITE(3,*) "AFTER RESTART: ITERNO=",ITERNO
       WRITE(3,*) "---------------------------"
       WRITE(*,*) "AFTER RESTART: ITERNO=",ITERNO
       WRITE(*,*) "---------------------------"
!=================================================
!==== same as the first iteration ============================
                               call norm &
      (N_start,N_end,g,g_inner,g_norm)

        write(*,*) "CURRENT GRADIENT NORM sqrt(g*g): g_norm=",g_norm
        write(3,*) "CURRENT GRADIENT NORM sqrt(g*g): g_norm=",g_norm

               CALL  PREP_RESTRT &
      (SGITOT,g_inner,ITERNO,GAMATOT,GAMANG)

                               call STEEP &
      (N_start,N_end,g,d)

!===============================================
                END IF
      END IF
!===============================================
!========END OF RESTART TEST==================
!===============================================

  return
  end subroutine ANGLE

!====================================================================
                             subroutine PREP_RESTRT &
     (SGITOT,YGRAD,ITERNO,GAMATOT,GAMANG)

!-----------------------------------------------------------------------
!   THIS PROGRAM PREPARES THE MINIMIZATION RESTART PARAMETERS
!-----------------------------------------------------------------------
      real,parameter :: PIQ=3.1415927
      real,parameter :: CSTR=180./PIQ
      complex Y,Y1,Y2,Z
!====================================================================
!     prepare for minimization restart 

               XGRAD=YGRAD*YGRAD

!===================
               IF(XGRAD.LE.0.) THEN
!===================
               SGITOT=0.
               GAMATOT=0.
!===================
               ELSE
!===================
               SGITOT=SGITOT+1./XGRAD

!----------calculate 'gama' needed for restart--------------------C

              GATOT=XGRAD*SGITOT

      write(3,*) "XGRAD=",XGRAD,"  SGITOT=",SGITOT
      write(*,*) "XGRAD=",XGRAD,"  SGITOT=",SGITOT
!=============================================================
!      FIRST, find the angle of the Fletcher-Reeves C-G method

      WRITE(3,*) " GATOT (GREATHER THAN ONE) = ",GATOT
      WRITE(*,*) " GATOT (GREATHER THAN ONE) = ",GATOT
         XGAT=1./GATOT
         XUSE=MIN(1.,SQRT(XGAT))
!!       XUSE=AMIN1(1.,SQRT(XGAT))
         GAMANG =ACOS(XUSE)*CSTR
       WRITE(3,7777) GAMANG
       WRITE(*,7777) GAMANG
 7777  FORMAT(' ','TOTAL ANGLE FOR FLETCHER-REEVES C-G:  ',F12.5)
!=============================================================
!-----Penalize for large angle ( > CRITANG) in early iterations-----
!
!     Use the third-order Butterworth digital filter response
!     with the bilinear transform 
!
!     (Roberts, R.A. and C.T. Mullis, 1987: Digital Signal Processing)
!

       WRITE(3,*) "CSTR GATOT XGRAD SGITOT ",CSTR,GATOT,XGRAD,SGITOT
       WRITE(*,*) "CSTR GATOT XGRAD SGITOT ",CSTR,GATOT,XGRAD,SGITOT

         CRITANG=75.
         A1=COS(CRITANG/CSTR)
         A2=1./GATOT 
         IF(GAMANG.LT.89.9) THEN
         A12=A1*A1/A2
         ELSE
         A12=1.E+10
         ENDIF
!===============
!     NN   = cut-off iteration number
!
!     NMAX = total number of iterations with "strict" test

            NN=13
            NMAX=2*NN+1
            XN=FLOAT(NN)
            XNMAX=FLOAT(NMAX-1)

          IF(ITERNO.LE.NMAX) THEN
             XIT=FLOAT(ITERNO-1)
          ELSE
             XIT=XNMAX
          ENDIF
             XTETA=PIQ*XIT/XNMAX
             Z=CMPLX(0.,XTETA)

             Y1=(1.+CEXP(-Z))**3
             Y2=2.*(3.+CEXP(-2.*Z))
             Y=Y1/Y2
             XX=CABS(Y)

       WRITE(3,*) "ITERNO XTETA  XX ",ITERNO,XTETA,XX
       WRITE(*,*) "ITERNO XTETA  XX ",ITERNO,XTETA,XX

!=================================================
!      DD = approximately the inverse condition number
!           of the (preconditioned) Hessian at the minimum

         DD=5.0E-01
         TAU=A12*XX + DD

              WRITE(3,*) "TAU= ",TAU
              WRITE(*,*) "TAU= ",TAU

!=============================================================
             IF(GATOT.LE.0.) THEN
             GAMATOT=0.
             ELSE
             GAMATOT=TAU/GATOT
             ENDIF
!===================
               ENDIF
!===================

  return
  end subroutine PREP_RESTRT

!===================================================================
                subroutine  READ_Ro_Y_S & 
      (K,N_start,N_end,Y,S,Ro)
!-----------------------------------------------------------------------
      integer,parameter :: igrad=40
      integer,parameter :: idesc=41
      integer,parameter :: ialfa=42

      real 	Ro,ALFA,SUMRo
      integer 	K
      integer 	N_start,N_end
      real,dimension(N_start:N_end) :: Y,S
      character*12 gradname
      character*11 descname
      character*8 alfaname
      real,allocatable,dimension(:) :: g,g_old,d
!-----------------------------------------------

      WRITE(3,*) "IN READ_Ro_Y_S: K= ",K

 100  format('gradient_',i2.2)
 200  format('descent_',i2.2)
 300  format('alfa_',i2.2)

!--- new gradient
      K1=K+1
      WRITE(3,*) "READ_Ro_Y_S: read gradient: K+1= ",K1
      WRITE(*,*) "READ_Ro_Y_S: read gradient: K+1= ",K1

      allocate(g(N_start:N_end))
      write(gradname,100) K1

      close(igrad)
      open(UNIT=igrad,FILE=gradname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' igrad OPEN UNIT ERROR IER=',IER
      read(igrad) (g(i),i=N_start,N_end)
      CLOSE(igrad,status='KEEP')

!--- old gradient
       WRITE(3,*) "                    "
       WRITE(3,*) "READ_Ro_Y_S: read gradient:  K= ",K
       WRITE(*,*) "READ_Ro_Y_S: read gradient:  K= ",K

      allocate(g_old(N_start:N_end))
      write(gradname,100) K

      close(igrad)
      open(UNIT=igrad,FILE=gradname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' igrad OPEN UNIT ERROR IER=',IER
      read(igrad) (g_old(i),i=N_start,N_end)
      CLOSE(igrad,status='KEEP')

!---------------y(k)=g(k+1)-g(k)----------------------------

             Y(:)=g(:)-g_old(:)

      deallocate(g) ; deallocate(g_old)

!--- descent direction
      WRITE(3,*) "READ_Ro_Y_S: read descent: K= ",K
      WRITE(*,*) "READ_Ro_Y_S: read descent: K= ",K

      allocate(d(N_start:N_end))
      write(descname,200) K

      close(idesc)
      open(UNIT=idesc,FILE=descname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' idesc OPEN UNIT ERROR IER=',IER
      read(idesc) (d(i),i=N_start,N_end)
      CLOSE(idesc,status='KEEP')

!--------------------------------------------------
      write(alfaname,300) K

      close(ialfa)
      open(UNIT=ialfa,FILE=alfaname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' ialfa OPEN UNIT ERROR IER=',IER
      read(ialfa) ALFA
      CLOSE(ialfa,status='KEEP')

      WRITE(3,*)"READ_Ro_Y_S:  ALFA= ",ALFA

!-----------get alfa*d(k)-----------------------

             S(:)=ALFA*d(:)

      deallocate(d)
!---------------------------------------------
!   calculate Ro=1./(y,s)

       SUMRo=0.0
         do i=N_start,N_end
         SUMRo=SUMRo+Y(i)*S(i)
         end do

      if(SUMRo.gt.0.) then
       Ro=1./SUMRo
      else
       Ro=1.
       write(*,*) "PROBLEM: negative SUMRo ",SUMRo
      end if

!-----------------------------------------------
  return
  end subroutine READ_Ro_Y_S

!========================================================================
                           subroutine norm &
      (N_start,N_end,x,x_inner,x_norm)

       integer :: N_start,N_end
       real :: x(N_start:N_end)
       real :: x_inner,x_norm

        x_inner=0.0
         do i=N_start,N_end
           x_inner=x_inner+x(i)*x(i)
          end do

        x_norm=sqrt(x_inner)

      return
      end subroutine norm