diff --git a/Lobatto.f90 b/Lobatto.f90
index 9f34793..7ee079f 100644
--- a/Lobatto.f90
+++ b/Lobatto.f90
@@ -14,8 +14,10 @@ module Lobatto
 
  subroutine AllLobattoKE(KE,npnt2)
   implicit none
-  real(rk) :: KE(1:npnt2,1:npnt2),a0,rmin
-  integer(ik)::       bra,ket,Ntot,Ntotp2,npnt2
+  integer(ik),intent(in) :: npnt2
+  real(rk),intent(out) :: KE(1:npnt2,1:npnt2)
+  integer(ik)::       bra,ket,Ntot,Ntotp2
+  real(rk) :: a0,rmin
   real(rk),allocatable,dimension(:) ::      rpt, w
   real(rk), allocatable, dimension(:,:) :: derLobMat
   !
@@ -50,10 +52,10 @@ end subroutine AllLobattoKE
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  !
  !
-  subroutine OlLobattoDVR(result,bra,ket,Ntot)
+  subroutine OlLobattoDVR(result,bra,ket)
     implicit none
     real(rk), intent(out)::   result
-    integer(ik), intent(in)::   bra, ket, Ntot
+    integer(ik), intent(in)::   bra, ket
     !
     !
      if (bra == ket) then
@@ -71,18 +73,19 @@ end subroutine OlLobattoDVR
 
 subroutine KeLobattoDVR(result,bra,ket,Ntot,rpt,w,derLobMat)
     implicit none
-    real(rk), intent(out)::   result,rpt(0:Ntot+1),w(0:Ntot+1)
     integer(ik), intent(in)::   bra,ket,Ntot
+    real(rk), intent(out)::   result,rpt(0:Ntot+1),w(0:Ntot+1)
     real(rk), intent(in)::    derLobMat(0:Ntot+1,0:Ntot+1)
     real(rk)::                pf, intresult
-    integer(ik)::               k
-
-
-pf=1d0!/sqrt(w(bra)*w(ket))
-intresult=0d0;
-do k=0,Ntot+1
-    intresult=intresult+w(k)*derLobMat(bra,k)*derLobMat(ket,k)
-end do
+    integer(ik)::             k
+    !
+    rpt = 0
+    w = 0
+    pf=1d0!/sqrt(w(bra)*w(ket))
+    intresult=0d0;
+    do k=0,Ntot+1
+      intresult=intresult+w(k)*derLobMat(bra,k)*derLobMat(ket,k)
+    end do
 
 
 result=intresult*pf;
@@ -117,21 +120,21 @@ end subroutine rminus2LobattoDVR
 
 subroutine derLobattoMat(result,Ntot,rpt,w)
         implicit none
-        real(rk), intent(out)     ::      result(0:Ntot+1,0:Ntot+1)
         integer(ik), intent(in)     ::      Ntot
+        real(rk), intent(out)     ::      result(0:Ntot+1,0:Ntot+1)
         real(rk), intent(in)      ::      rpt(0:Ntot+1),w(0:Ntot+1)
         real(rk)                ::      elementres
         integer(ik)                 ::      n,eta,k
 
-	do n=0,Ntot+1
-        	do eta=n,Ntot+1
-        	        call derLobatto(result(n,eta),n,eta,Ntot,rpt,w)
-			if(n.ne.eta) result(eta,n) = (-1.0d0) * (w(eta)/w(n)) * result(n,eta) 
-			result(n,eta) = result(n,eta) * 1.0d0/sqrt(w(n))
-			if(n.ne.eta) result(eta,n) = result(eta,n) * 1.0d0/sqrt(w(eta))
-        	end do
-	end do
-	return
+        do n=0,Ntot+1
+          do eta=n,Ntot+1
+             call derLobatto(result(n,eta),n,eta,Ntot,rpt,w)
+             if(n.ne.eta) result(eta,n) = (-1.0d0) * (w(eta)/w(n)) * result(n,eta) 
+             result(n,eta) = result(n,eta) * 1.0d0/sqrt(w(n))
+             if(n.ne.eta) result(eta,n) = result(eta,n) * 1.0d0/sqrt(w(eta))
+          end do
+        end do
+        return
 end subroutine derLobattoMat
 
 
@@ -142,41 +145,41 @@ subroutine derLobatto(result,n,eta,Ntot,rpt,w)
 !    Calculates u'_n(r_eta)
 !!!!!!
         implicit none
-	real(rk), intent(out)::   result
+    real(rk), intent(out)::   result
         integer(ik), intent(in)::   n, Ntot,eta
-	real(rk), intent(in)::    w(0:Ntot+1), rpt(0:Ntot+1)
-	real(rk)::                pf, intresult,intpf, factor(0:Ntot+1),signrecorder
+    real(rk), intent(in)::    w(0:Ntot+1), rpt(0:Ntot+1)
+    real(rk)::                pf, intresult,intpf, factor(0:Ntot+1),signrecorder
         integer(ik)::               j,m,Nhalf
 
-	Nhalf=(Ntot+2)/2
-	if (n==eta) then
-        	if (n==0) then
-        	        result=-1d0/(2d0*w(n))
-        	else if (n==Ntot+1) then
-        	        result=1d0/(2d0*w(n))
-        	else
-        	        result=0d0;
-        	end if
-	else
-        	result=1d0/((rpt(n)-rpt(eta)))
-		signrecorder = 1.0
-        	do j=0,Ntot+1
-        	        if ((j .ne. n) .and. (j .ne. eta)) then
-				factor(j) = (rpt(eta)-rpt(j))/(rpt(n)-rpt(j))
-				signrecorder = signrecorder * (factor(j) / abs(factor(j)))
-				factor(j) = abs(factor(j))
-			else
-				 factor(j) = 1.0d0
-        	        end if
-        	end do  
-		call QuickSort(factor,1,Ntot+2)
-		call Riffle(factor,Ntot+2,Nhalf)
-		do j = 0,Ntot+1
-			result = result * factor(j)
-		end do
-		result = result * signrecorder
-	end if
-	return
+    Nhalf=(Ntot+2)/2
+    if (n==eta) then
+            if (n==0) then
+                    result=-1d0/(2d0*w(n))
+            else if (n==Ntot+1) then
+                    result=1d0/(2d0*w(n))
+            else
+                    result=0d0;
+            end if
+    else
+            result=1d0/((rpt(n)-rpt(eta)))
+        signrecorder = 1.0
+            do j=0,Ntot+1
+                    if ((j .ne. n) .and. (j .ne. eta)) then
+               factor(j) = (rpt(eta)-rpt(j))/(rpt(n)-rpt(j))
+               signrecorder = signrecorder * (factor(j) / abs(factor(j)))
+               factor(j) = abs(factor(j))
+           else
+                factor(j) = 1.0d0
+                    end if
+            end do  
+        call QuickSort(factor,1,Ntot+2)
+        call Riffle(factor,Ntot+2,Nhalf)
+        do j = 0,Ntot+1
+           result = result * factor(j)
+        end do
+        result = result * signrecorder
+    end if
+    return
 end subroutine derLobatto
 
 
@@ -232,57 +235,58 @@ end subroutine LobattoAbsWeights
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 recursive subroutine QuickSort(a, first, last)
-	implicit none
-	real(rk) a(*), x, t
-	integer(ik) first, last
-	integer(ik) i,j
-
-	x = a( (first + last) / 2)
-	i = first
-	j = last
-	do
-		do while (a(i) < x)
-			i = i+1
-		end do
-		do while (x < a(j))
-		j = j-1
-		end do
-		if (i >= j) exit
-		t = a(i); a(i) = a(j); a(j) = t
-		i = i+1
-		j = j-1
-	end do
-	if (first < i-1) call quicksort(a,first,i-1)
-	if (j+1 < last) call quicksort(a,j+1,last)
+    implicit none
+    real(rk) a(*), x, t
+    integer(ik) first, last
+    integer(ik) i,j
+
+    x = a( (first + last) / 2)
+    i = first
+    j = last
+    do
+        do while (a(i) < x)
+           i = i+1
+        end do
+        do while (x < a(j))
+        j = j-1
+        end do
+        if (i >= j) exit
+        t = a(i); a(i) = a(j); a(j) = t
+        i = i+1
+        j = j-1
+    end do
+    if (first < i-1) call quicksort(a,first,i-1)
+    if (j+1 < last) call quicksort(a,j+1,last)
 end subroutine quicksort
 
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 subroutine Riffle(a,NAbsc,Nhalf)
-	implicit none
-	real(rk) ::  a(1:NAbsc),b(1:Nhalf),c(1:Nhalf)
-	integer(ik) NAbsc,i,Nhalf,Nhalfplusone
-
-	Nhalfplusone = Nhalf+1
-	if(mod(NAbsc,2).eq.0) then
-		do i=1,Nhalf
-			b(i) = a(i)
-			c(i) = a(NAbsc-i+1)
-		enddo	
-		do i = 1,Nhalf
-			a(2*i) = c(i)
-			a(2*i-1) = b(i)
-		enddo
-	else
-		do i=1,Nhalfplusone
-			b(i) = a(i)
-			c(i) = a(NAbsc-i+1)
-		enddo	
-		do i = 1,Nhalfplusone
-			if(i.ne.Nhalfplusone) a(2*i) = c(i)
-			a(2*i-1) = b(i)
-		enddo
-	endif	
+    implicit none
+    integer(ik),intent(in) :: NAbsc,Nhalf
+    real(rk) ::  a(1:NAbsc),b(1:Nhalf),c(1:Nhalf)
+    integer(ik) i,Nhalfplusone
+
+    Nhalfplusone = Nhalf+1
+    if(mod(NAbsc,2).eq.0) then
+        do i=1,Nhalf
+           b(i) = a(i)
+           c(i) = a(NAbsc-i+1)
+        enddo    
+        do i = 1,Nhalf
+           a(2*i) = c(i)
+           a(2*i-1) = b(i)
+        enddo
+    else
+        do i=1,Nhalfplusone
+           b(i) = a(i)
+           c(i) = a(NAbsc-i+1)
+        enddo    
+        do i = 1,Nhalfplusone
+           if(i.ne.Nhalfplusone) a(2*i) = c(i)
+           a(2*i-1) = b(i)
+        enddo
+    endif
 end subroutine Riffle
 
 end module Lobatto
\ No newline at end of file
diff --git a/MAKEFILES/makefile_g95 b/MAKEFILES/makefile_g95
index 47a25db..2c53f43 100644
--- a/MAKEFILES/makefile_g95
+++ b/MAKEFILES/makefile_g95
@@ -18,7 +18,7 @@ LIB     =   $(LAPACK)
 
 ###############################################################################
 
-OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o  symmetry.o dipole.o header_info.o
+OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o  symmetry.o dipole.o header_info.o Lobatto.o
 #compilation_details.o 
 
 duo.x:	$(OBJ) duo.o
@@ -27,10 +27,10 @@ duo.x:	$(OBJ) duo.o
 duo.o:	duo.f90 $(OBJ) 
 	$(FOR) -c duo.f90 $(FFLAGS)
 
-grids.o:	grids.f90 accuracy.o input.o
+grids.o:	grids.f90 accuracy.o input.o Lobatto.o
 	$(FOR) -c grids.f90 $(FFLAGS)
 
-diatom.o:	diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o
+diatom.o:	diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o Lobatto.o
 	$(FOR) -c diatom.f90 $(FFLAGS)
 
 refinement.o:	refinement.f90 accuracy.o input.o lapack.o diatom.o
@@ -39,7 +39,7 @@ refinement.o:	refinement.f90 accuracy.o input.o lapack.o diatom.o
 functions.o:	functions.f90 accuracy.o input.o lapack.o
 	$(FOR) -c functions.f90 $(FFLAGS)
 
-dipole.o:	dipole.f90 accuracy.o input.o lapack.o diatom.o
+dipole.o:	dipole.f90 accuracy.o input.o lapack.o diatom.o 
 	$(FOR) -c dipole.f90 $(FFLAGS)
 
 accuracy.o:  accuracy.f90
@@ -66,5 +66,8 @@ header_info.o:  accuracy.o
 atomic_and_nuclear_data.o:  atomic_and_nuclear_data.f90
 	$(FOR) -c atomic_and_nuclear_data.f90 $(FFLAGS)
 
+Lobatto.o: Lobatto.f90 timer.o
+	$(FOR) -c Lobatto.f90 $(FFLAGS)
+
 clean:
 	rm -f $(OBJ) *.mod *__genmod.f90 duo.o
diff --git a/MAKEFILES/makefile_gfortran b/MAKEFILES/makefile_gfortran
index 7aeea7e..b97cfbd 100644
--- a/MAKEFILES/makefile_gfortran
+++ b/MAKEFILES/makefile_gfortran
@@ -20,7 +20,7 @@ LIB     =   $(LAPACK)
 
 ###############################################################################
 
-OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o  symmetry.o dipole.o header_info.o
+OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o  symmetry.o dipole.o header_info.o Lobatto.o
 #compilation_details.o 
 
 duo.x:	$(OBJ) duo.o
@@ -29,10 +29,10 @@ duo.x:	$(OBJ) duo.o
 duo.o:	duo.f90 $(OBJ) 
 	$(FOR) -c duo.f90 $(FFLAGS)
 
-grids.o:	grids.f90 accuracy.o input.o
+grids.o:	grids.f90 accuracy.o input.o Lobatto.o
 	$(FOR) -c grids.f90 $(FFLAGS)
 
-diatom.o:	diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o
+diatom.o:	diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o Lobatto.o
 	$(FOR) -c diatom.f90 $(FFLAGS)
 
 refinement.o:	refinement.f90 accuracy.o input.o lapack.o diatom.o
@@ -41,7 +41,7 @@ refinement.o:	refinement.f90 accuracy.o input.o lapack.o diatom.o
 functions.o:	functions.f90 accuracy.o input.o lapack.o
 	$(FOR) -c functions.f90 $(FFLAGS)
 
-dipole.o:	dipole.f90 accuracy.o input.o lapack.o diatom.o
+dipole.o:	dipole.f90 accuracy.o input.o lapack.o diatom.o 
 	$(FOR) -c dipole.f90 $(FFLAGS)
 
 accuracy.o:  accuracy.f90
@@ -67,6 +67,9 @@ header_info.o:  accuracy.o
 
 atomic_and_nuclear_data.o:  atomic_and_nuclear_data.f90
 	$(FOR) -c atomic_and_nuclear_data.f90 $(FFLAGS)
-    
+
+Lobatto.o: Lobatto.f90 timer.o
+	$(FOR) -c Lobatto.f90 $(FFLAGS)
+
 clean:
 	rm -f $(OBJ) *.mod *__genmod.f90 duo.o
diff --git a/MAKEFILES/makefile_nag b/MAKEFILES/makefile_nag
index 446eef7..76eb649 100644
--- a/MAKEFILES/makefile_nag
+++ b/MAKEFILES/makefile_nag
@@ -19,19 +19,19 @@ LIB     =   $(LAPACK)
 
 ###############################################################################
 
-OBJ = grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o  symmetry.o dipole.o header_info.o
+OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o  symmetry.o dipole.o header_info.o Lobatto.o
 #compilation_details.o 
 
 duo.x:	$(OBJ) duo.o
-	$(FOR) -o duo$(PLAT).x $(OBJ) $(FFLAGS) duo.o $(LIB)
+	$(FOR) -o j-duo$(PLAT).x $(OBJ) $(FFLAGS) duo.o $(LIB)
 
 duo.o:	duo.f90 $(OBJ) 
 	$(FOR) -c duo.f90 $(FFLAGS)
 
-grids.o:	grids.f90 accuracy.o input.o
+grids.o:	grids.f90 accuracy.o input.o Lobatto.o
 	$(FOR) -c grids.f90 $(FFLAGS)
 
-diatom.o:	diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o
+diatom.o:	diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o Lobatto.o
 	$(FOR) -c diatom.f90 $(FFLAGS)
 
 refinement.o:	refinement.f90 accuracy.o input.o lapack.o diatom.o
@@ -40,7 +40,7 @@ refinement.o:	refinement.f90 accuracy.o input.o lapack.o diatom.o
 functions.o:	functions.f90 accuracy.o input.o lapack.o
 	$(FOR) -c functions.f90 $(FFLAGS)
 
-dipole.o:	dipole.f90 accuracy.o input.o lapack.o diatom.o
+dipole.o:	dipole.f90 accuracy.o input.o lapack.o diatom.o 
 	$(FOR) -c dipole.f90 $(FFLAGS)
 
 accuracy.o:  accuracy.f90
@@ -64,5 +64,11 @@ input.o:  input.f90
 header_info.o:  accuracy.o
 	$(FOR) -c header_info.f90 $(FFLAGS)
 
+atomic_and_nuclear_data.o:  atomic_and_nuclear_data.f90
+	$(FOR) -c atomic_and_nuclear_data.f90 $(FFLAGS)
+
+Lobatto.o: Lobatto.f90 timer.o
+	$(FOR) -c Lobatto.f90 $(FFLAGS)
+
 clean:
 	rm -f $(OBJ) *.mod *__genmod.f90 duo.o
diff --git a/MAKEFILES/makefile_pgi b/MAKEFILES/makefile_pgi
index 90332fe..677435d 100644
--- a/MAKEFILES/makefile_pgi
+++ b/MAKEFILES/makefile_pgi
@@ -19,19 +19,19 @@ LIB     =   $(LAPACK)
 
 ###############################################################################
 
-OBJ = grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o  symmetry.o dipole.o header_info.o
+OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o  symmetry.o dipole.o header_info.o Lobatto.o
 #compilation_details.o 
 
 duo.x:	$(OBJ) duo.o
-	$(FOR) -o duo$(PLAT).x $(OBJ) $(FFLAGS) duo.o $(LIB)
+	$(FOR) -o j-duo$(PLAT).x $(OBJ) $(FFLAGS) duo.o $(LIB)
 
 duo.o:	duo.f90 $(OBJ) 
 	$(FOR) -c duo.f90 $(FFLAGS)
 
-grids.o:	grids.f90 accuracy.o input.o
+grids.o:	grids.f90 accuracy.o input.o Lobatto.o
 	$(FOR) -c grids.f90 $(FFLAGS)
 
-diatom.o:	diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o
+diatom.o:	diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o Lobatto.o
 	$(FOR) -c diatom.f90 $(FFLAGS)
 
 refinement.o:	refinement.f90 accuracy.o input.o lapack.o diatom.o
@@ -40,7 +40,7 @@ refinement.o:	refinement.f90 accuracy.o input.o lapack.o diatom.o
 functions.o:	functions.f90 accuracy.o input.o lapack.o
 	$(FOR) -c functions.f90 $(FFLAGS)
 
-dipole.o:	dipole.f90 accuracy.o input.o lapack.o diatom.o
+dipole.o:	dipole.f90 accuracy.o input.o lapack.o diatom.o 
 	$(FOR) -c dipole.f90 $(FFLAGS)
 
 accuracy.o:  accuracy.f90
@@ -64,5 +64,11 @@ input.o:  input.f90
 header_info.o:  accuracy.o
 	$(FOR) -c header_info.f90 $(FFLAGS)
 
+atomic_and_nuclear_data.o:  atomic_and_nuclear_data.f90
+	$(FOR) -c atomic_and_nuclear_data.f90 $(FFLAGS)
+
+Lobatto.o: Lobatto.f90 timer.o
+	$(FOR) -c Lobatto.f90 $(FFLAGS)
+
 clean:
 	rm -f $(OBJ) *.mod *__genmod.f90 duo.o
diff --git a/diatom.f90 b/diatom.f90
index 13ef75b..967be27 100644
--- a/diatom.f90
+++ b/diatom.f90
@@ -3476,7 +3476,7 @@ subroutine ReadInput
       stop "Illegal number of states: ipo/=nestates"
     endif
     !
-    if (trim(solution_method)=="LOBATTO".xor.grid%nsub == 6) then 
+    if ((trim(solution_method)=="LOBATTO".and.grid%nsub /= 6).or.(trim(solution_method)/="LOBATTO".and.grid%nsub == 6)) then 
         write(out,"('Error: The grid 6 can be used without LOBATTO  (key word SOLUTIONMETHOD)')")
         write(out,"('solution_method, grid = ',a,i8)") trim(solution_method),grid%nsub 
     endif
@@ -5704,12 +5704,30 @@ subroutine duo_j0(iverbose_,J_list_,enerout,quantaout,nenerout)
                 stop
               endif
               vibmat(igrid,igrid) = vibmat(igrid,igrid) +(12._rk)* PI**2 / 3._rk
-
-               do jgrid =igrid+1, ngrid
-                 vibmat(igrid,jgrid) = +(12._rk)*2._rk* real( (-1)**(igrid+jgrid), rk) / real(igrid - jgrid, rk)**2
+              !
+              do jgrid =igrid+1, ngrid
+                vibmat(igrid,jgrid) = +(12._rk)*2._rk* real( (-1)**(igrid+jgrid), rk) / real(igrid - jgrid, rk)**2
+                vibmat(jgrid,igrid) = vibmat(igrid,jgrid)
+              enddo
+              !
+            case("LOBATTO") ! Implements a DVR method based on Lobatto quadrature
+                            ! Requires the Lobatto nonuniform grid to work
+              if(grid%nsub /= 6) then
+                write(out, '(A)') 'The Lobatto DVR method only works with the'
+                write(out, '(A)') 'Lobatto grid (grid type 6).'
+                stop
+              endif
+              !
+              do jgrid = igrid,ngrid
+                 do kgrid=1,ngrid
+                    vibTmat(igrid,jgrid) = vibTmat(igrid,jgrid) + (12._rk)*(hstep**2)*(LobWeights(kgrid))*& 
+                                           LobDerivs(igrid,kgrid)*LobDerivs(jgrid,kgrid)
+                 enddo
+                 vibTmat(jgrid,igrid) = vibTmat(igrid,jgrid)
+                 vibmat(igrid,jgrid) = vibmat(igrid,jgrid) + vibTmat(igrid,jgrid)
                  vibmat(jgrid,igrid) = vibmat(igrid,jgrid)
-               enddo
-
+              enddo
+              !
             case default
              write(out, '(A)') 'Error: unrecognized solution method' // trim(solution_method)
              write(out, '(A)') 'Possible options are: '
@@ -5717,7 +5735,7 @@ subroutine duo_j0(iverbose_,J_list_,enerout,quantaout,nenerout)
              write(out, '(A)') '                      SINC'
              write(out, '(A)') '                      LOBATTO'
             end select method_choice
-
+            !
        enddo
        !$omp end parallel do
        !
diff --git a/grids.f90 b/grids.f90
index 66e280b..d98bada 100644
--- a/grids.f90
+++ b/grids.f90
@@ -1,5 +1,7 @@
     subroutine gridred(nsub,alpha,re,mes,rmin,rmax,h,r,z,add,iverbose)
     use accuracy, only : rk, out, bohr
+    use Lobatto,   only : LobattoAbsWeights
+    
     implicit none
 !____________________________________________________________
 !
@@ -18,7 +20,7 @@ subroutine gridred(nsub,alpha,re,mes,rmin,rmax,h,r,z,add,iverbose)
 !
      real(kind=rk) :: h      ! uniform step of the z-coordinate treated
      real(kind=rk) :: r(mes) ! nonuniform grid points of the r-coordinate
-     real(kind=rk) :: z(mes),zmin,zmax,y,s,t
+     real(kind=rk) :: z(mes),zmin,zmax,y,s,t,w(mes)
      real(kind=rk) :: add(mes)
 !      intent(out) zmin,zmax,h,r !output variables
 ! -----------------------------------------------------------------
@@ -110,6 +112,13 @@ subroutine gridred(nsub,alpha,re,mes,rmin,rmax,h,r,z,add,iverbose)
               z(j)    = r(j)/t/re/(1._rk-y)
               add(j)  = (1._rk-(re**2)*(1._rk+t**2))/(2._rk*r(j))**2
            enddo
+!----------------------------------------------------------------------
+      case (6)         !  |   Lobatto quadrature grid points |
+!----------------------------------------------------------------------
+           CALL LobattoAbsWeights(r,w,mes,rmin,rmax)
+           zmin = rmin 
+           zmax = rmax 
+           h    = ( zmax - zmin )/real( mes - 1, rk)
 !----------------------------------------------------------------------
      case (0)      !   uniform grid points: z=r |
 !----------------------------------------------------------------------
diff --git a/timer.f90 b/timer.f90
index ecdbb0c..bb2c7be 100644
--- a/timer.f90
+++ b/timer.f90
@@ -681,7 +681,7 @@ subroutine ArrayMinus(name,isize,ikind,hsize)
       size_ = (ikind*real(hsize_))/real(1024**3) ! size in GByte
       !
       memory_now   = memory_now - size_
-      t%size = max(t%size-size_,0.0_hik)
+      t%size = max(t%size-size_,0.0_rk)
       mem = memory_now
       !
       !if (t%size<=0) then