diff --git a/.idea/misc.xml b/.idea/misc.xml
index b301bbd0cd67f8a7b66c93e1322ecf1bdeb1ebc4..ec201b41503a14c4ae2cdc34db70a1511394fe1d 100644
--- a/.idea/misc.xml
+++ b/.idea/misc.xml
@@ -3,6 +3,9 @@
   <component name="Black">
     <option name="sdkName" value="Python 3.8 (untitled)" />
   </component>
+  <component name="CMakePythonSetting">
+    <option name="pythonIntegrationState" value="YES" />
+  </component>
   <component name="CMakeWorkspace" PROJECT_DIR="$PROJECT_DIR$" />
   <component name="JsBuildToolPackageJson" sorting="DEFINITION_ORDER" />
   <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8 (untitled)" project-jdk-type="Python SDK" />
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6b6076e1cc668bdec98c13c5cdb1d3359fb63e64..35da8e83ce52be1ab2837e7e8e715312ee9155e8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -16,13 +16,12 @@ set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${dialect}")
 
 
 include_directories(src)
-set(plt_io
+set(src_io
         src/scm_io_plt.f90
-        #src/test_io_plt.f90
-        src/netcdf_io_module.f90)
-        #src/test_netcdf.f90)
+        src/io.f90)
 
-add_library(test_io ${plt_io})
+
+add_library(test_io ${src_io})
 
 target_include_directories(test_io PRIVATE ${netCDF_INCLUDE_DIR})
 
diff --git a/src/io.f90 b/src/io.f90
new file mode 100644
index 0000000000000000000000000000000000000000..4f1ee9137727671833dc50c69c2edc61ee26f2bb
--- /dev/null
+++ b/src/io.f90
@@ -0,0 +1,3389 @@
+! io.f90
+module io
+  use netcdf
+  implicit none
+
+
+  ! variable_metadata type:
+  !   This type contains metadata for a variable: the name of the variable,
+  !   the names of the dimensions, and character, integer, and real metadata.
+  !   The metadata can be written to a NetCDF file or an ASCII file.
+  type :: variable_metadata
+    character(len=32) :: name                       ! variable name
+    character(len=32), dimension(4) :: dim_names    ! dimension names
+    character(len=32) :: char_name(10) = ''         ! character metadata names
+    character(len=32) :: char_value(10) = ''        ! character metadata values
+    character(len=32) :: int_name(10) = ''          ! integer metadata names
+    integer :: int_value(10) = 0                    ! integer metadata values
+    character(len=32) :: real_name(10) = ''         ! real metadata names
+    real :: real_value(10) = 0.0                    ! real metadata values
+  end type variable_metadata
+
+  !interface write_var 
+  !  module procedure write_1d_real_nc, write_2d_real_nc, write_3d_real_nc, write_4d_real_nc
+  !  module procedure write_1d_double_nc, write_2d_double_nc, write_3d_double_nc, write_4d_double_nc
+  !  module procedure write_1d_int_nc, write_2d_int_nc, write_3d_int_nc, write_4d_int_nc
+  !  module procedure write_1d_logical_nc, write_2d_logical_nc, write_3d_logical_nc, write_4d_logical_nc
+  !end interface
+
+contains
+
+  ! Subroutine for writing a 1D real array to a NetCDF file
+  subroutine write_1d_real_nc(my_array, filepath, metadata)
+    real, dimension(:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(1)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_real, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_1d_real_nc
+
+
+
+  ! Subroutine for writing a 2D real array to a NetCDF file
+  subroutine write_2d_real_nc(my_array, filepath, metadata)
+    ! Arguments
+    real, dimension(:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(2)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_real, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_2d_real_nc
+
+
+
+  ! Subroutine for writing a 3D real array to a NetCDF file
+  subroutine write_3d_real_nc(my_array, filepath, metadata)
+    ! Arguments
+    real, dimension(:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(3)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+    call define_nc_dim(ncid, metadata%dim_names(3), size(my_array, 3), dimids(3))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_real, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_3d_real_nc
+
+
+
+  ! Subroutine for writing a 4D real array to a NetCDF file
+  subroutine write_4d_real_nc(my_array, filepath, metadata)
+    ! Arguments
+    real, dimension(:,:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(4)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+    call define_nc_dim(ncid, metadata%dim_names(3), size(my_array, 3), dimids(3))
+    call define_nc_dim(ncid, metadata%dim_names(4), size(my_array, 4), dimids(4))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_real, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_4d_real_nc
+
+
+
+  ! Subroutine for writing a 1D double array to a NetCDF file
+  subroutine write_1d_double_nc(my_array, filepath, metadata)
+    double precision, dimension(:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(1)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_double, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_1d_double_nc
+
+
+
+  ! Subroutine for writing a 2D double array to a NetCDF file
+  subroutine write_2d_double_nc(my_array, filepath, metadata)
+    ! Arguments
+    double precision, dimension(:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(2)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_double, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_2d_double_nc
+
+
+
+  ! Subroutine for writing a 3D double array to a NetCDF file
+  subroutine write_3d_double_nc(my_array, filepath, metadata)
+    ! Arguments
+    double precision, dimension(:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(3)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+    call define_nc_dim(ncid, metadata%dim_names(3), size(my_array, 3), dimids(3))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_double, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_3d_double_nc
+
+
+
+  ! Subroutine for writing a 4D double array to a NetCDF file
+  subroutine write_4d_double_nc(my_array, filepath, metadata)
+    ! Arguments
+    double precision, dimension(:,:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(4)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+    call define_nc_dim(ncid, metadata%dim_names(3), size(my_array, 3), dimids(3))
+    call define_nc_dim(ncid, metadata%dim_names(4), size(my_array, 4), dimids(4))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_double, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_4d_double_nc
+
+
+
+  ! Subroutine for writing a 1D integer array to a NetCDF file
+  subroutine write_1d_int_nc(my_array, filepath, metadata)
+    integer, dimension(:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(1)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_int, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_1d_int_nc
+
+
+
+  ! Subroutine for writing a 2D integer array to a NetCDF file
+  subroutine write_2d_int_nc(my_array, filepath, metadata)
+    ! Arguments
+    integer, dimension(:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(2)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_int, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_2d_int_nc
+
+
+
+  ! Subroutine for writing a 3D integer array to a NetCDF file
+  subroutine write_3d_int_nc(my_array, filepath, metadata)
+    ! Arguments
+    integer, dimension(:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(3)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+    call define_nc_dim(ncid, metadata%dim_names(3), size(my_array, 3), dimids(3))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_int, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_3d_int_nc
+
+
+
+  ! Subroutine for writing a 4D integer array to a NetCDF file
+  subroutine write_4d_int_nc(my_array, filepath, metadata)
+    ! Arguments
+    integer, dimension(:,:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(4)
+
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+    call define_nc_dim(ncid, metadata%dim_names(3), size(my_array, 3), dimids(3))
+    call define_nc_dim(ncid, metadata%dim_names(4), size(my_array, 4), dimids(4))
+
+    ! Defining the variable
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_int, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+
+    ! Writing the array to the file
+    status = nf90_put_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_4d_int_nc
+
+
+
+  ! Subroutine for writing a 1D logical array to a NetCDF file
+  subroutine write_1d_logical_nc(my_array, filepath, metadata)
+    ! Arguments
+    logical, dimension(:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(1)
+    integer, dimension(size(my_array)) :: int_array
+  
+    ! Convert logical array to integer array (1 for true, 0 for false)
+    int_array = merge(1, 0, my_array)
+  
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+  
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+  
+    ! Defining the variable (as integer)
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_int, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+  
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+  
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+  
+    ! Writing the integer array to the file
+    status = nf90_put_var(ncid, varid, int_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+  
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_1d_logical_nc  
+
+
+
+  ! Subroutine for writing a 2D logical array to a NetCDF file
+  subroutine write_2d_logical_nc(my_array, filepath, metadata)
+    ! Arguments
+    logical, dimension(:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(2)
+    integer, dimension(size(my_array, 1), size(my_array, 2)) :: int_array
+  
+    ! Convert logical array to integer array (1 for true, 0 for false)
+    int_array = merge(1, 0, my_array)
+  
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+  
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+  
+    ! Defining the variable (as integer)
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_int, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+  
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+  
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+  
+    ! Writing the integer array to the file
+    status = nf90_put_var(ncid, varid, int_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+  
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_2d_logical_nc
+
+
+
+  ! Subroutine for writing a 3D logical array to a NetCDF file
+  subroutine write_3d_logical_nc(my_array, filepath, metadata)
+    ! Arguments
+    logical, dimension(:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(3)
+    integer, dimension(size(my_array, 1), size(my_array, 2), size(my_array, 3)) :: int_array
+  
+    ! Convert logical array to integer array (1 for true, 0 for false)
+    int_array = merge(1, 0, my_array)
+  
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+  
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+    call define_nc_dim(ncid, metadata%dim_names(3), size(my_array, 3), dimids(3))
+  
+    ! Defining the variable (as integer)
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_int, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+  
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+  
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+  
+    ! Writing the integer array to the file
+    status = nf90_put_var(ncid, varid, int_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+  
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_3d_logical_nc  
+
+
+
+  ! Subroutine for writing a 4D logical array to a NetCDF file
+  subroutine write_4d_logical_nc(my_array, filepath, metadata)
+    ! Arguments
+    logical, dimension(:,:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(4)
+    integer, dimension(size(my_array, 1), size(my_array, 2), size(my_array, 3), size(my_array, 4)) :: int_array
+  
+    ! Convert logical array to integer array (1 for true, 0 for false)
+    int_array = merge(1, 0, my_array)
+  
+    ! Opening or creating the NetCDF file
+    call open_nc_file(filepath, ncid)
+  
+    ! Checking and defining dimensions
+    call define_nc_dim(ncid, metadata%dim_names(1), size(my_array, 1), dimids(1))
+    call define_nc_dim(ncid, metadata%dim_names(2), size(my_array, 2), dimids(2))
+    call define_nc_dim(ncid, metadata%dim_names(3), size(my_array, 3), dimids(3))
+    call define_nc_dim(ncid, metadata%dim_names(4), size(my_array, 4), dimids(4))
+  
+    ! Defining the variable (as integer)
+    status = nf90_def_var(ncid, trim(metadata%name), nf90_int, dimids, varid)
+    call handle_nc_error(status, 'Error defining variable: ' // trim(metadata%name))
+  
+    ! Writing variable attributes
+    call write_nc_var_attributes(ncid, varid, metadata)
+  
+    ! Ending define mode
+    status = nf90_enddef(ncid)
+    call handle_nc_error(status, 'Error ending NetCDF definition mode')
+  
+    ! Writing the integer array to the file
+    status = nf90_put_var(ncid, varid, int_array)
+    call handle_nc_error(status, 'Error writing array to NetCDF file')
+  
+    ! Closing the file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine write_4d_logical_nc
+
+
+
+  ! Subroutine for opening or creating a NetCDF file
+  subroutine open_nc_file(filepath, ncid)
+    character(len=*), intent(in) :: filepath
+    integer, intent(out) :: ncid
+    integer :: status
+
+    ! Checking for file existence and opening it
+    status = nf90_open(trim(filepath), nf90_write, ncid)
+    if (status == nf90_noerr) then
+      status = nf90_redef(ncid)
+      call handle_nc_error(status, 'Error switching to define mode')
+    else
+      ! If the file does not exist, create it
+      status = nf90_create(trim(filepath), nf90_clobber, ncid)
+      call handle_nc_error(status, 'Error creating NetCDF file')
+    end if
+  end subroutine open_nc_file
+
+
+
+  ! Subroutine for checking or defining dimensions
+  subroutine define_nc_dim(ncid, dim_name, dim_size, dimid)
+    integer, intent(in) :: ncid, dim_size
+    character(len=*), intent(in) :: dim_name
+    integer, intent(out) :: dimid
+    integer :: status, dim_len
+
+    ! Checking for the existence of the dimension
+    status = nf90_inq_dimid(ncid, trim(dim_name), dimid)
+    if (status == nf90_noerr) then
+      ! The dimension already exists, checking its size
+      status = nf90_inquire_dimension(ncid, dimid, len=dim_len)
+      if (dim_len /= dim_size) then
+        print *, 'Error: Dimension ', trim(dim_name), ' exists but has a different size.'
+        status = nf90_close(ncid)
+        return
+      end if
+    else
+      ! If the dimension does not exist, create it
+      status = nf90_def_dim(ncid, trim(dim_name), dim_size, dimid)
+      call handle_nc_error(status, 'Error defining dimension: ' // trim(dim_name))
+    end if
+  end subroutine define_nc_dim
+
+
+
+
+  subroutine write_nc_var_attributes(ncid, varid, metadata)
+    integer, intent(in) :: ncid, varid
+    type(variable_metadata), intent(in) :: metadata
+    integer :: status, i
+
+    ! Writing character metadata attributes
+    i = 1
+    do while (trim(metadata%char_name(i)) /= '')
+      status = nf90_put_att(ncid, varid, trim(metadata%char_name(i)), trim(metadata%char_value(i)))
+      call handle_nc_error(status, 'Error writing character attribute: '//trim(metadata%char_name(i)))
+      i = i + 1
+      if (i > size(metadata%char_name)) exit
+    end do
+
+    ! Writing integer metadata attributes
+    i = 1
+    do while (trim(metadata%int_name(i)) /= '')
+      status = nf90_put_att(ncid, varid, trim(metadata%int_name(i)), metadata%int_value(i))
+      call handle_nc_error(status, 'Error writing integer attribute: '//trim(metadata%int_name(i)))
+      i = i + 1
+      if (i > size(metadata%int_name)) exit
+    end do
+
+    ! Writing real metadata attributes
+    i = 1
+    do while (trim(metadata%real_name(i)) /= '')
+      status = nf90_put_att(ncid, varid, trim(metadata%real_name(i)), metadata%real_value(i))
+      call handle_nc_error(status, 'Error writing real attribute: '//trim(metadata%real_name(i)))
+      i = i + 1
+      if (i > size(metadata%real_name)) exit
+    end do
+  end subroutine write_nc_var_attributes
+
+
+
+  ! Subroutine for handling errors
+  subroutine handle_nc_error(status, errmsg)
+    integer, intent(in) :: status
+    character(len=*), intent(in) :: errmsg
+    character(len=256) :: netcdf_errmsg
+
+    if (status /= nf90_noerr) then
+      netcdf_errmsg = nf90_strerror(status)
+      print *, trim(errmsg), ': ', trim(netcdf_errmsg)
+      stop
+    end if
+  end subroutine handle_nc_error
+
+
+
+  ! Subroutine for writing a 1D real array to a text file
+  subroutine write_1d_real_ascii(my_array, filepath, metadata)
+    ! Arguments
+    real, dimension(:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, dim_size(1)
+    integer :: unit_number
+    character(len=8) :: format_string
+  
+    ! Get the size of the array dimension
+    dim_size(1) = size(my_array, 1)
+  
+    ! Use the specified format for precision
+    if (kind(my_array) == kind(1.0)) then
+      format_string = '(F24.7)'   ! Single precision (real)
+    else if (kind(my_array) == kind(1.0d0)) then
+      format_string = '(F24.16)'  ! Double precision (real(kind=8))
+    end if
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, format_string)
+  
+    ! Write dimension name and size on the third line
+    if (trim(metadata%dim_names(1)) /= '') then
+      write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(1)), ': ', dim_size(1)
+    end if
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 1D real array using the selected format
+    do i = 1, dim_size(1)
+      write(unit_number, format_string) my_array(i)
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_1d_real_ascii
+
+
+
+  ! Subroutine for writing a 2D real array to a text file 
+  subroutine write_2d_real_ascii(my_array, filepath, metadata)
+    ! Arguments
+    real, dimension(:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, dim_size(2)
+    integer :: unit_number
+    character(len=8) :: format_string
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+  
+    ! Use the specified format for precision
+    if (kind(my_array) == kind(1.0)) then
+      format_string = '(F24.7)'   ! Single precision (real)
+    else if (kind(my_array) == kind(1.0d0)) then
+      format_string = '(F24.16)'  ! Double precision (real(kind=8))
+    end if
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, format_string)
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 2
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 2D real array in a linearized fashion using the selected format
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        write(unit_number, format_string) my_array(i, j)
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_2d_real_ascii
+
+
+
+
+  ! Subroutine for writing a 3D real array to a text file
+  subroutine write_3d_real_ascii(my_array, filepath, metadata)
+    ! Arguments
+    real, dimension(:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, k, dim_size(3)
+    integer :: unit_number
+    character(len=8) :: format_string
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    dim_size(3) = size(my_array, 3)
+  
+    ! Use the specified format for precision
+    if (kind(my_array) == kind(1.0)) then
+      format_string = '(F24.7)'   ! Single precision (real)
+    else if (kind(my_array) == kind(1.0d0)) then
+      format_string = '(F24.16)'  ! Double precision (real(kind=8))
+    end if
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, format_string)
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 3
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 3D real array in a linearized fashion using the selected format
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          write(unit_number, format_string, advance='no') my_array(i, j, k)
+          write(unit_number, *)
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_3d_real_ascii
+
+
+
+  ! Subroutine for writing a 4D real array to a text file
+  subroutine write_4d_real_ascii(my_array, filepath, metadata)
+    ! Arguments
+    real, dimension(:,:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, k, l, dim_size(4)
+    integer :: unit_number
+    character(len=8) :: format_string
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    dim_size(3) = size(my_array, 3)
+    dim_size(4) = size(my_array, 4)
+  
+    ! Use the specified format for precision
+    if (kind(my_array) == kind(1.0)) then
+      format_string = '(F24.7)'   ! Single precision (real)
+    else if (kind(my_array) == kind(1.0d0)) then
+      format_string = '(F24.16)'  ! Double precision (real(kind=8))
+    end if
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, format_string)
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 4
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 4D real array in a linearized fashion using the selected format
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          do l = 1, dim_size(4)
+            write(unit_number, format_string) my_array(i, j, k, l)
+          end do
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_4d_real_ascii
+
+
+
+  ! Subroutine for writing a 1D double array to a text file
+  subroutine write_1d_double_ascii(my_array, filepath, metadata)
+    ! Arguments
+    real(kind=8), dimension(:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+    
+    ! Local variables
+    integer :: i, dim_size(1)
+    integer :: unit_number
+    character(len=12) :: format_string
+    
+    ! Get the size of the array dimension
+    dim_size(1) = size(my_array, 1)
+    
+    ! Set format for double precision
+    format_string = '(F24.16)'
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+    
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+    
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, format_string)
+    
+    ! Write dimension name and size on the third line
+    if (trim(metadata%dim_names(1)) /= '') then
+      write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(1)), ': ', dim_size(1)
+    end if
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+    
+    ! Write the 1D double precision array
+    do i = 1, dim_size(1)
+      write(unit_number, format_string) my_array(i)
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_1d_double_ascii
+
+
+
+  ! Subroutine for writing a 2D double array to a text file
+  subroutine write_2d_double_ascii(my_array, filepath, metadata)
+    ! Arguments
+    real(kind=8), dimension(:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+    
+    ! Local variables
+    integer :: i, j, dim_size(2)
+    integer :: unit_number
+    character(len=12) :: format_string
+    
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    
+    ! Set format for double precision
+    format_string = '(F24.16)'
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+    
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+    
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, format_string)
+    
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 2
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+    
+    ! Write the 2D double precision array
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        write(unit_number, format_string) my_array(i, j)
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_2d_double_ascii
+
+
+
+  ! Subroutine for writing a 3D double array to a text file
+  subroutine write_3d_double_ascii(my_array, filepath, metadata)
+    ! Arguments
+    real(kind=8), dimension(:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+    
+    ! Local variables
+    integer :: i, j, k, dim_size(3)
+    integer :: unit_number
+    character(len=12) :: format_string
+    
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    dim_size(3) = size(my_array, 3)
+    
+    ! Set format for double precision
+    format_string = '(F24.16)'
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+    
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+    
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, format_string)
+    
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 3
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+    
+    ! Write the 3D double precision array
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          write(unit_number, format_string) my_array(i, j, k)
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_3d_double_ascii
+
+
+  ! Subroutine for writing a 4D double array to a text file
+  subroutine write_4d_double_ascii(my_array, filepath, metadata)
+    ! Arguments
+    real(kind=8), dimension(:,:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+    
+    ! Local variables
+    integer :: i, j, k, l, dim_size(4)
+    integer :: unit_number
+    character(len=12) :: format_string
+    
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    dim_size(3) = size(my_array, 3)
+    dim_size(4) = size(my_array, 4)
+    
+    ! Set format for double precision
+    format_string = '(F24.16)'
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+    
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+    
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, format_string)
+    
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 4
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+    
+    ! Write the 4D double precision array
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          do l = 1, dim_size(4)
+            write(unit_number, format_string) my_array(i, j, k, l)
+          end do
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_4d_double_ascii
+
+
+
+  ! Subroutine for writing a 1D integer array to a text file
+  subroutine write_1d_int_ascii(my_array, filepath, metadata)
+    ! Arguments
+    integer, dimension(:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, dim_size(1)
+    integer :: unit_number
+  
+    ! Get the size of the array dimension
+    dim_size(1) = size(my_array, 1)
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, '(F24.7)')
+  
+    ! Write dimension name and size on the third line
+    if (trim(metadata%dim_names(1)) /= '') then
+      write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(1)), ': ', dim_size(1)
+    end if
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 1D integer array using the selected format
+    do i = 1, dim_size(1)
+      write(unit_number, '(I0)') my_array(i)
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_1d_int_ascii
+
+
+
+  ! Subroutine for writing a 2D integer array to a text file 
+  subroutine write_2d_int_ascii(my_array, filepath, metadata)
+    ! Arguments
+    integer, dimension(:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, dim_size(2)
+    integer :: unit_number
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, '(F24.7)')
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 2
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 2D integer array in a linearized fashion using the selected format
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        write(unit_number, '(I0)') my_array(i, j)
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_2d_int_ascii
+
+
+
+  ! Subroutine for writing a 3D integer array to a text file
+  subroutine write_3d_int_ascii(my_array, filepath, metadata)
+    ! Arguments
+    integer, dimension(:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, k, dim_size(3)
+    integer :: unit_number
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    dim_size(3) = size(my_array, 3)
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, '(F24.7)')
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 3
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 3D integer array in a linearized fashion using the selected format
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          write(unit_number, '(I0)') my_array(i, j, k)
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_3d_int_ascii
+
+
+
+  ! Subroutine for writing a 4D integer array to a text file
+  subroutine write_4d_int_ascii(my_array, filepath, metadata)
+    ! Arguments
+    integer, dimension(:,:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, k, l, dim_size(4)
+    integer :: unit_number
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    dim_size(3) = size(my_array, 3)
+    dim_size(4) = size(my_array, 4)
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, '(F24.7)')
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 4
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 4D integer array in a linearized fashion using the selected format
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          do l = 1, dim_size(4)
+            write(unit_number, '(I0)') my_array(i, j, k, l)
+          end do
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_4d_int_ascii
+
+
+
+  ! Subroutine for writing a 1D logical array to a text file as integers (1 for true, 0 for false)
+  subroutine write_1d_logical_ascii(my_array, filepath, metadata)
+    ! Arguments
+    logical, dimension(:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, dim_size(1)
+    integer :: unit_number
+    integer, dimension(size(my_array)) :: int_array
+  
+    ! Convert logical array to integer array (1 for true, 0 for false)
+    int_array = merge(1, 0, my_array)
+  
+    ! Get the size of the array dimension
+    dim_size(1) = size(my_array, 1)
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, '(F24.7)')
+  
+    ! Write dimension name and size on the third line
+    if (trim(metadata%dim_names(1)) /= '') then
+      write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(1)), ': ', dim_size(1)
+    end if
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 1D logical array (as integer)
+    do i = 1, dim_size(1)
+      write(unit_number, '(I0)') int_array(i)
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_1d_logical_ascii
+
+
+
+  ! Subroutine for writing a 2D logical array to a text file as integers (1 for true, 0 for false)
+  subroutine write_2d_logical_ascii(my_array, filepath, metadata)
+    ! Arguments
+    logical, dimension(:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, dim_size(2)
+    integer :: unit_number
+    integer, dimension(size(my_array, 1), size(my_array, 2)) :: int_array
+  
+    ! Convert logical array to integer array (1 for true, 0 for false)
+    int_array = merge(1, 0, my_array)
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, '(F24.7)')
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 2
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 2D logical array (as integer)
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        write(unit_number, '(I0)') int_array(i, j)
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_2d_logical_ascii
+
+
+
+  ! Subroutine for writing a 3D logical array to a text file as integers (1 for true, 0 for false)
+  subroutine write_3d_logical_ascii(my_array, filepath, metadata)
+    ! Arguments
+    logical, dimension(:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, k, dim_size(3)
+    integer :: unit_number
+    integer, dimension(size(my_array, 1), size(my_array, 2), size(my_array, 3)) :: int_array
+  
+    ! Convert logical array to integer array (1 for true, 0 for false)
+    int_array = merge(1, 0, my_array)
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    dim_size(3) = size(my_array, 3)
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, '(F24.7)')
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 3
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 3D logical array (as integer)
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          write(unit_number, '(I0)') int_array(i, j, k)
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_3d_logical_ascii
+
+
+
+  ! Subroutine for writing a 4D logical array to a text file as integers (1 for true, 0 for false)
+  subroutine write_4d_logical_ascii(my_array, filepath, metadata)
+    ! Arguments
+    logical, dimension(:,:,:,:), intent(in) :: my_array
+    character(len=*), intent(in) :: filepath
+    type(variable_metadata), intent(in) :: metadata
+  
+    ! Local variables
+    integer :: i, j, k, l, dim_size(4)
+    integer :: unit_number
+    integer, dimension(size(my_array, 1), size(my_array, 2), size(my_array, 3), size(my_array, 4)) :: int_array
+  
+    ! Convert logical array to integer array (1 for true, 0 for false)
+    int_array = merge(1, 0, my_array)
+  
+    ! Get the size of the array dimensions
+    dim_size(1) = size(my_array, 1)
+    dim_size(2) = size(my_array, 2)
+    dim_size(3) = size(my_array, 3)
+    dim_size(4) = size(my_array, 4)
+  
+    ! Open the file for writing
+    open(unit=unit_number, file=trim(filepath), status='replace', action='write')
+  
+    ! Write the metadata%name as the first line
+    write(unit_number, '(A)') trim(metadata%name)
+  
+    ! Write character, integer, and real metadata on the second line
+    call write_ascii_var_attributes(unit_number, metadata, '(F24.7)')
+  
+    ! Write dimension names and sizes on the third line with ';' separator
+    do i = 1, 4
+      if (trim(metadata%dim_names(i)) /= '') then
+        if (i > 1) write(unit_number, '(A)', advance='no') '; '  ! Add a separator
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%dim_names(i)), ': ', dim_size(i)
+      end if
+    end do
+    write(unit_number, *)  ! Move to the next line after writing dimensions
+  
+    ! Write the 4D logical array (as integer)
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          do l = 1, dim_size(4)
+            write(unit_number, '(I0)') int_array(i, j, k, l)
+          end do
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine write_4d_logical_ascii
+
+
+
+  ! Subroutine for writing variable attributes to a text file
+  subroutine write_ascii_var_attributes(unit_number, metadata, format_string)
+    ! Arguments
+    integer, intent(in) :: unit_number              ! File unit number
+    type(variable_metadata), intent(in) :: metadata ! Variable metadata
+    character(len=*), intent(in) :: format_string   ! Recording format for real metadata
+    integer :: i
+    character(len=32) :: real_value_str
+  
+    ! Writing character metadata
+    do i = 1, size(metadata%char_name)
+      if (trim(metadata%char_name(i)) /= '') then
+        write(unit_number, '(A, A, A)', advance='no') trim(metadata%char_name(i)), ': ', trim(metadata%char_value(i))
+        if (i < size(metadata%char_name)) write(unit_number, '(A)', advance='no') '; '
+      end if
+    end do
+  
+    ! Writing integer metadata
+    do i = 1, size(metadata%int_name)
+      if (trim(metadata%int_name(i)) /= '') then
+        write(unit_number, '(A, A, I0)', advance='no') trim(metadata%int_name(i)), ': ', metadata%int_value(i)
+        if (i < size(metadata%int_name)) write(unit_number, '(A)', advance='no') '; '
+      end if
+    end do
+  
+    ! Writing real metadata
+    do i = 1, size(metadata%real_name)
+      if (trim(metadata%real_name(i)) /= '') then
+        ! Convert the real metadata value to a string with the specified format
+        write(real_value_str, format_string) metadata%real_value(i)
+        ! Use trim and adjustl to remove trailing and leading spaces
+        write(unit_number, '(A, A, A)', advance='no') trim(metadata%real_name(i)), ': ', trim(adjustl(real_value_str))
+        if (i < size(metadata%real_name)) write(unit_number, '(A)', advance='no') '; '
+      end if
+    end do
+  
+    ! Move to the next line after writing all metadata
+    write(unit_number, *)
+  end subroutine write_ascii_var_attributes
+
+
+
+! Subroutine for reading 1D data from a NetCDF file
+  subroutine read_1d_real_nc(my_array, filepath, var_name)
+    ! Arguments
+    real, dimension(:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+  
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(1)
+    integer :: dim_size
+  
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+  
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+  
+    ! Get the dimension size
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+  
+    ! Retrieve dimension length
+    status = nf90_inquire_dimension(ncid, dimids(1), len=dim_size)
+    call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+  
+    ! Allocate the output array if needed
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size))
+    end if
+  
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+  
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  
+  end subroutine read_1d_real_nc
+
+
+
+  ! Subroutine for reading a 2D real array from a NetCDF file
+  subroutine read_2d_real_nc(my_array, filepath, var_name)
+    ! Arguments
+    real, dimension(:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+  
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(2)
+    integer :: dim_sizes(2)
+  
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+  
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+  
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+  
+    ! Retrieve dimension lengths
+    do i = 1, 2
+      status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+      call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+  
+    ! Allocate the output array if needed
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_sizes(1), dim_sizes(2)))
+    end if
+  
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+  
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  
+  end subroutine read_2d_real_nc
+
+
+
+  ! Subroutine for reading a 3D real array from a NetCDF file
+  subroutine read_3d_real_nc(my_array, filepath, var_name)
+    ! Arguments
+    real, dimension(:,:,:), allocatable,intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    character(len=32) :: my_array_dim_name
+  
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(3)
+    integer :: dim_sizes(3)
+  
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+  
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+  
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+  
+    ! Retrieve dimension lengths
+    do i = 1, 3
+      status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+      call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+  
+    ! Allocate the output array if needed
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_sizes(1), dim_sizes(2), dim_sizes(3)))
+    end if
+  
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+  
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  
+  end subroutine read_3d_real_nc
+
+
+
+  ! Subroutine for reading a 4D real array from a NetCDF file
+  subroutine read_4d_real_nc(my_array, filepath, var_name)
+    ! Arguments
+    real, dimension(:,:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+  
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(4)
+    integer :: dim_sizes(4)
+  
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+  
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+  
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+  
+    ! Retrieve dimension lengths
+    do i = 1, 4
+      status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+      call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+  
+    ! Allocate the output array if needed
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_sizes(1), dim_sizes(2), dim_sizes(3), dim_sizes(4)))
+    end if
+  
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+  
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  
+  end subroutine read_4d_real_nc
+
+
+
+! Subroutine for reading 1D double precision data from a NetCDF file
+  subroutine read_1d_double_nc(my_array, filepath, var_name)
+    ! Arguments
+    double precision, dimension(:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+  
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(1)
+    integer :: dim_size
+  
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+  
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+  
+    ! Get the dimension size
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+  
+    ! Retrieve dimension length
+    status = nf90_inquire_dimension(ncid, dimids(1), len=dim_size)
+    call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+  
+    ! Allocate the output array if needed
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size))
+    end if
+  
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+  
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  
+  end subroutine read_1d_double_nc
+
+
+
+  ! Subroutine for reading a 2D double array from a NetCDF file
+  subroutine read_2d_double_nc(my_array, filepath, var_name)
+    ! Arguments
+    double precision, dimension(:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+  
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(2)
+    integer :: dim_sizes(2)
+  
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+  
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+  
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+  
+    ! Retrieve dimension lengths
+    do i = 1, 2
+      status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+      call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+  
+    ! Allocate the output array if needed
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_sizes(1), dim_sizes(2)))
+    end if
+  
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+  
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  
+  end subroutine read_2d_double_nc
+
+
+
+  ! Subroutine for reading a 3D double array from a NetCDF file
+  subroutine read_3d_double_nc(my_array, filepath, var_name)
+    ! Arguments
+    double precision, dimension(:,:,:), allocatable,intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    character(len=32) :: my_array_dim_name
+  
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(3)
+    integer :: dim_sizes(3)
+  
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+  
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+  
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+  
+    ! Retrieve dimension lengths
+    do i = 1, 3
+      status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+      call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+  
+    ! Allocate the output array if needed
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_sizes(1), dim_sizes(2), dim_sizes(3)))
+    end if
+  
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+  
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  
+  end subroutine read_3d_double_nc
+
+
+
+  ! Subroutine for reading a 4D double array from a NetCDF file
+  subroutine read_4d_double_nc(my_array, filepath, var_name)
+    ! Arguments
+    double precision, dimension(:,:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+  
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(4)
+    integer :: dim_sizes(4)
+  
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+  
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+  
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+  
+    ! Retrieve dimension lengths
+    do i = 1, 4
+      status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+      call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+  
+    ! Allocate the output array if needed
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_sizes(1), dim_sizes(2), dim_sizes(3), dim_sizes(4)))
+    end if
+  
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+  
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  
+  end subroutine read_4d_double_nc
+
+
+  ! Subroutine for reading a 1D integer array from a NetCDF file
+  subroutine read_1d_int_nc(my_array, filepath, var_name)
+    integer, dimension(:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(1)
+    integer :: dim_size
+
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+
+    ! Get the dimension size
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+
+    ! Retrieve dimension length
+    status = nf90_inquire_dimension(ncid, dimids(1), len=dim_size)
+    call handle_nc_error(status, 'Error inquiring dimension size')
+
+    ! Allocate the output array
+    if (.not. allocated(my_array)) then
+        allocate(my_array(dim_size))
+    end if
+
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable for: ' // trim(var_name))
+
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file')
+  end subroutine read_1d_int_nc
+
+
+
+  ! Subroutine for reading a 2D integer array from a NetCDF file
+  subroutine read_2d_int_nc(my_array, filepath, var_name)
+    integer, dimension(:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(2)
+    integer :: dim_sizes(2)
+
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+
+    do i = 1, 2
+        status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+        call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+
+    ! Allocate the output array
+    if (.not. allocated(my_array)) then
+        allocate(my_array(dim_sizes(1), dim_sizes(2)))
+    end if
+
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file for: ' // trim(var_name))
+  end subroutine read_2d_int_nc
+
+
+
+  ! Subroutine for reading a 3D integer array from a NetCDF file
+  subroutine read_3d_int_nc(my_array, filepath, var_name)
+    integer, dimension(:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(3)
+    integer :: dim_sizes(3)
+
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+
+    do i = 1, 3
+        status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+        call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+
+    ! Allocate the output array
+    if (.not. allocated(my_array)) then
+        allocate(my_array(dim_sizes(1), dim_sizes(2), dim_sizes(3)))
+    end if
+
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file for: ' // trim(var_name))
+  end subroutine read_3d_int_nc
+
+
+
+  ! Subroutine for reading a 4D integer array from a NetCDF file
+  subroutine read_4d_int_nc(my_array, filepath, var_name)
+    integer, dimension(:,:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(4)
+    integer :: dim_sizes(4)
+
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+
+    do i = 1, 4
+        status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+        call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+
+    ! Allocate the output array
+    if (.not. allocated(my_array)) then
+        allocate(my_array(dim_sizes(1), dim_sizes(2), dim_sizes(3), dim_sizes(4)))
+    end if
+
+    ! Read the data into the array
+    status = nf90_get_var(ncid, varid, my_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file for: ' // trim(var_name))
+  end subroutine read_4d_int_nc
+
+
+
+  ! Subroutine for reading a 1D logical array from a NetCDF file
+  subroutine read_1d_logical_nc(my_array, filepath, var_name)
+    logical, dimension(:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+
+    ! Local variables
+    integer :: ncid, varid, status
+    integer :: dimids(1)
+    integer :: dim_size
+    integer, dimension(:), allocatable :: int_array
+
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+
+    ! Get the dimension size
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+
+    status = nf90_inquire_dimension(ncid, dimids(1), len=dim_size)
+    call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+
+    ! Allocate the integer array and logical array
+    allocate(int_array(dim_size))
+    allocate(my_array(dim_size))
+
+    ! Read the data into the integer array
+    status = nf90_get_var(ncid, varid, int_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+
+    ! Convert integer array to logical array
+    my_array = merge(.true., .false., int_array /= 0)
+
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file for: ' // trim(var_name))
+  end subroutine read_1d_logical_nc
+
+
+
+  ! Subroutine for reading a 2D logical array from a NetCDF file
+  subroutine read_2d_logical_nc(my_array, filepath, var_name)
+    logical, dimension(:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(2)
+    integer :: dim_sizes(2)
+    integer, dimension(:,:), allocatable :: int_array
+
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+
+    do i = 1, 2
+        status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+        call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+
+    ! Allocate the integer array and logical array
+    allocate(int_array(dim_sizes(1), dim_sizes(2)))
+    allocate(my_array(dim_sizes(1), dim_sizes(2)))
+
+    ! Read the data into the integer array
+    status = nf90_get_var(ncid, varid, int_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+
+    ! Convert integer array to logical array
+    my_array = merge(.true., .false., int_array /= 0)
+
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file for: ' // trim(var_name))
+  end subroutine read_2d_logical_nc
+
+
+
+  ! Subroutine for reading a 3D logical array from a NetCDF file
+  subroutine read_3d_logical_nc(my_array, filepath, var_name)
+    logical, dimension(:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(3)
+    integer :: dim_sizes(3)
+    integer, dimension(:,:,:), allocatable :: int_array
+
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+
+    do i = 1, 3
+        status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+        call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+
+    ! Allocate the integer array and logical array
+    allocate(int_array(dim_sizes(1), dim_sizes(2), dim_sizes(3)))
+    allocate(my_array(dim_sizes(1), dim_sizes(2), dim_sizes(3)))
+
+    ! Read the data into the integer array
+    status = nf90_get_var(ncid, varid, int_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+
+    ! Convert integer array to logical array
+    my_array = merge(.true., .false., int_array /= 0)
+
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file for: ' // trim(var_name))
+  end subroutine read_3d_logical_nc
+
+
+
+  ! Subroutine for reading a 4D logical array from a NetCDF file
+  subroutine read_4d_logical_nc(my_array, filepath, var_name)
+    logical, dimension(:,:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+
+    ! Local variables
+    integer :: ncid, varid, status, i
+    integer :: dimids(4)
+    integer :: dim_sizes(4)
+    integer, dimension(:,:,:,:), allocatable :: int_array
+
+    ! Open the NetCDF file for reading
+    status = nf90_open(trim(filepath), nf90_nowrite, ncid)
+    call handle_nc_error(status, 'Error opening NetCDF file: ' // trim(filepath))
+
+    ! Get the variable ID
+    status = nf90_inq_varid(ncid, trim(var_name), varid)
+    call handle_nc_error(status, 'Error getting variable ID for: ' // trim(var_name))
+
+    ! Get the dimension sizes
+    status = nf90_inquire_variable(ncid, varid, dimids=dimids)
+    call handle_nc_error(status, 'Error inquiring variable dimensions for: ' // trim(var_name))
+
+    do i = 1, 4
+        status = nf90_inquire_dimension(ncid, dimids(i), len=dim_sizes(i))
+        call handle_nc_error(status, 'Error inquiring dimension size for: ' // trim(var_name))
+    end do
+
+    ! Allocate the integer array and logical array
+    allocate(int_array(dim_sizes(1), dim_sizes(2), dim_sizes(3), dim_sizes(4)))
+    allocate(my_array(dim_sizes(1), dim_sizes(2), dim_sizes(3), dim_sizes(4)))
+
+    ! Read the data into the integer array
+    status = nf90_get_var(ncid, varid, int_array)
+    call handle_nc_error(status, 'Error reading data for variable: ' // trim(var_name))
+
+    ! Convert integer array to logical array
+    my_array = merge(.true., .false., int_array /= 0)
+
+    ! Close the NetCDF file
+    status = nf90_close(ncid)
+    call handle_nc_error(status, 'Error closing NetCDF file for: ' // trim(var_name))
+  end subroutine read_4d_logical_nc
+
+
+
+  ! Subroutine for reading a 1D real array from a text file
+  subroutine read_1d_real_ascii(my_array, filepath, var_name, start_line, dim_line)
+    real, dimension(:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+  
+    integer :: i, dim_size(1)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+  
+    data_start = 4
+    dims_line = 3
+  
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+  
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+  
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+  
+    read(unit_number, '(A)') line
+    call parse_dim_line_1d(line, dim_size)
+  
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+  
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1)))
+    end if
+  
+    do i = 1, dim_size(1)
+      read(unit_number, *) my_array(i)
+    end do
+  
+    close(unit_number)
+  end subroutine read_1d_real_ascii
+
+
+
+  ! Subroutine for reading a 2D real array from a text file
+  subroutine read_2d_real_ascii(my_array, filepath, var_name, start_line, dim_line)
+    real, dimension(:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+  
+    integer :: i, j, dim_size(2)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+  
+    data_start = 4
+    dims_line = 3
+  
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+  
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+  
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+  
+    read(unit_number, '(A)') line
+    call parse_dim_line_2d(line, dim_size)
+  
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+  
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2)))
+    end if
+  
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        read(unit_number, *) my_array(i, j)
+      end do
+    end do
+  
+    close(unit_number)
+  end subroutine read_2d_real_ascii
+
+
+
+  ! Subroutine for reading a 3D real array from a text file
+  subroutine read_3d_real_ascii(my_array, filepath, var_name, start_line, dim_line)
+    ! Arguments
+    real, dimension(:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+  
+    ! Local variables
+    integer :: i, j, k, dim_size(3)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+  
+    ! Set default values for optional arguments
+    data_start = 4   ! Default: data starts from 4th line
+    dims_line = 3    ! Default: dimensions are on the 3rd line
+  
+    ! If optional arguments are provided, override the defaults
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+  
+    ! Open the text file for reading
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+  
+    ! Skip lines until we reach the line with dimensions (dim_line)
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+  
+    ! Read the dimensions from the dim_line
+    read(unit_number, '(A)') line
+    call parse_dim_line_3d(line, dim_size)
+  
+    ! Check if data_start is the next line after dims_line
+    if (data_start > dims_line + 1) then
+      ! Skip to the line with the data (start_line)
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+  
+    ! Allocate the 3D real array
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2), dim_size(3)))
+    end if
+  
+    ! Read the array data from the file
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          read(unit_number, *) my_array(i, j, k)
+        end do
+      end do
+    end do
+  
+    ! Close the file
+    close(unit_number)
+  end subroutine read_3d_real_ascii
+
+
+
+  ! Subroutine for reading a 4D real array from a text file
+  subroutine read_4d_real_ascii(my_array, filepath, var_name, start_line, dim_line)
+    real, dimension(:,:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+    
+    integer :: i, j, k, l, dim_size(4)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+    
+    data_start = 4
+    dims_line = 3
+    
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+    
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+    
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+    
+    read(unit_number, '(A)') line
+    call parse_dim_line_4d(line, dim_size)
+    
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+    
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2), dim_size(3), dim_size(4)))
+    end if
+    
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          do l = 1, dim_size(4)
+            read(unit_number, *) my_array(i, j, k, l)
+          end do
+        end do
+      end do
+    end do
+    
+    close(unit_number)
+  end subroutine read_4d_real_ascii
+
+
+
+  ! Subroutine for reading a 1D double array from a text file
+  subroutine read_1d_double_ascii(my_array, filepath, var_name, start_line, dim_line)
+    double precision, dimension(:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+  
+    integer :: i, dim_size(1)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+  
+    data_start = 4
+    dims_line = 3
+  
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+  
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+  
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+  
+    read(unit_number, '(A)') line
+    call parse_dim_line_1d(line, dim_size)
+  
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+  
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1)))
+    end if
+  
+    do i = 1, dim_size(1)
+      read(unit_number, *) my_array(i)
+    end do
+  
+    close(unit_number)
+  end subroutine read_1d_double_ascii
+
+
+
+  ! Subroutine for reading a 2D double array from a text file
+  subroutine read_2d_double_ascii(my_array, filepath, var_name, start_line, dim_line)
+    double precision, dimension(:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+  
+    integer :: i, j, dim_size(2)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+  
+    data_start = 4
+    dims_line = 3
+  
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+  
+    open(unit_number, file=trim(filepath), status='old', action='read')
+  
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+  
+    read(unit_number, '(A)') line
+    call parse_dim_line_2d(line, dim_size)
+  
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+  
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2)))
+    end if
+  
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        read(unit_number, *) my_array(i, j)
+      end do
+    end do
+  
+    close(unit_number)
+  end subroutine read_2d_double_ascii
+
+
+
+  ! Subroutine for reading a 3D double array from a text file
+  subroutine read_3d_double_ascii(my_array, filepath, var_name, start_line, dim_line)
+    double precision, dimension(:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+  
+    integer :: i, j, k, dim_size(3)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+  
+    data_start = 4
+    dims_line = 3
+  
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+  
+    open(unit_number, file=trim(filepath), status='old', action='read')
+  
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+  
+    read(unit_number, '(A)') line
+    call parse_dim_line_3d(line, dim_size)
+  
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+  
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2), dim_size(3)))
+    end if
+  
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          read(unit_number, *) my_array(i, j, k)
+        end do
+      end do
+    end do
+  
+    close(unit_number)
+  end subroutine read_3d_double_ascii
+
+
+
+  ! Subroutine for reading a 4D double array from a text file
+  subroutine read_4d_double_ascii(my_array, filepath, var_name, start_line, dim_line)
+    double precision, dimension(:,:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+  
+    integer :: i, j, k, l, dim_size(4)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+  
+    data_start = 4
+    dims_line = 3
+  
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+  
+    open(unit_number, file=trim(filepath), status='old', action='read')
+  
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+  
+    read(unit_number, '(A)') line
+    call parse_dim_line_4d(line, dim_size)
+  
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+  
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2), dim_size(3), dim_size(4)))
+    end if
+  
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          do l = 1, dim_size(4)
+            read(unit_number, *) my_array(i, j, k, l)
+          end do
+        end do
+      end do
+    end do
+  
+    close(unit_number)
+  end subroutine read_4d_double_ascii
+
+
+
+  ! Subroutine for reading a 1D integer array from a text file
+  subroutine read_1d_int_ascii(my_array, filepath, var_name, start_line, dim_line)
+    integer, dimension(:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+    ! Local variables
+    integer :: i, dim_size(1)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+
+    ! Set default values for optional arguments
+    data_start = 4   ! Default: data starts from 4th line
+    dims_line = 3    ! Default: dimensions are on the 3rd line
+
+    ! If optional arguments are provided, override the defaults
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+
+    ! Open the text file for reading
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+
+    ! Skip lines until we reach the line with dimensions (dim_line)
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+
+    ! Read the dimensions from the dim_line
+    read(unit_number, '(A)') line
+    call parse_dim_line_1d(line, dim_size)
+
+    ! Skip lines until we reach the line with data (data_start)
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1)))
+    end if
+
+    ! Read the data from the data lines
+    do i = 1, dim_size(1)
+      read(unit_number, *) my_array(i)
+    end do
+
+    close(unit_number)
+  end subroutine read_1d_int_ascii
+
+
+
+  ! Subroutine for reading a 2D integer array from a text file
+  subroutine read_2d_int_ascii(my_array, filepath, var_name, start_line, dim_line)
+    integer, dimension(:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+    ! Local variables
+    integer :: i, j, dim_size(2)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+
+    ! Set default values for optional arguments
+    data_start = 4   ! Default: data starts from 4th line
+    dims_line = 3    ! Default: dimensions are on the 3rd line
+
+    ! If optional arguments are provided, override the defaults
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+
+    ! Open the text file for reading
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+
+    ! Skip lines until we reach the line with dimensions (dim_line)
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+
+    ! Read the dimensions from the dim_line
+    read(unit_number, '(A)') line
+    call parse_dim_line_2d(line, dim_size)
+
+    ! Skip lines until we reach the line with data (data_start)
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2)))
+    end if
+
+    ! Read the data from the data lines
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        read(unit_number, *) my_array(i, j)
+      end do
+    end do
+
+    close(unit_number)
+  end subroutine read_2d_int_ascii
+
+
+
+  ! Subroutine for reading a 3D integer array from a text file
+  subroutine read_3d_int_ascii(my_array, filepath, var_name, start_line, dim_line)
+    integer, dimension(:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+    ! Local variables
+    integer :: i, j, k, dim_size(3)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+
+    ! Set default values for optional arguments
+    data_start = 4   ! Default: data starts from 4th line
+    dims_line = 3    ! Default: dimensions are on the 3rd line
+
+    ! If optional arguments are provided, it override the defaults
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+
+    ! Open the text file for reading
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+
+    ! Skip lines until we reach the line with dimensions (dim_line)
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+
+    ! Read the dimensions from the dim_line
+    read(unit_number, '(A)') line
+    call parse_dim_line_3d(line, dim_size)
+
+    ! Skip lines until we reach the line with data (data_start)
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2), dim_size(3)))
+    end if
+
+    ! Read the data from the data lines
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          read(unit_number, *) my_array(i, j, k)
+        end do
+      end do
+    end do
+
+    close(unit_number)
+  end subroutine read_3d_int_ascii
+
+
+
+  ! Subroutine for reading a 4D integer array from a text file
+  subroutine read_4d_int_ascii(my_array, filepath, var_name, start_line, dim_line)
+    integer, dimension(:,:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+
+    integer :: i, j, k, l, dim_size(4)
+    integer :: unit_number, line_num, data_start, dims_line
+    character(len=256) :: line
+
+    ! Set default values for optional arguments
+    data_start = 4   ! Default: data starts from 4th line
+    dims_line = 3    ! Default: dimensions are on the 3rd line
+
+    ! If optional arguments are provided, it override the defaults
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+
+    ! Open the text file for reading
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+
+    ! Skip lines until we reach the line with dimensions (dim_line)
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+
+    ! Read the dimensions from the dim_line
+    read(unit_number, '(A)') line
+    call parse_dim_line_4d(line, dim_size)
+
+    ! Skip lines until we reach the line with data (data_start)
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2), dim_size(3), dim_size(4)))
+    end if
+
+    ! Read the data from the data lines
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          do l = 1, dim_size(4)
+            read(unit_number, *) my_array(i, j, k, l)
+          end do
+        end do
+      end do
+    end do
+
+    close(unit_number)
+  end subroutine read_4d_int_ascii
+
+
+
+  ! Subroutine for reading a 1D logical array from a text file
+  subroutine read_1d_logical_ascii(my_array, filepath, var_name, start_line, dim_line)
+    logical, dimension(:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+
+    integer :: i, dim_size(1)
+    integer :: unit_number, line_num, data_start, dims_line
+    integer, dimension(:), allocatable :: int_array
+    character(len=256) :: line
+
+    data_start = 4
+    dims_line = 3
+
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+
+    read(unit_number, '(A)') line
+    call parse_dim_line_1d(line, dim_size)
+
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+
+    if (.not. allocated(int_array)) then
+      allocate(int_array(dim_size(1)))
+    end if
+
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1)))
+    end if
+
+    do i = 1, dim_size(1)
+      read(unit_number, *) int_array(i)
+    end do
+
+    my_array = merge(.true., .false., int_array /= 0)
+
+    close(unit_number)
+  end subroutine read_1d_logical_ascii
+
+
+
+  ! Subroutine for reading a 2D logical array from a text file
+  subroutine read_2d_logical_ascii(my_array, filepath, var_name, start_line, dim_line)
+    logical, dimension(:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+
+    integer :: i, j, dim_size(2)
+    integer :: unit_number, line_num, data_start, dims_line
+    integer, dimension(:,:), allocatable :: int_array
+    character(len=256) :: line
+
+    data_start = 4
+    dims_line = 3
+
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+
+    read(unit_number, '(A)') line
+    call parse_dim_line_2d(line, dim_size)
+
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+
+    if (.not. allocated(int_array)) then
+      allocate(int_array(dim_size(1), dim_size(2)))
+    end if
+
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2)))
+    end if
+
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        read(unit_number, *) int_array(i, j)
+      end do
+    end do
+
+    my_array = merge(.true., .false., int_array /= 0)
+
+    close(unit_number)
+  end subroutine read_2d_logical_ascii
+
+
+
+  ! Subroutine for reading a 3D logical array from a text file
+  subroutine read_3d_logical_ascii(my_array, filepath, var_name, start_line, dim_line)
+    logical, dimension(:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+
+    integer :: i, j, k, dim_size(3)
+    integer :: unit_number, line_num, data_start, dims_line
+    integer, dimension(:,:,:), allocatable :: int_array
+    character(len=256) :: line
+
+    data_start = 4
+    dims_line = 3
+
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+
+    read(unit_number, '(A)') line
+    call parse_dim_line_3d(line, dim_size)
+
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+
+    if (.not. allocated(int_array)) then
+      allocate(int_array(dim_size(1), dim_size(2), dim_size(3)))
+    end if
+
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2), dim_size(3)))
+    end if
+
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          read(unit_number, *) int_array(i, j, k)
+        end do
+      end do
+    end do
+
+    my_array = merge(.true., .false., int_array /= 0)
+
+    close(unit_number)
+  end subroutine read_3d_logical_ascii
+
+
+
+  ! Subroutine for reading a 4D logical array from a text file
+  subroutine read_4d_logical_ascii(my_array, filepath, var_name, start_line, dim_line)
+    logical, dimension(:,:,:,:), allocatable, intent(out) :: my_array
+    character(len=*), intent(in) :: filepath
+    character(len=*), intent(in) :: var_name
+    integer, intent(in), optional :: start_line, dim_line
+
+    integer :: i, j, k, l, dim_size(4)
+    integer :: unit_number, line_num, data_start, dims_line
+    integer, dimension(:,:,:,:), allocatable :: int_array
+    character(len=256) :: line
+
+    data_start = 4
+    dims_line = 3
+
+    if (present(start_line)) data_start = start_line
+    if (present(dim_line)) dims_line = dim_line
+
+    open(unit=unit_number, file=trim(filepath), status='old', action='read')
+
+    do line_num = 1, dims_line-1
+      read(unit_number, '(A)')
+    end do
+
+    read(unit_number, '(A)') line
+    call parse_dim_line_4d(line, dim_size)
+
+    if (data_start > dims_line + 1) then
+      do line_num = dims_line+1, data_start-1
+        read(unit_number, '(A)')
+      end do
+    end if
+
+    if (.not. allocated(int_array)) then
+      allocate(int_array(dim_size(1), dim_size(2), dim_size(3), dim_size(4)))
+    end if
+
+    if (.not. allocated(my_array)) then
+      allocate(my_array(dim_size(1), dim_size(2), dim_size(3), dim_size(4)))
+    end if
+
+    do i = 1, dim_size(1)
+      do j = 1, dim_size(2)
+        do k = 1, dim_size(3)
+          do l = 1, dim_size(4)
+            read(unit_number, *) int_array(i, j, k, l)
+          end do
+        end do
+      end do
+    end do
+
+    my_array = merge(.true., .false., int_array /= 0)
+
+    close(unit_number)
+  end subroutine read_4d_logical_ascii
+
+
+
+  ! Subroutine for parsing 1 dimension line
+  subroutine parse_dim_line_1d(line, dim_size)
+    character(len=*), intent(inout) :: line
+    integer, dimension(1), intent(out) :: dim_size
+    integer :: io_status
+    character(len=32) :: token, key, value
+    character(len=1) :: delimiter
+    delimiter = ';'
+  
+    dim_size = 0
+    call get_next_token(line, token, delimiter)
+    if (len_trim(token) > 0) then
+      call split_key_value(token, key, value)
+      if (len_trim(value) > 0) then
+        read(value, *, iostat=io_status) dim_size(1)
+        if (io_status /= 0) then
+          print *, 'Error: Failed to read dimension size for ', key
+          stop
+        end if
+      else
+        print *, 'Error: Missing value for dimension ', key
+        stop
+      end if
+    else
+      print *, 'Error: Missing dimension information in the input file.'
+      stop
+    end if
+  end subroutine parse_dim_line_1d
+
+
+
+  ! Subroutine for parsing 2 dimension line
+  subroutine parse_dim_line_2d(line, dim_size)
+    character(len=*), intent(inout) :: line
+    integer, dimension(2), intent(out) :: dim_size
+    integer :: i, io_status
+    character(len=32) :: token, key, value
+    character(len=1) :: delimiter
+    delimiter = ';'
+  
+    dim_size = 0
+    do i = 1, 2
+      call get_next_token(line, token, delimiter)
+      if (len_trim(token) > 0) then
+        call split_key_value(token, key, value)
+        if (len_trim(value) > 0) then
+          read(value, *, iostat=io_status) dim_size(i)
+          if (io_status /= 0) then
+            print *, 'Error: Failed to read dimension size for ', key
+            stop
+          end if
+        else
+          print *, 'Error: Missing value for dimension ', key
+          stop
+        end if
+      else
+        print *, 'Error: Missing dimension information in the input file.'
+        stop
+      end if
+    end do
+  end subroutine parse_dim_line_2d
+
+
+
+  ! Subroutine for parsing the 3 dimension line
+  subroutine parse_dim_line_3d(line, dim_size)
+    character(len=*), intent(inout) :: line
+    integer, dimension(3), intent(out) :: dim_size
+    integer :: i, io_status
+    character(len=32) :: token, key, value
+    character(len=1) :: delimiter
+    delimiter = ';'
+  
+    ! Initialize dim_size to zero to catch uninitialized values
+    dim_size = 0
+  
+    ! Split the line by the delimiter and extract dimension sizes
+    do i = 1, 3
+      ! Get the next token (key:value pair)
+      call get_next_token(line, token, delimiter)
+  
+      ! Check if token is not empty
+      if (len_trim(token) > 0) then
+        call split_key_value(token, key, value)
+        
+        ! Ensure the value part is present and is a valid number
+        if (len_trim(value) > 0) then
+          read(value, *, iostat=io_status) dim_size(i)
+          if (io_status /= 0) then
+            print *, 'Error: Failed to read dimension size for ', key
+            stop
+          end if
+        else
+          print *, 'Error: Missing value for dimension ', key
+          stop
+        end if
+      else
+        print *, 'Error: Missing dimension information in the input file.'
+        stop
+      end if
+    end do
+  end subroutine parse_dim_line_3d
+
+
+
+  ! Subroutine for parsing the 4 dimension line
+  subroutine parse_dim_line_4d(line, dim_size)
+    character(len=*), intent(inout) :: line
+    integer, dimension(4), intent(out) :: dim_size
+    integer :: i, io_status
+    character(len=32) :: token, key, value
+    character(len=1) :: delimiter
+    delimiter = ';'
+    
+    dim_size = 0
+    do i = 1, 4
+      call get_next_token(line, token, delimiter)
+      if (len_trim(token) > 0) then
+        call split_key_value(token, key, value)
+        if (len_trim(value) > 0) then
+          read(value, *, iostat=io_status) dim_size(i)
+          if (io_status /= 0) then
+            print *, 'Error: Failed to read dimension size for ', key
+            stop
+          end if
+        else
+          print *, 'Error: Missing value for dimension ', key
+          stop
+        end if
+      else
+        print *, 'Error: Missing dimension information in the input file.'
+        stop
+      end if
+    end do
+  end subroutine parse_dim_line_4d
+
+
+
+  ! Helper subroutine for splitting "key: value" into key and value
+  subroutine split_key_value(token, key, value)
+    character(len=*), intent(in) :: token
+    character(len=*), intent(out) :: key, value
+    integer :: pos
+  
+    ! Find the position of the colon (:) and split the token into key and value
+    pos = index(token, ':')
+    if (pos > 0) then
+      key = trim(token(1:pos-1))
+      value = trim(token(pos+1:))
+    else
+      key = trim(token)
+      value = ''
+    end if
+  end subroutine split_key_value
+
+
+
+  ! Helper subroutine for token extraction
+  subroutine get_next_token(line, token, delimiter)
+    character(len=*), intent(inout) :: line
+    character(len=*), intent(out) :: token
+    character(len=1), intent(in) :: delimiter
+    integer :: pos
+  
+    ! Find the position of the delimiter and extract the token
+    pos = index(line, delimiter)
+    if (pos > 0) then
+      token = line(1:pos-1)
+      line = line(pos+len(delimiter):)
+    else
+      token = trim(line)
+      line = ''
+    end if
+  end subroutine get_next_token  
+
+end module io
diff --git a/src/netcdf_io_module.f90 b/src/netcdf_io_module.f90
deleted file mode 100644
index 9d18197e94a9c3cb4a5dee07c22f4d8f626ca2cc..0000000000000000000000000000000000000000
--- a/src/netcdf_io_module.f90
+++ /dev/null
@@ -1,211 +0,0 @@
-module netcdf_io_module
-    use netcdf
-    implicit none
-
-    type :: io_struct
-        integer :: ncid
-        integer :: time_dimid, height_dimid, lat_dimid, lon_dimid
-        integer :: varid_time, varid_series, varid_timescan, varid_profile
-        integer :: varid_2Dtime, varid_3Dtime, varid_2D, varid_3D
-        logical :: is_open = .false.
-        character(len=128) :: series_name = 'timeseries of variable'
-        character(len=128) :: series_long_name = 'timeseries of some variable'
-        character(len=128) :: series_units = 'unit'
-        character(len=128) :: timescan_name = 'vertical timescan of variable'
-        character(len=128) :: timescan_long_name = 'vertical timescan of some variable'
-        character(len=128) :: timescan_units = 'unit'
-        character(len=128) :: profile_name = 'vertical profile of variable'
-        character(len=128) :: profile_long_name = 'vertical profile of some variable'
-        character(len=128) :: profile_units = 'unit'
-        character(len=128) :: d2time_name = '2D variable with time'
-        character(len=128) :: d2time_long_name = '2D variable with time'
-        character(len=128) :: d2time_units = 'unit'
-        character(len=128) :: d3time_name = '3D variable with time'
-        character(len=128) :: d3time_long_name = '3D variable with time'
-        character(len=128) :: d3time_units = 'unit'
-        character(len=128) :: d2_name = '2D variable'
-        character(len=128) :: d2_long_name = '2D variable'
-        character(len=128) :: d2_units = 'unit'
-        character(len=128) :: d3_name = '3D variable'
-        character(len=128) :: d3_long_name = '3D variable'
-        character(len=128) :: d3_units = 'unit'
-    end type io_struct
-
-contains
-    subroutine open_netcdf(filename, ios, time_len, height_len, lat_len, lon_len)
-        character(len=*), intent(in) :: filename
-        type(io_struct), intent(out) :: ios
-        integer, intent(in) :: time_len, height_len, lat_len, lon_len
-        integer :: ierr
-
-        ierr = nf90_create(filename, nf90_clobber, ios%ncid)
-        if (ierr == nf90_noerr) then
-            ios%is_open = .true.
-            ierr = nf90_def_dim(ios%ncid, 'time', time_len, ios%time_dimid)
-            ierr = nf90_def_dim(ios%ncid, 'height', height_len, ios%height_dimid)
-            ierr = nf90_def_dim(ios%ncid, 'lat', lat_len, ios%lat_dimid)
-            ierr = nf90_def_dim(ios%ncid, 'lon', lon_len, ios%lon_dimid)
-            ierr = nf90_def_var(ios%ncid, 'time', nf90_float, ios%time_dimid, ios%varid_time)
-            ierr = nf90_def_var(ios%ncid, 'series', nf90_float, ios%time_dimid, ios%varid_series)
-            ierr = nf90_def_var(ios%ncid, 'timescan', nf90_float, &
-                    [ios%time_dimid, ios%height_dimid], ios%varid_timescan)
-            ierr = nf90_def_var(ios%ncid, 'profile', nf90_float, ios%height_dimid, ios%varid_profile)
-            ierr = nf90_def_var(ios%ncid, '2Dtime', nf90_float, &
-                    [ios%time_dimid, ios%lat_dimid, ios%lon_dimid], ios%varid_2Dtime)
-            ierr = nf90_def_var(ios%ncid, '3Dtime', nf90_float, &
-                    [ios%time_dimid, ios%lat_dimid, ios%lon_dimid, ios%height_dimid], ios%varid_3Dtime)
-            ierr = nf90_def_var(ios%ncid, '2D', nf90_float, &
-                    [ios%lat_dimid, ios%lon_dimid], ios%varid_2D)
-            ierr = nf90_def_var(ios%ncid, '3D', nf90_float, &
-                    [ios%lat_dimid, ios%lon_dimid, ios%height_dimid], ios%varid_3D)
-
-            ierr = nf90_put_att(ios%ncid, ios%varid_series, 'standard_name', ios%series_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_series, 'long_name', ios%series_long_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_series, 'units', ios%series_units)
-            ierr = nf90_put_att(ios%ncid, ios%varid_timescan, 'standard_name', ios%timescan_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_timescan, 'long_name', ios%timescan_long_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_timescan, 'units', ios%timescan_units)
-            ierr = nf90_put_att(ios%ncid, ios%varid_profile, 'standard_name', ios%profile_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_profile, 'long_name', ios%profile_long_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_profile, 'units', ios%profile_units)
-            ierr = nf90_put_att(ios%ncid, ios%varid_2Dtime, 'standard_name', ios%d2time_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_2Dtime, 'long_name', ios%d2time_long_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_2Dtime, 'units', ios%d2time_units)
-            ierr = nf90_put_att(ios%ncid, ios%varid_3Dtime, 'standard_name', ios%d3time_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_3Dtime, 'long_name', ios%d3time_long_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_3Dtime, 'units', ios%d3time_units)
-            ierr = nf90_put_att(ios%ncid, ios%varid_2D, 'standard_name', ios%d2_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_2D, 'long_name', ios%d2_long_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_2D, 'units', ios%d2_units)
-            ierr = nf90_put_att(ios%ncid, ios%varid_3D, 'standard_name', ios%d3_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_3D, 'long_name', ios%d3_long_name)
-            ierr = nf90_put_att(ios%ncid, ios%varid_3D, 'units', ios%d3_units)
-
-            ierr = nf90_enddef(ios%ncid)
-        else
-            print *, 'Error opening file: ', ierr
-            ios%is_open = .false.
-        endif
-    end subroutine open_netcdf
-
-    subroutine write_series(ios, time_data, series_data)
-        type(io_struct), intent(in) :: ios
-        real, dimension(:), intent(in) :: time_data, series_data
-        integer :: ierr
-
-        if (.not. ios%is_open) then
-            print *, "File is not open."
-            return
-        endif
-
-        ierr = nf90_put_var(ios%ncid, ios%varid_time, time_data)
-        ierr = nf90_put_var(ios%ncid, ios%varid_series, series_data)
-        if (ierr /= nf90_noerr) then
-            print *, 'Error writing series data:', ierr
-        endif
-    end subroutine write_series
-
-    subroutine write_timescan(ios, timescan_data)
-        type(io_struct), intent(in) :: ios
-        real, dimension(:,:), intent(in) :: timescan_data
-        integer :: ierr
-
-        if (.not. ios%is_open) then
-            print *, "File is not open."
-            return
-        endif
-
-        ierr = nf90_put_var(ios%ncid, ios%varid_timescan, timescan_data)
-        if (ierr /= nf90_noerr) then
-            print *, 'Error writing timescan data:', ierr
-        endif
-    end subroutine write_timescan
-
-    subroutine write_profile(ios, profile_data)
-        type(io_struct), intent(in) :: ios
-        real, dimension(:), intent(in) :: profile_data
-        integer :: ierr
-
-        if (.not. ios%is_open) then
-            print *, "File is not open."
-            return
-        endif
-
-        ierr = nf90_put_var(ios%ncid, ios%varid_profile, profile_data)
-        if (ierr /= nf90_noerr) then
-            print *, 'Error writing profile data:', ierr
-        endif
-    end subroutine write_profile
-
-    subroutine write_2Dtime(ios, d2time_data)
-        type(io_struct), intent(in) :: ios
-        real, dimension(:,:,:) :: d2time_data
-        integer :: ierr
-
-        if (.not. ios%is_open) then
-            print *, "File is not open."
-            return
-        endif
-
-        ierr = nf90_put_var(ios%ncid, ios%varid_2Dtime, d2time_data)
-        if (ierr /= nf90_noerr) then
-            print *, 'Error writing 2Dtime data:', ierr
-        endif
-    end subroutine write_2Dtime
-
-    subroutine write_3Dtime(ios, d3time_data)
-        type(io_struct), intent(in) :: ios
-        real, dimension(:,:,:,:) :: d3time_data
-        integer :: ierr
-
-        if (.not. ios%is_open) then
-            print *, "File is not open."
-            return
-        endif
-
-        ierr = nf90_put_var(ios%ncid, ios%varid_3Dtime, d3time_data)
-        if (ierr /= nf90_noerr) then
-            print *, "Error writing 3Dtime data:", ierr
-        endif
-    end subroutine write_3Dtime
-
-    subroutine write_2D(ios, d2_data)
-        type(io_struct), intent(in) :: ios
-        real, dimension(:,:) :: d2_data
-        integer :: ierr
-
-        if (.not. ios%is_open) then
-            print *, "File is not open."
-            return
-        endif
-
-        ierr = nf90_put_var(ios%ncid, ios%varid_2D, d2_data)
-        if (ierr /= nf90_noerr) then
-            print *, 'Error writing 2D data:', ierr
-        endif
-    end subroutine write_2D
-
-    subroutine write_3D(ios, d3_data)
-        type(io_struct), intent(in) :: ios
-        real, dimension(:,:,:) :: d3_data
-        integer :: ierr
-
-        if (.not. ios%is_open) then
-            print *, "File is not open."
-            return
-        endif
-
-        ierr = nf90_put_var(ios%ncid, ios%varid_3D, d3_data)
-        if (ierr /= nf90_noerr) then
-            print *, 'Error writing 3D data:', ierr
-        endif
-    end subroutine write_3D
-
-    subroutine close_netcdf(ios)
-        type(io_struct), intent(inout) :: ios
-        integer :: ierr
-
-        ierr = nf90_close(ios%ncid)
-        ios%is_open = .false.
-    end subroutine close_netcdf
-end module netcdf_io_module