diff --git a/examples/MHD/input_mhdchan.i3d b/examples/MHD/input_mhdchan.i3d
index aaa20393e..82431ca6f 100644
--- a/examples/MHD/input_mhdchan.i3d
+++ b/examples/MHD/input_mhdchan.i3d
@@ -121,8 +121,8 @@ nclzSn = 0
 &MHDParam
 !====================
 
-mhd_equation = 'potential'
-hartmann = 6.d0
+mhd_equation = 'potential' ! MHD type of equation (potential or induction)
+hartmann = 6.d0            ! Hartmann number (Ha) ratio of electromagnetic force to the viscous force (BLsqrt(sigma/nu)) 
 
 /End
 
diff --git a/examples/MHD/input_mhdchan_test.i3d b/examples/MHD/input_mhdchan_test.i3d
index e400b1f7b..999f9c8fa 100644
--- a/examples/MHD/input_mhdchan_test.i3d
+++ b/examples/MHD/input_mhdchan_test.i3d
@@ -121,8 +121,8 @@ nclzSn = 0
 &MHDParam
 !====================
 
-mhd_equation = 'potential'
-hartmann = 6.d0
+mhd_equation = 'potential' ! MHD type of equation (potential or induction)
+hartmann = 6.d0            ! Hartmann number (Ha) ratio of electromagnetic force to the viscous force (BLsqrt(sigma/nu)) 
 
 /End
 
diff --git a/examples/MHD/input_otv.i3d b/examples/MHD/input_otv.i3d
index 5d475869a..74ffb77e8 100644
--- a/examples/MHD/input_otv.i3d
+++ b/examples/MHD/input_otv.i3d
@@ -8,8 +8,8 @@
 itype = 2
 
 ! Domain decomposition
-p_row=2                ! Row partition
-p_col=2                ! Column partition
+p_row=0                ! Row partition
+p_col=0                ! Column partition
  
 ! Mesh 
 nx=128                 ! X-direction nodes
@@ -97,9 +97,9 @@ nstat = 1             ! Size arrays for statistic collection
 &MHDParam
 !====================
 
-mhd_equation = 'induction'
-Stuart = 17.54386 
-Rem    = 17.54386
+mhd_equation = 'induction' ! MHD type of equation (potential or induction)
+Stuart = 17.54386          ! Stuart number ratio of electromagnetic to inertial forces (Ha^2/Rem)
+Rem    = 17.54386          ! Magnetic Reynolds number (UL/eta)
 
 /End
 
diff --git a/examples/MHD/input_otv_test.i3d b/examples/MHD/input_otv_test.i3d
index e2d327acc..59c1276b6 100644
--- a/examples/MHD/input_otv_test.i3d
+++ b/examples/MHD/input_otv_test.i3d
@@ -97,9 +97,9 @@ nstat = 1             ! Size arrays for statistic collection
 &MHDParam
 !====================
 
-mhd_equation = 'induction'
-Stuart = 50.d0
-Rem    = 50.d0
+mhd_equation = 'induction' ! MHD type of equation (potential or induction)
+Stuart = 50.0              ! Stuart number ratio of electromagnetic to inertial forces (Ha^2/Rem)
+Rem    = 50.0              ! Magnetic Reynolds number (UL/eta)
 
 /End
 
diff --git a/src/Case-ABL.f90 b/src/Case-ABL.f90
index 9e729d6d0..f59b1cccc 100644
--- a/src/Case-ABL.f90
+++ b/src/Case-ABL.f90
@@ -119,21 +119,6 @@ subroutine init_abl(ux1,uy1,uz1,ep1,phi1)
           endif
         endif
       enddo
-
-      !! Add random noise
-      !do j=1,xsize(2)
-      !  if (istret==0) y=real(j + xstart(2)-1-1,mytype)*dy
-      !  if (istret/=0) y=yp(j+xstart(2)-1)
-      !  !if (y.lt.50) then 
-      !  !  do k=1,xsize(3)
-      !  !  do i=1,xsize(1)
-      !  !    call random_number(phinoise)
-      !  !    phinoise=0.1*(phinoise*2.-1.)
-      !  !    phi1(i,j,k,1)=phi1(i,j,k,1)+phinoise
-      !  !  enddo
-      !  !  enddo
-      !  !endif
-      !enddo
     endif
 
     return
diff --git a/src/Case-Channel.f90 b/src/Case-Channel.f90
index 9d8d30482..078c476fd 100644
--- a/src/Case-Channel.f90
+++ b/src/Case-Channel.f90
@@ -29,7 +29,7 @@ subroutine init_channel (ux1,uy1,uz1,ep1,phi1)
     use variables
     use param
     use MPI
-    use mhd, only : mhd_active, mhd_equation,Bm,Bmean
+    use mhd, only : mhd_equation,Bm,Bmean
 
     implicit none
 
@@ -93,9 +93,9 @@ subroutine init_channel (ux1,uy1,uz1,ep1,phi1)
                 if (idir_stream == 1) then
                    ux1(i,j,k)=one-y*y
                    uy1(i,j,k)=zero
-                   uz1(i,j,k)=zero !sin(real(i-1,mytype)*dx)+cos(real(k-1,mytype)*dz)
+                   uz1(i,j,k)=sin(real(i-1,mytype)*dx)+cos(real(k+xstart(3)-2,mytype)*dz)
                 else
-                        print *,'test'
+                   ! TODO: check if ux1 needs sin and cos 
                    uz1(i,j,k)=one-y*y
                    uy1(i,j,k)=zero
                    ux1(i,j,k)=zero
@@ -155,7 +155,7 @@ subroutine boundary_conditions_channel (ux,uy,uz,phi)
     use param
     use var, only : di2
     use variables
-    use mhd, only : Bm, mhd_active, mhd_equation 
+    use mhd, only : Bm, mhd_equation 
 
     implicit none
 
@@ -276,100 +276,13 @@ subroutine postprocess_channel(ux1,uy1,uz1,pp3,phi1,ep1)
     real(mytype), intent(in), dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi1
     real(mytype), intent(in), dimension(ph1%zst(1):ph1%zen(1),ph1%zst(2):ph1%zen(2),nzmsize,npress) :: pp3
 
-    real(mytype) :: ufriction,umax,resdux1,resduy1,resduz1,temp1,temp2
-    integer :: nxc, nyc, nzc, xsize1, xsize2, xsize3
-    integer :: i,j,k
-    integer :: code
-
-    if (nclx1==1 .and. xend(1)==nx) then
-       xsize1=xsize(1)-1
-    else
-       xsize1=xsize(1)
-    endif
-    if (ncly1==1 .and. xend(2)==ny) then
-       xsize2=xsize(2)-1
-    else
-       xsize2=xsize(2)
-    endif
-    if (nclz1==1 .and. xend(3)==nz) then
-       xsize3=xsize(3)-1
-    else
-       xsize3=xsize(3)
-    endif
-    if (nclx1==1) then
-       nxc=nxm
-    else
-       nxc=nx
-    endif
-    if (ncly1==1) then
-       nyc=nym
-    else
-       nyc=ny
-    endif
-    if (nclz1==1) then
-       nzc=nzm
-    else
-       nzc=nz
-    endif
-
-    !we only collect statistics every 10 time steps to save computational time
-    if (mod(itime, 10) == 0) then
-
-       temp1=zero
-
-       if (ncly1==2 .and. xstart(2)==1) then
-
-         ! bottom wall, one node away from wall
-         j=2
-         do k=1,xsize(3)
-             do i=1,xsize(1)
-              temp1=temp1+ux1(i,j,k)*ppy(j-1)
-             enddo
-         enddo
-
-       endif
-       !
-       if (nclyn==2 .and. xend(2)==ny) then
-
-         ! upper wall, one node away from wall
-         j=xsize(2)-1
-         do k=1,xsize(3)
-             do i=1,xsize(1)
-              temp1=temp1-ux1(i,j,k)*ppy(j+1)
-             enddo
-         enddo
-
-       endif
-       !
-       temp2=zero
-       do k=1,xsize3
-          do j=1,xsize2
-             do i=1,xsize1
-              temp2=max(temp2,ux1(i,j,k))
-             enddo
-          enddo
-       enddo
-       
-       call MPI_ALLREDUCE(temp1,ufriction,1,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-       call MPI_ALLREDUCE(temp2,     umax,1,real_type,MPI_MAX,MPI_COMM_WORLD,code)
-
-       ufriction=xnu*1.5_mytype*ufriction/(nxc*nzc) ! ub=2/3
-   
-       if (nrank==0) then
-          write(52,'(3(E20.12))') (itime-1)*dt,ufriction,umax
-          flush(52)
-       endif
-
-    endif
-
-
   end subroutine postprocess_channel
 
   subroutine visu_channel_init(visu_initialised)
 
     use decomp_2d_io, only : decomp_2d_register_variable
     use visu, only : io_name, output2D
-    use mhd, only : mhd_active
+    use param, only : mhd_active
     
     implicit none
 
@@ -403,7 +316,8 @@ subroutine visu_channel(ux1, uy1, uz1, pp3, phi1, ep1, num)
     use var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3
     use var, ONLY : nzmsize
     use visu, only : write_field
-    use mhd, only : mhd_active,mhd_equation,Je, Bm
+    use param, only : mhd_active
+    use mhd, only : mhd_equation,Je, Bm
     
     use ibm_param, only : ubcx,ubcy,ubcz
 
diff --git a/src/Case-TGV.f90 b/src/Case-TGV.f90
index bff26aa8c..f55014377 100644
--- a/src/Case-TGV.f90
+++ b/src/Case-TGV.f90
@@ -28,7 +28,7 @@ subroutine init_tgv (ux1,uy1,uz1,ep1,phi1)
     use variables
     use param
     use MPI
-    use mhd,  only: mhd_active,Bm,Bmean
+    use mhd,  only: Bm,Bmean
     use visu, only: write_snapshot, end_snapshot
 
     implicit none
@@ -196,7 +196,7 @@ subroutine postprocess_tgv(ux1,uy1,uz1,phi1,ep1)
     USE var, only : ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1,di1
     USE var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3
     USE ibm_param
-    use mhd, only : mhd_active,Bm,Je,Rem
+    use mhd, only : Bm,Je,Rem
 
     real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3)) :: ux1, uy1, uz1,ep1
     real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi1
@@ -485,7 +485,6 @@ subroutine visu_tgv_init (visu_initialised)
 
     use decomp_2d_io, only : decomp_2d_register_variable
     use visu, only : io_name, output2D
-    use mhd, only : mhd_active
     
     implicit none
 
@@ -521,7 +520,7 @@ subroutine visu_tgv(ux1, uy1, uz1, num)
     USE var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3
     use var, ONLY : nxmsize, nymsize, nzmsize
     use visu, only : write_field
-    use mhd, only : mhd_active,Bm,Je,del_cross_prod,rem
+    use mhd, only : Bm,Je,del_cross_prod,rem
     use ibm_param, only : ubcx,ubcy,ubcz
 
     implicit none
diff --git a/src/case.f90 b/src/case.f90
index 9118b9d03..245ccb268 100644
--- a/src/case.f90
+++ b/src/case.f90
@@ -493,7 +493,7 @@ end subroutine visu_case
   !##################################################################
   subroutine momentum_forcing(dux1, duy1, duz1, rho1, ux1, uy1, uz1, phi1)
 
-    use mhd, only: mhd_active,momentum_forcing_mhd
+    use mhd, only: momentum_forcing_mhd
 
     implicit none
 
diff --git a/src/mhd.f90 b/src/mhd.f90
index 824de83d0..0bdc47dc0 100644
--- a/src/mhd.f90
+++ b/src/mhd.f90
@@ -11,11 +11,9 @@ module mhd
   !
   implicit none
   !
-  logical :: mhd_active
   character(len=9) :: mhd_equation
   real(8) :: hartmann,stuart,rem
   !+------------+------------------------------------------------------+
-  !|  mhd_active| the swith to activate the mhd module.                |
   !|    hartmann| hartmann number, the ratio of lorentz force to       |
   !|            | viscous force                                        |
   !|      stuart| Stuart number, magnetic interaction parameter, ratio |
@@ -48,7 +46,9 @@ subroutine mhd_init
     !
     use param, only: re,ntime,itype,itype_channel,itype_tgv,zlz,dz
     use param, only: zero, one
+    use decomp_2d_mpi, only : decomp_2d_abort
     !
+    ! TODO fix these if
     ! stuart=hartmann**2/re
     if(stuart<=1.d-15) then
       stuart=hartmann**2/re
@@ -78,16 +78,15 @@ subroutine mhd_init
     !
     allocate(Bmean(xsize(1),xsize(2),xsize(3),1:3))
     !
-    if(nrank==0) print*,'** MHD fields allocated'
-    !
     if(mhd_equation == 'induction') then
-      
+      ! TODO Check what is going to happen for other flow types  
       if(itype.eq.itype_tgv) then
         Bmean(:,:,:,1)=zero
         Bmean(:,:,:,2)=zero
         Bmean(:,:,:,3)=zero
       endif
     elseif(mhd_equation == 'potential') then
+      ! TODO Check what is going to happen for other flow types  
       if(itype.eq.itype_channel) then
         Bmean=zero
         !
@@ -95,10 +94,8 @@ subroutine mhd_init
         Bm(:,:,:,2)=one
         Bm(:,:,:,3)=zero
       endif
-      !
-      if(nrank==0) print*,'** MHD fields initlised'
     else
-      stop ' mhd_equation not induction or potential'
+      call decomp_2d_abort(1, "Error in setup mhd equation induction or potential not selected")
     endif
     !
   end subroutine mhd_init
@@ -183,6 +180,10 @@ end function vortcal
   !+-------------------------------------------------------------------+
   subroutine momentum_forcing_mhd(dux1,duy1,duz1,ux1,uy1,uz1)
     !
+    USE decomp_2d_mpi, only : decomp_2d_abort
+    
+    implicit none
+    
     ! arguments
     real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3)) ::   &
                                                 ux1,uy1,uz1
@@ -198,7 +199,7 @@ subroutine momentum_forcing_mhd(dux1,duy1,duz1,ux1,uy1,uz1)
     elseif(mhd_equation == 'potential') then
       Je=solve_mhd_potential_poisson(ux1,uy1,uz1)
     else
-      stop ' mhd_equation error '
+      call decomp_2d_abort(1, "Error in setup mhd equation induction or potential not selected")
     endif
     !
     do k = 1, xsize(3)
@@ -431,7 +432,8 @@ subroutine mhd_rhs_eq(dB,B,B0,ux1,uy1,uz1)
     integer :: i,j,k,is
     !
 
-    ! local 
+    ! local
+    ! TODO Need to check is all these local variable are necessary -> potential issue with INTEL  
     real(mytype) :: rrem
     real(mytype), save, allocatable, dimension(:,:,:) :: tx1,ty1,tz1,tx2,ty2,tz2, &
                                                          tx3,ty3,tz3,bx2,by2,bz2, &
@@ -696,6 +698,8 @@ subroutine mhd_rhs_eq(dB,B,B0,ux1,uy1,uz1)
 
   end subroutine mhd_rhs_eq
   !
+  ! TODO Check if the normal divergence can be used  
+  ! TODO Check if already allocated arrays can be re-used
   subroutine solve_poisson_mhd
     !
     use decomp_2d, only : zsize, ph1
@@ -749,6 +753,7 @@ function divergence_scalar(vec,nlock) result(pp3)
 
     nvect3=(ph1%zen(1)-ph1%zst(1)+1)*(ph1%zen(2)-ph1%zst(2)+1)*nzmsize
 
+    ! TODO If we properly use tmp var this should not be used 
     ta1(:,:,:) = vec(:,:,:,1)
     tb1(:,:,:) = vec(:,:,:,2)
     tc1(:,:,:) = vec(:,:,:,3)
diff --git a/src/module_param.f90 b/src/module_param.f90
index db5085da4..a3f629b37 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -327,6 +327,7 @@ module param
   real(mytype) :: C_filter
   character(len=512) :: inflowpath
   logical :: validation_restart
+  logical :: mhd_active
 
   ! Logical, true when synchronization is needed
   logical, save :: sync_vel_needed = .true.
diff --git a/src/navier.f90 b/src/navier.f90
index 76c587c52..2b54d467e 100644
--- a/src/navier.f90
+++ b/src/navier.f90
@@ -623,13 +623,23 @@ subroutine pre_correc(ux,uy,uz,ep)
                 dpdzy1(i,k)=dpdzy1(i,k)*gdt(itr)
              enddo
           enddo
-          do k=1,xsize(3)
-             do i=1,xsize(1)
-                ux(i,1,k)=byx1(i,k) !+dpdxy1(i,k)
-                uy(i,1,k)=byy1(i,k)
-                uz(i,1,k)=byz1(i,k) !+dpdzy1(i,k)
+          if (mhd_active) then
+             do k=1,xsize(3)
+                do i=1,xsize(1)
+                   ux(i,1,k)=byx1(i,k) 
+                   uy(i,1,k)=byy1(i,k)
+                   uz(i,1,k)=byz1(i,k) 
+                enddo
              enddo
-          enddo
+          else
+             do k=1,xsize(3)
+                do i=1,xsize(1)
+                   ux(i,1,k)=byx1(i,k) + dpdxy1(i,k)
+                   uy(i,1,k)=byy1(i,k)
+                   uz(i,1,k)=byz1(i,k) + dpdzy1(i,k)
+                enddo
+             enddo
+          endif
        endif
     endif
 
@@ -643,21 +653,41 @@ subroutine pre_correc(ux,uy,uz,ep)
           enddo
        endif
        if (dims(1)==1) then
-          do k=1,xsize(3)
-             do i=1,xsize(1)
-                ux(i,xsize(2),k)=byxn(i,k) !+dpdxyn(i,k)
-                uy(i,xsize(2),k)=byyn(i,k)
-                uz(i,xsize(2),k)=byzn(i,k) !+dpdzyn(i,k)
+          if (mhd_active) then
+             do k=1,xsize(3)
+                do i=1,xsize(1)
+                   ux(i,xsize(2),k)=byxn(i,k) 
+                   uy(i,xsize(2),k)=byyn(i,k)
+                   uz(i,xsize(2),k)=byzn(i,k) 
+                enddo
              enddo
-          enddo
+          else
+             do k=1,xsize(3)
+                do i=1,xsize(1)
+                   ux(i,xsize(2),k)=byxn(i,k) + dpdxyn(i,k)
+                   uy(i,xsize(2),k)=byyn(i,k)
+                   uz(i,xsize(2),k)=byzn(i,k) + dpdzyn(i,k)
+                enddo
+             enddo
+          endif
        elseif (ny - (nym / dims(1)) == xstart(2)) then
-          do k=1,xsize(3)
-             do i=1,xsize(1)
-                ux(i,xsize(2),k)=byxn(i,k) !+dpdxyn(i,k)
-                uy(i,xsize(2),k)=byyn(i,k)
-                uz(i,xsize(2),k)=byzn(i,k) !+dpdzyn(i,k)
+          if (mhd_active) then
+             do k=1,xsize(3)
+                do i=1,xsize(1)
+                   ux(i,xsize(2),k)=byxn(i,k) 
+                   uy(i,xsize(2),k)=byyn(i,k)
+                   uz(i,xsize(2),k)=byzn(i,k) 
+                enddo
              enddo
-          enddo
+          else
+             do k=1,xsize(3)
+                do i=1,xsize(1)
+                   ux(i,xsize(2),k)=byxn(i,k) + dpdxyn(i,k)
+                   uy(i,xsize(2),k)=byyn(i,k)
+                   uz(i,xsize(2),k)=byzn(i,k) + dpdzyn(i,k)
+                enddo
+             enddo
+          endif
        endif
     endif
 
diff --git a/src/parameters.f90 b/src/parameters.f90
index 44c28e434..1aa2614fc 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -31,7 +31,7 @@ subroutine parameter(input_i3d)
   use visu, only : output2D
   use forces, only : iforces, nvol, setup_forces
 
-  use mhd, only: mhd_active,mhd_equation,hartmann,stuart,rem
+  use mhd, only: mhd_equation,hartmann,stuart,rem
 
   implicit none
 
@@ -599,7 +599,7 @@ subroutine parameter_defaults()
   use visu, only : output2D
   use forces, only : iforces, nvol
 
-  use mhd, only: mhd_active, mhd_equation, rem, stuart, hartmann 
+  use mhd, only: mhd_equation, rem, stuart, hartmann 
 
   implicit none
 
diff --git a/src/schemes.f90 b/src/schemes.f90
index b09754277..f8145bcfc 100644
--- a/src/schemes.f90
+++ b/src/schemes.f90
@@ -16,7 +16,7 @@ subroutine schemes()
   USE variables
   USE var
   USE ydiff_implicit, only : init_implicit, implicit_schemes
-  use mhd, only: mhd_active, mhd_equation
+  use mhd, only: mhd_equation
 
   implicit none
 
diff --git a/src/time_integrators.f90 b/src/time_integrators.f90
index 26f390051..eafdf74e1 100644
--- a/src/time_integrators.f90
+++ b/src/time_integrators.f90
@@ -203,9 +203,10 @@ SUBROUTINE int_time(rho1, ux1, uy1, uz1, phi1, drho1, dux1, duy1, duz1, dphi1)
     use param, only : iimplicit, sc_even
     use param, only : primary_species, massfrac
     use param, only : scalar_lbound, scalar_ubound
+    use param, only : mhd_active
     use variables, only : numscalar,nu0nu
     use var, only : ta1, tb1
-    use mhd,     only : mhd_active,mhd_equation,int_time_magnet
+    use mhd,     only : mhd_equation,int_time_magnet
     use MPI
 
 
diff --git a/src/tools.f90 b/src/tools.f90
index 702b19c6d..0e9dd9c09 100644
--- a/src/tools.f90
+++ b/src/tools.f90
@@ -209,7 +209,7 @@ subroutine restart(ux1,uy1,uz1,dux1,duy1,duz1,ep1,pp3,phi1,dphi1,px1,py1,pz1,rho
     use param
     use MPI
     use navier, only : gradp
-    use mhd, only : mhd_active,mhd_equation,Bm,dBm
+    use mhd, only : mhd_equation,Bm,dBm
 
     implicit none
 
@@ -507,9 +507,9 @@ subroutine init_restart_adios2()
 
     use decomp_2d_io, only : decomp_2d_register_variable, decomp_2d_init_io
     use variables, only : numscalar
-    use param, only : ilmn, nrhotime, ntime
+    use param, only : ilmn, nrhotime, ntime, mhd_active
     use var, only : itimescheme, iibm
-    use mhd, only : mhd_active, mhd_equation
+    use mhd, only : mhd_equation
     
     implicit none
 
@@ -793,8 +793,9 @@ end subroutine write_outflow
   subroutine compute_cfldiff()
      use param, only : xnu,dt,dx,dy,dz,istret
      use param, only : cfl_diff_sum, cfl_diff_x, cfl_diff_y, cfl_diff_z
+     use param, only : mhd_active
      use variables, only : dyp
-     use mhd, only: mhd_active, mhd_equation,rem
+     use mhd, only: mhd_equation,rem
 
      implicit none
 
diff --git a/src/transeq.f90 b/src/transeq.f90
index 14686f659..26340950d 100644
--- a/src/transeq.f90
+++ b/src/transeq.f90
@@ -21,8 +21,8 @@ subroutine calculate_transeq_rhs(drho1,dux1,duy1,duz1,dphi1,rho1,ux1,uy1,uz1,ep1
 
     use decomp_2d, only : xsize, zsize
     use variables, only : numscalar
-    use param, only : ntime, ilmn, nrhotime, ilmn_solve_temp
-    use mhd,   only : mhd_active,mhd_equation,calculate_mhd_transeq_rhs
+    use param, only : ntime, ilmn, nrhotime, ilmn_solve_temp, mhd_active
+    use mhd,   only : mhd_equation,calculate_mhd_transeq_rhs
 
     implicit none
 
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index d02ac04d3..5ecbcef77 100644
--- a/src/xcompact3d.f90
+++ b/src/xcompact3d.f90
@@ -16,8 +16,9 @@ program xcompact3d
   use ibm_param
   use ibm, only : body
   use genepsi, only : genepsi3d
-  use mhd,    only : Bm,mhd_active,mhd_equation,test_magnetic, &
+  use mhd,    only : Bm,mhd_equation,test_magnetic, &
                      solve_poisson_mhd
+  use param, only : mhd_active
 
   implicit none
 
@@ -119,7 +120,7 @@ subroutine init_xcompact3d()
        init_inflow_outflow, read_inflow
 
   use param, only : ilesmod, jles,itype
-  use param, only : irestart
+  use param, only : irestart, mhd_active
 
   use variables, only : nx, ny, nz, nxm, nym, nzm
   use variables, only : p_row, p_col
@@ -135,7 +136,7 @@ subroutine init_xcompact3d()
 
   use probes, only : init_probes
 
-  use mhd, only: mhd_active,mhd_init
+  use mhd, only: mhd_init
 
   implicit none
 
@@ -278,17 +279,9 @@ subroutine init_xcompact3d()
      endif
   endif
 
-  if (itype==3) then
-     if (nrank.eq.0)then
-        open(52,file='chan_time_evol.dat',form='formatted')
-     endif
-  endif
-  
-  if (itype==5) then
-     if (iforces == 1) then
-        if(nrank.eq.0)then
-           open(38,file='forces.dat',form='formatted')
-        endif
+  if (iforces == 1) then
+     if(nrank.eq.0)then
+        open(38,file='forces.dat',form='formatted')
      endif
   endif
   
@@ -309,11 +302,11 @@ subroutine finalise_xcompact3d()
   use decomp_2d_io, only : decomp_2d_io_finalise
 
   use tools, only : simu_stats
-  use param, only : itype, jles, ilesmod
+  use param, only : itype, jles, ilesmod, mhd_active
   use probes, only : finalize_probes
   use visu, only : visu_finalise
   use les, only: finalise_explicit_les
-  use mhd, only: mhd_active, mhd_fin
+  use mhd, only: mhd_fin
   use case, only: visu_case_finalise
   use forces, only: iforces
 
@@ -330,17 +323,9 @@ subroutine finalise_xcompact3d()
      endif
   endif
 
-  if (itype==3) then
+  if (iforces == 1) then
      if(nrank.eq.0)then
-        close(52)
-     endif
-  endif
-
-  if (itype==5) then
-     if (iforces == 1) then
-        if(nrank.eq.0)then
-           close(38)
-        endif
+        close(38)
      endif
   endif