diff --git a/diag.f90 b/diag.F90
similarity index 98%
rename from diag.f90
rename to diag.F90
index 722c0d1..d6b559e 100644
--- a/diag.f90
+++ b/diag.F90
@@ -1,10 +1,10 @@
 module diag
- 
-!dec$ define arpack_ = 0
-!dec$ define blacs_  = 0
-!dec$ define mpi_    = 0
-!dec$ define omparpack_ = 0
-!dec$ define propack_ = 0
+
+#define arpack_    0
+#define blacs_     0
+#define mpi_       0
+#define omparpack_ 0
+#define propack_   0
 
 !
 !  Simplistic type-agnostic LAPACK interface
@@ -2088,17 +2088,17 @@ subroutine diag_dseupd(n,bterm,nroots,factor,maxitr_,tol,h,e)
         ! has been exceeded.                          
         !
 
-        !dec$ if (arpack_ > 0)
+#if arpack_ > 0
            !
            call dsaupd ( ido, bmat, n, which, nev, tol, resid, &
                          ncv, v, ldv, iparam, ipntr, workd, workl, &
                          lworkl, info )
-        !dec$ else 
+#else 
             !
             write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
             stop 'Arpack was not activated'
             !
-        !dec$ end if
+#endif
         !
         !if (verbose>=4.and.iparam(5)>0) then 
         !  !
@@ -2134,7 +2134,7 @@ subroutine diag_dseupd(n,bterm,nroots,factor,maxitr_,tol,h,e)
       !             
       rvec = .true.
       !
-      !dec$ if (arpack_ > 0)
+#if arpack_ > 0
         !
         if (verbose>=5) write(out,"(/'Arpack: dseupd')") 
         !
@@ -2148,12 +2148,12 @@ subroutine diag_dseupd(n,bterm,nroots,factor,maxitr_,tol,h,e)
         !
         if (verbose>=5) write(out,"(/'Arpack: done!')") 
         !
-      !dec$ else 
+#else 
           !
           write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
           stop 'Arpack was not activated'
           !
-      !dec$ end if
+#endif
       !
       if ( ierr < 0 ) then
          write(out,"(/'Error with_seupd, info = ',i8)") ierr
@@ -2439,17 +2439,17 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
 	  myid = 1
 	  nprow = 1
       !
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         call BLACS_GRIDINFO( comm, nprow, npcol, myprow, mypcol )
         myid= myprow
-      !dec$ elseif (mpi_ > 0)
+#elif mpi_ > 0
         call MPI_COMM_RANK( comm, myid, ierr )
         call MPI_COMM_SIZE( comm, nprocs_, ierr )
         if (nprocs_/=nprocs) then
           write(out,"('matvec_p: inconsistent number of  nprocs  = ',2i)") nprocs_,nprocs
           stop 'matvec_p: inconsistent number of  nprocs s'
         endif
-      !dec$ end if
+#endif
       !
       kend = kstart(myid) + nloc-1
       !
@@ -2479,12 +2479,12 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
            !
            nx = iend-istart+1
            !
-           !dec$ if (blacs_ > 0)
+#if blacs_ > 0
              call dgerv2d( comm, nx, 1, mv_buf(istart:iend), nx, iprev, mypcol )
-           !dec$ end if
-           !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
              call mpi_recv(mv_buf(istart),nx,MPI_DOUBLE_PRECISION,iprev,myid,comm,ierr)
-           !dec$ end if
+#endif
            !
            istart = max(istart,bterm(k,1))
            !
@@ -2498,12 +2498,12 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
            !
            nx = iend-istart+1
            !
-           !dec$ if (blacs_ > 0)
+#if blacs_ > 0
              call dgesd2d( comm, nx, 1, z(istart:iend), nx, iprev, mypcol )
-           !dec$ end if
-           !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
              call mpi_send(mv_buf(istart),nx,MPI_DOUBLE_PRECISION,iprev,iprev,comm,ierr)
-           !dec$ end if
+#endif
            !
          enddo
          !
@@ -2520,12 +2520,12 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
            !
            nx = iend-istart+1
            !
-           !dec$ if (blacs_ > 0)
+#if blacs_ > 0
              call dgerv2d( comm, nx, 1, mv_buf(istart:iend), nx, inext, mypcol )
-           !dec$ end if
-           !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
              call mpi_recv(mv_buf(istart),nx,MPI_DOUBLE_PRECISION,inext,myid,comm,ierr)
-           !dec$ end if
+#endif
            !
            iend = min(bterm(k,2),iend)
            !
@@ -2539,12 +2539,12 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
            !
            nx = iend-istart+1
            !
-           !dec$ if (blacs_ > 0)
+#if blacs_ > 0
              call dgesd2d( comm, nx, 1, z(istart:iend), nx, inext, mypcol )
-           !dec$ end if
-           !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
              call mpi_send(mv_buf(istart),nx,MPI_DOUBLE_PRECISION,inext,inext,comm,ierr)
-           !dec$ end if
+#endif
            !
          enddo
          !
@@ -2674,19 +2674,19 @@ subroutine diag_dseupd_p(n,bterm,nroots,factor,maxitr_,tol,h,e)
       write(out,"('Start PARPACK-diagonalization')")
 
 
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         write(out,"('BLAS-PINFO-start')")
         call BLACS_PINFO( iam, nprocs )
         print *,nprocs
         blacs_or_mpi = 'BLACS'
-      !dec$ end if
+#endif
       !
       write(out,"('BLAS-PINFO-done')")
       !
       !call BLACS_PINFO( iam, nprocs )
       !print *,nprocs
 
-      !dec$ if (mpi_ > 0)
+#if mpi_ > 0
         call MPI_INIT( ierr )
         comm = MPI_COMM_WORLD
         call MPI_COMM_RANK( comm, myid, ierr )
@@ -2701,7 +2701,7 @@ subroutine diag_dseupd_p(n,bterm,nroots,factor,maxitr_,tol,h,e)
         !
         blacs_or_mpi = 'MPI'
         !
-      !dec$ end if
+#endif
       !
 
 
@@ -2710,9 +2710,9 @@ subroutine diag_dseupd_p(n,bterm,nroots,factor,maxitr_,tol,h,e)
 !
       if (nprocs .lt. 1) then
          nprocs = 1
-         !dec$ if (blacs_ > 0)
+#if blacs_ > 0
            call BLACS_SETUP( iam, nprocs )
-         !dec$ end if
+#endif
 	 !
 	 print *,nprocs
 	 !
@@ -2801,14 +2801,14 @@ subroutine diag_dseupd_p(n,bterm,nroots,factor,maxitr_,tol,h,e)
       !
       myprow = 1 ; mypcol = 1 ; myid = 1
       !
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         call BLACS_GET( 0, 0, comm )
         call BLACS_GRIDINIT( comm, 'Row', nprow, npcol )
         call BLACS_GRIDINFO( comm, nprow, npcol, myprow, mypcol )
         !
         write(out,"(' myprow, nprow, mypcol, npcol = ',4i8)") myprow, nprow, mypcol, npcol
         !
-      !dec$ end if
+#endif
       !
       if (verbose>=2.and.trim(blacs_or_mpi)=='BLACS') write(out,"('myprow, nprow, mypcol, npcol, nprocs = ',5i8)") myprow, nprow, mypcol, npcol, nprocs
       if (verbose>=2.and.trim(blacs_or_mpi)=='MPI') write(out,"('myid,nproc = ',2i8)") myid,nprocs
@@ -2874,19 +2874,19 @@ subroutine diag_dseupd_p(n,bterm,nroots,factor,maxitr_,tol,h,e)
         ! has been exceeded.                          
         !
 
-        !dec$ if (blacs_ > 0.or.mpi_ > 0)
+#if (blacs_ > 0) || (mpi_ > 0)
            !
            call pdsaupd ( comm, ido, bmat, nloc, which, nev, tol, resid, &
                           ncv, v, ldv, iparam, ipntr, workd, workl, &
                           lworkl, info )
         !
-        !dec$ else 
+#else 
             !
            call dsaupd ( ido, bmat, n, which, nev, tol, resid, &
                           ncv, v, ldv, iparam, ipntr, workd, workl, &
                           lworkl, info )
             !
-        !dec$ end if
+#endif
         !
         !if (verbose>=4.and.iparam(5)>0) then 
         !  !
@@ -2920,7 +2920,7 @@ subroutine diag_dseupd_p(n,bterm,nroots,factor,maxitr_,tol,h,e)
       !             
       rvec = .true.
       !
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         !
         if (verbose>=5) write(out,"(/'Arpack: dseupd')") 
         !
@@ -2932,12 +2932,12 @@ subroutine diag_dseupd_p(n,bterm,nroots,factor,maxitr_,tol,h,e)
         if (verbose>=5) write(out,"(/'Arpack: done!')") 
 
         !                    
-      !dec$ else 
+#else 
           !
           write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
           stop 'Arpack was not activated'
           !
-      !dec$ end if
+#endif
       !
       if ( ierr < 0 ) then
          write(out,"(/'Error with_seupd, info = ',i8)") ierr
@@ -2985,13 +2985,13 @@ subroutine diag_dseupd_p(n,bterm,nroots,factor,maxitr_,tol,h,e)
 !
  9000 continue
 !
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         call BLACS_GRIDEXIT ( comm )
         call BLACS_EXIT(0)
-      !dec$ end if
-      !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
          call MPI_FINALIZE(rc)
-      !dec$ end if
+#endif
 
 
       deallocate(v,workl,workd,d,resid,select,mv_buf)
@@ -3136,17 +3136,17 @@ subroutine dseupd_omp_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
         ! has been exceeded.                          
         !
 
-        !dec$ if (omparpack_ > 0)
+#if omparpack_ > 0
            !
            call dsaupd_ ( ido, bmat, n, which, nev, tol, resid, &
                          ncv, v, ldv, iparam, ipntr, workd, workl, &
                          lworkl, info )
-        !dec$ else 
+#else 
             !
             write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
             stop 'Arpack was not activated'
             !
-        !dec$ end if
+#endif
         !
         !if (verbose>=4.and.iparam(5)>0) then 
         !  !
@@ -3182,7 +3182,7 @@ subroutine dseupd_omp_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
       !             
       rvec = .true.
       !
-      !dec$ if (omparpack_ > 0)
+#if omparpack_ > 0
         !
         if (verbose>=5) write(out,"(/'Arpack: dseupd')") 
         !
@@ -3196,12 +3196,12 @@ subroutine dseupd_omp_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
         !
         if (verbose>=5) write(out,"(/'Arpack: done!')") 
         !
-      !dec$ else 
+#else 
           !
           write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
           stop 'Arpack was not activated'
           !
-      !dec$ end if
+#endif
       !
       if ( ierr < 0 ) then
          write(out,"(/'Error with_seupd, info = ',i8)") ierr
@@ -3294,7 +3294,7 @@ subroutine diag_propack(n,bterm,nroots,factor,maxiter,iverbose,tol,h,e)
       if (iverbose>=2) call TimerStart('diag_propack: diagonalization')
       !
       !
-      !dec$ if (propack_ < 1)
+#if propack_ < 1
         !
         write(out,'("PROPACK is not activated!")')
         stop 'PROPACK is not activated!'
@@ -3302,7 +3302,7 @@ subroutine diag_propack(n,bterm,nroots,factor,maxiter,iverbose,tol,h,e)
         nroots = 0
         return
         !
-      !dec$ endif
+#endif
       !
       nev = nroots
       !
diff --git a/dipole.f90 b/dipole.F90
similarity index 99%
rename from dipole.f90
rename to dipole.F90
index 6d588ee..e858d05 100644
--- a/dipole.f90
+++ b/dipole.F90
@@ -1,10 +1,11 @@
 module dipole
 
-!dec$ define dipole_debug = 0 ! set dipole_debug > 2 with small expansions only
-                              ! 0 - none
-                              ! 1 - some checkings only (no printing)
-                              ! 2 - minimal printing
-                              ! 3 - extendent printing
+#define dipole_debug 0
+! 0 - none
+! 1 - some checkings only (no printing)
+! 2 - minimal printing
+! 3 - extendent printing
+! set dipole_debug > 2 with small expansions only
 
  use accuracy,     only : hik, ik, rk, ark, cl, wl, out, vellgt, planck, avogno, boltz, pi, small_, rad
  use fields,       only : manifold,job,analysis,bset
@@ -848,10 +849,10 @@ subroutine dm_intensity_symmvec(Jval)
       !
     enddo
     !
-    !dec$ if (dipole_debug >= 0)
+#if (dipole_debug >= 0)
        write(out,"(' Total number of lower states = ',i8)") nlevelI
        write(out,"(' Total number of transitions  = ',i10)") Ntransit
-    !dec$ end if
+#endif
     !
     !
     ! In order to speed up the line strength evaluation, 
@@ -3485,10 +3486,10 @@ subroutine dm_intensity(Jval)
     !
     write(out,"('Number of states for each symm = ',<Nrepresen>i8)") nlevelsG(:)
     !
-    !dec$ if (dipole_debug >= 0)
+#if (dipole_debug >= 0)
        write(out,"(' Total number of lower states = ',i8)") nlevelI
        write(out,"(' Total number of transitions  = ',i8)") Ntransit
-    !dec$ end if
+#endif
     !
     !
     ! In order to speed up the line strength evaluation, 
@@ -3676,11 +3677,11 @@ subroutine dm_intensity(Jval)
               !
               !if (ndegI>1) ideg_range = mod(idegI,2)+1
               !
-              !dec$ if (dipole_debug >=3)
+#if (dipole_debug >=3)
               !  !
               !  ideg_range(1) = 1; ideg_range(2) = ndegF 
               !  !
-              !dec$ end if
+#endif
               !
               if (job%rotsym_do) then 
                 !
@@ -3872,11 +3873,11 @@ subroutine dm_intensity(Jval)
                 !!
                 !if (ndegI>1) ideg_range = mod(idegI,2)+1
                 !!
-                !!dec$ if (dipole_debug >=3)
+                !#if (dipole_debug >=3)
                 !  !
                 !  ideg_range(1) = 1; ideg_range(2) = ndegF 
                 !  !
-                !!dec$ end if
+                !#endif
                 !
                 ! swtich to a reduced C3v/D3h case by commenting the next loop
                 !
@@ -5822,4 +5823,4 @@ subroutine convert_symvector_to_contrvector(indI,dimenI,igammaI,idegI,ijterm,sym
      end subroutine convert_symvector_to_contrvector
 
 
-end module dipole
\ No newline at end of file
+end module dipole
diff --git a/lapack.f90 b/lapack.F90
similarity index 98%
rename from lapack.f90
rename to lapack.F90
index e936d36..a9c0453 100644
--- a/lapack.f90
+++ b/lapack.F90
@@ -1,9 +1,9 @@
 module lapack
  
-!dec$ define arpack_ = 0
-!dec$ define blacs_  = 0
-!dec$ define mpi_    = 0
-!dec$ define omparpack_ = 0
+#define arpack_ 0
+#define blacs_  0
+#define mpi_    0
+#define omparpack_ 0
 
 !
 !  Simplistic type-agnostic LAPACK interface
@@ -1299,18 +1299,18 @@ subroutine dseupd_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e,IO)
         nev2 = nev
         rnorm = 1e-5
         !
-        !dec$ if (arpack_ > 0)
+#if arpack_ > 0
         !
         call dsaupd( ido, bmat, n, which, nev, tol, resid, &
                       ncv, v, ldv, iparam, ipntr, workd, workl, &
                       lworkl, info) !,np, rnorm, nconv, nev2 )
         !
-        !dec$ else 
+#else 
             !
             write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
             stop 'Arpack was not activated'
             !
-        !dec$ end if
+#endif
         !
         ido = -1
         !
@@ -1328,17 +1328,17 @@ subroutine dseupd_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e,IO)
         ! has been exceeded.                          
         !
 
-        !dec$ if (arpack_ > 0)
+#if arpack_ > 0
            !
            call dsaupd ( ido, bmat, n, which, nev, tol, resid, &
                          ncv, v, ldv, iparam, ipntr, workd, workl, &
                          lworkl, info) !,np, rnorm, nconv, nev2)
-        !dec$ else 
+#else 
             !
             write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
             stop 'Arpack was not activated'
             !
-        !dec$ end if
+#endif
         !
         !if (verbose>=4.and.iparam(5)>0) then 
         !  !
@@ -1399,7 +1399,7 @@ subroutine dseupd_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e,IO)
       !             
       rvec = .true.
       !
-      !dec$ if (arpack_ > 0)
+#if arpack_ > 0
         !
         if (verbose>=5) write(out,"(/'Arpack: dseupd')") 
         !
@@ -1413,12 +1413,12 @@ subroutine dseupd_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e,IO)
         !
         if (verbose>=5) write(out,"(/'Arpack: done!')") 
         !
-      !dec$ else 
+#else 
           !
           write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
           stop 'Arpack was not activated'
           !
-      !dec$ end if
+#endif
       !
       if ( ierr < 0 ) then
          write(out,"(/'Error with_seupd, info = ',i8)") ierr
@@ -1537,17 +1537,17 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
 	  myid = 1
 	  nprow = 1
       !
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         call BLACS_GRIDINFO( comm, nprow, npcol, myprow, mypcol )
         myid= myprow
-      !dec$ elseif (mpi_ > 0)
+#elif mpi_ > 0
         call MPI_COMM_RANK( comm, myid, ierr )
         call MPI_COMM_SIZE( comm, nprocs_, ierr )
         if (nprocs_/=nprocs) then
           write(out,"('matvec_p: inconsistent number of  nprocs  = ',2i)") nprocs_,nprocs
           stop 'matvec_p: inconsistent number of  nprocs s'
         endif
-      !dec$ end if
+#endif
       !
       kend = kstart(myid) + nloc-1
       !
@@ -1577,12 +1577,12 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
            !
            nx = iend-istart+1
            !
-           !dec$ if (blacs_ > 0)
+#if blacs_ > 0
              call dgerv2d( comm, nx, 1, mv_buf(istart:iend), nx, iprev, mypcol )
-           !dec$ end if
-           !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
              call mpi_recv(mv_buf(istart),nx,MPI_DOUBLE_PRECISION,iprev,myid,comm,ierr)
-           !dec$ end if
+#endif
            !
            istart = max(istart,bterm(k,1))
            !
@@ -1596,12 +1596,12 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
            !
            nx = iend-istart+1
            !
-           !dec$ if (blacs_ > 0)
+#if blacs_ > 0
              call dgesd2d( comm, nx, 1, z(istart:iend), nx, iprev, mypcol )
-           !dec$ end if
-           !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
              call mpi_send(mv_buf(istart),nx,MPI_DOUBLE_PRECISION,iprev,iprev,comm,ierr)
-           !dec$ end if
+#endif
            !
          enddo
          !
@@ -1618,12 +1618,12 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
            !
            nx = iend-istart+1
            !
-           !dec$ if (blacs_ > 0)
+#if blacs_ > 0
              call dgerv2d( comm, nx, 1, mv_buf(istart:iend), nx, inext, mypcol )
-           !dec$ end if
-           !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
              call mpi_recv(mv_buf(istart),nx,MPI_DOUBLE_PRECISION,inext,myid,comm,ierr)
-           !dec$ end if
+#endif
            !
            iend = min(bterm(k,2),iend)
            !
@@ -1637,12 +1637,12 @@ subroutine matvec_p(comm,n,nloc,nprocs,kstart,dx,bterm,h,mv_buf,z,w)
            !
            nx = iend-istart+1
            !
-           !dec$ if (blacs_ > 0)
+#if blacs_ > 0
              call dgesd2d( comm, nx, 1, z(istart:iend), nx, inext, mypcol )
-           !dec$ end if
-           !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
              call mpi_send(mv_buf(istart),nx,MPI_DOUBLE_PRECISION,inext,inext,comm,ierr)
-           !dec$ end if
+#endif
            !
          enddo
          !
@@ -1768,12 +1768,12 @@ subroutine dseupd_p_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
 !     %-----------------------%
 !
 
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         call BLACS_PINFO( iam, nprocs )
         blacs_or_mpi = 'BLACS'
-      !dec$ end if
+#endif
 
-      !dec$ if (mpi_ > 0)
+#if mpi_ > 0
         call MPI_INIT( ierr )
         comm = MPI_COMM_WORLD
         call MPI_COMM_RANK( comm, myid, ierr )
@@ -1788,7 +1788,7 @@ subroutine dseupd_p_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
         !
         blacs_or_mpi = 'MPI'
         !
-      !dec$ end if
+#endif
       !
 
 
@@ -1797,9 +1797,9 @@ subroutine dseupd_p_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
 !
       if (nprocs .lt. 1) then
          nprocs = 1
-         !dec$ if (blacs_ > 0)
+#if blacs_ > 0
            call BLACS_SETUP( iam, nprocs )
-         !dec$ end if
+#endif
       endif
       if (nprocs >maxnprocs) stop 'nprocs > maxnprocs'
       !
@@ -1884,11 +1884,11 @@ subroutine dseupd_p_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
 !
       !
 	  myprow = 1 ; mypcol = 1 ; myid = 1
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         call BLACS_GET( 0, 0, comm )
         call BLACS_GRIDINIT( comm, 'Row', nprow, npcol )
         call BLACS_GRIDINFO( comm, nprow, npcol, myprow, mypcol )
-      !dec$ end if
+#endif
       !
       if (verbose>=2.and.trim(blacs_or_mpi)=='BLACS') write(out,"('myprow, nprow, mypcol, npcol, nprocs = ',5i8)") myprow, nprow, mypcol, npcol, nprocs
       if (verbose>=2.and.trim(blacs_or_mpi)=='MPI') write(out,"('myid,nproc = ',2i8)") myid,nprocs
@@ -1954,19 +1954,19 @@ subroutine dseupd_p_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
         ! has been exceeded.                          
         !
 
-        !dec$ if (blacs_ > 0.or.mpi_ > 0)
+#if (blacs_ > 0) || (mpi_ > 0)
            !
            call pdsaupd ( comm, ido, bmat, nloc, which, nev, tol, resid, &
                           ncv, v, ldv, iparam, ipntr, workd, workl, &
                           lworkl, info )
         !
-        !dec$ else 
+#else 
             !
            call dsaupd ( ido, bmat, n, which, nev, tol, resid, &
                           ncv, v, ldv, iparam, ipntr, workd, workl, &
                           lworkl, info )
             !
-        !dec$ end if
+#endif
         !
         !if (verbose>=4.and.iparam(5)>0) then 
         !  !
@@ -2000,7 +2000,7 @@ subroutine dseupd_p_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
       !             
       rvec = .true.
       !
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         !
         if (verbose>=5) write(out,"(/'Arpack: dseupd')") 
         !
@@ -2012,12 +2012,12 @@ subroutine dseupd_p_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
         if (verbose>=5) write(out,"(/'Arpack: done!')") 
 
         !                    
-      !dec$ else 
+#else 
           !
           write(out,"(/'Arpack was not activated yet. Please uncomment dsaupd and  dseupd bellow')")
           stop 'Arpack was not activated'
           !
-      !dec$ end if
+#endif
       !
       if ( ierr < 0 ) then
          write(out,"(/'Error with_seupd, info = ',i8)") ierr
@@ -2065,13 +2065,13 @@ subroutine dseupd_p_arpack(n,bterm,nroots,factor,maxitr_,tol,h,e)
 !
  9000 continue
 !
-      !dec$ if (blacs_ > 0)
+#if blacs_ > 0
         call BLACS_GRIDEXIT ( comm )
         call BLACS_EXIT(0)
-      !dec$ end if
-      !dec$ if (mpi_ > 0)
+#endif
+#if mpi_ > 0
          call MPI_FINALIZE(rc)
-      !dec$ end if
+#endif
 
 
       deallocate(v,workl,workd,d,resid,select,mv_buf)
diff --git a/makefile b/makefile
index 7fb460c..28ef091 100644
--- a/makefile
+++ b/makefile
@@ -26,6 +26,8 @@ LIB     =   $(LAPACK)
 
 %.o : %.f90
 	$(FOR) -c $(FFLAGS) $<
+%.o : %.F90
+	$(FOR) -c $(FFLAGS) $<
 
 ###############################################################################
 
diff --git a/plasma.f90 b/plasma.F90
similarity index 97%
rename from plasma.f90
rename to plasma.F90
index d311830..cf5615c 100644
--- a/plasma.f90
+++ b/plasma.F90
@@ -1,6 +1,6 @@
 module plasma
- 
-!dec$ define plasma_ = 0
+
+#define plasma_ 0
 
 !
 !  Simplistic type-agnostic PLASMA interface
@@ -21,13 +21,13 @@ module plasma
 
   subroutine plasma_diag_dsytrdx(n,a,w,nroots,Ethres_)
     !
-    !dec$ if (plasma_ > 0) 
+#if plasma_ > 0 
       INCLUDE "plasmaf.h"
       integer(ik), parameter :: VEC = PlasmaVec
       integer(ik), parameter :: UPLO = PlasmaLower
       EXTERNAL PLASMA_DSYTRDX
       INTEGER PLASMA_DSYTRDX
-    !dec$ end if
+#endif
     !
     integer         , intent(in)    :: n
     double precision, intent(inout) :: a(n,n)  ! In:  symmetric matrix to be diagonalized
@@ -69,15 +69,15 @@ subroutine plasma_diag_dsytrdx(n,a,w,nroots,Ethres_)
     !
     if (present(Ethres_)) Ethres = Ethres_
     !
-    !dec$ if (plasma_ == 0) 
+#if plasma_ == 0 
        write(out,"('Plasma is not activated, in plasma.f90 please set plasma_ to 1')")
-    !dec$ end if
+#endif
     !
-    !dec$ if (plasma_ > 0)
+#if plasma_ > 0
       !
       !call getsize(N,LDA,nprocs_)
       !
-    !dec$ end if
+#endif
     !  
     !$omp parallel private(tid)
       if (tid==0) then
@@ -184,7 +184,7 @@ subroutine plasma_diag_dsytrdx(n,a,w,nroots,Ethres_)
     ! 
     ! set up my PLASMA_DSYTRDX environmenta, then call the eigensolver
     !
-    !dec$ if (plasma_ > 0)
+#if plasma_ > 0
         !
         call resetcore(corea,corec)
         CALL PRINTARGS(VEC, UPLO, N, LDA, LDQ, COREA, COREB, COREC, NB, IB)
@@ -202,7 +202,7 @@ subroutine plasma_diag_dsytrdx(n,a,w,nroots,Ethres_)
         END IF
         CALL USETPLASMAENV()
         !
-    !dec$ end if
+#endif
     !
     real_end = get_real_time()
     cpu_end  = get_cpu_time ()
diff --git a/refinement.f90 b/refinement.F90
similarity index 99%
rename from refinement.f90
rename to refinement.F90
index 1f3c4ae..510b2ca 100644
--- a/refinement.f90
+++ b/refinement.F90
@@ -1,6 +1,6 @@
 module refinement
 
-!dec$ define fit_debug = 1
+#define fit_debug 1
 
 
  use accuracy,     only : ik, hik, rk, ark, cl, out, small_
@@ -488,9 +488,9 @@ subroutine sf_fitting(Jval)
        !
        do i=1,pot_npts
          !
-         !dec$ if (fit_debug > 6)
+#if (fit_debug > 6)
            write (out,"('i = ',i)") i
-         !dec$ end if
+#endif
          !
          read (potunit,*) ar_t(1:molec%ncoords),pot_values(i),wtall(en_npts+i)
          local(:,i) = ar_t(:)
@@ -726,9 +726,9 @@ subroutine sf_fitting(Jval)
                        !
                        mat(jentry,ientry) = mat(ientry,jentry)
                        !
-                       !dec$ if (fit_debug > 3)
+#if (fit_debug > 3)
                          write (out,"('mat (',i,',',i,')= ',es14.7)") ientry,jentry,mat(ientry,jentry)
-                       !dec$ end if
+#endif
                        !
                      enddo
                      !
@@ -749,11 +749,11 @@ subroutine sf_fitting(Jval)
                     endif
                  enddo
                  !
-                 !dec$ if (fit_debug > 2)
+#if (fit_debug > 2)
                     !
                     write (out,"(/'Smallest diag. value of mat = ',es14.7,'at k = ',i8,' upper_range = ',es14.7/)") mat(k,k),k,job%upper_ener+mat(k,k)
                     !
-                 !dec$ end if
+#endif
                  !
                  if (allocated(energy_)) deallocate(energy_)
                  !
@@ -1218,18 +1218,18 @@ subroutine sf_fitting(Jval)
                  do icolumn=1,irow    
                    al(irow,icolumn)=sum(rjacob(1:npts,icolumn)*rjacob(1:npts,irow)*wtall(1:npts))
                    al(icolumn,irow)=al(irow,icolumn)
-                   !dec$ if (fit_debug > 2)
+#if (fit_debug > 2)
                      write (out,"('al (',i,',',i,')= ',es14.7)") irow,icolumn,al(irow,icolumn)
-                   !dec$ end if
+#endif
                  enddo
                enddo
                !
                ! form B matrix 
                do irow=1,numpar      
                  bl(irow)=sum(eps(1:npts)*rjacob(1:npts,irow)*wtall(1:npts))
-                 !dec$ if (fit_debug > 2)
+#if (fit_debug > 2)
                    write (out,"('bl (',i,')= ',es14.7)") irow,bl(irow)
-                 !dec$ end if
+#endif
                enddo  
                !
                ! Two types of the linear solver are availible: 
@@ -1438,8 +1438,8 @@ subroutine sf_fitting(Jval)
             !
             ! Print the potential energy points into a separate unit. 
             !
-            !dec$ if (fit_debug > 1)
-            !dec$ end if
+#if (fit_debug > 1)
+#endif
             !
             if (job%verbose>=6) call TimerReport
             !
@@ -2034,9 +2034,9 @@ subroutine bandcentres_fitting(Jval)
                        !
                        mat(jentry,ientry) = mat(ientry,jentry)
                        !
-                       !dec$ if (fit_debug > 3)
+#if (fit_debug > 3)
                          write (out,"('mat (',i,',',i,')= ',es14.7)") ientry,jentry,mat(ientry,jentry)
-                       !dec$ end if
+#endif
                        !
                      enddo
                      !
@@ -2058,11 +2058,11 @@ subroutine bandcentres_fitting(Jval)
                     endif
                  enddo
                  !
-                 !dec$ if (fit_debug > 2)
+#if (fit_debug > 2)
                     !
                     write (out,"(/'Smallest diag. value of mat = ',es14.7,'at k = ',i8,' upper_range = ',es14.7/)") mat(k,k),k,job%upper_ener+mat(k,k)
                     !
-                 !dec$ end if
+#endif
                  !
                  if (allocated(energy_)) deallocate(energy_)
                  !
@@ -2412,18 +2412,18 @@ subroutine bandcentres_fitting(Jval)
                  do icolumn=1,irow    
                    al(irow,icolumn)=sum(rjacob(1:npts,icolumn)*rjacob(1:npts,irow)*wtall(1:npts))
                    al(icolumn,irow)=al(irow,icolumn)
-                   !dec$ if (fit_debug > 2)
+#if (fit_debug > 2)
                      write (out,"('al (',i,',',i,')= ',es14.7)") irow,icolumn,al(irow,icolumn)
-                   !dec$ end if
+#endif
                  enddo
                enddo
                !
                ! form B matrix 
                do irow=1,numpar      
                  bl(irow)=sum(eps(1:npts)*rjacob(1:npts,irow)*wtall(1:npts))
-                 !dec$ if (fit_debug > 2)
+#if (fit_debug > 2)
                    write (out,"('bl (',i,')= ',es14.7)") irow,bl(irow)
-                 !dec$ end if
+#endif
                enddo  
                !
                ! Two types of the linear solver are availible: 
@@ -2667,8 +2667,8 @@ subroutine bandcentres_fitting(Jval)
             !
             ! Print the potential energy points into a separate unit. 
             !
-            !dec$ if (fit_debug > 1)
-            !dec$ end if
+#if (fit_debug > 1)
+#endif
             !
             if (job%verbose>=6) call TimerReport
             !
@@ -2792,9 +2792,9 @@ subroutine restore_vib_matrix_elements
    integer(hik)       :: rootsize,rootsize2
 
 
-   !dec$ if (fit_debug > 1)
+#if (fit_debug > 1)
       write(out, '(/a, 1x, a)') 'read vibrational contracted matrix elements from file', trim(job%extFmat_file)
-   !dec$ end if
+#endif
    !
    job_is ='external field contracted matrix elements for J=0'
    call IOStart(trim(job_is),chkptIO)
@@ -2819,9 +2819,9 @@ subroutine restore_vib_matrix_elements
    rootsize = int(ncontr_t*(ncontr_t+1)/2,hik)
    rootsize2= int(ncontr_t*ncontr_t,hik)
    !
-   !dec$ if (fit_debug > 2)
+#if (fit_debug > 2)
       write(out,"(/'restore_vib_matrix_elements...: Number of elements: ',i8)") ncontr_t
-   !dec$ end if
+#endif
    !
    allocate(poten_me(ncontr_t,ncontr_t,extF%rank),stat=alloc)
    call ArrayStart('poten_me',alloc,1,kind(poten_me),rootsize2)
@@ -2847,9 +2847,9 @@ subroutine restore_vib_matrix_elements
    !
    close(chkptIO,status='keep')
 
-   !dec$ if (fit_debug > 1)
+#if (fit_debug > 1)
       write(out, '(/a)') 'done'
-   !dec$ end if
+#endif
 
    !
  end subroutine restore_vib_matrix_elements
@@ -2922,9 +2922,9 @@ subroutine Hamiltonian_vector(ientry,indJ,nentries,tmat,cdimen,poten,mat)
     end do
     !$omp end parallel do
     !
-    !dec$ if (fit_debug > 3)
+#if (fit_debug > 3)
       write (out,"('ientry = ',i,'; dimen = ',i)") ientry,dimen
-    !dec$ end if
+#endif
     !
     !loop over final states
     !
diff --git a/tran.f90 b/tran.F90
similarity index 98%
rename from tran.f90
rename to tran.F90
index 9c233d8..5eb9723 100644
--- a/tran.f90
+++ b/tran.F90
@@ -4,7 +4,8 @@
 !
 module tran
 
-!dec$ define tran_debug = 1 ! set tran_debug > 2 with small vibrational bases and small expansions only
+#define tran_debug 1
+! set tran_debug > 2 with small vibrational bases and small expansions only
 
  use accuracy,     only : ik, rk, hik, ark, cl, out, small_
  use timer,        only : IOstart,IOstop,arraystart,arraystop,arrayminus,Timerstart,Timerstop,TimerReport,MemoryReport
@@ -95,10 +96,10 @@ subroutine read_contrind(njval, jval)
        filename = trim(job%eigenfile%filebase)//'_quanta'//trim(adjustl(jchar))//'.chk'
 
 
-       !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
           write(out, '(/a, 1x, i2, 2(1x, a))') 'read contraction indexes for J =', jval(jind), &
                                                'from file', trim(filename)
-       !dec$ end if
+#endif
 
 
        write(ioname, '(a, i4)') 'contraction indexes for J=', jval(jind)
@@ -120,11 +121,11 @@ subroutine read_contrind(njval, jval)
        bset_contr(jind)%nclasses     = nclasses
 
 
-       !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
           write(out, '(3(/1x, a, 1x, i5))') 'ncases ', ncases, 'ndegmax', nlambdas, 'ncontr ', ncontr
           write(out, '(/1x, a, 1x, a, 4x, a/40x, 1x, a, 1x, a, 6x, a)') 'icase', 'ndeg', 'ilevel(0:nclasses)', 'ideg', &
                                                                         'iroot', 'ideg(0:nclasses)'
-       !dec$ end if
+#endif
        !
        allocate(bset_contr(jind)%index_deg(ncases),bset_contr(jind)%contractive_space(0:nclasses, ncases),stat = info)
        call ArrayStart('bset_contr',info,size(bset_contr(jind)%contractive_space),kind(bset_contr(jind)%contractive_space))
@@ -148,10 +149,10 @@ subroutine read_contrind(njval, jval)
           !
           bset_contr(jind)%index_deg(icase)%size1 = nlambdas
           !
-          !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
              write(out, '(1x, i5, 2x, i3, 1x, <nclasses>(1x, i3), 1x, i3)')                                            &
              icase, nlambdas, bset_contr(jind)%contractive_space(0:nclasses, icase)
-          !dec$ end if
+#endif
           !
           allocate(bset_contr(jind)%index_deg(icase)%icoeffs(0:nclasses, nlambdas))
           call ArrayStart('bset_contr',info,size(bset_contr(jind)%index_deg(icase)%icoeffs),kind(bset_contr(jind)%index_deg(icase)%icoeffs))
@@ -161,10 +162,10 @@ subroutine read_contrind(njval, jval)
              read(iounit, '(i8, i6, <nclasses>i6)') iroot, bset_contr(jind)%index_deg(icase)%icoeffs(0:nclasses, ilambda)
 
 
-             !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
                 write(out, '(42x, i3, 1x, i5, 3x, <nclasses>(1x, i3), 1x, i3)')                                        &
                 ilambda, iroot, bset_contr(jind)%index_deg(icase)%icoeffs(0:nclasses, ilambda)
-             !dec$ end if
+#endif
 
 
              icontr = icontr + 1
@@ -188,10 +189,10 @@ subroutine read_contrind(njval, jval)
        read(iounit, '(2i8)') ncases, nlambdas
 
 
-       !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
           write(out, '(2(/1x, a, 1x, i2))') 'nrot      ', ncases, 'nrotdegmax', nlambdas
           write(out, '(/1x, a, 1x, a, 8x, a, 3x, a, 1x, a)') 'irot', 'irotdeg', 'J', 'K', 'Tau'
-       !dec$ end if
+#endif
 
 
        allocate(bset_contr(jind)%rot_index(ncases, nlambdas), stat = info)
@@ -205,11 +206,11 @@ subroutine read_contrind(njval, jval)
                                        bset_contr(jind)%rot_index(icase, ilambda)%k,                                       &
                                        bset_contr(jind)%rot_index(icase, ilambda)%tau
 
-          !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
              write(out, '(2x, i3, 5x, i3, 5x, 3(1x, i3))') icase, ilambda, bset_contr(jind)%rot_index(icase, ilambda)%j,   &
                                                                            bset_contr(jind)%rot_index(icase, ilambda)%k,   &
                                                                            bset_contr(jind)%rot_index(icase, ilambda)%tau
-          !dec$ end if
+#endif
 
        end do
 
@@ -307,19 +308,19 @@ subroutine index_correlation(njval, jval)
        allocate(cnu_i(1:nclasses),cnu_j(1:nclasses),stat = info)
        if (info /= 0) stop 'index_correlation: cnu_i allocation error'
        !
-       !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
           write(out, '(/a, 1x, i2)') 'find correlation between contraction indexes for J = 0 and J =', jval(jind)
-       !dec$ end if
+#endif
        !
        if (jind > size(bset_contr)) stop 'index_correlation error: jind > size(bset_contr)'
        !
        allocate(bset_contr(jind)%icontr_correlat_j0(bset_contr(jind)%Maxsymcoeffs,bset_contr(jind)%max_deg_size), stat = info)
        call ArrayStart('bset_contr',info,size(bset_contr(jind)%icontr_correlat_j0),kind(bset_contr(jind)%icontr_correlat_j0))
        !
-       !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
           write(out, '(/7x, a, i3, 17x, a/1x, a, 1x, a, 1x, a, 4x, a, 1x, a, 1x, a)')                                  &
           'J =', jval(jind), 'J =  0', 'icase', 'ilambda', 'iroot', 'icase', 'ilambda', 'iroot'
-       !dec$ end if
+#endif
        !
        l_icase : do icase = 1, bset_contr(jind)%Maxsymcoeffs
           !
@@ -357,13 +358,13 @@ subroutine index_correlation(njval, jval)
              !
              bset_contr(jind)%icontr_correlat_j0(icase, ilambda) = jcontr
              !
-             !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
                 write(out, '(i6, 6x, i2, i6, 3x, i6, 6x, i2, i6)')   &
                 !                      case,                      lambda, contr
                 icase,ilambda,bset_contr(jind)%icase2icontr(icase,ilambda), &
                 bset_contr(1)%icontr2icase(jcontr,1), bset_contr(1)%icontr2icase(jcontr,2), jcontr
                 !
-             !dec$ end if
+#endif
              !
           end do l_ilambda
           !
@@ -401,15 +402,15 @@ subroutine index_correlation(njval, jval)
             bset_contr(jind)%k(iroot)    = ideg
           endif 
           !
-          !dec$ if (tran_debug >= 3)
+#if (tran_debug >= 3)
             write (out,"('iroot,icase,ilambda,jcontr = ',4i7)") iroot,icase,ilambda,jcontr 
-          !dec$ end if
+#endif
           !
        enddo
        !
-       !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
           write(out, '(/a)') 'done'
-       !dec$ end if
+#endif
        !
        deallocate(cnu_i,cnu_j)
        !
@@ -518,9 +519,9 @@ subroutine read_eigenval(njval, jval)
           !
           filename = trim(job%eigenfile%filebase)//'_descr'//trim(adjustl(jchar))//'_'//trim(adjustl(gchar))//'.chk'
           !
-          !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
              write(out, '(/a, 1x, i2, 2(1x, a))') 'read eigenvalues for J =', jval(jind), ', gamma = ',i2,' from file', trim(filename)
-          !dec$ end if
+#endif
           !
           if (jind > size(bset_contr)) stop 'read_eigenval error: jind > size(bset_contr)'
           !
@@ -838,7 +839,7 @@ subroutine read_eigenval(njval, jval)
           !
           !
           !print energies
-          !dec$ if (tran_debug > 2)
+#if (tran_debug > 2)
              write(out, '(/1x, a, 2x, i8/1x, a, 1x, i8)') 'number of roots', nroots, 'number of levels', nlevels
              write(out, '(/1x, a, 11x, a, 1x, a, 1x, a, 8x, a, <nmodes>(2x), 1x, a)') 'ilevel', 'energy', 'ndeg',       &
              'igamma', 'nu(0:nmodes)', 'irec'
@@ -853,7 +854,7 @@ subroutine read_eigenval(njval, jval)
                 !
              end do
              write(out, '(a)') '...done!'
-          !dec$ end if
+#endif
           !
        enddo
        !