PROGRAM space
        USE nrtype
        USE error
        USE space_pars, ONLY: read_space_pars,sn_ratio,sn_flag,&
         &TGM_force,error_est,flag_norm,flag_ABD_loop,&
         &flag_move_user
        USE share
        USE utils
        USE make_model
        USE read_GCOG
        USE func_poly
!        USE quick_est
        USE interfaces, ONLY: fit_cont,write_res
        USE minimize
        USE uncertains2
        IMPLICIT NONE
        INTEGER(I2B) :: i, j
        INTEGER(I1B) :: AllocateStatus
        REAL(DP) :: sn_new
        REAL(DP),DIMENSION(5) :: temp_gridS,logg_gridS,met_gridS
        REAL(DP), DIMENSION(3) :: TGM_old, residTGM
        REAL(DP), DIMENSION(4) :: pars
        !variables for the routine mrqmin
        REAL(DP), DIMENSION(5,5) :: covarTGM,alpha_covarTGM
        REAL(DP), DIMENSION(:), ALLOCATABLE :: sig_noise
        REAL(DP), DIMENSION(:,:), ALLOCATABLE :: covarABD,alpha_covarABD
        REAL(DP) :: chisq,chisqTGM,chisqABD
        LOGICAL :: flag_lim,flag_limS,flag_move
        CHARACTER(LEN=5),DIMENSION(3) :: TGMc

        !initialize variables
        TGM_prox=REAL((/-99.,-99.,-99./),SP)
        flag_lim=.FALSE.
        flag_limS=.FALSE.
        flag_move=.FALSE.
        convTGM=.TRUE.
        convABD=.TRUE.
        !read the parameters file
        CALL read_space_pars('space.par')
        !load the spectrum, the line lists
        CALL read_files(obs_sp_file,llist,llist_rej)!
        !----------------------
        !allocate the mask for the llist that will be use for the measures
        ALLOCATE(select_ll_mask(dim_ll), STAT = AllocateStatus)
        IF (AllocateStatus/=0) CALL stop_msg('Not enough memory to allocate select_ll_mask')

        sn=sn_ratio
        !allocate the array sig_noise for the mrqmin routine
        ALLOCATE(sig_noise(dimsp), STAT = AllocateStatus)
        IF (AllocateStatus/=0) CALL stop_msg('Not enough memory to allocate sig_noise')

        ! allocate the array ele2meas which contains the list of elements
        ! present in the line list. It also allocate ABD, ADB_old,residABD
        ! ABD_mask. It outputs dim_ele
        CALL alloc_ABD(ele_ll,dim_ll,dim_ele)

        !allocate the arrays covarABD and alpha_covarABD 
        ALLOCATE(covarABD(dim_ele,dim_ele), STAT = AllocateStatus)
        IF (AllocateStatus/=0) CALL stop_msg('Not enough memory to allocate covar')
        ALLOCATE(alpha_covarABD(dim_ele,dim_ele), STAT = AllocateStatus)
        IF (AllocateStatus/=0) CALL stop_msg('Not enough memory to allocate alpha_covar')

        !initialize TGM_mask
        TGM_mask=.TRUE.
        !initialize covar
        covarTGM=0.
        covarABD=0.
        !initialize alpha_covar
        alpha_covarTGM=0.
        alpha_covarABD=0.

        !give a starting point
        TGM(1:3)=(/5000.,3.0,-0.4/)
        !initialize TGM to the user's TGM, if any
        IF (TGM_force(1)>-9.) THEN
          TGM(1)=TGM_force(1)
          TGM_mask(1)=.FALSE.
        END IF  
        IF (TGM_force(2)>-9.) THEN
          TGM(2)=TGM_force(2)
          TGM_mask(2)=.FALSE.
        END IF  
        !initialize other TGM variables
        TGM(4)=sigma
        TGM(5)=0.

        !initialize the internal normalized spectrum
        f_sp_norm=f_sp
        !first TGM estimation
        CALL load_GCOG_4deg_quick
        weights=1.0_dp
        sig_noise=(1.0_dp/sn)/weights
        select_ll_mask=.TRUE.

        CALL minimTGM_Q(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
        &,covarTGM,alpha_covarTGM,chisqTGM)
        write(TGMc(1),fmt='(I5)') INT(TGM(1))
        write(TGMc(2),fmt='(F5.2)') TGM(2)
        write(TGMc(3),fmt='(F5.2)') TGM(3)
        write(*,*) 'TGM first est.        ',' Teff=',TGMc(1),' log g=',TGMc(2),' [M/H]=',TGMc(3)

        !find the closest point in the grid
        CALL find_prox(TGM,TGM_prox)
        !make the sub-grids, upload the GCOG,
        !select the linelist, set the sig_noise
        !all this is done by the routine update_gridS_ll
        CALL update_gridS_ll

        !now re-normalize the spectrum
        CALL make_model_TGM(f_model,TGM)
        IF(flag_norm) THEN
         CALL fit_cont(f_sp,f_model,cont);
         f_sp_norm=f_sp/cont
        END IF

        !compute the new S/N
        IF(sn_flag) THEN
         CALL new_sn(sn)
        END IF

!#########################
!###### start iterations
!#########################
        outer: DO i=1,30
!        write(*,*) 'outer loop TGM_prox',TGM_prox
        !##############
        !search for TGM
        !##############
!        write(*,*) 'start TGM minim, outer loop=',i
        CALL minimTGM(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
        &,covarTGM,alpha_covarTGM,chisqTGM)
        !x=w_sp;y=f_sp_norm;sig=sig_noise;a=TGM;maska=TGM_mask;covar=covar;alpha=alpha_covar;
        !chisq=chisq;
        write(TGMc(1),fmt='(I5)') INT(TGM(1))
        write(TGMc(2),fmt='(F5.2)') TGM(2)
        write(TGMc(3),fmt='(F5.2)') TGM(3)
        write(*,*) 'TGM outer loop',i,' Teff=',TGMc(1),' log g=',TGMc(2),' [M/H]=',TGMc(3)
!        write(*,*) 'TGM outer loop',i,'Teff=',TGM(1),'log g=',TGM(2),'[M/H]=',TGM(3),'fwhm=',2.35*TGM(4),'Rv=',TGM(5)

        !check the new sn
        IF(sn_flag) THEN
         sn_new=sn
         CALL new_sn(sn_new)
         IF(sn<0) CALL error_msg('impossible to derive S/N (too few pixels?)')
         IF(abs(sn_new-sn)>(sn*0.05_dp)) THEN
          sn=sn_new
          CALL update_gridS_ll
          CYCLE
         END IF
        END IF
        !update the value sigma
        sigma=ABS(TGM(4))
        !find the closest point in the sub-grid and verify if TGM is inside the grid
        !write(*,*) temp_gridS,logg_gridS,met_gridS
        CALL find_proxS(TGM,temp_gridS,logg_gridS,met_gridS,TGM_prox,&
        &flag_lim,flag_limS,flag_move)
        !if flag_lim=.TRUE. means that TGM is out of the grid, then exit
!        write(*,*) 'after outer loop TGM_prox',TGM(1:3),TGM_prox,flag_lim,flag_limS,flag_move
        IF(flag_lim) THEN
           CALL error_msg('stellar parameters out of the grid')
        END IF
        !if flag_limS=.TRUE. means that TGM is out of the 
        !gridS, then update the GCOG, the linelist and 
        !cycle the outer loop
        IF(flag_limS) THEN
           CALL update_gridS_ll
           TGM_old=TGM(1:3)
           CYCLE
        END IF
        IF(flag_move) THEN
           CALL update_gridS_ll
           TGM_old=TGM(1:3)
           CYCLE
        END IF



         !###############################
         !### start the inner loop
         !### which normalize the spectrum iteratively and 
         !### estimate TGM and ABD
         inner: DO j=1,30
         !###############################
          !if SP_Ace arrives here means that
          !the outer loop was successful! then, after re-normalization
          !TGM will be re-estimated. 

          !Now store TGM in TGM_old
          TGM_old=TGM(1:3)
          !now re-normalize the spectrum
          CALL make_model_TGM(f_model,TGM)
          
          IF(flag_norm) THEN
           CALL fit_cont(f_sp,f_model,cont);
           f_sp_norm=f_sp/cont
          END IF

          !##############
          !search for TGM
          !##############

          !write(*,*) 'start TGM minim, inner loop=',j,TGM(1:3)
          CALL minimTGM(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
          &,covarTGM,alpha_covarTGM,chisqTGM)
          !x=w_sp;y=f_sp_norm;sig=sig_noise;a=TGM;maska=TGM_mask;covar=covar;alpha=alpha_covar;
          !chisq=chisq;
          write(TGMc(1),fmt='(I5)') INT(TGM(1))
          write(TGMc(2),fmt='(F5.2)') TGM(2)
          write(TGMc(3),fmt='(F5.2)') TGM(3)
          write(*,*) 'TGM inner loop',j,' Teff=',TGMc(1),' log g=',TGMc(2),' [M/H]=',TGMc(3)
!         write(*,*) 'TGM inner loop',j,'Teff=',TGM(1),'log g=',TGM(2),'[M/H]=',TGM(3),'fwhm=',2.35*TGM(4),'Rv=',TGM(5),sn
          !check the new sn
          IF(sn_flag) THEN
           sn_new=sn
           CALL new_sn(sn_new)
           IF(sn<0) CALL error_msg('impossible to derive S/N (too few pixels?)')
           IF(abs(sn_new-sn)>(sn*0.05_dp)) THEN
            sn=sn_new
            CALL update_gridS_ll
            CYCLE
           END IF
          END IF

          !find the closest point in the sub-grid and verify if TGM is inside the grid
          CALL find_proxS(TGM,temp_gridS,logg_gridS,met_gridS,TGM_prox,&
          &flag_lim,flag_limS,flag_move)
!         write(*,*) 'after inner loop TGM_prox',TGM(1:3),TGM_prox,flag_lim,flag_limS,flag_move
          !if flag_lim=.TRUE. means that TGM is out of the grid, then exit
          !the program
          IF(flag_lim) THEN
           CALL error_msg('stellar parameters out of the grid')
          END IF
          !if flag_limS=.TRUE. means that TGM is out of the 
          !gridS, then update the GCOG, the linelist, 
          !exit the inner loop and continue with the outer loop
          IF(flag_limS) THEN
           CALL update_gridS_ll
           TGM_old=TGM(1:3)
           CYCLE
          END IF
         !if flag_move=.TRUE. means that the closest gridS point to TGM is changed.
         !then update the GCOG, the linelist and cycle the inner loop
         IF(flag_move.AND.flag_move_user) THEN
           CALL update_gridS_ll
           TGM_old=TGM(1:3)
           CYCLE
         END IF

         !update some variables
         sigma=ABS(TGM(4))
         residTGM=ABS(TGM_old-TGM(1:3))
         TGM_old=TGM(1:3)

         IF(residTGM(1)>20..OR.ANY(residTGM(2:3)>0.05).AND.j<30) THEN
          CYCLE
         ELSE

         !check if some of the parameters are meaningful.
         !if not, quit.
         IF((TGM(4)*2.35)>10.) THEN
          CALL error_msg('FWHM does not converge')
         END IF
         !IF(ABS(TGM(5))*w_sp(1)/300000.)>TGM(4)) THEN
         IF(ABS(TGM(5))>100.) THEN
          CALL error_msg('RV too large (|Rv|>100km/sec)')
         END IF
         IF(ALL(PACK((/5000.,3.0/),TGM_mask(1:2))==PACK(TGM(1:2),TGM_mask(1:2))).AND.ALL(TGM_mask)) THEN
          CALL error_msg('No convergence point found in the parameter space')
         END IF


          !##################
          !now search for ABD
          !##################
!          write(*,*) 'start ABD minim, inner loop=',j

          CALL minimABD(w_sp,f_sp_norm,sig_noise,ABD,ABD_mask&
          &,covarABD,alpha_covarABD,chisqABD)
          !x=w_sp;y=f_sp_norm;sig=sig_noise;a=TGM;maska=TGM_mask;covar=covar;alpha=alpha_covar;
          !chisq=chisq;
        !write(*,*) 'ele2write',(PACK(ele2meas,write_ABD_mask))
        !write(*,*) 'ABD',(PACK(ABD,write_ABD_mask))
        !write(*,*) ele2meas

          IF(flag_ABD_loop) THEN
           !write(*,*) 'ABD_mask',(PACK(ABD_mask,write_ABD_mask))
           !write(*,*) 'ABD_mask',(PACK(ABD_mask,write_ABD_mask))
           !write(*,*) 'ABD_write',ele2write
           !write(*,*) 'ABD',(PACK(ABD,write_ABD_mask))

           !compute the differences between the old and new ABD estimations
           residABD=ABS(ABD-ABD_old)
      !write(*,*) 'residABD',(PACK(residABD,write_ABD_mask))
           !if the differences are large, then re-do the inner loop 
           IF(ANY((PACK(residABD,write_ABD_mask))>0.05).AND.j<30) THEN
           !Now substitute the Fe abundance (if the internal loop cycled at least one time) in TMG(3)  
           !if this is the first loop for the internal loop, then ABD=0 and TGM does not change
           TGM(3)=TGM(3)+ABD(1)
           ABD=ABD-ABD(1)
           !do not allow ABD to go beyond the grid
           WHERE(ABD>0.8.AND.ABD_mask) ABD=0.8_dp
           WHERE(ABD<-0.6.AND.ABD_mask) ABD=-0.6_dp
           WHERE((ABD>0.79.OR.ABD<-0.59).AND.ABD_mask) ABD_mask=.FALSE.
           WHERE(.NOT.ABD_mask) ABD=0.0

           !compute the differences between the old and new ABD estimations
           residABD=ABS(ABD-ABD_old)
           !update the ABD_old
           ABD_old=ABD
            CYCLE
           END IF
          END IF

          EXIT
          END IF

         !##################
         !### end inner loop
         END DO inner
         !#################

        IF(convTGM.AND.convABD) THEN
         chisq=chisqABD
         EXIT
        ELSE IF (convTGM.AND..NOT.convABD) THEN
         chisq=chisqTGM
         EXIT
        ELSE IF (.NOT.convTGM) THEN
         CALL error_msg('minimization routines do not converge. ')
         EXIT
        END IF

        END DO outer
!######################
        CALL make_model_TGM(f_model,TGM)
        chisq=SUM(((f_sp_norm-f_model)*sn_var*weights)**2)


        IF(error_est) THEN
         write(*,*) 'errors estimation in progress'
         CALL TGM_errors(chisq)
         ELSE
         CALL TGM_errors_null
        END IF

        !write the results and exit
        CALL make_model_TGM(f_model,TGM)
        CALL write_res(error_est,chisq)

        CONTAINS

!#######################################
        SUBROUTINE update_gridS_ll
        IMPLICIT NONE
        LOGICAL,DIMENSION(dimsp) :: mask1

        !make a new subgrids
        CALL make_gridS(TGM_prox,temp_gridS,logg_gridS,met_gridS)

        ! allocate and upload the table coeff_4deg which contains the coefficient
        ! to construct the polynomial GCOG
        CALL load_GCOG_4deg(TGM_prox(1),TGM_prox(2),TGM_prox(3))

        !remove the cosmic rays with a 4sigma clipping
        CALL cosmic_rej

        IF(sn_flag) THEN
         CALL find_sn_var(f_sp,f_model,TGM(4),sn_var)
        ELSE
         sn_var=sn
        END IF
        !select the lines strong enough to be measurable in this
        !region of the parameters space
        pars=(/TGM_prox(1),TGM_prox(2),TGM_prox(3),0._dp/)
        CALL select_lines(ew_poly4,coeff_4deg,pars,sigma)

        sig_noise=(1.0_dp/sn_var)/weights
!        write(*,*) 'new TGM_prox', TGM_prox

!   OPEN(unit=9,file='butta1',action='WRITE')
!   write(9,fmt="(F9.5,TR1,F9.5,TR1,F9.5)") (sig_noise(j),weights(j),sn, j=1,dimsp)
!   CLOSE(9)

        !give to n_weig the number of pixels having
        !weight>0.01
        mask1=.TRUE.
        WHERE(weights<0.01) mask1=.FALSE.
        n_weig=INT(COUNT(mask1),I2B)
        IF(n_weig<1) THEN
         n_weig=dimsp
        END IF

        END SUBROUTINE update_gridS_ll
!########################################
END PROGRAM space