diff --git a/Source/ccib.f90 b/Source/ccib.f90
index be466134231..eafbe5b069e 100644
--- a/Source/ccib.f90
+++ b/Source/ccib.f90
@@ -9,7 +9,6 @@ MODULE CC_SCALARS
 USE COMPLEX_GEOMETRY
 USE PRECISION_PARAMETERS
 USE GLOBAL_CONSTANTS
-USE MESH_VARIABLES
 USE MESH_POINTERS
 USE COMP_FUNCTIONS, ONLY: CURRENT_TIME, GET_FILE_NUMBER
 USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE
diff --git a/Source/dump.f90 b/Source/dump.f90
index bb969ccb1bc..6084a9e28c6 100644
--- a/Source/dump.f90
+++ b/Source/dump.f90
@@ -7,7 +7,6 @@
 MODULE DUMP
 
 USE PRECISION_PARAMETERS
-USE MESH_VARIABLES
 USE GLOBAL_CONSTANTS
 USE OUTPUT_CLOCKS
 USE MESH_POINTERS
diff --git a/Source/geom.f90 b/Source/geom.f90
index ef5b7996546..db30b109da9 100644
--- a/Source/geom.f90
+++ b/Source/geom.f90
@@ -4,7 +4,6 @@ MODULE COMPLEX_GEOMETRY
 
 USE PRECISION_PARAMETERS
 USE GLOBAL_CONSTANTS
-USE MESH_VARIABLES
 USE MESH_POINTERS
 USE COMP_FUNCTIONS, ONLY: CHECKREAD,CHECK_XB,GET_FILE_NUMBER,SHUTDOWN,CURRENT_TIME
 USE MEMORY_FUNCTIONS, ONLY: ChkMemErr
diff --git a/Source/part.f90 b/Source/part.f90
index 4d974a7a903..a0af3bd19db 100644
--- a/Source/part.f90
+++ b/Source/part.f90
@@ -7,8 +7,8 @@ MODULE PART
 USE MESH_POINTERS
 USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
 USE COMPLEX_GEOMETRY, ONLY: CC_IDCC,CC_IDCF,CC_CGSC,CC_SOLID
-IMPLICIT NONE (TYPE,EXTERNAL)
 
+IMPLICIT NONE (TYPE,EXTERNAL)
 PRIVATE
 
 PUBLIC INSERT_ALL_PARTICLES,MOVE_PARTICLES,PARTICLE_MASS_ENERGY_TRANSFER,REMOVE_PARTICLES,&
@@ -1594,6 +1594,7 @@ END SUBROUTINE VOLUME_INIT_PARTICLE
 !> \brief Set up the properties of a single, newly inserted particle
 
 SUBROUTINE INITIALIZE_SINGLE_PARTICLE
+
 USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
 REAL(EB) :: AREA,SCALE_FACTOR,RADIUS,LP_VOLUME
 INTEGER :: N
diff --git a/Source/radi.f90 b/Source/radi.f90
index 093fe41cbb9..8107692fdfd 100644
--- a/Source/radi.f90
+++ b/Source/radi.f90
@@ -1,7 +1,9 @@
 MODULE SPECDATA
 
 USE PRECISION_PARAMETERS
+
 IMPLICIT NONE (TYPE,EXTERNAL)
+
 INTEGER ::  I,J
 INTEGER, PARAMETER :: NWATERK=183
 REAL(EB) :: CPLXREF_WATER(NWATERK,2)
@@ -315,15 +317,17 @@ MODULE SPECDATA
 
 END MODULE SPECDATA
 
+
 MODULE WSGG_ARRAYS
-USE PRECISION_PARAMETERS
 
+USE PRECISION_PARAMETERS
 REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: WSGG_B1_ARRAY,WSGG_B2_ARRAY,WSGG_D_ARRAY
 REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: WSGG_C_ARRAY
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: WSGG_KAPPAP1_ARRAY,WSGG_KAPPAP2_ARRAY
 
 END MODULE WSGG_ARRAYS
 
+
 MODULE MIEV
 
 USE PRECISION_PARAMETERS
@@ -2666,14 +2670,17 @@ END FUNCTION WRONG
 END SUBROUTINE TESTMI
 
 
+!> \brief Print out a warning or error message;  ABORT IF ERROR AFTER MAKING SYMBOLIC DUMP (MACHINE-SPECIFIC)
+
 SUBROUTINE ERRMSG( MESSAG, FATAL )
-! PRINT OUT A WARNING OR ERROR MESSAGE;  ABORT IF ERROR AFTER MAKING SYMBOLIC DUMP (MACHINE-SPECIFIC)
+
 USE COMP_FUNCTIONS, ONLY: SHUTDOWN
 CHARACTER :: MESSAG*(*)
 LOGICAL ::  FATAL,MSGLIM
 INTEGER ::  MAXMSG, NUMMSG
 SAVE      MAXMSG, NUMMSG, MSGLIM
 DATA      NUMMSG / 0 /,  MAXMSG / 100 /,  MSGLIM /.FALSE./
+
 IF (FATAL) CALL SHUTDOWN(MESSAG)
 NUMMSG = NUMMSG + 1
 IF( MSGLIM ) RETURN
@@ -2688,19 +2695,21 @@ SUBROUTINE ERRMSG( MESSAG, FATAL )
 END SUBROUTINE ERRMSG
 
 
+!> \brief Write names of erroneous variables and return T
+!> \param VARNAM Name of erroneous variable to be written
+
 LOGICAL FUNCTION WRTBAD( VARNAM )
-!          WRITE NAMES OF ERRONEOUS VARIABLES AND RETURN 'TRUE'
-!      INPUT :   VARNAM = NAME OF ERRONEOUS VARIABLE TO BE WRITTEN
-!                         ( CHARACTER, ANY LENGTH )
+
 CHARACTER VARNAM*(*)
 INTEGER   MAXMSG, NUMMSG
 SAVE      NUMMSG, MAXMSG
 DATA      NUMMSG / 0 /, MAXMSG / 50 /
+
 WRTBAD = .TRUE.
 NUMMSG = NUMMSG + 1
 WRITE(LU_ERR, '(3A)' ) ' ****  INPUT VARIABLE  ', VARNAM,'  IN ERROR  ****'
 IF( NUMMSG==MAXMSG ) CALL ERRMSG( 'TOO MANY INPUT ERRORS.  ABORTING...',.TRUE.)
-RETURN
+
 END FUNCTION WRTBAD
 
 
@@ -2717,7 +2726,7 @@ LOGICAL FUNCTION WRTDIM( DIMNAM, MINVAL )
 
 WRITE(LU_ERR, '(3A,I7)' ) ' ****  SYMBOLIC DIMENSION  ', DIMNAM,'  SHOULD BE INCREASED TO AT LEAST ', MINVAL
 WRTDIM = .TRUE.
-RETURN
+
 END FUNCTION WRTDIM
 
 
@@ -2739,15 +2748,15 @@ END FUNCTION SQ
 
 END MODULE MIEV
 
-MODULE RAD
 
-! Radiation heat transfer
+!> \brief Radiation heat transfer
+
+MODULE RAD
 
 USE PRECISION_PARAMETERS
 USE GLOBAL_CONSTANTS
 USE OUTPUT_CLOCKS, ONLY: RADF_CLOCK,RADF_COUNTER
 USE MESH_POINTERS
-USE MESH_VARIABLES
 USE RADCONS
 USE RADCAL_VAR
 
@@ -2762,7 +2771,7 @@ MODULE RAD
 
 
 !> \brief Initialize radiation arrays.
-!>
+
 SUBROUTINE INIT_RADIATION
 
 ! Meanings of some variables defined here:
@@ -4710,12 +4719,11 @@ REAL(EB) FUNCTION GET_KAPPA(Z_IN,TMP,IBND)
 END FUNCTION GET_KAPPA
 
 
-!==================================================================================
-!Function to compute the absorption coefficient according to Bordbar et al. (2014)
-!==================================================================================
+!> \brief Function to compute the absorption coefficient according to Bordbar et al. (2014)
+
 REAL(EB) FUNCTION KAPPA_WSGG(X_H2O, X_CO2, MOL_RATIO,PARTIAL_PRESSURE,JWSGG)
-USE WSGG_ARRAYS
 
+USE WSGG_ARRAYS
 INTEGER, INTENT(IN) :: JWSGG
 INTEGER :: NN
 REAL(EB), INTENT(IN) :: X_H2O, X_CO2, MOL_RATIO,PARTIAL_PRESSURE
@@ -4754,26 +4762,21 @@ REAL(EB) FUNCTION KAPPA_WSGG(X_H2O, X_CO2, MOL_RATIO,PARTIAL_PRESSURE,JWSGG)
 
 END FUNCTION KAPPA_WSGG
 
-!===================================================================================
-!Function to compute the temperature coefficient according to Bordbar et al. (2014)
-!===================================================================================
-REAL(EB) RECURSIVE FUNCTION A_WSGG(TTMP,MOL_RATIO,JWSGG) &
-   RESULT(A_FUNC_RES)
-USE WSGG_ARRAYS
 
+!> \brief Function to compute the temperature coefficient according to Bordbar et al. (2014)
+
+REAL(EB) RECURSIVE FUNCTION A_WSGG(TTMP,MOL_RATIO,JWSGG) RESULT(A_FUNC_RES)
+
+USE WSGG_ARRAYS
 INTEGER,INTENT(IN) :: JWSGG
 INTEGER :: MM,NN
 REAL(EB),INTENT(IN) :: MOL_RATIO,TTMP
 REAL(EB) :: TREF,SUM_A,SUM_B,SUM_C
 
-!------------------------
-!Parameters of the model
-!------------------------
 TREF = 1200._EB ! Reference temperature
 
-!------------------------------------------------------------------
-!Computing the temperature coefficient for the transparent windows
-!------------------------------------------------------------------
+! Computing the temperature coefficient for the transparent windows
+
 IF (JWSGG==5) THEN
    SUM_A = 0._EB
    DO MM=1,4
@@ -4782,19 +4785,17 @@ REAL(EB) RECURSIVE FUNCTION A_WSGG(TTMP,MOL_RATIO,JWSGG) &
    A_FUNC_RES = 1._EB - SUM_A
 
 ELSE
-!---------------------------------------------------------
-!Computing the temperature coefficient for the gray gases
-!(within three intervals of molar ratio)
-!---------------------------------------------------------
+   ! Computing the temperature coefficient for the gray gases (within three intervals of molar ratio)
+
    IF (MOL_RATIO < 0.01_EB) THEN               !Only CO2
-      !Computing the polynomial
+      ! Computing the polynomial
       SUM_B = 0._EB
       DO MM=0,4
          SUM_B = SUM_B + WSGG_B1_ARRAY(JWSGG,MM)*(TTMP/TREF)**(REAL(MM,EB))
       ENDDO
 
    ELSEIF (MOL_RATIO > 4._EB) THEN               !Only H2O
-      !Computing the polynomial
+      ! Computing the polynomial
       SUM_B = 0._EB
       DO MM=0,4
          SUM_B = SUM_B + WSGG_B2_ARRAY(JWSGG,MM)*(TTMP/TREF)**(REAL(MM,EB))
@@ -4818,22 +4819,16 @@ REAL(EB) RECURSIVE FUNCTION A_WSGG(TTMP,MOL_RATIO,JWSGG) &
 END FUNCTION A_WSGG
 
 
-!====================================================
-!Function to compute the gray absorption coefficient
-!of soot (same function as the one used by Fluent)
-!====================================================
+!> \brief Function to compute the gray absorption coefficient of soot (same function as the one used by Fluent)
+
 REAL(EB) FUNCTION KAPPA_SOOT(SOOT_MASS_CONCENTRATION,TTMP)
 
 REAL(EB),INTENT(IN) :: SOOT_MASS_CONCENTRATION,TTMP
 
-   !KAPPA_SOOT = 1232.4_EB*SOOT_MASS_CONCENTRATION*(1._EB+4.8E-4_EB*(TTMP-2000._EB))
-   KAPPA_SOOT = SOOT_MASS_CONCENTRATION*(1232.4_EB+0.591552_EB*(TTMP-2000._EB))
+KAPPA_SOOT = SOOT_MASS_CONCENTRATION*(1232.4_EB+0.591552_EB*(TTMP-2000._EB))
 
 END FUNCTION KAPPA_SOOT
 
 
 END MODULE RAD
 
-
-
-
diff --git a/Source/read.f90 b/Source/read.f90
index 4144f103acb..ed39065b688 100644
--- a/Source/read.f90
+++ b/Source/read.f90
@@ -7,7 +7,6 @@
 MODULE READ_INPUT
 
 USE PRECISION_PARAMETERS
-USE MESH_VARIABLES
 USE GLOBAL_CONSTANTS
 USE TRAN
 USE MESH_POINTERS
diff --git a/Source/turb.f90 b/Source/turb.f90
index c91046d3bab..9fb8accb7d4 100644
--- a/Source/turb.f90
+++ b/Source/turb.f90
@@ -6,7 +6,6 @@ MODULE TURBULENCE
 USE PRECISION_PARAMETERS
 USE GLOBAL_CONSTANTS
 USE MESH_POINTERS
-USE MESH_VARIABLES
 USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER,SHUTDOWN
 
 IMPLICIT NONE (TYPE,EXTERNAL)
@@ -3120,7 +3119,6 @@ MODULE MANUFACTURED_SOLUTIONS
 
 SUBROUTINE SAAD_MMS_1(NM)
 USE MESH_POINTERS
-USE MESH_VARIABLES
 USE GLOBAL_CONSTANTS, ONLY: P_INF,TMPA,N_TRACKED_SPECIES
 USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT
 INTEGER, INTENT(IN) :: NM
@@ -3152,7 +3150,6 @@ END SUBROUTINE SAAD_MMS_1
 
 SUBROUTINE SHUNN_MMS_3(DT,NM)
 USE MESH_POINTERS
-USE MESH_VARIABLES
 USE GLOBAL_CONSTANTS, ONLY: P_INF,TMPA,N_TRACKED_SPECIES
 USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT
 INTEGER, INTENT(IN) :: NM