Clean up the code using gfortran and removing not used variables

Lobatto.f90

@@ -147,8 +147,8 @@ subroutine derLobatto(result,n,eta,Ntot,rpt,w)
  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):: intpf, factor(0:Ntot+1),signrecorder
- integer(ik):: j,m,Nhalf
+ real(rk):: factor(0:Ntot+1),signrecorder
+ integer(ik):: j,Nhalf
 
  Nhalf=(Ntot+2)/2
  if (n==eta) then

diatom.f90

@@ -5,7 +5,7 @@ module diatom_module
  use functions, only : analytical_fieldT
  use symmetry, only : sym,SymmetryInitialize
  use Lobatto, only : LobattoAbsWeights,derLobattoMat
- use me_numer, only : ME_numerov
+ !use me_numer, only : ME_numerov
  !
  implicit none
  ! by Lorenzo Lodi
@@ -6323,8 +6323,12 @@ subroutine duo_j0(iverbose_,J_list_,enerout,quantaout,nenerout)
  !
  mu_rr = 2.0_rk*b_rot
  !
- call ME_numerov(nroots,(/grid%rmin,grid%rmax/),ngrid-1,ngrid-1,r,poten(istate)%gridvalue,mu_rr,1,0,&
- job%vibenermax(istate),iverbose,vibener,vibmat)
+ !
+ write(out,"('Error: ME_numerov is not implemented yet')")
+ stop 'Error: ME_numerov is not implemented yet'
+ !
+ !call ME_numerov(nroots,(/grid%rmin,grid%rmax/),ngrid-1,ngrid-1,r,poten(istate)%gridvalue,mu_rr,1,0,&
+ ! job%vibenermax(istate),iverbose,vibener,vibmat)
  !
  deallocate(mu_rr)
  call ArrayStop('mu_rr')

dipole.f90

@@ -348,7 +348,7 @@ subroutine dm_intensity(Jval,iverbose)
  real(rk), allocatable :: vecI(:), vecF(:)
  real(rk),allocatable :: half_linestr(:)
  !
- integer(ik) :: jind,nlevels,j1,j2
+ integer(ik) :: jind,nlevels
  !
  integer(ik) :: iroot,NlevelsI,NlevelsF,nlower,k,k_,iLF,iflag_rich
  !
@@ -374,7 +374,7 @@ subroutine dm_intensity(Jval,iverbose)
  !
  integer(ik) :: alloc_p
  !
- integer(ik) :: Jmax_,ID_J,nMs
+ integer(ik) :: Jmax_,ID_J
  real(rk) :: J_
  character(len=12) :: char_Jf,char_Ji,char_LF
  integer(ik),allocatable :: richunit(:,:)
@@ -1370,7 +1370,7 @@ subroutine find_igamma_pair(igamma_pair)
  !
  enddo
  !
- if ( intensity%gns(igammaI)/=intensity%gns(igamma_pair(igammaI)) ) then 
+ if ( nint(intensity%gns(igammaI)-intensity%gns(igamma_pair(igammaI)))/=0 ) then 
  !
  write(out,"('dm_intensity: selection rules do not agree with Gns')")
  stop 'dm_intensity: selection rules do not agree with Gns!'
@@ -1511,7 +1511,7 @@ subroutine intens_filter(jI,jF,energyI,energyF,isymI,isymF,igamma_pair,passed)
  !
  passed = passed.and. &
  !
- (jF/=intensity%J(1).or.jI/=intensity%J(1).or.nint(jI+jF)==1).and.  &
+ (nint(jF-intensity%J(1))/=0.or.nint(jI-intensity%J(1))/=0.or.nint(jI+jF)==1).and. &
  !
  !( ( nint(jF-intensity%J(1))/=0.or.nint(jI-intensity%J(1))/=0 ).and.intensity%J(1)>0 ).and. &
  ( intensity%J(1)+intensity%J(2)>0 ).and. &
@@ -1584,7 +1584,7 @@ subroutine matelem_filter(jI,jF,energyI,energyF,isymI,isymF,igamma_pair,passed)
  !
  passed = passed.and. &
  !
- (jF/=intensity%J(1).or.jI/=intensity%J(1).or.nint(jI+jF)==1).and. &
+ (nint(jF-intensity%J(1))/=0.or.nint(jI-intensity%J(1))/=0.or.nint(jI+jF)==1).and. &
  !
  !( ( nint(jF-intensity%J(1))/=0.or.nint(jI-intensity%J(1))/=0 ).and.intensity%J(1)>0 ).and. &
  ( intensity%J(1)+intensity%J(2)>0 ).and. &
@@ -2007,7 +2007,7 @@ end subroutine do_1st_half_tm
 
 
 
- subroutine  (iLF,iunit,jI,jF,icount)
+ subroutine do_LF_matrix_elements(iLF,iunit,jI,jF,icount)
  implicit none 
  real(rk),intent(in) :: jI,jF
  integer(ik),intent(in) :: iLF,iunit
@@ -2129,7 +2129,8 @@ function three_j0(a,b,c,al,be,ga)
  if(c.lt.abs(a-b)) return
  if(a.lt.0.or.b.lt.0.or.c.lt.0) return
  if(a.lt.abs(al).or.b.lt.abs(be).or.c.lt.abs(ga)) return
- if(-ga.ne.al+be) return
+ !if(-ga.ne.al+be) return
+ if(nint(ga+al+be).ne.0) return
 !
 !
 ! compute delta(abc)

polarizability.f90

@@ -301,15 +301,13 @@ subroutine pol_intensity(Jval, iverbose)
  real(rk), allocatable :: vecI(:), vecF(:)
  real(rk), allocatable :: half_linestr(:), half_linestr_zero(:)
  !
- integer(ik) :: jind, nlevels, j1, j2
+ integer(ik) :: jind, nlevels
  !
- integer(ik) :: iroot, NlevelsI, NlevelsF, nlower, k, k_, iLF, iflag_rich
+ integer(ik) :: iroot, NlevelsI, NlevelsF, nlower,iLF, iflag_rich
  !
  integer(ik) :: igamma_pair(sym%Nrepresen), igamma, istateI, istateF,&
  ivibI, ivibF, ivI, ivF, ilambdaI, ilambdaF, iparityI, itau
- integer(ik) :: ivF_, ilambdaF_
- real(rk) :: spinI, spinF, omegaI, omegaF, sigmaI, sigmaF,& 
- sigmaF_, omegaF_, spinF_
+ real(rk) :: spinI, spinF, omegaI, omegaF, sigmaI, sigmaF
  integer(hik) :: matsize
  !
  character(len=1) :: branch, ef, pm
@@ -319,7 +317,7 @@ subroutine pol_intensity(Jval, iverbose)
  type(quantaT),pointer :: quantaI, quantaF
  !
  real(rk) :: boltz_fc, beta, intens_cm_mol, emcoef,&
- A_coef_s_1, A_einst, absorption_int, lande
+ A_coef_s_1, A_einst, absorption_int
  !
  character(len=130) :: my_fmt !format for I/O specification
  !
@@ -331,7 +329,7 @@ subroutine pol_intensity(Jval, iverbose)
  !
  integer(ik) :: alloc_p
  !
- integer(ik) :: Jmax_, ID_J, nMs
+ integer(ik) :: Jmax_, ID_J
  real(rk) :: J_
  !
  character(len=12) :: char_Jf,char_Ji,char_LF

refinement.f90

@@ -1446,7 +1446,7 @@ subroutine sf_fitting
  !
  if (fitting%robust>small_) then
  !
- call robust_fit(numpar,a_wats,sigma(1:npts),eps(1:npts),wtall(1:npts))
+ call robust_fit(a_wats,sigma(1:npts),eps(1:npts),wtall(1:npts))
  !
  ssq=sum(eps(1:npts)*eps(1:npts)*wtall(1:npts))
  !
@@ -2023,9 +2023,8 @@ subroutine MLinvmat(al,ai,dimen,ierr)
  end subroutine MLinvmat
 
 
- subroutine Robust_fit(numpar,a_wats,sigma,eps,wt)
+ subroutine Robust_fit(a_wats,sigma,eps,wt)
 
- integer(ik),intent(in) :: numpar
  real(rk),intent(inout) :: a_wats
  real(rk),intent(in) :: sigma(:),eps(:)
  real(rk),intent(inout) :: wt(:)
@@ -2044,26 +2043,6 @@ subroutine Robust_fit(numpar,a_wats,sigma,eps,wt)
  !
  !Watson alpha-parameter
  ! 
- ! Comment by Lorenzo Lodi: the following loop is never executed!
- !do i = 1,-1
- ! !
- ! da1 = 0
- ! da2 = 0
- ! !
- ! do nrow=1,npts
- ! if (wt(nrow)>small_) then 
- ! da1 = da1+eps(nrow)**2/( sigma(nrow)**2+a_wats*eps(nrow)**2 )
- ! da2 = da2+eps(nrow)**4/( sigma(nrow)**2+a_wats*eps(nrow)**2 )**2
- ! endif 
- ! enddo 
- ! !
- ! da =( da1 - real(nused-numpar,rk) )/da2
- ! a_wats = a_wats + da
- ! !
- ! if (a_wats<sqrt(small_)) a_wats = 1e-3+real(i,rk)*1e-2
- ! !
- !enddo
- !
  a_wats = 0.001_rk
  !
  if (verbose>=4) write(out,"('Watson parameter =',f18.8)") a_wats