subroutine comp_f0(smv,smn,buf1,fre,nf,n_win,jmin,jmax,
     #                    f1,f0,f2,ierr)
c
c  Find an estimate of Fo according to Chatelain & Guillier method
c  1.25 < rf < 1.50;     f0/rf < window_width < f0*rf
c 
c a.tento 06/03  v.1.0.
c a.tento 02/04  v.1.1. --> geometric mean of Fo
c                       --> Fo = -9 if it cannot be estimated
c
        implicit none
c
        real  smv, smn, buf1, fre, f1, f0, f2
        integer ierr, jmin, jmax, nf, n_win
        dimension  buf1(nf), smv(nf,n_win), smn(nf,n_win), fre(nf)
c
        integer i, j, jpeak, j1, j2, kwin, kf0 
        real  rf, fmin, fmax, ff1, ff2, sd, bb
        parameter ( fmin = 0.2, fmax = 20.0 )
c
c
        ierr = 0
c
        jpeak = kf0(buf1,nf,1,1,jmin,jmax)
c
        f0 = fre(jpeak)
        if ( jpeak  .eq.  jmin   .or.   jpeak  .eq.  jmax ) then
                write(*,*) ' WARNING : Fo does not appear in h/v'
                f0 = -9.0
                f1 = f0
                f2 = f0
                ierr = 2
                return
        endif
c
c
        rf = 1.5 - 0.25*(f0 - fmin) / (fmax - fmin)
        if ( rf  .lt.  1.25 ) rf = 1.25
        if ( rf  .gt.  1.5  ) rf = 1.5
        ff1 = f0 / rf
        ff2 = f0 * rf
c
c
        if ( ff1  .lt.  fre(jmin) ) then
        write(*,100) fre(jmin), ff1
        ierr = 2
100     format(' WARNING : Available lower frequency for searching of',
     #  ' Fo is:',f8.3,/,' instead of: ', f8.3)
        j1 = jmin
        else
        do 24 i = jmin, jmax
        j1 = i 
        if ( fre(i)  .ge.  ff1 ) goto 25
24      continue  
25      continue
        endif
c         
c
        if ( ff2  .gt.  fre(jmax) ) then
        write(*,110) fre(jmax), ff2
        ierr = 2
110     format(' WARNING : Available upper frequency for searching of',
     #  ' Fo is:',f8.3,/,' instead of: ', f8.3)
        j2 = jmax
        else
        do 26 i = jmin, jmax
        j2 = i 
        if ( fre(i)  .ge.  ff2 ) goto 27
26      continue  
27      continue
        endif
c
c
        kwin = 0
        do 3 j = 1, n_win 
        jpeak = kf0(smv,nf,n_win,j,j1,j2)
        if ( jpeak  .eq.  j1   .or.   jpeak  .eq.  j2 ) goto 3
        kwin = kwin + 1
        smn(1,kwin) = fre(jpeak)
c        smn(2,kwin) = alog(smv(jpeak,j))
3       continue
c
c       
        if ( kwin  .eq. 0 ) then
            write(*,*) ' WARNING : it is not possible to compute Fo'
            f0 = -9.0
            f1 = f0
            f2 = f0
            ierr = 2
            return
        endif
        f0 = 0.0
        do 4 i = 1, kwin
        f0 = f0 + alog(smn(1,i))
4       continue
        f0 = f0 / float(kwin)
c
        sd = 0.0
        if ( kwin   .ge.   2 ) then
        do 6 i = 1, kwin
        bb = alog(smn(1,i)) - f0
        sd = sd + bb*bb
6       continue
        sd = sqrt(sd / float(kwin - 1))
        endif
c
        f1 = f0 - sd
        f2 = f0 + sd
c
        f0 = exp(f0)
        f1 = exp(f1)
        f2 = exp(f2)
c
c
        return
c
        end