Index: palm/trunk/UTIL/inifor/src/inifor_io.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/inifor_io.f90	(revision 3537)
+++ palm/trunk/UTIL/inifor/src/inifor_io.f90	(revision 3557)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
+! 
+! 
+! 3537 2018-11-20 10:53:14Z eckhard
 ! New routine get_netcdf_dim_vector()
 ! 
@@ -73,5 +77,5 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -92,4 +96,12 @@
     IMPLICIT NONE
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
+!> variable 'in_var % name' from the file 'in_file'. The netCDF data array is
+!> stored in the 'buffer' array and metadata added to the respective members of
+!> the given 'in_var'.
+!------------------------------------------------------------------------------!
     INTERFACE get_netcdf_variable
         MODULE PROCEDURE get_netcdf_variable_int
@@ -101,4 +113,10 @@
  CONTAINS
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_variable_int() implements the integer variant for the
+!> get_netcdf_variable interface.
+!------------------------------------------------------------------------------!
     SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
 
@@ -139,4 +157,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_variable_real() implements the real variant for the
+!> get_netcdf_variable interface.
+!------------------------------------------------------------------------------!
     SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
 
@@ -177,9 +201,15 @@
 
 
-    SUBROUTINE get_netcdf_dim_vector(filename, varname, array)
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
+!> netCDF file 'filename'.
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
 
        CHARACTER(LEN=*), INTENT(IN)         ::  filename
-       CHARACTER(LEN=*), INTENT(IN)         ::  varname
-       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  array(:)
+       CHARACTER(LEN=*), INTENT(IN)         ::  coordname
+       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  coords(:)
 
        INTEGER ::  ncid, varid, dimlen
@@ -187,15 +217,15 @@
 
        IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
-            nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR )  THEN
+            nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR )  THEN
 
           CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
           CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
 
-          ALLOCATE(array(dimlen))
-          CALL check(nf90_get_var( ncid, varid, array ))
+          ALLOCATE(coords(dimlen))
+          CALL check(nf90_get_var( ncid, varid, coords))
 
        ELSE
 
-          message = "Failed to read '" // TRIM(varname) // &
+          message = "Failed to read '" // TRIM(coordname) // &
              "' from file '" // TRIM(filename) // "'."
           CALL abort('get_netcdf_dim_vector', message)
@@ -206,4 +236,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
+!> by 'in_var % name' into 'in_var' data structure.
+!------------------------------------------------------------------------------!
     SUBROUTINE get_input_dimensions(in_var, ncid)
 
@@ -232,4 +268,11 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> get_netcdf_start_and_count() gets the start position and element counts for
+!> the given netCDF file. This information is used in get_netcdf_variable_int()
+!> and _real() for reading input variables..
+!------------------------------------------------------------------------------!
     SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
 
@@ -251,8 +294,10 @@
        start = (/ 1, 1, 1 /)
        IF ( TRIM(in_var % name) .EQ. 'T_SO' )  THEN
-          ! Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
+!
+!--       Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
           in_var % dimlen(3) = in_var % dimlen(3) - 1
 
-          ! Start reading from second level, e.g. depth = 0.005 instead of 0.0
+!
+!--       Start reading from second level, e.g. depth = 0.005 instead of 0.0
           start(3) = 2
        END IF
@@ -269,4 +314,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
+!> output.
+!------------------------------------------------------------------------------!
     SUBROUTINE netcdf_define_variable(var, ncid)
 
@@ -286,4 +337,11 @@
     
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
+!> them in the given the 'var' data structure. This information is used later 
+!> for writing output variables in update_output().
+!------------------------------------------------------------------------------!
     SUBROUTINE netcdf_get_dimensions(var, ncid)
 
@@ -293,5 +351,4 @@
         CHARACTER(SNAME)            ::  null
 
-        ! Remember dimension lenghts for NetCDF writing routine
         DO i = 1, var % ndim
            CALL check(nf90_inquire_dimension(ncid, var % dimids(i), &
@@ -306,8 +363,6 @@
 ! Description:
 ! ------------
-!> This routine initializes Inifor. This includes parsing command-line
-!> arguments, setting the names of the input and output file names as well as
-!> the name of the input namelist and, subsequently, reading in and setting grid
-!> parameters for the PALM-4U computational grid.
+!> This routine parses and interpretes the command-line options and stores the
+!> resulting settings in the 'cfg' data structure.
 !------------------------------------------------------------------------------!
     SUBROUTINE parse_command_line_arguments( cfg )
@@ -330,10 +385,4 @@
        IF (arg_count .GT. 0)  THEN
 
-          ! Every option should have an argument.
-          !IF ( MOD(arg_count, 2) .NE. 0 )  THEN
-          !   message = "Syntax error in command line."
-          !   CALL abort('parse_command_line_arguments', message)
-          !END IF
-          
           message = "The -clon and -clat command line options are depricated. " // &
              "Please remove them form your inifor command and specify the " // &
@@ -360,10 +409,8 @@
                 cfg % debug = .TRUE.
 
-             ! Elevation of the PALM-4U domain above sea level
              CASE( '-z0', '-z', '--elevation' )
                 CALL get_option_argument( i, arg )
                 READ(arg, *) cfg % z0
 
-             ! surface pressure, at z0
              CASE( '-p0', '-r', '--surface-pressure' )
                 cfg % p0_is_set = .TRUE.
@@ -371,5 +418,4 @@
                 READ(arg, *) cfg % p0
 
-             ! geostrophic wind in x direction
              CASE( '-ug', '-u', '--geostrophic-u' )
                 cfg % ug_is_set = .TRUE.
@@ -377,5 +423,4 @@
                 READ(arg, *) cfg % ug
 
-             ! geostrophic wind in y direction
              CASE( '-vg', '-v', '--geostrophic-v' )
                 cfg % vg_is_set = .TRUE.
@@ -383,11 +428,6 @@
                 READ(arg, *) cfg % vg
 
-             ! domain centre geographical longitude and latitude
              CASE( '-clon', '-clat' )
                 CALL abort('parse_command_line_arguments', message)          
-                !READ(arg, *) lambda_cg
-                !lambda_cg = lambda_cg * TO_RADIANS
-                !READ(arg, *) phi_cg
-                !phi_cg = phi_cg * TO_RADIANS
 
              CASE( '-path', '-p', '--path' )
@@ -444,10 +484,8 @@
                 cfg % namelist_file = TRIM(arg)
 
-             ! initial condition mode: 'profile' / 'volume'
              CASE( '-mode', '-i', '--init-mode' )
                 CALL get_option_argument( i, arg )
                 cfg % ic_mode = TRIM(arg)
 
-             ! boundary conditions / forcing mode: 'ideal' / 'real'
              CASE( '-f', '--forcing-mode' )
                 CALL get_option_argument( i, arg )
@@ -484,4 +522,10 @@
 
    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Get the argument of the i'th command line option, which is at the location
+!> i+1 of the argument list.
+!------------------------------------------------------------------------------!
    SUBROUTINE get_option_argument(i, arg)
       CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
@@ -494,4 +538,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Checks the INIFOR configuration 'cfg' for plausibility.
+!------------------------------------------------------------------------------!
    SUBROUTINE validate_config(cfg)
       TYPE(inifor_config), INTENT(IN) ::  cfg
@@ -503,5 +552,6 @@
       all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file, 'SOILTYP')
 
-      ! Only check optional static driver file name, if it has been given.
+!
+!--   Only check optional static driver file name, if it has been given.
       IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
          all_files_present = all_files_present .AND. file_present(cfg % static_driver_file, 'static driver')
@@ -557,4 +607,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check whether the given file is present on the filesystem.
+!------------------------------------------------------------------------------!
    LOGICAL FUNCTION file_present(filename, kind)
       CHARACTER(LEN=PATH), INTENT(IN) ::  filename
@@ -585,9 +640,9 @@
 ! Description:
 ! ------------
-!> This routine initializes the Inifor output file, i.e. the PALM-4U
-!> initializing and forcing data as a NetCDF file.
+!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
+!> file.
 !> 
 !> Besides writing metadata, such as global attributes, coordinates, variables,
-!> in the NetCDF file, various NetCDF IDs are saved for later, when Inifor
+!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
 !> writes the actual data.
 !------------------------------------------------------------------------------!
@@ -611,5 +666,6 @@
        CALL report('setup_netcdf_dimensions', message)
 
-       ! Create the NetCDF file. NF90_CLOBBER selects overwrite mode.
+!
+!--    Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
 #if defined( __netcdf4 )
        CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
@@ -646,7 +702,8 @@
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
-       ! FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
-       ! FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
-       ! FIXME: Standard v1.9., origin_z)
+!
+!--    FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
+!--    FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
+!--    FIXME: Standard v1.9., origin_z)
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z',       TRIM(real_to_str(z0, '(F18.13)'))))
        CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
@@ -658,57 +715,79 @@
 !- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
 !------------------------------------------------------------------------------
-       dimids = (/0, 0, 0/) ! reset dimids
-          CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
-          CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
-          CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
-          output_file % dimids_scl = dimids ! save dimids for later
-
-       dimvarids = (/0, 0, 0/) ! reset dimvarids
-          CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
-          CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
-          CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
-
-          CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
-          CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
-          CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
-
-          CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
-          CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
-          CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
-       output_file % dimvarids_scl = dimvarids ! save dimvarids for later
-
-       ! overwrite third dimid with the one of depth
+!
+!--    reset dimids first
+       dimids = (/0, 0, 0/)
+       CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
+       CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
+       CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
+!
+!--    save dimids for later
+       output_file % dimids_scl = dimids 
+
+!
+!--    reset dimvarids first
+       dimvarids = (/0, 0, 0/)
+       CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
+       CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
+       CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
+
+       CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
+       CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
+       CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
+
+       CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
+       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
+       CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
+!
+!--    save dimvarids for later
+       output_file % dimvarids_scl = dimvarids
+
+!
+!--    overwrite third dimid with the one of depth
        CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) )
-       output_file % dimids_soil = dimids ! save dimids for later
-
-       ! overwrite third dimvarid with the one of depth
+!
+!--    save dimids for later
+       output_file % dimids_soil = dimids
+
+!
+!--    overwrite third dimvarid with the one of depth
        CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
        CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
        CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
        CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
-       output_file % dimvarids_soil = dimvarids ! save dimvarids for later
+!
+!--    save dimvarids for later
+       output_file % dimvarids_soil = dimvarids
 !
 !------------------------------------------------------------------------------
 !- Section 2b: Define dimensions for cell faces/velocities
 !------------------------------------------------------------------------------
-       dimids = (/0, 0, 0/) ! reset dimids
-          CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
-          CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
-          CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
-       output_file % dimids_vel = dimids ! save dimids for later
-
-       dimvarids = (/0, 0, 0/) ! reset dimvarids
-          CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
-          CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
-          CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
-
-          CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
-          CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
-          CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
-
-          CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
-          CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
-          CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
-       output_file % dimvarids_vel = dimvarids ! save dimvarids for later
+!
+!--    reset dimids first
+       dimids = (/0, 0, 0/)
+       CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
+       CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
+       CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
+!
+!--    save dimids for later
+       output_file % dimids_vel = dimids
+
+!
+!--    reset dimvarids first
+       dimvarids = (/0, 0, 0/)
+       CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
+       CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
+       CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
+
+       CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
+       CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
+       CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
+
+       CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
+       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
+       CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
+!
+!--    save dimvarids for later
+       output_file % dimvarids_vel = dimvarids
 
 !
@@ -739,11 +818,14 @@
        CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(3), palm_grid%zw))
 
-       ! TODO Read in soil depths from input file before this.
+!
+!--    TODO Read in soil depths from input file before this.
        CALL check(nf90_put_var(ncid, output_file % dimvarids_soil(3), palm_grid%depths))
 
-       ! Write time vector
+!
+!--    Write time vector
        CALL check(nf90_put_var(ncid, output_file % dimvarid_time, output_file % time))
 
-       ! Close the file
+!
+!--    Close the file
        CALL check(nf90_close(ncid))
 
@@ -751,4 +833,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Defines the netCDF variables to be written to the dynamic driver file
+!------------------------------------------------------------------------------!
     SUBROUTINE setup_netcdf_variables(filename, output_variable_table, debug)
 
@@ -824,6 +911,8 @@
 !- Section 1: Load input buffers for accumulated variables
 !------------------------------------------------------------------------------
+!
+!--    radiation budgets, precipitation
        IF (group % kind == 'running average' .OR.                              &
-           group % kind == 'accumulated')  THEN ! radiation budgets, precipitation
+           group % kind == 'accumulated')  THEN
 
           IF (SIZE(group % in_var_list) .GT. 1 ) THEN
@@ -835,10 +924,14 @@
           END IF
 
-          ! use two buffer arrays
+!
+!--       use two buffer arrays
           nbuffers = 2
           IF ( .NOT. ALLOCATED( buffer ) )  ALLOCATE( buffer(nbuffers) )
 
-          ! chose correct buffer array
-          hour = iter - 1! hour of the day
+!
+!--       hour of the day
+          hour = iter - 1
+!
+!--       chose correct buffer array
           buf_id = select_buffer(hour)
 
@@ -859,20 +952,23 @@
        ELSE
 
-          ! Allocate one input buffer per input_variable. If more quantities
-          ! have to be computed than input variables exist in this group,
-          ! allocate more buffers. For instance, in the thermodynamics group,
-          ! there are three input variabels (PP, T, Qv) and four quantities
-          ! necessart (P, Theta, Rho, qv) for the corresponding output fields
-          ! (p0, Theta, qv, ug, and vg)
+!
+!--       Allocate one input buffer per input_variable. If more quantities
+!--       have to be computed than input variables exist in this group,
+!--       allocate more buffers. For instance, in the thermodynamics group,
+!--       there are three input variabels (PP, T, Qv) and four quantities
+!--       necessart (P, Theta, Rho, qv) for the corresponding output fields
+!--       (p0, Theta, qv, ug, and vg)
           nbuffers = MAX( group % n_inputs, group % n_output_quantities )
           ALLOCATE( buffer(nbuffers) )
  CALL run_control('time', 'alloc')
           
-          ! Read in all input variables, leave extra buffers-if any-untouched.
+!
+!--       Read in all input variables, leave extra buffers-if any-untouched.
           DO ivar = 1, group % n_inputs
 
              input_var => group % in_var_list(ivar)
 
-             ! Check wheather P or PP is present in input file
+!
+!            Check wheather P or PP is present in input file
              IF (input_var % name == 'P')  THEN
                 input_var % name = TRIM( get_pressure_varname(input_file) )
@@ -891,4 +987,10 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Select the appropriate buffer ID for accumulated COSMO input variables
+!> depending on the current hour.
+!------------------------------------------------------------------------------!
     INTEGER FUNCTION select_buffer(hour)
        INTEGER, INTENT(IN) ::  hour
@@ -943,4 +1045,9 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read the given global attribute form the given netCDF file.
+!------------------------------------------------------------------------------!
     FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
 
@@ -948,5 +1055,5 @@
        REAL(dp)                     ::  attribute_value
 
-       INTEGER                      :: ncid
+       INTEGER                      ::  ncid
 
        IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
@@ -966,4 +1073,11 @@
 
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Updates the dynamic driver with the interpolated field of the current
+!> variable at the current time step.
+!------------------------------------------------------------------------------!
     SUBROUTINE update_output(var, array, iter, output_file, cfg)
        TYPE(nc_var), INTENT(IN)  ::  var
@@ -980,5 +1094,6 @@
        )
 
-       ! Skip time dimension for output
+!
+!--    Skip time dimension for output
        ndim = var % ndim
        IF ( var_is_time_dependent )  ndim = var % ndim - 1
@@ -990,5 +1105,6 @@
        CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
 
-       ! Reduce dimension of output array according to variable kind
+!
+!--    Reduce dimension of output array according to variable kind
        SELECT CASE (TRIM(var % kind))
        
@@ -1086,46 +1202,9 @@
 
 
-    SUBROUTINE write_netcdf_variable_2d(var, iter, output_file, buffer)
-       TYPE(nc_var), INTENT(IN)          ::  var
-       INTEGER, INTENT(IN)               ::  iter
-       TYPE(nc_file), INTENT(IN)         ::  output_file
-       REAL(dp), INTENT(IN)              ::  buffer(:,:,:)
-
-       INTEGER ::  ncid, ndim_out, start(4), count(4)
-       LOGICAL ::  last_dimension_is_time
-
-       ndim_out = var % ndim
-
-       last_dimension_is_time = var % dimids( var % ndim ) == output_file % dimid_time
-       IF ( last_dimension_is_time )  THEN
-          ndim_out = ndim_out - 1
-       END IF
-
-       start(:)      = (/1,1,1,iter/)
-       count(1:ndim_out) = var%dimlen(1:ndim_out)
-
-       CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
-
-       IF (TRIM(var % kind) .EQ. 'left/right scalar')  THEN
-
-          CALL check(nf90_put_var( ncid, var%varid, buffer(1,:,:) ) )
-
-       ELSE IF (TRIM(var % kind) .EQ. 'north/south scalar')  THEN
-
-          CALL check(nf90_put_var( ncid, var%varid, buffer(:,1,:) ) )
-
-       ELSE IF (TRIM(var % kind) .EQ. 'top scalar')  THEN
-
-          CALL check(nf90_put_var( ncid, var%varid, buffer(:,:,1) ) )
-       ELSE
-
-          CALL check(nf90_put_var( ncid, var%varid, buffer ) )
-
-       END IF
-       CALL check(nf90_close(ncid))
-
-    END SUBROUTINE write_netcdf_variable_2d
-
-
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Checks the status of a netCDF API call and aborts if an error occured
+!------------------------------------------------------------------------------!
     SUBROUTINE check(status)
 
