SUBROUTINE fit_cont(f_sp,f_model,cont)
        USE nrtype
        USE stats
        USE share, ONLY: dimsp,weights,rad_pix
        IMPLICIT NONE

        REAL(DP) :: avg_f,avg_r,avg_m,var,sig,sig_r,avg_mm,sig_mm,set,weig
        REAL(DP),DIMENSION(dimsp), INTENT(IN) :: f_sp,f_model
        REAL(DP),DIMENSION(dimsp) :: resid
        REAL(DP),DIMENSION(dimsp), INTENT(OUT) :: cont
        INTEGER(I2B) :: i,iinf,isup,dimsp_loop
        LOGICAL,DIMENSION(dimsp) :: mask

        resid=f_sp-f_model
 

        IF (rad_pix(1)<dimsp) THEN
         dimsp_loop=dimsp
        ELSE
         dimsp_loop=1
        END IF

       DO i=1,dimsp_loop

        IF(dimsp_loop>1) THEN
        weig=weights(i)
        ELSE
        weig=1.
        END IF
        
        IF(weig>1e-3) THEN
           mask=.FALSE.

           isup=MIN(i+rad_pix(i),dimsp);
           iinf=MAX(i-rad_pix(i),1_I2B);
           WHERE(weights(iinf:isup)>1e-3)
            mask(iinf:isup)=.TRUE.
           END WHERE

           IF(ANY(mask(iinf:isup))) THEN
            CALL avg_std_(f_sp(iinf:isup),mask(iinf:isup),avg_f,var,sig)
            CALL avg_std_(resid(iinf:isup),mask(iinf:isup),avg_r,var,sig_r)
            !pseudo-continuum level of the model
            CALL avg_std_(f_model(iinf:isup),mask(iinf:isup),avg_mm,var,sig_mm)
            !clip 2 sigma low
            WHERE(f_sp(iinf:isup)<=(avg_f-2.0_dp*sig_r)) mask(iinf:isup)=.FALSE. 
            !clip 4 sigma high
            WHERE(f_sp(iinf:isup)>=(1.0+4.0_dp*sig_r)) mask(iinf:isup)=.FALSE.
           END IF !IF(ANY(mask(iinf:isup))) THEN

           IF(COUNT(mask)>1) THEN
            CALL avg_std_(f_sp(iinf:isup),mask(iinf:isup),avg_f,var,sig)
            CALL avg_std_(f_model(iinf:isup),mask(iinf:isup),avg_m,var,sig)
            CALL avg_std_(resid(iinf:isup),mask(iinf:isup),avg_r,var,sig_r)
           ELSE
            avg_f=1.0_dp
            avg_m=1.0_dp
            sig_r=0.0_dp
           END IF !IF(COUNT(mask)>1) THEN

          set=4.*(1.-exp((1.-avg_mm)**3))
!          write(*,*) avg_mm,set
           cont(i)=1.0_dp+(avg_f-avg_m-set);

        ELSE !IF(weig>1e-3) THEN

           IF(i==1) THEN
            cont(i)=1.0_dp
           ELSE
            cont(i)=cont(i-1)
           END IF
          
        END IF !IF(weig>1e-3) THEN

       END DO

        !if the radius is bigger than the extention of
        ! the spectrum, then cont(1) is valid for
        !the whole spectrum.
        IF (rad_pix(1)>=dimsp) THEN
         cont=cont(1)
        END IF



        END SUBROUTINE fit_cont