subroutine thermo_length
c Integrating |PMF| to claculate thermodynamic length
      implicit none
      include 'DIMENSIONS.PMF'
      include 'COMMON.IOUNITS'
      include 'COMMON.MD'
      include 'COMMON.PMF'
      double precision aL(0:maxHdim,maxT_h)
      double precision aLdelta,qdelta,qtemp,
      integer nT,tmin,tmax,i,j,t,tlow,tup,nW
      qdelta=(qmax-qmin)/(nW-1)
      do j=1,nT
        do t=tmin,tmax
          aL(t,j)=0.0d0
        enddo
        dH=hfin(tmin+1,j)-hfin(tmin,j)
        print *,"tmin",tmin,hfin(tmin+1,j),hfin(tmin,j)," dH",dH
        do t=tmin-1,0,-1
          hfin(t,j)=hfin(t+1,j)-dH
        enddo
        aL(0,j)=0.5d0*(dabs(hfin(0,j))+dabs(hfin(1,j)))
        do t=1,tmax-1
          aL(t,j)=aL(t-1,j)+0.5d0*(dabs(hfin(t,j))+dabs(hfin(t+1,j)))
        enddo
        aL(tmax,j)=aL(tmax-1,j)
        do t=0,tmax
         aL(t,j)=(aL(t,j)-0.5d0*(dabs(hfin(0,j))+dabs(hfin(1,j))))*delta
        enddo
        write (*,*) "T=",1.0d0/(1.987D-3*beta_h(j))," PMF and L"
        do t=0,tmax 
          write (*,'(f5.2,2f10.5)') t*delta,hfin(t,j),aL(t,j)
        enddo
c Determining window centers
        tlow=(qmin+1.0d-4)/delta
        tup=(qmax+1.0d-4)/delta
        aLspan=aL(tup,j)-aL(tlow,j)
        aLdelta=aLspan/(nW-1)
        print *,"tlow",tlow," tup",tup," aLspan",aLspan,
     &   " aLdelta",aLdelta
        q0(1,j)=qmin
        q0(nW,j)=qmax
        do i=2,nW
          aLw=aLdelta*(i-1)
          print *,"i",i," aLw",aLw
          do t=tlow,tup-1
            if (aLw.ge.aL(t,j).and.aLw.le.aL(t+1,j)) then
              print *,"i",i," t",t," aLw",aLw,aL(t,j),aL(t+1,j)
              q0(i,j)=t*delta+delta*(aLw-aL(t,j))/(aL(t+1,j)-aL(t,j))
              exit
            endif
          enddo
        enddo
        write (*,*) "q0"
        write (*,'(20f10.5)') (q0(i,j),i=1,nW)
        return
        end
c Force constants
        do t=1,tmax-1
          Wprim=(hfin(t+1,j)-hfin(t-1,j))/(2*delta)
          Wbis=(hfin(t+1,j)+hfin(t-1,j)-2*hfin(t,j))/delta**2
          keff=0.5d0*Wprim**2+dabs(Wbis)
     &    +dsqrt((0.5d0*Wprim**2+dabs(Wbis))**2-Wbis**2)
          write (*,'(f5.2,3f10.1)') t*delta,Wprim,Wbis,keff
        enddo
      enddo ! j
      end