Index: palm/trunk/SOURCE/read_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 4497)
+++ palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 4498)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! argument removed from rd_mpi_io_open
+! 
+! 4496 2020-04-15 08:37:26Z raasch
 ! bugfix: MPI barrier removed, coupling character added to input filename
 ! 
@@ -845,5 +848,5 @@
 !
 !--    Open the MPI-IO restart file.
-       CALL rd_mpi_io_open( 'read', 'BININ' // TRIM( coupling_char ), only_global = .TRUE. )
+       CALL rd_mpi_io_open( 'read', 'BININ' // TRIM( coupling_char ) )
 
 !
Index: palm/trunk/SOURCE/restart_data_mpi_io_mod.f90
===================================================================
--- palm/trunk/SOURCE/restart_data_mpi_io_mod.f90	(revision 4497)
+++ palm/trunk/SOURCE/restart_data_mpi_io_mod.f90	(revision 4498)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! bugfix for creation of filetypes, argument removed from rd_mpi_io_open
+! 
+! 4497 2020-04-15 10:20:51Z raasch
 ! last bugfix deactivated because of compile problems
 ! 
@@ -244,5 +247,5 @@
 !> Open restart file for read or write with MPI-IO
 !--------------------------------------------------------------------------------------------------!
- SUBROUTINE rd_mpi_io_open( action, file_name, only_global )
+ SUBROUTINE rd_mpi_io_open( action, file_name )
 
     IMPLICIT NONE
@@ -251,7 +254,4 @@
     CHARACTER(LEN=*), INTENT(IN)  ::  file_name                        !<
 
-    LOGICAL, INTENT(IN), OPTIONAL ::  only_global                      !<
-    LOGICAL                       ::  set_filetype                     !<
-
     INTEGER(iwp)                  ::  i                                !<
     INTEGER(iwp)                  ::  gh_size                          !<
@@ -276,14 +276,9 @@
 !-- Create subarrays and file types
     filetypes_created = .FALSE.
-    set_filetype      = .TRUE.
-
-    IF ( PRESENT( only_global ) )  THEN
-       IF ( only_global )  set_filetype = .FALSE.
-    ENDIF
 
 !
 !-- In case of read it is not known yet if data include total domain. Filetypes will be created
 !-- further below.
-    IF ( set_filetype  .AND.  wr_flag)  THEN
+    IF ( wr_flag)  THEN
        CALL rs_mpi_io_create_filetypes
        filetypes_created = .TRUE.
@@ -368,9 +363,8 @@
 
 !
-!--    File types deoend on if boundaries of the total domain is included in data
-       IF ( set_filetype )  THEN
-          CALL rs_mpi_io_create_filetypes
-          filetypes_created = .TRUE.
-       ENDIF
+!--    File types depend on if boundaries of the total domain is included in data. This has been
+!--    checked with the previous statement.
+       CALL rs_mpi_io_create_filetypes
+       filetypes_created = .TRUE.
 
 #if defined( __parallel )
Index: palm/trunk/SOURCE/user_data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/user_data_output_mask.f90	(revision 4497)
+++ palm/trunk/SOURCE/user_data_output_mask.f90	(revision 4498)
@@ -1,23 +1,23 @@
 !> @file user_data_output_mask.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
 !
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
 !
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
 !
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
 !
 ! Copyright 1997-2020 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+!
 !
 ! Current revisions:
-! ------------------
+! -----------------
 ! 
 ! 
@@ -25,17 +25,21 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+!
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
-! 
+!
 ! 4168 2019-08-16 13:50:17Z suehring
-! Remove dependency on surface_mod + example for terrain-following output 
+! Remove dependency on surface_mod + example for terrain-following output
 ! adjusted
-! 
+!
 ! 4069 2019-07-01 14:05:51Z Giersch
-! Masked output running index mid has been introduced as a local variable to 
+! Masked output running index mid has been introduced as a local variable to
 ! avoid runtime error (Loop variable has been modified) in time_integration
-! 
+!
 ! 3768 2019-02-27 14:35:58Z raasch
 ! variables commented + statement added to avoid compiler warnings about unused variables
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! Add terrain-following output
@@ -45,34 +49,32 @@
 ! Description:
 ! ------------
-!> Resorts the user-defined output quantity with indices (k,j,i) to a
-!> temporary array with indices (i,j,k) for masked data output.
-!------------------------------------------------------------------------------!
+!> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with
+!> indices (i,j,k) for masked data output.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE user_data_output_mask( av, variable, found, local_pf, mid )
- 
+
 
     USE control_parameters
-        
+
     USE indices
-    
+
     USE kinds
-    
+
     USE user
 
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable  !< 
+    CHARACTER(LEN=*) ::  variable  !<
 
-    INTEGER(iwp) ::  av             !< 
-    INTEGER(iwp) ::  mid            !< masked output running index
-!    INTEGER(iwp) ::  i              !<
-!    INTEGER(iwp) ::  j              !<
-!    INTEGER(iwp) ::  k              !<
-!    INTEGER(iwp) ::  topo_top_index !< k index of highest horizontal surface
+    INTEGER(iwp) ::  av              !<
+    INTEGER(iwp) ::  mid             !< masked output running index
+!    INTEGER(iwp) ::  i               !<
+!    INTEGER(iwp) ::  j               !<
+!    INTEGER(iwp) ::  k               !<
+!    INTEGER(iwp) ::  topo_top_index  !< k index of highest horizontal surface
 
-    LOGICAL ::  found               !< 
+    LOGICAL ::  found  !<
 
-    REAL(wp),                                                                  &
-       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
-          local_pf   !< 
+    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf  !<
 
 !
@@ -87,6 +89,6 @@
 
 !--    Uncomment and extend the following lines, if necessary.
-!--    The arrays for storing the user defined quantities (here u2 and u2_av)
-!--    have to be declared and defined by the user!
+!--    The arrays for storing the user defined quantities (here u2 and u2_av) have to be declared
+!--    and defined by the user!
 !--    Sample for user-defined output:
 !       CASE ( 'u2' )
@@ -98,6 +100,6 @@
 !                   DO  j = 1, mask_size_l(mid,2)
 !                      DO  k = 1, mask_size_l(mid,3)
-!                         local_pf(i,j,k) = u2(mask_k(mid,k),                  &
-!                                              mask_j(mid,j),                  &
+!                         local_pf(i,j,k) = u2(mask_k(mid,k),                                       &
+!                                              mask_j(mid,j),                                       &
 !                                              mask_i(mid,i))
 !                      ENDDO
@@ -111,15 +113,10 @@
 !!
 !!--                   Get k index of highest horizontal surface
-!                      topo_top_index = topo_top_ind( &
-!                                        mask_j(mid,j), &
-!                                        mask_i(mid,i), &
-!                                        1          )
+!                      topo_top_index = topo_top_ind( mask_j(mid,j), mask_i(mid,i), 1 )
 !!
 !!--                   Save output array
 !                      DO  k = 1, mask_size_l(mid,3)
-!                         local_pf(i,j,k) = u2(MIN( topo_top_index+mask_k(mid,k),&
-!                                                   nzt+1 ),                     &
-!                                              mask_j(mid,j),                    &
-!                                              mask_i(mid,i)                   )
+!                         local_pf(i,j,k) = u2(MIN( topo_top_index + mask_k(mid,k), nzt+1 ),        &
+!                                              mask_j(mid,j), mask_i(mid,i) )
 !                      ENDDO
 !                   ENDDO
@@ -133,7 +130,5 @@
 !                   DO  j = 1, mask_size_l(mid,2)
 !                      DO  k = 1, mask_size_l(mid,3)
-!                          local_pf(i,j,k) = u2_av(mask_k(mid,k),              &
-!                                                  mask_j(mid,j),              &
-!                                                  mask_i(mid,i) )
+!                          local_pf(i,j,k) = u2_av(mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )
 !                       ENDDO
 !                    ENDDO
@@ -146,16 +141,10 @@
 !!
 !!--                   Get k index of highest horizontal surface
-!                      topo_top_index = topo_top_ind(   &
-!                                        mask_j(mid,j), &
-!                                        mask_i(mid,i), &
-!                                        1 )
+!                      topo_top_index = topo_top_ind( mask_j(mid,j), mask_i(mid,i), 1 )
 !!
 !!--                   Save output array
 !                      DO  k = 1, mask_size_l(mid,3)
-!                         local_pf(i,j,k) = u2_av(                               &
-!                                              MIN( topo_top_index+mask_k(mid,k),&
-!                                                   nzt+1 ),                     &
-!                                              mask_j(mid,j),                    &
-!                                              mask_i(mid,i)                   )
+!                         local_pf(i,j,k) = u2_av( MIN( topo_top_index+mask_k(mid,k), nzt+1 ),      &
+!                                                  mask_j(mid,j), mask_i(mid,i) )
 !                      ENDDO
 !                   ENDDO
Index: palm/trunk/SOURCE/user_flight.f90
===================================================================
--- palm/trunk/SOURCE/user_flight.f90	(revision 4497)
+++ palm/trunk/SOURCE/user_flight.f90	(revision 4498)
@@ -1,23 +1,23 @@
 !> @file user_flight.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
 !
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
 !
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
 !
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
 !
 ! Copyright 1997-2020 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+!
 !
 ! Current revisions:
-! ------------------
+! -----------------
 ! 
 ! 
@@ -25,12 +25,16 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+!
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
-! 
+!
 ! 3768 2019-02-27 14:35:58Z raasch
 ! unused variables commented out + statement added to avoid compiler warnings
-! 
+!
 ! 3684 2019-01-20 20:20:58Z knoop
 ! Corrected "Former revisions" section
-! 
+!
 ! 1957 2016-07-07 10:43:48Z suehring
 ! Initial revision
@@ -38,11 +42,10 @@
 ! Description:
 ! ------------
-!> Calculation of user-defined output quantity for flight measurements after 
-!> each timestep. 
-!------------------------------------------------------------------------------!
+!> Calculation of user-defined output quantity for flight measurements after each timestep.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE user_flight( var, id )
 
     USE control_parameters
-    
+
     USE grid_variables
 
@@ -60,7 +63,7 @@
 !    INTEGER(iwp) ::  j  !< index along y
 !    INTEGER(iwp) ::  k  !< index along z
-    INTEGER(iwp) ::  id !< variable identifyer, according to the settings in user_init_flight
-        
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !< treated variable
+    INTEGER(iwp) ::  id  !< variable identifyer, according to the settings in user_init_flight
+
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< treated variable
 
 !
@@ -69,31 +72,29 @@
 
 !
-!-- Here, the respective variable is calculated. There is no call of 
-!-- exchange_horiz necessary. 
-!-- The variable identifyer (id) must be set according to the settings in 
-!-- user_init_flight.
+!-- Here, the respective variable is calculated. There is no call of exchange_horiz necessary.
+!-- The variable identifyer (id) must be set according to the settings in user_init_flight.
 !-- Please note, so far, variable must be located at the center of a grid box.
 !     var = 0.0_wp
 
 !     SELECT CASE ( id )
-! 
+!
 !        CASE ( 1 )
-!           DO i = nxl-1, nxr+1
-!              DO j = nys-1, nyn+1
-!                 DO k = nzb, nzt
+!           DO  i = nxl-1, nxr+1
+!              DO  j = nys-1, nyn+1
+!                 DO  k = nzb, nzt
 !                    var(k,j,i) = ABS( u(k,j,i )
 !                 ENDDO
 !              ENDDO
 !           ENDDO
-!           
+!
 !        CASE ( 2 )
-!           DO i = nxl-1, nxr+1
-!              DO j = nys-1, nyn+1
-!                 DO k = nzb, nzt
+!           DO  i = nxl-1, nxr+1
+!              DO  j = nys-1, nyn+1
+!                 DO  k = nzb, nzt
 !                    var(k,j,i) = ABS( v(k,j,i) )
 !                 ENDDO
 !              ENDDO
 !           ENDDO
-! 
+!
 !     END SELECT
 
Index: palm/trunk/SOURCE/user_init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/user_init_3d_model.f90	(revision 4497)
+++ palm/trunk/SOURCE/user_init_3d_model.f90	(revision 4498)
@@ -1,20 +1,20 @@
 !> @file user_init_3d_model.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
 !
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
 !
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
 !
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
 !
 ! Copyright 1997-2020 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+!
 !
 ! Current revisions:
@@ -25,9 +25,13 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+!
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
-! 
+!
 ! 3768 2019-02-27 14:35:58Z raasch
 ! variables commented out to avoid compiler warnings about unused variables
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! Corrected "Former revisions" section
@@ -42,29 +46,29 @@
 !> @attention The user is responsible to set at least all those quantities which
 !>            are normally set within init_3d_model!
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE user_init_3d_model
- 
+
 
     USE arrays_3d
-    
+
     USE control_parameters
-    
+
     USE indices
-    
+
     USE kinds
 
     USE surface_mod
-    
+
     USE user
 
     IMPLICIT NONE
 
-!    INTEGER(iwp) ::  l !< running index surface orientation
-!    INTEGER(iwp) ::  m !< running index surface elements
+!    INTEGER(iwp) ::  l  !< running index surface orientation
+!    INTEGER(iwp) ::  m  !< running index surface elements
 
 !
 !-- Initialization of surface-related quantities.
-!-- The following example shows required initialization of surface quantitites
-!-- at default-type upward-facing surfaces.  
+!-- The following example shows required initialization of surface quantitites at default-type
+!-- upward-facing surfaces.
 !   DO  m = 1, surf_def_h(0)%ns
 !      surf_def_h(0)%ol(m)   = ...    ! Obukhov length
@@ -87,15 +91,15 @@
 !         surf_def_h(0)%ssws(m) = ... ! surface latent heat flux
 !      ENDIF
-!   ENDDO 
+!   ENDDO
 !
 !-- Same for natural and urban type surfaces
 !   DO  m = 1, surf_lsm_h%ns
 !      ...
-!   ENDDO 
+!   ENDDO
 !   DO  m = 1, surf_usm_h%ns
 !      ...
 !   ENDDO
 !
-!-- Also care for vertically aligned surfaces (default-, natural-, and 
+!-- Also care for vertically aligned surfaces (default-, natural-, and
 !-- urban-type).
 !   DO  l = 0, 3
Index: palm/trunk/SOURCE/user_init_flight.f90
===================================================================
--- palm/trunk/SOURCE/user_init_flight.f90	(revision 4497)
+++ palm/trunk/SOURCE/user_init_flight.f90	(revision 4498)
@@ -1,20 +1,20 @@
 !> @file user_init_flight.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
 !
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
 !
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
 !
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
 !
 ! Copyright 1997-2020 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+!
 !
 ! Current revisions:
@@ -25,12 +25,16 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+!
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
-! 
+!
 ! 3768 2019-02-27 14:35:58Z raasch
 ! statements commented or added to avoid compiler warnings about unused variables
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! Corrected "Former revisions" section
-! 
+!
 ! 1957 2016-07-07 10:43:48Z suehring
 ! Initial revision
@@ -39,27 +43,28 @@
 ! ------------
 !> Execution of user-defined initialization for flight measurements.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE user_init_flight( init, k, id, label_leg )
- 
+
 
     USE control_parameters
-    
+
     USE indices
-    
+
     USE kinds
-    
-!    USE netcdf_interface,                                                      &
-!        ONLY: dofl_label, dofl_unit
-    
+
+!    USE netcdf_interface,                                                                          &
+!        ONLY: dofl_label,                                                                          &
+!              dofl_unit
+
     USE user
 
     IMPLICIT NONE
-    
+
     CHARACTER(LEN=10), OPTIONAL ::  label_leg     !< label of the respective leg
-   
+
     INTEGER(iwp), OPTIONAL                ::  id  !< variable index
     INTEGER(iwp), OPTIONAL, INTENT(INOUT) ::  k   !< index for respective variable and leg
-    
-    LOGICAL ::  init  !< variable to recognize initial call 
+
+    LOGICAL ::  init  !< variable to recognize initial call
 
 !
@@ -71,7 +76,6 @@
 !
 !-- Sample for user-defined flight-time series.
-!-- For each quantity you have to give a label and a unit, which will be used
-!-- for the output into NetCDF file. They must not contain more than
-!-- twenty characters. 
+!-- For each quantity you have to give a label and a unit, which will be used for the output into
+!-- NetCDF file. They must not contain more than twenty characters.
 
 
@@ -80,30 +84,30 @@
 !--    The number of user-defined quantity has to be increased appropriately.
 !--    In the following example, 2 user-defined quantities are added.
-!        num_var_fl_user = num_var_fl_user + 2 
+!        num_var_fl_user = num_var_fl_user + 2
 
        init = .FALSE.
-    
+
     ELSE
-    
+
 !
 !--    Please add the respective number of new variables as following:
-      
+
 !        SELECT CASE ( id )
-!        
+!
 !           CASE ( 1 )
 !              dofl_label(k)   = TRIM(label_leg) // '_' // 'abs_u'
 !              dofl_unit(k)    = 'm/s'
 !              k               = k + 1
-!              
+!
 !           CASE ( 2 )
-!     
+!
 !              dofl_label(k)   = TRIM(label_leg) // '_' // 'abs_v'
 !              dofl_unit(k)    = 'm/s'
 !              k               = k + 1
-!              
+!
 !        END SELECT
 
     ENDIF
-       
+
  END SUBROUTINE user_init_flight
 
Index: palm/trunk/SOURCE/user_init_grid.f90
===================================================================
--- palm/trunk/SOURCE/user_init_grid.f90	(revision 4497)
+++ palm/trunk/SOURCE/user_init_grid.f90	(revision 4498)
@@ -1,20 +1,20 @@
 !> @file user_init_grid.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
 !
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
 !
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
 !
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
 !
 ! Copyright 1997-2020 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+!
 !
 ! Current revisions:
@@ -25,9 +25,13 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+!
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
-! 
+!
 ! 3768 2019-02-27 14:35:58Z raasch
 ! variables commented + statement added to avoid compiler warnings about unused variables
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! dz was replaced by dz(1)
@@ -40,20 +44,20 @@
 ! ------------
 !> Execution of user-defined grid initializing actions
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE user_init_grid( topo_3d )
- 
+
 
     USE control_parameters
-    
+
     USE indices
-    
+
     USE kinds
-    
+
     USE user
 
     IMPLICIT NONE
 
-!    INTEGER(iwp)                                           ::  k_topo      !< topography top index
-    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d     !< 3D topography field
+!    INTEGER(iwp)                                           ::  k_topo   !< topography top index
+    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d  !< 3D topography field
 
 !    REAL(wp) ::  h_topo !< user-defined topography height
@@ -78,31 +82,28 @@
        CASE ( 'user_defined_topography_1' )
 !
-!--       Here the user can define his own topography.
+!--       Here the user can define their own topography.
 !--       After definition, please remove the following three lines!
           message_string = 'topography "' // topography // '" not available yet'
           CALL message( 'user_init_grid', 'UI0005', 1, 2, 0, 6, 0 )
 !
-!--       The user is allowed to set surface-mounted as well as non-surface
-!--       mounted topography (e.g. overhanging structures). For both, use 
-!--       3D array topo_3d and set bit 0. The convention is: bit is zero inside 
-!--       topography, bit is 1 for atmospheric grid point. 
-!--       The following example shows how to prescribe sine-like topography 
-!--       along x-direction with amplitude of 10 * dz(1) and wavelength 10 * dy.
+!--       The user is allowed to set surface-mounted as well as non-surface mounted topography
+!--       (e.g. overhanging structures). For both, use 3D array topo_3d and set bit 0. The
+!--       convention is: bit is zero inside topography, bit is 1 for atmospheric grid point.
+!--       The following example shows how to prescribe sine-like topography along x-direction with
+!--       amplitude of 10 * dz(1) and wavelength 10 * dy.
 !           DO  i = nxlg, nxrg
-!              h_topo = 10.0_wp * dz(1) * (SIN(3.14_wp*0.5_wp)*i*dx / ( 5.0_wp * dy ) )**2
-! 
+!              h_topo = 10.0_wp * dz(1) * ( SIN( 3.14_wp * 0.5_wp) * i * dx / ( 5.0_wp * dy ) )**2
+!
 !              k_topo = MINLOC( ABS( zw - h_topo ), 1 ) - 1
-! 
-!              topo_3d(k_topo+1:nzt+1,:,i) =                                     &
-!                                          IBSET( topo_3d(k_topo+1:nzt+1,:,i), 0 ) 
-!           ENDDO 
-! 
+!
+!              topo_3d(k_topo+1:nzt+1,:,i) = IBSET( topo_3d(k_topo+1:nzt+1,:,i), 0 )
+!           ENDDO
+!
 !           CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
 
        CASE DEFAULT
 !
-!--       The DEFAULT case is reached if the parameter topography contains a
-!--       wrong character string that is neither recognized in init_grid nor 
-!--       here in user_init_grid.
+!--       The DEFAULT case is reached if the parameter topography contains a wrong character string
+!--       that is neither recognized in init_grid nor here in user_init_grid.
           message_string = 'unknown topography "' // topography // '"'
           CALL message( 'user_init_grid', 'UI0006', 1, 2, 0, 6, 0 )
@@ -110,7 +111,4 @@
     END SELECT
 
-
-
-
  END SUBROUTINE user_init_grid
 
Index: palm/trunk/SOURCE/virtual_measurement_mod.f90
===================================================================
--- palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 4497)
+++ palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 4498)
@@ -1,20 +1,20 @@
 !> @virtual_measurement_mod.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
 !
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017-2020 iz Universitaet Hannover
-!------------------------------------------------------------------------------!
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2020 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------------------------!
+!
 !
 ! Current revisions:
@@ -25,17 +25,21 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+!
+! 4481 2020-03-31 18:55:54Z maronga
 ! bugfix: cpp-directives for serial mode added
-! 
+!
 ! 4438 2020-03-03 20:49:28Z suehring
 ! Add cpu-log points
-! 
+!
 ! 4422 2020-02-24 22:45:13Z suehring
 ! Missing trim()
-! 
+!
 ! 4408 2020-02-14 10:04:39Z gronemeier
 ! - Output of character string station_name after DOM has been enabled to
 !   output character variables
 ! - Bugfix, missing coupling_char statement when opening the input file
-! 
+!
 ! 4408 2020-02-14 10:04:39Z gronemeier
 ! write fill_value attribute
@@ -116,126 +120,124 @@
 ! Description:
 ! ------------
-!> The module acts as an interface between 'real-world' observations and
-!> model simulations. Virtual measurements will be taken in the model at the
-!> coordinates representative for the 'real-world' observation coordinates.
-!> More precisely, coordinates and measured quanties will be read from a
-!> NetCDF file which contains all required information. In the model,
-!> the same quantities (as long as all the required components are switched-on)
-!> will be sampled at the respective positions and output into an extra file,
-!> which allows for straight-forward comparison of model results with
-!> observations.
-!------------------------------------------------------------------------------!
+!> The module acts as an interface between 'real-world' observations and model simulations.
+!> Virtual measurements will be taken in the model at the coordinates representative for the
+!> 'real-world' observation coordinates. More precisely, coordinates and measured quanties will be
+!> read from a NetCDF file which contains all required information. In the model, the same
+!> quantities (as long as all the required components are switched-on) will be sampled at the
+!> respective positions and output into an extra file, which allows for straight-forward comparison
+!> of model results with observations.
+!--------------------------------------------------------------------------------------------------!
  MODULE virtual_measurement_mod
 
-    USE arrays_3d,                                                             &
-        ONLY:  dzw,                                                            &
-               exner,                                                          &
-               hyp,                                                            &
-               q,                                                              &
-               ql,                                                             &
-               pt,                                                             &
-               rho_air,                                                        &
-               u,                                                              &
-               v,                                                              &
-               w,                                                              &
-               zu,                                                             &
+    USE arrays_3d,                                                                                 &
+        ONLY:  dzw,                                                                                &
+               exner,                                                                              &
+               hyp,                                                                                &
+               q,                                                                                  &
+               ql,                                                                                 &
+               pt,                                                                                 &
+               rho_air,                                                                            &
+               u,                                                                                  &
+               v,                                                                                  &
+               w,                                                                                  &
+               zu,                                                                                 &
                zw
 
-    USE basic_constants_and_equations_mod,                                     &
-        ONLY:  convert_utm_to_geographic,                                      &
-               degc_to_k,                                                      &
-               magnus,                                                         &
-               pi,                                                             &
+    USE basic_constants_and_equations_mod,                                                         &
+        ONLY:  convert_utm_to_geographic,                                                          &
+               degc_to_k,                                                                          &
+               magnus,                                                                             &
+               pi,                                                                                 &
                rd_d_rv
 
-    USE chem_gasphase_mod,                                                     &
+    USE chem_gasphase_mod,                                                                         &
         ONLY:  nvar
 
-    USE chem_modules,                                                          &
+    USE chem_modules,                                                                              &
         ONLY:  chem_species
 
-    USE control_parameters,                                                    &
-        ONLY:  air_chemistry,                                                  &
-               coupling_char,                                                  &
-               dz,                                                             &
-               end_time,                                                       &
-               humidity,                                                       &
-               message_string,                                                 &
-               neutral,                                                        &
-               origin_date_time,                                               &
-               rho_surface,                                                    &
-               surface_pressure,                                               &
-               time_since_reference_point,                                     &
+    USE control_parameters,                                                                        &
+        ONLY:  air_chemistry,                                                                      &
+               coupling_char,                                                                      &
+               dz,                                                                                 &
+               end_time,                                                                           &
+               humidity,                                                                           &
+               message_string,                                                                     &
+               neutral,                                                                            &
+               origin_date_time,                                                                   &
+               rho_surface,                                                                        &
+               surface_pressure,                                                                   &
+               time_since_reference_point,                                                         &
                virtual_measurement
 
-    USE cpulog,                                                                &
-        ONLY:  cpu_log,                                                        &
+    USE cpulog,                                                                                    &
+        ONLY:  cpu_log,                                                                            &
                log_point_s
 
     USE data_output_module
 
-    USE grid_variables,                                                        &
-        ONLY:  ddx,                                                            &
-               ddy,                                                            &
-               dx,                                                             &
+    USE grid_variables,                                                                            &
+        ONLY:  ddx,                                                                                &
+               ddy,                                                                                &
+               dx,                                                                                 &
                dy
 
-    USE indices,                                                               &
-        ONLY:  nbgp,                                                           &
-               nzb,                                                            &
-               nzt,                                                            &
-               nxl,                                                            &
-               nxlg,                                                           &
-               nxr,                                                            &
-               nxrg,                                                           &
-               nys,                                                            &
-               nysg,                                                           &
-               nyn,                                                            &
-               nyng,                                                           &
-               topo_top_ind,                                                   &
+    USE indices,                                                                                   &
+        ONLY:  nbgp,                                                                               &
+               nzb,                                                                                &
+               nzt,                                                                                &
+               nxl,                                                                                &
+               nxlg,                                                                               &
+               nxr,                                                                                &
+               nxrg,                                                                               &
+               nys,                                                                                &
+               nysg,                                                                               &
+               nyn,                                                                                &
+               nyng,                                                                               &
+               topo_top_ind,                                                                       &
                wall_flags_total_0
 
     USE kinds
 
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  close_input_file,                                               &
-               coord_ref_sys,                                                  &
-               crs_list,                                                       &
-               get_attribute,                                                  &
-               get_dimension_length,                                           &
-               get_variable,                                                   &
-               init_model,                                                     &
-               input_file_atts,                                                &
-               input_file_vm,                                                  &
-               input_pids_static,                                              &
-               input_pids_vm,                                                  &
-               inquire_fill_value,                                             &
-               open_read_file,                                                 &
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  close_input_file,                                                                   &
+               coord_ref_sys,                                                                      &
+               crs_list,                                                                           &
+               get_attribute,                                                                      &
+               get_dimension_length,                                                               &
+               get_variable,                                                                       &
+               init_model,                                                                         &
+               input_file_atts,                                                                    &
+               input_file_vm,                                                                      &
+               input_pids_static,                                                                  &
+               input_pids_vm,                                                                      &
+               inquire_fill_value,                                                                 &
+               open_read_file,                                                                     &
                pids_id
 
     USE pegrid
 
-    USE surface_mod,                                                           &
-        ONLY:  surf_lsm_h,                                                     &
+    USE surface_mod,                                                                               &
+        ONLY:  surf_lsm_h,                                                                         &
                surf_usm_h
 
-    USE land_surface_model_mod,                                                &
-        ONLY:  m_soil_h,                                                       &
-               nzb_soil,                                                       &
-               nzt_soil,                                                       &
-               t_soil_h,                                                       &
+    USE land_surface_model_mod,                                                                    &
+        ONLY:  m_soil_h,                                                                           &
+               nzb_soil,                                                                           &
+               nzt_soil,                                                                           &
+               t_soil_h,                                                                           &
                zs
 
-    USE radiation_model_mod,                                                   &
-        ONLY:  rad_lw_in,                                                      &
-               rad_lw_out,                                                     &
-               rad_sw_in,                                                      &
-               rad_sw_in_diff,                                                 &
-               rad_sw_out,                                                     &
+    USE radiation_model_mod,                                                                       &
+        ONLY:  rad_lw_in,                                                                          &
+               rad_lw_out,                                                                         &
+               rad_sw_in,                                                                          &
+               rad_sw_in_diff,                                                                     &
+               rad_sw_out,                                                                         &
                radiation_scheme
 
-    USE urban_surface_mod,                                                     &
-        ONLY:  nzb_wall,                                                       &
-               nzt_wall,                                                       &
+    USE urban_surface_mod,                                                                         &
+        ONLY:  nzb_wall,                                                                           &
+               nzt_wall,                                                                           &
                t_wall_h
 
@@ -244,42 +246,41 @@
 
     TYPE virt_general
-       INTEGER(iwp) ::  nvm = 0   !< number of virtual measurements
+       INTEGER(iwp) ::  nvm = 0  !< number of virtual measurements
     END TYPE virt_general
 
     TYPE virt_var_atts
-       CHARACTER(LEN=100) ::  coordinates          !< defined longname of the variable
-       CHARACTER(LEN=100) ::  grid_mapping         !< defined longname of the variable
-       CHARACTER(LEN=100) ::  long_name            !< defined longname of the variable
-       CHARACTER(LEN=100) ::  name                 !< variable name
-       CHARACTER(LEN=100) ::  standard_name        !< defined standard name of the variable
-       CHARACTER(LEN=100) ::  units                !< unit of the output variable
-
-       REAL(wp)           ::  fill_value = -9999.0 !< _FillValue attribute
+       CHARACTER(LEN=100) ::  coordinates           !< defined longname of the variable
+       CHARACTER(LEN=100) ::  grid_mapping          !< defined longname of the variable
+       CHARACTER(LEN=100) ::  long_name             !< defined longname of the variable
+       CHARACTER(LEN=100) ::  name                  !< variable name
+       CHARACTER(LEN=100) ::  standard_name         !< defined standard name of the variable
+       CHARACTER(LEN=100) ::  units                 !< unit of the output variable
+
+       REAL(wp)           ::  fill_value = -9999.0  !< _FillValue attribute
     END TYPE virt_var_atts
 
     TYPE virt_mea
-
-       CHARACTER(LEN=100)  ::  feature_type                     !< type of the real-world measurement
-       CHARACTER(LEN=100)  ::  feature_type_out = 'timeSeries'  !< type of the virtual measurement
-                                                                !< (all will be timeSeries, even trajectories)
-       CHARACTER(LEN=100)  ::  nc_filename                      !< name of the NetCDF output file for the station
-       CHARACTER(LEN=100)  ::  site                             !< name of the measurement site
-
-       CHARACTER(LEN=1000) ::  data_content = REPEAT(' ', 1000) !< string of measured variables (data output only)
-
-       INTEGER(iwp) ::  end_coord_a = 0     !< end coordinate in NetCDF file for local atmosphere observations
-       INTEGER(iwp) ::  end_coord_s = 0     !< end coordinate in NetCDF file for local soil observations
-       INTEGER(iwp) ::  file_time_index = 0 !< time index in NetCDF output file
-       INTEGER(iwp) ::  ns = 0              !< number of observation coordinates on subdomain, for atmospheric measurements
-       INTEGER(iwp) ::  ns_tot = 0          !< total number of observation coordinates, for atmospheric measurements
-       INTEGER(iwp) ::  n_tr_st             !< number of trajectories / station of a measurement
-       INTEGER(iwp) ::  nmeas               !< number of measured variables (atmosphere + soil)
-       INTEGER(iwp) ::  ns_soil = 0         !< number of observation coordinates on subdomain, for soil measurements
-       INTEGER(iwp) ::  ns_soil_tot = 0     !< total number of observation coordinates, for soil measurements
-       INTEGER(iwp) ::  start_coord_a = 0   !< start coordinate in NetCDF file for local atmosphere observations
-       INTEGER(iwp) ::  start_coord_s = 0   !< start coordinate in NetCDF file for local soil observations
-
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory
-                                                         !< or station that are no _FillValues
+       CHARACTER(LEN=100)  ::  feature_type                      !< type of the real-world measurement
+       CHARACTER(LEN=100)  ::  feature_type_out = 'timeSeries'   !< type of the virtual measurement
+                                                                 !< (all will be timeSeries, even trajectories)
+       CHARACTER(LEN=100)  ::  nc_filename                       !< name of the NetCDF output file for the station
+       CHARACTER(LEN=100)  ::  site                              !< name of the measurement site
+
+       CHARACTER(LEN=1000) ::  data_content = REPEAT(' ', 1000)  !< string of measured variables (data output only)
+
+       INTEGER(iwp) ::  end_coord_a     = 0  !< end coordinate in NetCDF file for local atmosphere observations
+       INTEGER(iwp) ::  end_coord_s     = 0  !< end coordinate in NetCDF file for local soil observations
+       INTEGER(iwp) ::  file_time_index = 0  !< time index in NetCDF output file
+       INTEGER(iwp) ::  ns              = 0  !< number of observation coordinates on subdomain, for atmospheric measurements
+       INTEGER(iwp) ::  ns_tot          = 0  !< total number of observation coordinates, for atmospheric measurements
+       INTEGER(iwp) ::  n_tr_st              !< number of trajectories / station of a measurement
+       INTEGER(iwp) ::  nmeas                !< number of measured variables (atmosphere + soil)
+       INTEGER(iwp) ::  ns_soil         = 0  !< number of observation coordinates on subdomain, for soil measurements
+       INTEGER(iwp) ::  ns_soil_tot     = 0  !< total number of observation coordinates, for soil measurements
+       INTEGER(iwp) ::  start_coord_a   = 0  !< start coordinate in NetCDF file for local atmosphere observations
+       INTEGER(iwp) ::  start_coord_s   = 0  !< start coordinate in NetCDF file for local soil observations
+
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t  !< number observations individual for each trajectory
+                                                          !< or station that are no _FillValues
 
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< grid index for measurement position in x-direction
@@ -291,19 +292,18 @@
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k_soil  !< grid index for measurement position in k-direction
 
-       LOGICAL ::  trajectory         = .FALSE. !< flag indicating that the observation is a mobile observation
-       LOGICAL ::  timseries          = .FALSE. !< flag indicating that the observation is a stationary point measurement
-       LOGICAL ::  timseries_profile  = .FALSE. !< flag indicating that the observation is a stationary profile measurement
-       LOGICAL ::  soil_sampling      = .FALSE. !< flag indicating that soil state variables were sampled
-
-       REAL(wp) ::  fill_eutm                            !< fill value for UTM coordinates in case of missing values
-       REAL(wp) ::  fill_nutm                            !< fill value for UTM coordinates in case of missing values
-       REAL(wp) ::  fill_zar                             !< fill value for heigth coordinates in case of missing values
-       REAL(wp) ::  fillout = -9999.0                    !< fill value for output in case a observation is taken
-                                                         !< e.g. from inside a building
-       REAL(wp) ::  origin_x_obs                         !< origin of the observation in UTM coordiates in x-direction
-       REAL(wp) ::  origin_y_obs                         !< origin of the observation in UTM coordiates in y-direction
-
+       LOGICAL ::  soil_sampling      = .FALSE.  !< flag indicating that soil state variables were sampled
+       LOGICAL ::  trajectory         = .FALSE.  !< flag indicating that the observation is a mobile observation
+       LOGICAL ::  timseries          = .FALSE.  !< flag indicating that the observation is a stationary point measurement
+       LOGICAL ::  timseries_profile  = .FALSE.  !< flag indicating that the observation is a stationary profile measurement
+
+       REAL(wp) ::  fill_eutm          !< fill value for UTM coordinates in case of missing values
+       REAL(wp) ::  fill_nutm          !< fill value for UTM coordinates in case of missing values
+       REAL(wp) ::  fill_zar           !< fill value for heigth coordinates in case of missing values
+       REAL(wp) ::  fillout = -9999.0  !< fill value for output in case an observation is taken e.g. from inside a building
+       REAL(wp) ::  origin_x_obs       !< origin of the observation in UTM coordiates in x-direction
+       REAL(wp) ::  origin_y_obs       !< origin of the observation in UTM coordiates in y-direction
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth         !< measurement depth in soil
        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zar           !< measurement height above ground level
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth         !< measurement depth in soil
 
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars       !< measured variables
@@ -311,16 +311,12 @@
 
        TYPE( virt_var_atts ), DIMENSION(:), ALLOCATABLE ::  var_atts !< variable attributes
-
     END TYPE virt_mea
 
-    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"                      !< dimension name for UTM coordinate easting
-    CHARACTER(LEN=11) ::  char_feature = "featureType"             !< attribute name for feature type
+    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"            !< dimension name for UTM coordinate easting
+    CHARACTER(LEN=11) ::  char_feature = "featureType"   !< attribute name for feature type
 
     ! This need to be generalized
     CHARACTER(LEN=10) ::  char_fill = '_FillValue'                 !< attribute name for fill value
     CHARACTER(LEN=9)  ::  char_long = 'long_name'                  !< attribute name for long_name
-    CHARACTER(LEN=13) ::  char_standard = 'standard_name'          !< attribute name for standard_name
-    CHARACTER(LEN=5)  ::  char_unit = 'units'                      !< attribute name for standard_name
-    CHARACTER(LEN=11) ::  char_soil = "soil_sample"                !< attribute name for soil sampling indication
     CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
     CHARACTER(LEN=5)  ::  char_nutm = "N_UTM"                      !< dimension name for UTM coordinate northing
@@ -329,5 +325,8 @@
     CHARACTER(LEN=8)  ::  char_origy = "origin_y"                  !< attribute name for station coordinate in y
     CHARACTER(LEN=4)  ::  char_site = "site"                       !< attribute name for site name
+    CHARACTER(LEN=11) ::  char_soil = "soil_sample"                !< attribute name for soil sampling indication
+    CHARACTER(LEN=13) ::  char_standard = 'standard_name'          !< attribute name for standard_name
     CHARACTER(LEN=9)  ::  char_station_h = "station_h"             !< variable name indicating height of the site
+    CHARACTER(LEN=5)  ::  char_unit = 'units'                      !< attribute name for standard_name
     CHARACTER(LEN=1)  ::  char_zar = "z"                           !< attribute name indicating height above reference level
     CHARACTER(LEN=10) ::  type_ts   = 'timeSeries'                 !< name of stationary point measurements
@@ -335,41 +334,38 @@
     CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
 
-    CHARACTER(LEN=6), DIMENSION(1:5) ::  soil_vars       = (/                  & !< list of soil variables
-                            't_soil',                                          &
-                            'm_soil',                                          &
-                            'lwc   ',                                          &
-                            'lwcs  ',                                          &
-                            'smp   '                       /)
-
-    CHARACTER(LEN=10), DIMENSION(0:1,1:8) ::  chem_vars = RESHAPE( (/          &
-                                              'mcpm1     ', 'PM1       ',      &
-                                              'mcpm2p5   ', 'PM2.5     ',      &
-                                              'mcpm10    ', 'PM10      ',      &
-                                              'mfno2     ', 'NO2       ',      &
-                                              'mfno      ', 'NO        ',      &
-                                              'mcno2     ', 'NO2       ',      &
-                                              'mcno      ', 'NO        ',      &
-                                              'tro3      ', 'O3        '       &
-                                                                   /), (/ 2, 8 /) )
-
-    LOGICAL ::  global_attribute = .TRUE.           !< flag indicating a global attribute
-    LOGICAL ::  initial_write_coordinates = .FALSE. !< flag indicating a global attribute
-    LOGICAL ::  use_virtual_measurement = .FALSE.   !< Namelist parameter
-
-    INTEGER(iwp) ::  maximum_name_length = 32 !< maximum name length of station names
-    INTEGER(iwp) ::  ntimesteps               !< number of timesteps defined in NetCDF output file
-    INTEGER(iwp) ::  off_pr = 1               !< number neighboring grid points (in each direction) where virtual profile
-                                              !< measurements shall be taken, in addition to the given coordinates in the driver
-    INTEGER(iwp) ::  off_ts = 1               !< number neighboring grid points (in each direction) where virtual timeseries
-                                              !< measurements shall be taken, in addition to the given coordinates in the driver
-    INTEGER(iwp) ::  off_tr = 1               !< number neighboring grid points (in each direction) where virtual trajectory
-                                              !< measurements shall be taken, in addition to the given coordinates in the driver
-
-    REAL(wp) ::  dt_virtual_measurement = 0.0_wp    !< sampling interval
+    CHARACTER(LEN=6), DIMENSION(1:5) ::  soil_vars = (/ 't_soil', & !< list of soil variables
+                                                        'm_soil',                                  &
+                                                        'lwc   ',                                  &
+                                                        'lwcs  ',                                  &
+                                                        'smp   ' /)
+
+    CHARACTER(LEN=10), DIMENSION(0:1,1:8) ::  chem_vars = RESHAPE( (/ 'mcpm1     ', 'PM1       ',  &
+                                                                      'mcpm2p5   ', 'PM2.5     ',  &
+                                                                      'mcpm10    ', 'PM10      ',  &
+                                                                      'mfno2     ', 'NO2       ',  &
+                                                                      'mfno      ', 'NO        ',  &
+                                                                      'mcno2     ', 'NO2       ',  &
+                                                                      'mcno      ', 'NO        ',  &
+                                                                      'tro3      ', 'O3        '   &
+                                                                    /), (/ 2, 8 /) )
+
+    INTEGER(iwp) ::  maximum_name_length = 32  !< maximum name length of station names
+    INTEGER(iwp) ::  ntimesteps                !< number of timesteps defined in NetCDF output file
+    INTEGER(iwp) ::  off_pr              = 1   !< number of neighboring grid points (in each direction) where virtual profile
+                                               !< measurements shall be taken, in addition to the given coordinates in the driver
+    INTEGER(iwp) ::  off_ts              = 1   !< number of neighboring grid points (in each direction) where virtual timeseries
+                                               !< measurements shall be taken, in addition to the given coordinates in the driver
+    INTEGER(iwp) ::  off_tr              = 1   !< number of neighboring grid points (in each direction) where virtual trajectory
+                                               !< measurements shall be taken, in addition to the given coordinates in the driver
+    LOGICAL ::  global_attribute          = .TRUE.   !< flag indicating a global attribute
+    LOGICAL ::  initial_write_coordinates = .FALSE.  !< flag indicating a global attribute
+    LOGICAL ::  use_virtual_measurement   = .FALSE.  !< Namelist parameter
+
+    REAL(wp) ::  dt_virtual_measurement   = 0.0_wp  !< sampling interval
     REAL(wp) ::  time_virtual_measurement = 0.0_wp  !< time since last sampling
-    REAL(wp) ::  vm_time_start = 0.0                !< time after which sampling shall start
-
-    TYPE( virt_general )                        ::  vmea_general !< data structure which encompass global variables
-    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea         !< data structure contain station-specific variables
+    REAL(wp) ::  vm_time_start            = 0.0     !< time after which sampling shall start
+
+    TYPE( virt_general )                        ::  vmea_general  !< data structure which encompasses global variables
+    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea          !< data structure containing station-specific variables
 
     INTERFACE vm_check_parameters
@@ -403,17 +399,17 @@
 !
 !-- Public interfaces
-    PUBLIC  vm_check_parameters,                                               &
-            vm_data_output,                                                    &
-            vm_init,                                                           &
-            vm_init_output,                                                    &
-            vm_parin,                                                          &
+    PUBLIC  vm_check_parameters,                                                                   &
+            vm_data_output,                                                                        &
+            vm_init,                                                                               &
+            vm_init_output,                                                                        &
+            vm_parin,                                                                              &
             vm_sampling
 
 !
 !-- Public variables
-    PUBLIC  dt_virtual_measurement,                                            &
-            time_virtual_measurement,                                          &
-            vmea,                                                              &
-            vmea_general,                                                      &
+    PUBLIC  dt_virtual_measurement,                                                                &
+            time_virtual_measurement,                                                              &
+            vmea,                                                                                  &
+            vmea_general,                                                                          &
             vm_time_start
 
@@ -421,9 +417,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check parameters for virtual measurement module
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE vm_check_parameters
 
@@ -432,5 +428,5 @@
 !-- Virtual measurements require a setup file.
     IF ( .NOT. input_pids_vm )  THEN
-       message_string = 'If virtual measurements are taken, a setup input ' // &
+       message_string = 'If virtual measurements are taken, a setup input ' //                     &
                         'file for the site locations is mandatory.'
        CALL message( 'vm_check_parameters', 'PA0533', 1, 2, 0, 6, 0 )
@@ -438,10 +434,9 @@
 !
 !-- In case virtual measurements are taken, a static input file is required.
-!-- This is because UTM coordinates for the PALM domain origin are required
-!-- for correct mapping of the measurements.
+!-- This is because UTM coordinates for the PALM domain origin are required for correct mapping of
+!-- the measurements.
 !-- ToDo: Revise this later and remove this requirement.
     IF ( .NOT. input_pids_static )  THEN
-       message_string = 'If virtual measurements are taken, a static input ' //&
-                        'file is mandatory.'
+       message_string = 'If virtual measurements are taken, a static input file is mandatory.'
        CALL message( 'vm_check_parameters', 'PA0534', 1, 2, 0, 6, 0 )
     ENDIF
@@ -451,17 +446,15 @@
 !-- In case of non-parallel NetCDF the virtual measurement output is not
 !-- working. This is only designed for parallel NetCDF.
-    message_string = 'If virtual measurements are taken, parallel ' //         &
-                     'NetCDF is required.'
+    message_string = 'If virtual measurements are taken, parallel NetCDF is required.'
     CALL message( 'vm_check_parameters', 'PA0708', 1, 2, 0, 6, 0 )
 #endif
 !
-!-- Check if the given number of neighboring grid points do not exceeds the number
+!-- Check if the given number of neighboring grid points do not exceed the number
 !-- of ghost points.
-    IF ( off_pr > nbgp - 1  .OR.  off_ts > nbgp - 1  .OR.  off_tr > nbgp - 1 ) &
-    THEN
-       WRITE(message_string,*)                                                 &
-                        'If virtual measurements are taken, the number ' //    &
-                        'of surrounding grid points must not be larger ' //    &
-                        'than the number of ghost points - 1, which is: ',     &
+    IF ( off_pr > nbgp - 1  .OR.  off_ts > nbgp - 1  .OR.  off_tr > nbgp - 1 )  THEN
+       WRITE(message_string,*)                                                                     &
+                        'If virtual measurements are taken, the number ' //                        &
+                        'of surrounding grid points must not be larger ' //                        &
+                        'than the number of ghost points - 1, which is: ',                         &
                         nbgp - 1
        CALL message( 'vm_check_parameters', 'PA0705', 1, 2, 0, 6, 0 )
@@ -475,332 +468,331 @@
  END SUBROUTINE vm_check_parameters
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Subroutine defines variable attributes according to UC2 standard. Note, later
-!> this list can be moved to the data-output module where it can be re-used also
-!> for other output.
-!------------------------------------------------------------------------------!
-  SUBROUTINE vm_set_attributes( output_variable )
-
-     TYPE( virt_var_atts ), INTENT(INOUT) ::  output_variable !< data structure with attributes that need to be set
-
-     output_variable%long_name     = 'none'
-     output_variable%standard_name = 'none'
-     output_variable%units         = 'none'
-     output_variable%coordinates   = 'lon lat E_UTM N_UTM x y z time station_name'
-     output_variable%grid_mapping  = 'crs'
-
-     SELECT CASE ( TRIM( output_variable%name ) )
-
-        CASE ( 'u' )
-           output_variable%long_name     = 'u wind component'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'ua' )
-           output_variable%long_name     = 'eastward wind'
-           output_variable%standard_name = 'eastward_wind'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'v' )
-           output_variable%long_name     = 'v wind component'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'va' )
-           output_variable%long_name     = 'northward wind'
-           output_variable%standard_name = 'northward_wind'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'w' )
-           output_variable%long_name     = 'w wind component'
-           output_variable%standard_name = 'upward_air_velocity'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'wspeed' )
-           output_variable%long_name     = 'wind speed'
-           output_variable%standard_name = 'wind_speed'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'wdir' )
-           output_variable%long_name     = 'wind from direction'
-           output_variable%standard_name = 'wind_from_direction'
-           output_variable%units         = 'degrees'
-
-        CASE ( 'theta' )
-           output_variable%long_name     = 'air potential temperature'
-           output_variable%standard_name = 'air_potential_temperature'
-           output_variable%units         = 'K'
-
-        CASE ( 'utheta' )
-           output_variable%long_name     = 'eastward kinematic sensible heat flux in air'
-           output_variable%units         = 'K m s-1'
-
-        CASE ( 'vtheta' )
-           output_variable%long_name     = 'northward kinematic sensible heat flux in air'
-           output_variable%units         = 'K m s-1'
-
-        CASE ( 'wtheta' )
-           output_variable%long_name     = 'upward kinematic sensible heat flux in air'
-           output_variable%units         = 'K m s-1'
-
-        CASE ( 'ta' )
-           output_variable%long_name     = 'air temperature'
-           output_variable%standard_name = 'air_temperature'
-           output_variable%units         = 'degree_C'
-
-        CASE ( 'tva' )
-           output_variable%long_name     = 'virtual acoustic temperature'
-           output_variable%units         = 'K'
-
-        CASE ( 'haa' )
-           output_variable%long_name     = 'absolute atmospheric humidity'
-           output_variable%units         = 'kg m-3'
-
-        CASE ( 'hus' )
-           output_variable%long_name     = 'specific humidity'
-           output_variable%standard_name = 'specific_humidity'
-           output_variable%units         = 'kg kg-1'
-
-        CASE ( 'hur' )
-           output_variable%long_name     = 'relative humidity'
-           output_variable%standard_name = 'relative_humidity'
-           output_variable%units         = '1'
-
-        CASE ( 'rlu' )
-           output_variable%long_name     = 'upwelling longwave flux in air'
-           output_variable%standard_name = 'upwelling_longwave_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'rlus' )
-           output_variable%long_name     = 'surface upwelling longwave flux in air'
-           output_variable%standard_name = 'surface_upwelling_longwave_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'rld' )
-           output_variable%long_name     = 'downwelling longwave flux in air'
-           output_variable%standard_name = 'downwelling_longwave_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'rsddif' )
-           output_variable%long_name     = 'diffuse downwelling shortwave flux in air'
-           output_variable%standard_name = 'diffuse_downwelling_shortwave_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'rsd' )
-           output_variable%long_name     = 'downwelling shortwave flux in air'
-           output_variable%standard_name = 'downwelling_shortwave_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'rnds' )
-           output_variable%long_name     = 'surface net downward radiative flux'
-           output_variable%standard_name = 'surface_net_downward_radiative_flux'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'rsu' )
-           output_variable%long_name     = 'upwelling shortwave flux in air'
-           output_variable%standard_name = 'upwelling_shortwave_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'rsus' )
-           output_variable%long_name     = 'surface upwelling shortwave flux in air'
-           output_variable%standard_name = 'surface_upwelling_shortwave_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'rsds' )
-           output_variable%long_name     = 'surface downwelling shortwave flux in air'
-           output_variable%standard_name = 'surface_downwelling_shortwave_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'hfss' )
-           output_variable%long_name     = 'surface upward sensible heat flux'
-           output_variable%standard_name = 'surface_upward_sensible_heat_flux'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'hfls' )
-           output_variable%long_name     = 'surface upward latent heat flux'
-           output_variable%standard_name = 'surface_upward_latent_heat_flux'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'ts' )
-           output_variable%long_name     = 'surface temperature'
-           output_variable%standard_name = 'surface_temperature'
-           output_variable%units         = 'K'
-
-        CASE ( 'thetas' )
-           output_variable%long_name     = 'surface layer temperature scale'
-           output_variable%units         = 'K'
-
-        CASE ( 'us' )
-           output_variable%long_name     = 'friction velocity'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'uw' )
-           output_variable%long_name     = 'upward eastward kinematic momentum flux in air'
-           output_variable%units         = 'm2 s-2'
-
-        CASE ( 'vw' )
-           output_variable%long_name     = 'upward northward kinematic momentum flux in air'
-           output_variable%units         = 'm2 s-2'
-
-        CASE ( 'uv' )
-           output_variable%long_name     = 'eastward northward kinematic momentum flux in air'
-           output_variable%units         = 'm2 s-2'
-
-        CASE ( 'plev' )
-           output_variable%long_name     = 'air pressure'
-           output_variable%standard_name = 'air_pressure'
-           output_variable%units         = 'Pa'
-
-        CASE ( 'm_soil' )
-           output_variable%long_name     = 'soil moisture volumetric'
-           output_variable%units         = 'm3 m-3'
-
-        CASE ( 't_soil' )
-           output_variable%long_name     = 'soil temperature'
-           output_variable%standard_name = 'soil_temperature'
-           output_variable%units         = 'degree_C'
-
-        CASE ( 'hfdg' )
-           output_variable%long_name     = 'downward heat flux at ground level in soil'
-           output_variable%standard_name = 'downward_heat_flux_at_ground_level_in_soil'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'hfds' )
-           output_variable%long_name     = 'downward heat flux in soil'
-           output_variable%standard_name = 'downward_heat_flux_in_soil'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'hfla' )
-           output_variable%long_name     = 'upward latent heat flux in air'
-           output_variable%standard_name = 'upward_latent_heat_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'hfsa' )
-           output_variable%long_name     = 'upward latent heat flux in air'
-           output_variable%standard_name = 'upward_sensible_heat_flux_in_air'
-           output_variable%units         = 'W m-2'
-
-        CASE ( 'jno2' )
-           output_variable%long_name     = 'photolysis rate of nitrogen dioxide'
-           output_variable%standard_name = 'photolysis_rate_of_nitrogen_dioxide'
-           output_variable%units         = 's-1'
-
-        CASE ( 'lwcs' )
-           output_variable%long_name     = 'liquid water content of soil layer'
-           output_variable%standard_name = 'liquid_water_content_of_soil_layer'
-           output_variable%units         = 'kg m-2'
-
-        CASE ( 'lwp' )
-           output_variable%long_name     = 'liquid water path'
-           output_variable%standard_name = 'atmosphere_mass_content_of_cloud_liquid_water'
-           output_variable%units         = 'kg m-2'
-
-        CASE ( 'ps' )
-           output_variable%long_name     = 'surface air pressure'
-           output_variable%standard_name = 'surface_air_pressure'
-           output_variable%units         = 'hPa'
-
-        CASE ( 'pswrtg' )
-           output_variable%long_name     = 'platform speed wrt ground'
-           output_variable%standard_name = 'platform_speed_wrt_ground'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'pswrta' )
-           output_variable%long_name     = 'platform speed wrt air'
-           output_variable%standard_name = 'platform_speed_wrt_air'
-           output_variable%units         = 'm s-1'
-
-        CASE ( 'pwv' )
-           output_variable%long_name     = 'water vapor partial pressure in air'
-           output_variable%standard_name = 'water_vapor_partial_pressure_in_air'
-           output_variable%units         = 'hPa'
-
-        CASE ( 'ssdu' )
-           output_variable%long_name     = 'duration of sunshine'
-           output_variable%standard_name = 'duration_of_sunshine'
-           output_variable%units         = 's'
-
-        CASE ( 't_lw' )
-           output_variable%long_name     = 'land water temperature'
-           output_variable%units         = 'degree_C'
-
-        CASE ( 'tb' )
-           output_variable%long_name     = 'brightness temperature'
-           output_variable%standard_name = 'brightness_temperature'
-           output_variable%units         = 'K'
-
-        CASE ( 'uqv' )
-           output_variable%long_name     = 'eastward kinematic latent heat flux in air'
-           output_variable%units         = 'g kg-1 m s-1'
-
-        CASE ( 'vqv' )
-           output_variable%long_name     = 'northward kinematic latent heat flux in air'
-           output_variable%units         = 'g kg-1 m s-1'
-
-        CASE ( 'wqv' )
-           output_variable%long_name     = 'upward kinematic latent heat flux in air'
-           output_variable%units         = 'g kg-1 m s-1'
-
-        CASE ( 'zcb' )
-           output_variable%long_name     = 'cloud base altitude'
-           output_variable%standard_name = 'cloud_base_altitude'
-           output_variable%units         = 'm'
-
-        CASE ( 'zmla' )
-           output_variable%long_name     = 'atmosphere boundary layer thickness'
-           output_variable%standard_name = 'atmosphere_boundary_layer_thickness'
-           output_variable%units         = 'm'
-
-        CASE ( 'mcpm1' )
-           output_variable%long_name     = 'mass concentration of pm1 ambient aerosol particles in air'
-           output_variable%standard_name = 'mass_concentration_of_pm1_ambient_aerosol_particles_in_air'
-           output_variable%units         = 'kg m-3'
-
-        CASE ( 'mcpm10' )
-           output_variable%long_name     = 'mass concentration of pm10 ambient aerosol particles in air'
-           output_variable%standard_name = 'mass_concentration_of_pm10_ambient_aerosol_particles_in_air'
-           output_variable%units         = 'kg m-3'
-
-        CASE ( 'mcpm2p5' )
-           output_variable%long_name     = 'mass concentration of pm2p5 ambient aerosol particles in air'
-           output_variable%standard_name = 'mass_concentration_of_pm2p5_ambient_aerosol_particles_in_air'
-           output_variable%units         = 'kg m-3'
-
-        CASE ( 'mfno', 'mcno'  )
-           output_variable%long_name     = 'mole fraction of nitrogen monoxide in air'
-           output_variable%standard_name = 'mole_fraction_of_nitrogen_monoxide_in_air'
-           output_variable%units         = 'ppm' !'mol mol-1'
-
-        CASE ( 'mfno2', 'mcno2'  )
-           output_variable%long_name     = 'mole fraction of nitrogen dioxide in air'
-           output_variable%standard_name = 'mole_fraction_of_nitrogen_dioxide_in_air'
-           output_variable%units         = 'ppm' !'mol mol-1'
-
-        CASE ( 'tro3'  )
-           output_variable%long_name     = 'mole fraction of ozone in air'
-           output_variable%standard_name = 'mole_fraction_of_ozone_in_air'
-           output_variable%units         = 'ppm' !'mol mol-1'
-
-        CASE DEFAULT
-
-     END SELECT
-
-  END SUBROUTINE vm_set_attributes
-
-
-!------------------------------------------------------------------------------!
+!> Subroutine defines variable attributes according to UC2 standard. Note, later  this list can be
+!> moved to the data-output module where it can be re-used also for other output.
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE vm_set_attributes( output_variable )
+
+    TYPE( virt_var_atts ), INTENT(INOUT) ::  output_variable !< data structure with attributes that need to be set
+
+    output_variable%long_name     = 'none'
+    output_variable%standard_name = 'none'
+    output_variable%units         = 'none'
+    output_variable%coordinates   = 'lon lat E_UTM N_UTM x y z time station_name'
+    output_variable%grid_mapping  = 'crs'
+
+    SELECT CASE ( TRIM( output_variable%name ) )
+
+       CASE ( 'u' )
+          output_variable%long_name     = 'u wind component'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'ua' )
+          output_variable%long_name     = 'eastward wind'
+          output_variable%standard_name = 'eastward_wind'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'v' )
+          output_variable%long_name     = 'v wind component'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'va' )
+          output_variable%long_name     = 'northward wind'
+          output_variable%standard_name = 'northward_wind'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'w' )
+          output_variable%long_name     = 'w wind component'
+          output_variable%standard_name = 'upward_air_velocity'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'wspeed' )
+          output_variable%long_name     = 'wind speed'
+          output_variable%standard_name = 'wind_speed'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'wdir' )
+          output_variable%long_name     = 'wind from direction'
+          output_variable%standard_name = 'wind_from_direction'
+          output_variable%units         = 'degrees'
+
+       CASE ( 'theta' )
+          output_variable%long_name     = 'air potential temperature'
+          output_variable%standard_name = 'air_potential_temperature'
+          output_variable%units         = 'K'
+
+       CASE ( 'utheta' )
+          output_variable%long_name     = 'eastward kinematic sensible heat flux in air'
+          output_variable%units         = 'K m s-1'
+
+       CASE ( 'vtheta' )
+          output_variable%long_name     = 'northward kinematic sensible heat flux in air'
+          output_variable%units         = 'K m s-1'
+
+       CASE ( 'wtheta' )
+          output_variable%long_name     = 'upward kinematic sensible heat flux in air'
+          output_variable%units         = 'K m s-1'
+
+       CASE ( 'ta' )
+          output_variable%long_name     = 'air temperature'
+          output_variable%standard_name = 'air_temperature'
+          output_variable%units         = 'degree_C'
+
+       CASE ( 'tva' )
+          output_variable%long_name     = 'virtual acoustic temperature'
+          output_variable%units         = 'K'
+
+       CASE ( 'haa' )
+          output_variable%long_name     = 'absolute atmospheric humidity'
+          output_variable%units         = 'kg m-3'
+
+       CASE ( 'hus' )
+          output_variable%long_name     = 'specific humidity'
+          output_variable%standard_name = 'specific_humidity'
+          output_variable%units         = 'kg kg-1'
+
+       CASE ( 'hur' )
+          output_variable%long_name     = 'relative humidity'
+          output_variable%standard_name = 'relative_humidity'
+          output_variable%units         = '1'
+
+       CASE ( 'rlu' )
+          output_variable%long_name     = 'upwelling longwave flux in air'
+          output_variable%standard_name = 'upwelling_longwave_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'rlus' )
+          output_variable%long_name     = 'surface upwelling longwave flux in air'
+          output_variable%standard_name = 'surface_upwelling_longwave_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'rld' )
+          output_variable%long_name     = 'downwelling longwave flux in air'
+          output_variable%standard_name = 'downwelling_longwave_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'rsddif' )
+          output_variable%long_name     = 'diffuse downwelling shortwave flux in air'
+          output_variable%standard_name = 'diffuse_downwelling_shortwave_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'rsd' )
+          output_variable%long_name     = 'downwelling shortwave flux in air'
+          output_variable%standard_name = 'downwelling_shortwave_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'rnds' )
+          output_variable%long_name     = 'surface net downward radiative flux'
+          output_variable%standard_name = 'surface_net_downward_radiative_flux'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'rsu' )
+          output_variable%long_name     = 'upwelling shortwave flux in air'
+          output_variable%standard_name = 'upwelling_shortwave_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'rsus' )
+          output_variable%long_name     = 'surface upwelling shortwave flux in air'
+          output_variable%standard_name = 'surface_upwelling_shortwave_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'rsds' )
+          output_variable%long_name     = 'surface downwelling shortwave flux in air'
+          output_variable%standard_name = 'surface_downwelling_shortwave_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'hfss' )
+          output_variable%long_name     = 'surface upward sensible heat flux'
+          output_variable%standard_name = 'surface_upward_sensible_heat_flux'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'hfls' )
+          output_variable%long_name     = 'surface upward latent heat flux'
+          output_variable%standard_name = 'surface_upward_latent_heat_flux'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'ts' )
+          output_variable%long_name     = 'surface temperature'
+          output_variable%standard_name = 'surface_temperature'
+          output_variable%units         = 'K'
+
+       CASE ( 'thetas' )
+          output_variable%long_name     = 'surface layer temperature scale'
+          output_variable%units         = 'K'
+
+       CASE ( 'us' )
+          output_variable%long_name     = 'friction velocity'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'uw' )
+          output_variable%long_name     = 'upward eastward kinematic momentum flux in air'
+          output_variable%units         = 'm2 s-2'
+
+       CASE ( 'vw' )
+          output_variable%long_name     = 'upward northward kinematic momentum flux in air'
+          output_variable%units         = 'm2 s-2'
+
+       CASE ( 'uv' )
+          output_variable%long_name     = 'eastward northward kinematic momentum flux in air'
+          output_variable%units         = 'm2 s-2'
+
+       CASE ( 'plev' )
+          output_variable%long_name     = 'air pressure'
+          output_variable%standard_name = 'air_pressure'
+          output_variable%units         = 'Pa'
+
+       CASE ( 'm_soil' )
+          output_variable%long_name     = 'soil moisture volumetric'
+          output_variable%units         = 'm3 m-3'
+
+       CASE ( 't_soil' )
+          output_variable%long_name     = 'soil temperature'
+          output_variable%standard_name = 'soil_temperature'
+          output_variable%units         = 'degree_C'
+
+       CASE ( 'hfdg' )
+          output_variable%long_name     = 'downward heat flux at ground level in soil'
+          output_variable%standard_name = 'downward_heat_flux_at_ground_level_in_soil'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'hfds' )
+          output_variable%long_name     = 'downward heat flux in soil'
+          output_variable%standard_name = 'downward_heat_flux_in_soil'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'hfla' )
+          output_variable%long_name     = 'upward latent heat flux in air'
+          output_variable%standard_name = 'upward_latent_heat_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'hfsa' )
+          output_variable%long_name     = 'upward latent heat flux in air'
+          output_variable%standard_name = 'upward_sensible_heat_flux_in_air'
+          output_variable%units         = 'W m-2'
+
+       CASE ( 'jno2' )
+          output_variable%long_name     = 'photolysis rate of nitrogen dioxide'
+          output_variable%standard_name = 'photolysis_rate_of_nitrogen_dioxide'
+          output_variable%units         = 's-1'
+
+       CASE ( 'lwcs' )
+          output_variable%long_name     = 'liquid water content of soil layer'
+          output_variable%standard_name = 'liquid_water_content_of_soil_layer'
+          output_variable%units         = 'kg m-2'
+
+       CASE ( 'lwp' )
+          output_variable%long_name     = 'liquid water path'
+          output_variable%standard_name = 'atmosphere_mass_content_of_cloud_liquid_water'
+          output_variable%units         = 'kg m-2'
+
+       CASE ( 'ps' )
+          output_variable%long_name     = 'surface air pressure'
+          output_variable%standard_name = 'surface_air_pressure'
+          output_variable%units         = 'hPa'
+
+       CASE ( 'pswrtg' )
+          output_variable%long_name     = 'platform speed wrt ground'
+          output_variable%standard_name = 'platform_speed_wrt_ground'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'pswrta' )
+          output_variable%long_name     = 'platform speed wrt air'
+          output_variable%standard_name = 'platform_speed_wrt_air'
+          output_variable%units         = 'm s-1'
+
+       CASE ( 'pwv' )
+          output_variable%long_name     = 'water vapor partial pressure in air'
+          output_variable%standard_name = 'water_vapor_partial_pressure_in_air'
+          output_variable%units         = 'hPa'
+
+       CASE ( 'ssdu' )
+          output_variable%long_name     = 'duration of sunshine'
+          output_variable%standard_name = 'duration_of_sunshine'
+          output_variable%units         = 's'
+
+       CASE ( 't_lw' )
+          output_variable%long_name     = 'land water temperature'
+          output_variable%units         = 'degree_C'
+
+       CASE ( 'tb' )
+          output_variable%long_name     = 'brightness temperature'
+          output_variable%standard_name = 'brightness_temperature'
+          output_variable%units         = 'K'
+
+       CASE ( 'uqv' )
+          output_variable%long_name     = 'eastward kinematic latent heat flux in air'
+          output_variable%units         = 'g kg-1 m s-1'
+
+       CASE ( 'vqv' )
+          output_variable%long_name     = 'northward kinematic latent heat flux in air'
+          output_variable%units         = 'g kg-1 m s-1'
+
+       CASE ( 'wqv' )
+          output_variable%long_name     = 'upward kinematic latent heat flux in air'
+          output_variable%units         = 'g kg-1 m s-1'
+
+       CASE ( 'zcb' )
+          output_variable%long_name     = 'cloud base altitude'
+          output_variable%standard_name = 'cloud_base_altitude'
+          output_variable%units         = 'm'
+
+       CASE ( 'zmla' )
+          output_variable%long_name     = 'atmosphere boundary layer thickness'
+          output_variable%standard_name = 'atmosphere_boundary_layer_thickness'
+          output_variable%units         = 'm'
+
+       CASE ( 'mcpm1' )
+          output_variable%long_name     = 'mass concentration of pm1 ambient aerosol particles in air'
+          output_variable%standard_name = 'mass_concentration_of_pm1_ambient_aerosol_particles_in_air'
+          output_variable%units         = 'kg m-3'
+
+       CASE ( 'mcpm10' )
+          output_variable%long_name     = 'mass concentration of pm10 ambient aerosol particles in air'
+          output_variable%standard_name = 'mass_concentration_of_pm10_ambient_aerosol_particles_in_air'
+          output_variable%units         = 'kg m-3'
+
+       CASE ( 'mcpm2p5' )
+          output_variable%long_name     = 'mass concentration of pm2p5 ambient aerosol particles in air'
+          output_variable%standard_name = 'mass_concentration_of_pm2p5_ambient_aerosol_particles_in_air'
+          output_variable%units         = 'kg m-3'
+
+       CASE ( 'mfno', 'mcno'  )
+          output_variable%long_name     = 'mole fraction of nitrogen monoxide in air'
+          output_variable%standard_name = 'mole_fraction_of_nitrogen_monoxide_in_air'
+          output_variable%units         = 'ppm' !'mol mol-1'
+
+       CASE ( 'mfno2', 'mcno2'  )
+          output_variable%long_name     = 'mole fraction of nitrogen dioxide in air'
+          output_variable%standard_name = 'mole_fraction_of_nitrogen_dioxide_in_air'
+          output_variable%units         = 'ppm' !'mol mol-1'
+
+       CASE ( 'tro3'  )
+          output_variable%long_name     = 'mole fraction of ozone in air'
+          output_variable%standard_name = 'mole_fraction_of_ozone_in_air'
+          output_variable%units         = 'ppm' !'mol mol-1'
+
+       CASE DEFAULT
+
+    END SELECT
+
+ END SUBROUTINE vm_set_attributes
+
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Read namelist for the virtual measurement module
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE vm_parin
 
-    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
-
-    NAMELIST /virtual_measurement_parameters/  dt_virtual_measurement,         &
-                                               off_ts,                         &
-                                               off_pr,                         &
-                                               off_tr,                         &
-                                               use_virtual_measurement,        &
+    CHARACTER(LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
+
+    NAMELIST /virtual_measurement_parameters/  dt_virtual_measurement,                             &
+                                               off_ts,                                             &
+                                               off_pr,                                             &
+                                               off_tr,                                             &
+                                               use_virtual_measurement,                            &
                                                vm_time_start
 
@@ -810,5 +802,5 @@
     REWIND ( 11 )
     line = ' '
-    DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
+    DO  WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
        READ ( 11, '(A)', END=20 )  line
     ENDDO
@@ -822,5 +814,4 @@
 !-- Set flag that indicates that the virtual measurement module is switched on
     IF ( use_virtual_measurement )  virtual_measurement = .TRUE.
-
     GOTO 20
 
@@ -834,49 +825,48 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Initialize virtual measurements: read coordiante arrays and measured
-!> variables, set indicies indicating the measurement points, read further
-!> attributes, etc..
-!------------------------------------------------------------------------------!
+!> Initialize virtual measurements: read coordiante arrays and measured variables, set indicies
+!> indicating the measurement points, read further attributes, etc..
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE vm_init
 
-    CHARACTER(LEN=5)                  ::  dum                          !< dummy string indicating station id
-    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
-    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed
-
-    INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate
-    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
-    INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
-    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
-    INTEGER(iwp) ::  js        !< grid index of real observation point of the respective station in y-direction
-    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
-    INTEGER(iwp) ::  kl        !< lower vertical index of surrounding grid points of an observation coordinate
-    INTEGER(iwp) ::  ks        !< grid index of real observation point of the respective station in z-direction
-    INTEGER(iwp) ::  ksurf     !< topography top index
-    INTEGER(iwp) ::  ku        !< upper vertical index of surrounding grid points of an observation coordinate
-    INTEGER(iwp) ::  l         !< running index over all stations
-    INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
-    INTEGER(iwp) ::  ll        !< running index over all measured variables in file
-    INTEGER(iwp) ::  m         !< running index for surface elements
-    INTEGER(iwp) ::  n         !< running index over trajectory coordinates
-    INTEGER(iwp) ::  nofill    !< dummy for nofill return value (not used)
-    INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
-    INTEGER(iwp) ::  off       !< number of surrounding grid points to be sampled
-    INTEGER(iwp) ::  t         !< running index over number of trajectories
-
-    INTEGER(KIND=1)                             ::  soil_dum !< dummy variable to input a soil flag
-
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE     ::  ns_all !< dummy array used to sum-up the number of observation coordinates
+    CHARACTER(LEN=5)                  ::  dum                           !< dummy string indicating station id
+    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables_file = ''  !< array with all measured variables read from NetCDF
+    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables      = ''  !< dummy array with all measured variables that are allowed
+
+    INTEGER(iwp) ::  dim_ntime  !< dimension size of time coordinate
+    INTEGER(iwp) ::  i          !< grid index of virtual observation point in x-direction
+    INTEGER(iwp) ::  is         !< grid index of real observation point of the respective station in x-direction
+    INTEGER(iwp) ::  j          !< grid index of observation point in x-direction
+    INTEGER(iwp) ::  js         !< grid index of real observation point of the respective station in y-direction
+    INTEGER(iwp) ::  k          !< grid index of observation point in x-direction
+    INTEGER(iwp) ::  kl         !< lower vertical index of surrounding grid points of an observation coordinate
+    INTEGER(iwp) ::  ks         !< grid index of real observation point of the respective station in z-direction
+    INTEGER(iwp) ::  ksurf      !< topography top index
+    INTEGER(iwp) ::  ku         !< upper vertical index of surrounding grid points of an observation coordinate
+    INTEGER(iwp) ::  l          !< running index over all stations
+    INTEGER(iwp) ::  len_char   !< character length of single measured variables without Null character
+    INTEGER(iwp) ::  ll         !< running index over all measured variables in file
+    INTEGER(iwp) ::  m          !< running index for surface elements
+    INTEGER(iwp) ::  n          !< running index over trajectory coordinates
+    INTEGER(iwp) ::  nofill     !< dummy for nofill return value (not used)
+    INTEGER(iwp) ::  ns         !< counter variable for number of observation points on subdomain
+    INTEGER(iwp) ::  off        !< number of surrounding grid points to be sampled
+    INTEGER(iwp) ::  t          !< running index over number of trajectories
+
+    INTEGER(KIND=1)                             ::  soil_dum  !< dummy variable to input a soil flag
+
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE     ::  ns_all  !< dummy array used to sum-up the number of observation coordinates
 
 #if defined( __parallel )
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_atmos !< number of observation points for each station on each mpi rank
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_soil  !< number of observation points for each station on each mpi rank
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_atmos  !< number of observation points for each station on each mpi rank
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_soil   !< number of observation points for each station on each mpi rank
 #endif
 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
-
-    LOGICAL  ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag  !< mask array indicating measurement positions
+
+    LOGICAL  ::  on_pe  !< flag indicating that the respective measurement coordinate is on subdomain
 
     REAL(wp) ::  fill_eutm !< _FillValue for coordinate array E_UTM
@@ -884,10 +874,10 @@
     REAL(wp) ::  fill_zar  !< _FillValue for height coordinate
 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm     !< easting UTM coordinate, temporary variable
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm     !< northing UTM coordinate, temporary variable
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm_tmp !< EUTM coordinate before rotation
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm_tmp !< NUTM coordinate before rotation
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  station_h !< station height above reference
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zar       !< observation height above reference
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm      !< easting UTM coordinate, temporary variable
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm_tmp  !< EUTM coordinate before rotation
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm      !< northing UTM coordinate, temporary variable
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm_tmp  !< NUTM coordinate before rotation
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  station_h  !< station height above reference
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zar        !< observation height above reference
 #if defined( __netcdf )
 !
@@ -896,18 +886,13 @@
 !
 !-- Obtain number of sites.
-    CALL get_attribute( pids_id,                                               &
-                        char_numstations,                                      &
-                        vmea_general%nvm,                                      &
-                        global_attribute )
-!
-!-- Allocate data structure which encompass all required information, such as
-!-- grid points indicies, absolute UTM coordinates, the measured quantities,
-!-- etc. .
+    CALL get_attribute( pids_id, char_numstations, vmea_general%nvm, global_attribute )
+!
+!-- Allocate data structure which encompasses all required information, such as  grid points indicies,
+!-- absolute UTM coordinates, the measured quantities, etc. .
     ALLOCATE( vmea(1:vmea_general%nvm) )
 !
-!-- Allocate flag array. This dummy array is used to identify grid points
-!-- where virtual measurements should be taken. Please note, in order to
-!-- include also the surrounding grid points of the original coordinate
-!-- ghost points are required.
+!-- Allocate flag array. This dummy array is used to identify grid points where virtual measurements
+!-- should be taken. Please note, in order to include also the surrounding grid points of the
+!-- original coordinate, ghost points are required.
     ALLOCATE( meas_flag(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     meas_flag = 0
@@ -916,6 +901,5 @@
     DO  l = 1, vmea_general%nvm
 !
-!--    Determine suffix which contains the ID, ordered according to the number
-!--    of measurements.
+!--    Determine suffix which contains the ID, ordered according to the number of measurements.
        IF( l < 10 )  THEN
           WRITE( dum, '(I1)')  l
@@ -931,25 +915,13 @@
 !
 !--    Read the origin site coordinates (UTM).
-       CALL get_attribute( pids_id,                                            &
-                           char_origx // TRIM( dum ),                          &
-                           vmea(l)%origin_x_obs,                               &
-                           global_attribute )
-       CALL get_attribute( pids_id,                                            &
-                           char_origy // TRIM( dum ),                          &
-                           vmea(l)%origin_y_obs,                               &
-                           global_attribute )
+       CALL get_attribute( pids_id, char_origx // TRIM( dum ), vmea(l)%origin_x_obs, global_attribute )
+       CALL get_attribute( pids_id, char_origy // TRIM( dum ), vmea(l)%origin_y_obs, global_attribute )
 !
 !--    Read site name.
-       CALL get_attribute( pids_id,                                            &
-                           char_site // TRIM( dum ),                           &
-                           vmea(l)%site,                                       &
-                           global_attribute )
-!
-!--    Read a flag which indicates that also soil quantities are take at the
-!--    respective site (is part of the virtual measurement driver).
-       CALL get_attribute( pids_id,                                            &
-                           char_soil // TRIM( dum ),                           &
-                           soil_dum,                                           &
-                           global_attribute )
+       CALL get_attribute( pids_id, char_site // TRIM( dum ), vmea(l)%site, global_attribute )
+!
+!--    Read a flag which indicates that also soil quantities are take at the respective site
+!--    (is part of the virtual measurement driver).
+       CALL get_attribute( pids_id, char_soil // TRIM( dum ), soil_dum, global_attribute )
 !
 !--    Set flag indicating soil-sampling.
@@ -957,8 +929,5 @@
 !
 !--    Read type of the measurement (trajectory, profile, timeseries).
-       CALL get_attribute( pids_id,                                            &
-                           char_feature // TRIM( dum ),                        &
-                           vmea(l)%feature_type,                               &
-                           global_attribute )
+       CALL get_attribute( pids_id, char_feature // TRIM( dum ), vmea(l)%feature_type, global_attribute )
 !
 !---   Set logicals depending on the type of the measurement
@@ -972,7 +941,5 @@
 !--    Give error message in case the type matches non of the pre-defined types.
        ELSE
-          message_string = 'Attribue featureType = ' //                        &
-                           TRIM( vmea(l)%feature_type ) //                     &
-                           ' is not allowed.'
+          message_string = 'Attribue featureType = ' // TRIM( vmea(l)%feature_type ) // ' is not allowed.'
           CALL message( 'vm_init', 'PA0535', 1, 2, 0, 6, 0 )
        ENDIF
@@ -980,28 +947,24 @@
 !--    Read string with all measured variables at this site.
        measured_variables_file = ''
-       CALL get_variable( pids_id,                                             &
-                          char_mv // TRIM( dum ),                              &
-                          measured_variables_file )
+       CALL get_variable( pids_id, char_mv // TRIM( dum ), measured_variables_file )
 !
 !--    Count the number of measured variables.
-!--    Please note, for some NetCDF interal reasons characters end with a NULL,
-!--    i.e. also empty characters contain a NULL. Therefore, check the strings
-!--    for a NULL to get the correct character length in order to compare
-!--    them with the list of allowed variables.
+!--    Please note, for some NetCDF interal reasons, characters end with a NULL, i.e. also empty
+!--    characters contain a NULL. Therefore, check the strings for a NULL to get the correct
+!--    character length in order to compare them with the list of allowed variables.
        vmea(l)%nmeas  = 1
-       DO ll = 1, SIZE( measured_variables_file )
-          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
+       DO  ll = 1, SIZE( measured_variables_file )
+          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.                                  &
                measured_variables_file(ll)(1:1) /= ' ')  THEN
 !
 !--          Obtain character length of the character
              len_char = 1
-             DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
-                 .AND.  measured_variables_file(ll)(len_char:len_char) /= ' ' )
+             DO  WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)  .AND.          &
+                 measured_variables_file(ll)(len_char:len_char) /= ' ' )
                 len_char = len_char + 1
              ENDDO
              len_char = len_char - 1
 
-             measured_variables(vmea(l)%nmeas) =                               &
-                                       measured_variables_file(ll)(1:len_char)
+             measured_variables(vmea(l)%nmeas) = measured_variables_file(ll)(1:len_char)
              vmea(l)%nmeas = vmea(l)%nmeas + 1
 
@@ -1010,57 +973,47 @@
        vmea(l)%nmeas = vmea(l)%nmeas - 1
 !
-!--    Allocate data-type array for the measured variables names and attributes
-!--    at the respective site.
+!--    Allocate data-type array for the measured variables names and attributes at the respective
+!--    site.
        ALLOCATE( vmea(l)%var_atts(1:vmea(l)%nmeas) )
 !
-!--    Store the variable names in a data structures, which assigns further
-!--    attributes to this name. Further, for data output reasons, create a
-!--    string of output variables, which will be written into the attribute
-!--    data_content.
+!--    Store the variable names in a data structure, which assigns further attributes to this name.
+!--    Further, for data output reasons, create a string of output variables, which will be written
+!--    into the attribute data_content.
        DO  ll = 1, vmea(l)%nmeas
           vmea(l)%var_atts(ll)%name = TRIM( measured_variables(ll) )
 
-          vmea(l)%data_content = TRIM( vmea(l)%data_content ) // " " //        &
+          vmea(l)%data_content = TRIM( vmea(l)%data_content ) // " " //                            &
                                  TRIM( vmea(l)%var_atts(ll)%name )
        ENDDO
 !
-!--    Read all the UTM coordinates for the site. Based on the coordinates,
-!--    define the grid-index space on each subdomain where virtual measurements
-!--    should be taken. Note, the entire coordinate array (on the entire model
-!--    domain) won't be stored as this would exceed memory requirements,
-!--    particularly for trajectories.
+!--    Read all the UTM coordinates for the site. Based on the coordinates, define the grid-index
+!--    space on each subdomain where virtual measurements should be taken. Note, the entire
+!--    coordinate array (on the entire model domain) won't be stored as this would exceed memory
+!--    requirements, particularly for trajectories.
        IF ( vmea(l)%nmeas > 0 )  THEN
 !
-!--       For stationary measurements UTM coordinates are just one value and
-!--       its dimension is "station", while for mobile measurements UTM
-!--       coordinates are arrays depending on the number of trajectories and
-!--       time, according to (UC)2 standard. First, inquire dimension length
-!--       of the UTM coordinates.
+!--       For stationary measurements UTM coordinates are just one value and its dimension is
+!--       "station", while for mobile measurements UTM coordinates are arrays depending on the
+!--       number of trajectories and time, according to (UC)2 standard. First, inquire dimension
+!--       length of the UTM coordinates.
           IF ( vmea(l)%trajectory )  THEN
 !
-!--          For non-stationary measurements read the number of trajectories
-!--          and the number of time coordinates.
-             CALL get_dimension_length( pids_id,                               &
-                                        vmea(l)%n_tr_st,                       &
-                                        "traj" // TRIM( dum ) )
-             CALL get_dimension_length( pids_id,                               &
-                                        dim_ntime,                             &
-                                        "ntime" // TRIM( dum ) )
-!
-!--       For stationary measurements the dimension for UTM is station
-!--       and for the time-coordinate it is one.
+!--          For non-stationary measurements read the number of trajectories and the number of time
+!--          coordinates.
+             CALL get_dimension_length( pids_id, vmea(l)%n_tr_st, "traj" // TRIM( dum ) )
+             CALL get_dimension_length( pids_id, dim_ntime, "ntime" // TRIM( dum ) )
+!
+!--       For stationary measurements the dimension for UTM is station and for the time-coordinate
+!--       it is one.
           ELSE
-             CALL get_dimension_length( pids_id,                               &
-                                        vmea(l)%n_tr_st,                       &
-                                        "station" // TRIM( dum ) )
+             CALL get_dimension_length( pids_id, vmea(l)%n_tr_st, "station" // TRIM( dum ) )
              dim_ntime = 1
           ENDIF
 !
-!-        Allocate array which defines individual time/space frame for each
-!--       trajectory or station.
+!-        Allocate array which defines individual time/space frame for each trajectory or station.
           ALLOCATE( vmea(l)%dim_t(1:vmea(l)%n_tr_st) )
 !
-!--       Allocate temporary arrays for UTM and height coordinates. Note,
-!--       on file UTM coordinates might be 1D or 2D variables
+!--       Allocate temporary arrays for UTM and height coordinates. Note, on file UTM coordinates
+!--       might be 1D or 2D variables
           ALLOCATE( e_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
           ALLOCATE( n_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
@@ -1075,56 +1028,28 @@
           ALLOCATE( n_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
 !
-!--       Read UTM and height coordinates coordinates for all trajectories and
-!--       times. Note, in case these obtain any missing values, replace them
-!--       with default _FillValues.
-          CALL inquire_fill_value( pids_id,                                    &
-                                   char_eutm // TRIM( dum ),                   &
-                                   nofill,                                     &
-                                   fill_eutm )
-          CALL inquire_fill_value( pids_id,                                    &
-                                   char_nutm // TRIM( dum ),                   &
-                                   nofill,                                     &
-                                   fill_nutm )
-          CALL inquire_fill_value( pids_id,                                    &
-                                   char_zar // TRIM( dum ),                    &
-                                   nofill,                                     &
-                                   fill_zar )
-!
-!--       Further line is just to avoid compiler warnings. nofill might be used
-!--       in future.
+!--       Read UTM and height coordinates for all trajectories and times. Note, in case
+!--       these obtain any missing values, replace them with default _FillValues.
+          CALL inquire_fill_value( pids_id, char_eutm // TRIM( dum ), nofill, fill_eutm )
+          CALL inquire_fill_value( pids_id, char_nutm // TRIM( dum ), nofill, fill_nutm )
+          CALL inquire_fill_value( pids_id, char_zar // TRIM( dum ), nofill, fill_zar )
+!
+!--       Further line is just to avoid compiler warnings. nofill might be used in future.
           IF ( nofill == 0  .OR.  nofill /= 0 )  CONTINUE
 !
-!--       Read observation coordinates. Please note, for trajectories the
-!--       observation height is stored directly in z, while for timeSeries
-!--       it is stored in z - station_h, according to UC2-standard.
+!--       Read observation coordinates. Please note, for trajectories the observation height is
+!--       stored directly in z, while for timeSeries it is stored in z - station_h, according to
+!--       UC2-standard.
           IF ( vmea(l)%trajectory )  THEN
-             CALL get_variable( pids_id,                                       &
-                                char_eutm // TRIM( dum ),                      &
-                                e_utm,                                         &
-                                0, dim_ntime-1,                                &
-                                0, vmea(l)%n_tr_st-1 )
-             CALL get_variable( pids_id,                                       &
-                                char_nutm // TRIM( dum ),                      &
-                                n_utm,                                         &
-                                0, dim_ntime-1,                                &
-                                0, vmea(l)%n_tr_st-1 )
-             CALL get_variable( pids_id,                                       &
-                                char_zar // TRIM( dum ),                       &
-                                zar,                                           &
-                                0, dim_ntime-1,                                &
-                                0, vmea(l)%n_tr_st-1 )
+             CALL get_variable( pids_id, char_eutm // TRIM( dum ), e_utm, 0, dim_ntime-1, 0,       &
+                                vmea(l)%n_tr_st-1 )
+             CALL get_variable( pids_id, char_nutm // TRIM( dum ), n_utm, 0, dim_ntime-1, 0,       &
+                                vmea(l)%n_tr_st-1 )
+             CALL get_variable( pids_id, char_zar // TRIM( dum ), zar, 0, dim_ntime-1, 0,          &
+                                vmea(l)%n_tr_st-1 )
           ELSE
-             CALL get_variable( pids_id,                                       &
-                                char_eutm // TRIM( dum ),                      &
-                                e_utm(:,1) )
-             CALL get_variable( pids_id,                                       &
-                                char_nutm // TRIM( dum ),                      &
-                                n_utm(:,1) )
-             CALL get_variable( pids_id,                                       &
-                                char_station_h // TRIM( dum ),                 &
-                                station_h(:,1) )
-             CALL get_variable( pids_id,                                       &
-                                char_zar // TRIM( dum ),                       &
-                                zar(:,1) )
+             CALL get_variable( pids_id, char_eutm // TRIM( dum ), e_utm(:,1) )
+             CALL get_variable( pids_id, char_nutm // TRIM( dum ), n_utm(:,1) )
+             CALL get_variable( pids_id, char_station_h // TRIM( dum ), station_h(:,1) )
+             CALL get_variable( pids_id, char_zar // TRIM( dum ), zar(:,1) )
           ENDIF
 
@@ -1136,45 +1061,39 @@
           zar  = zar - station_h
 !
-!--       Based on UTM coordinates, check if the measurement station or parts
-!--       of the trajectory are on subdomain. This case, setup grid index space
-!--       sample these quantities.
+!--       Based on UTM coordinates, check if the measurement station or parts of the trajectory are
+!--       on subdomain. This case, setup grid index space sample these quantities.
           meas_flag = 0
           DO  t = 1, vmea(l)%n_tr_st
 !
-!--          First, compute relative x- and y-coordinates with respect to the
-!--          lower-left origin of the model domain, which is the difference
-!--          between UTM coordinates. Note, if the origin is not correct, the
-!--          virtual sites will be misplaced. Further, in case of an rotated
-!--          model domain, the UTM coordinates must be also rotated.
+!--          First, compute relative x- and y-coordinates with respect to the lower-left origin of
+!--          the model domain, which is the difference between UTM coordinates. Note, if the origin
+!--          is not correct, the virtual sites will be misplaced. Further, in case of an rotated
+!--          model domain, the UTM coordinates must also be rotated.
              e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
              n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
-             e_utm(t,1:dim_ntime) = COS( init_model%rotation_angle * pi / 180.0_wp ) &
-                                    * e_utm_tmp(t,1:dim_ntime)                       &
-                                  - SIN( init_model%rotation_angle * pi / 180.0_wp ) &
+             e_utm(t,1:dim_ntime) = COS( init_model%rotation_angle * pi / 180.0_wp )               &
+                                    * e_utm_tmp(t,1:dim_ntime)                                     &
+                                  - SIN( init_model%rotation_angle * pi / 180.0_wp )               &
                                     * n_utm_tmp(t,1:dim_ntime)
-             n_utm(t,1:dim_ntime) = SIN( init_model%rotation_angle * pi / 180.0_wp ) &
-                                    * e_utm_tmp(t,1:dim_ntime)                       &
-                                  + COS( init_model%rotation_angle * pi / 180.0_wp ) &
+             n_utm(t,1:dim_ntime) = SIN( init_model%rotation_angle * pi / 180.0_wp )               &
+                                    * e_utm_tmp(t,1:dim_ntime)                                     &
+                                  + COS( init_model%rotation_angle * pi / 180.0_wp )               &
                                     * n_utm_tmp(t,1:dim_ntime)
 !
-!--          Determine the individual time coordinate length for each station and
-!--          trajectory. This is required as several stations and trajectories
-!--          are merged into one file but they do not have the same number of
-!--          points in time, hence, missing values may occur and cannot be
-!--          processed further. This is actually a work-around for the specific
-!--          (UC)2 dataset, but it won't harm anyway.
+!--          Determine the individual time coordinate length for each station and trajectory. This
+!--          is required as several stations and trajectories are merged into one file but they do
+!--          not have the same number of points in time, hence, missing values may occur and cannot
+!--          be processed further. This is actually a work-around for the specific (UC)2 dataset,
+!--          but it won't harm anyway.
              vmea(l)%dim_t(t) = 0
              DO  n = 1, dim_ntime
-                IF ( e_utm(t,n) /= fill_eutm  .AND.                            &
-                     n_utm(t,n) /= fill_nutm  .AND.                            &
+                IF ( e_utm(t,n) /= fill_eutm  .AND.  n_utm(t,n) /= fill_nutm  .AND.                &
                      zar(t,n)   /= fill_zar )  vmea(l)%dim_t(t) = n
              ENDDO
 !
-!--          Compute grid indices relative to origin and check if these are
-!--          on the subdomain. Note, virtual measurements will be taken also
-!--          at grid points surrounding the station, hence, check also for
-!--          these grid points.
-!--          The number of surrounding grid points is set according to the
-!--          featureType.
+!--          Compute grid indices relative to origin and check if these are on the subdomain. Note,
+!--          virtual measurements will be taken also at grid points surrounding the station, hence,
+!--          check also for these grid points. The number of surrounding grid points is set
+!--          according to the featureType.
              IF ( vmea(l)%timseries_profile )  THEN
                 off = off_pr
@@ -1186,22 +1105,20 @@
 
              DO  n = 1, vmea(l)%dim_t(t)
-                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
-                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
+                 is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
+                 js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
 !
 !--             Is the observation point on subdomain?
-                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.                   &
-                          js >= nys  .AND.  js <= nyn )
-!
-!--             Check if observation coordinate is on subdomain
+                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.  js >= nys  .AND.  js <= nyn )
+!
+!--             Check if observation coordinate is on subdomain.
                 IF ( on_pe )  THEN
 !
-!--                Determine vertical index which correspond to the observation
-!--                height.
+!--                Determine vertical index which corresponds to the observation height.
                    ksurf = topo_top_ind(js,is,0)
                    ks = MINLOC( ABS( zu - zw(ksurf) - zar(t,n) ), DIM = 1 ) - 1
 !
-!--                Set mask array at the observation coordinates. Also, flag the
-!--                surrounding coordinate points, but first check whether the
-!--                surrounding coordinate points are on the subdomain.
+!--                Set mask array at the observation coordinates. Also, flag the surrounding
+!--                coordinate points, but first check whether the surrounding coordinate points are
+!--                on the subdomain.
                    kl = MERGE( ks-off, ksurf, ks-off >= nzb  .AND. ks-off >= ksurf )
                    ku = MERGE( ks+off, nzt,   ks+off < nzt+1 )
@@ -1210,9 +1127,6 @@
                       DO  j = js-off, js+off
                          DO  k = kl, ku
-                            meas_flag(k,j,i) = MERGE(                           &
-                                          IBSET( meas_flag(k,j,i), 0 ),         &
-                                          0,                                    &
-                                          BTEST( wall_flags_total_0(k,j,i), 0 ) &
-                                                    )
+                            meas_flag(k,j,i) = MERGE( IBSET( meas_flag(k,j,i), 0 ), 0,             &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
                          ENDDO
                       ENDDO
@@ -1223,8 +1137,7 @@
           ENDDO
 !
-!--       Based on the flag array count the number of sampling coordinates.
-!--       Please note, sampling coordinates in atmosphere and soil may be
-!--       different, as within the soil all levels will be measured.
-!--       Hence, count individually. Start with atmoshere.
+!--       Based on the flag array, count the number of sampling coordinates. Please note, sampling
+!--       coordinates in atmosphere and soil may be different, as within the soil all levels will be
+!--       measured. Hence, count individually. Start with atmoshere.
           ns = 0
           DO  i = nxl-off, nxr+off
@@ -1237,6 +1150,6 @@
 
 !
-!--       Store number of observation points on subdomain and allocate index
-!--       arrays as well as array containing height information.
+!--       Store number of observation points on subdomain and allocate index arrays as well as array
+!--       containing height information.
           vmea(l)%ns = ns
 
@@ -1246,6 +1159,6 @@
           ALLOCATE( vmea(l)%zar(1:vmea(l)%ns) )
 !
-!--       Based on the flag array store the grid indices which correspond to
-!--       the observation coordinates.
+!--       Based on the flag array store the grid indices which correspond to the observation
+!--       coordinates.
           ns = 0
           DO  i = nxl-off, nxr+off
@@ -1263,21 +1176,16 @@
           ENDDO
 !
-!--       Same for the soil. Based on the flag array, count the number of
-!--       sampling coordinates in soil. Sample at all soil levels in this case.
-!--       Please note, soil variables can only be sampled on subdomains, not
-!--       on ghost layers.
+!--       Same for the soil. Based on the flag array, count the number of sampling coordinates in
+!--       soil. Sample at all soil levels in this case. Please note, soil variables can only be
+!--       sampled on subdomains, not on ghost layers.
           IF ( vmea(l)%soil_sampling )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
-                      IF ( surf_lsm_h%start_index(j,i) <=                      &
-                           surf_lsm_h%end_index(j,i) )  THEN
-                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
-                                                      nzt_soil - nzb_soil + 1
+                      IF ( surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) )  THEN
+                         vmea(l)%ns_soil = vmea(l)%ns_soil + nzt_soil - nzb_soil + 1
                       ENDIF
-                      IF ( surf_usm_h%start_index(j,i) <=                      &
-                           surf_usm_h%end_index(j,i) )  THEN
-                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
-                                                      nzt_wall - nzb_wall + 1
+                      IF ( surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) )  THEN
+                         vmea(l)%ns_soil = vmea(l)%ns_soil + nzt_wall - nzb_wall + 1
                       ENDIF
                    ENDIF
@@ -1286,6 +1194,5 @@
           ENDIF
 !
-!--       Allocate index arrays as well as array containing height information
-!--       for soil.
+!--       Allocate index arrays as well as array containing height information for soil.
           IF ( vmea(l)%soil_sampling )  THEN
              ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) )
@@ -1301,6 +1208,5 @@
                 DO  j = nys, nyn
                    IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
-                      IF ( surf_lsm_h%start_index(j,i) <=                      &
-                           surf_lsm_h%end_index(j,i) )  THEN
+                      IF ( surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) )  THEN
                          m = surf_lsm_h%start_index(j,i)
                          DO  k = nzb_soil, nzt_soil
@@ -1313,6 +1219,5 @@
                       ENDIF
 
-                      IF ( surf_usm_h%start_index(j,i) <=                      &
-                           surf_usm_h%end_index(j,i) )  THEN
+                      IF ( surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) )  THEN
                          m = surf_usm_h%start_index(j,i)
                          DO  k = nzb_wall, nzt_wall
@@ -1332,13 +1237,11 @@
           ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) )
 
-          IF ( vmea(l)%soil_sampling )                                         &
-             ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,           &
-                                                  1:vmea(l)%nmeas) )
+          IF ( vmea(l)%soil_sampling )                                                             &
+             ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil, 1:vmea(l)%nmeas) )
 !
 !--       Initialize with _FillValues
           vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) = vmea(l)%fillout
-          IF ( vmea(l)%soil_sampling )                                         &
-             vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nmeas) =   &
-                                                                vmea(l)%fillout
+          IF ( vmea(l)%soil_sampling )                                                             &
+             vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nmeas) = vmea(l)%fillout
 !
 !--       Deallocate temporary coordinate arrays
@@ -1366,6 +1269,6 @@
     ns_all = 0
 #if defined( __parallel )
-    CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER,  &
-                        MPI_SUM, comm2d, ierr )
+    CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm,                                   &
+                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
 #else
     ns_all(:) = vmea(:)%ns
@@ -1376,5 +1279,5 @@
     ns_all = 0
 #if defined( __parallel )
-    CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,          &
+    CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,                              &
                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
 #else
@@ -1385,6 +1288,6 @@
     DEALLOCATE( ns_all )
 !
-!-- In case of parallel NetCDF the start coordinate for each mpi rank needs to
-!-- be defined, so that each processor knows where to write the data.
+!-- In case of parallel NetCDF the start coordinate for each mpi rank needs to be defined, so that
+!-- each processor knows where to write the data.
 #if defined( __netcdf4_parallel )
     ALLOCATE( ns_atmos(0:numprocs-1,1:vmea_general%nvm) )
@@ -1399,7 +1302,7 @@
 
 #if defined( __parallel )
-    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_atmos, numprocs * vmea_general%nvm,   &
+    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_atmos, numprocs * vmea_general%nvm,                       &
                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
-    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_soil, numprocs * vmea_general%nvm,    &
+    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_soil, numprocs * vmea_general%nvm,                        &
                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
 #else
@@ -1409,8 +1312,7 @@
 
 !
-!-- Determine the start coordinate in NetCDF file for the local arrays.
-!-- Note, start coordinates are initialized with zero for sake of simplicity
-!-- in summation. However, in NetCDF the start coordinates must be >= 1,
-!-- so that a one needs to be added at the end.
+!-- Determine the start coordinate in NetCDF file for the local arrays. Note, start coordinates are
+!-- initialized with zero for sake of simplicity in summation. However, in NetCDF the start
+!-- coordinates must be >= 1, so that a one needs to be added at the end.
     DO  l = 1, vmea_general%nvm
        DO  n  = 0, myid - 1
@@ -1435,29 +1337,29 @@
 #endif
 
-  END SUBROUTINE vm_init
-
-
-!------------------------------------------------------------------------------!
+ END SUBROUTINE vm_init
+
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Initialize output using data-output module
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE vm_init_output
 
     CHARACTER(LEN=100) ::  variable_name  !< name of output variable
 
-    INTEGER(iwp) ::  l              !< loop index
-    INTEGER(iwp) ::  n              !< loop index
-    INTEGER      ::  return_value   !< returned status value of called function
+    INTEGER(iwp) ::  l             !< loop index
+    INTEGER(iwp) ::  n             !< loop index
+    INTEGER      ::  return_value  !< returned status value of called function
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ndim !< dummy to write dimension
 
-    REAL(wp) ::  dum_lat !< transformed geographical coordinate (latitude)
-    REAL(wp) ::  dum_lon !< transformed geographical coordinate (longitude)
+    REAL(wp) ::  dum_lat  !< transformed geographical coordinate (latitude)
+    REAL(wp) ::  dum_lon  !< transformed geographical coordinate (longitude)
 
 !
 !-- Determine the number of output timesteps.
-    ntimesteps = CEILING(                                                      &
-                  ( end_time - MAX( vm_time_start, time_since_reference_point )&
+    ntimesteps = CEILING(                                                                          &
+                  ( end_time - MAX( vm_time_start, time_since_reference_point )                    &
                   ) / dt_virtual_measurement )
 !
@@ -1472,7 +1374,6 @@
 !
 !--    Define output file.
-       WRITE( vmea(l)%nc_filename, '(A,I4.4)' ) 'VM_OUTPUT' //                 &
-                                                TRIM( coupling_char ) // '/' //&
-                                                'site', l
+       WRITE( vmea(l)%nc_filename, '(A,I4.4)' ) 'VM_OUTPUT' // TRIM( coupling_char ) // '/' //     &
+              'site', l
 
        return_value = dom_def_file( vmea(l)%nc_filename, 'netcdf4-parallel' )
@@ -1480,98 +1381,95 @@
 !--    Define global attributes.
 !--    Before, transform UTM into geographical coordinates.
-       CALL convert_utm_to_geographic( crs_list,                               &
-                                       vmea(l)%origin_x_obs,                   &
-                                       vmea(l)%origin_y_obs,                   &
-                                       dum_lon,                                &
-                                       dum_lat )
-
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'site',                    &
+       CALL convert_utm_to_geographic( crs_list, vmea(l)%origin_x_obs, vmea(l)%origin_y_obs,       &
+                                       dum_lon, dum_lat )
+
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'site',                                        &
                                    value = TRIM( vmea(l)%site ) )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'title',                   &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'title',                                       &
                                    value = 'Virtual measurement output')
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'source',                  &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'source',                                      &
                                    value = 'PALM-4U')
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'institution',             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'institution',                                 &
                                    value = input_file_atts%institution )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'acronym',                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'acronym',                                     &
                                    value = input_file_atts%acronym )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'author',                  &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'author',                                      &
                                    value = input_file_atts%author )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'contact_person',          &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'contact_person',                              &
                                    value = input_file_atts%contact_person )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'iop',                     &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'iop',                                         &
                                    value = input_file_atts%campaign )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'campaign',                &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'campaign',                                    &
                                    value = 'PALM-4U' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'origin_time ',            &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'origin_time ',                                &
                                    value = origin_date_time)
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'location',                &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'location',                                    &
                                    value = input_file_atts%location )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'origin_x',                &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'origin_x',                                    &
                                    value = vmea(l)%origin_x_obs )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'origin_y',                &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'origin_y',                                    &
                                    value = vmea(l)%origin_y_obs )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'origin_lon',              &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'origin_lon',                                  &
                                    value = dum_lon )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'origin_lat',              &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'origin_lat',                                  &
                                    value = dum_lat )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'origin_z',                &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'origin_z',                                    &
                                    value = 0.0 )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'rotation_angle',          &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'rotation_angle',                              &
                                    value = input_file_atts%rotation_angle )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'featureType',             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'featureType',                                 &
                                    value = TRIM( vmea(l)%feature_type_out ) )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'data_content',            &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'data_content',                                &
                                    value = TRIM( vmea(l)%data_content ) )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'creation_time',           &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'creation_time',                               &
                                    value = input_file_atts%creation_time )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'version',                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'version',                                     &
                                    value = 1 ) !input_file_atts%version )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'creation_time',           &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'creation_time',                               &
                                    value = TRIM( vmea(l)%site ) )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'Conventions',             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'Conventions',                                 &
                                    value = input_file_atts%conventions )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'dependencies',            &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'dependencies',                                &
                                    value = input_file_atts%dependencies )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'history',                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'history',                                     &
                                    value = input_file_atts%history )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'references',              &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'references',                                  &
                                    value = input_file_atts%references )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'comment',                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'comment',                                     &
                                    value = input_file_atts%comment )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'keywords',                &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'keywords',                                    &
                                    value = input_file_atts%keywords )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   attribute_name = 'licence',                 &
-                                   value = '[UC]2 Open Licence; see [UC]2 ' // &
-                                           'data policy available at ' //      &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   attribute_name = 'licence',                                     &
+                                   value = '[UC]2 Open Licence; see [UC]2 ' //                     &
+                                           'data policy available at ' //                          &
                                            'www.uc2-program.org/uc2_data_policy.pdf' )
 !
@@ -1582,8 +1480,6 @@
           ndim(n) = n
        ENDDO
-       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
-                                   dimension_name = 'station',                 &
-                                   output_type = 'int32',                      &
-                                   bounds = (/1_iwp, vmea(l)%ns_tot/),         &
+       return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'station',                &
+                                   output_type = 'int32', bounds = (/1_iwp, vmea(l)%ns_tot/),      &
                                    values_int32 = ndim )
        DEALLOCATE( ndim )
@@ -1595,8 +1491,6 @@
        ENDDO
 
-       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
-                                   dimension_name = 'ntime',                   &
-                                   output_type = 'int32',                      &
-                                   bounds = (/1_iwp, ntimesteps/),             &
+       return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'ntime',                  &
+                                   output_type = 'int32', bounds = (/1_iwp, ntimesteps/),          &
                                    values_int32 = ndim )
        DEALLOCATE( ndim )
@@ -1608,8 +1502,6 @@
        ENDDO
 
-       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
-                                   dimension_name = 'nv',                      &
-                                   output_type = 'int32',                      &
-                                   bounds = (/1_iwp, 2_iwp/),                  &
+       return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'nv',                     &
+                                   output_type = 'int32', bounds = (/1_iwp, 2_iwp/),               &
                                    values_int32 = ndim )
        DEALLOCATE( ndim )
@@ -1621,9 +1513,7 @@
        ENDDO
 
-       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
-                                   dimension_name = 'max_name_len',            &
-                                   output_type = 'int32',                      &
-                                   bounds = (/1_iwp, maximum_name_length /),   &
-                                   values_int32 = ndim )
+       return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'max_name_len',           &
+                                   output_type = 'int32',                                          &
+                                   bounds = (/1_iwp, maximum_name_length /), values_int32 = ndim )
        DEALLOCATE( ndim )
 !
@@ -1631,307 +1521,282 @@
 !--    time
        variable_name = 'time'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/ 'station  ',           &
-                                                        'ntime    '/),         &
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/ 'station  ', 'ntime    '/),                &
                                    output_type = 'real32' )
 !
 !--    station_name
        variable_name = 'station_name'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/ 'max_name_len',        &
-                                                        'station     ' /),     &
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/ 'max_name_len', 'station     ' /),         &
                                    output_type = 'char' )
 !
 !--    vrs (vertical reference system)
        variable_name = 'vrs'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/ 'station' /),          &
-                                   output_type = 'int8' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/ 'station' /), output_type = 'int8' )
 !
 !--    crs (coordinate reference system)
        variable_name = 'crs'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/ 'station' /),          &
-                                   output_type = 'int8' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/ 'station' /), output_type = 'int8' )
 !
 !--    z
        variable_name = 'z'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/'station'/),            &
-                                   output_type = 'real32' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/'station'/), output_type = 'real32' )
 !
 !--    station_h
        variable_name = 'station_h'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/'station'/),            &
-                                   output_type = 'real32' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/'station'/), output_type = 'real32' )
 !
 !--    x
        variable_name = 'x'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/'station'/),            &
-                                   output_type = 'real32' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/'station'/), output_type = 'real32' )
 !
 !--    y
        variable_name = 'y'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/'station'/),            &
-                                   output_type = 'real32' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/'station'/), output_type = 'real32' )
 !
 !--    E-UTM
        variable_name = 'E_UTM'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/'station'/),            &
-                                   output_type = 'real32' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/'station'/), output_type = 'real32' )
 !
 !--    N-UTM
        variable_name = 'N_UTM'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/'station'/),            &
-                                   output_type = 'real32' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/'station'/), output_type = 'real32' )
 !
 !--    latitude
        variable_name = 'lat'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/'station'/),            &
-                                   output_type = 'real32' )
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/'station'/), output_type = 'real32' )
 !
 !--    longitude
        variable_name = 'lon'
-       return_value = dom_def_var( vmea(l)%nc_filename,                        &
-                                   variable_name = variable_name,              &
-                                   dimension_names = (/'station'/),            &
-                                   output_type = 'real32' )
-!
-!--    Set attributes for the coordinate variables. Note, not all coordinates
-!--    have the same number of attributes.
+       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
+                                   dimension_names = (/'station'/), output_type = 'real32' )
+!
+!--    Set attributes for the coordinate variables. Note, not all coordinates have the same number
+!--    of attributes.
 !--    Units
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'time',                     &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'time',                                         &
+                                   attribute_name = char_unit,                                     &
                                    value = 'seconds since ' // origin_date_time )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'z',                        &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'z',                                            &
+                                   attribute_name = char_unit,                                     &
                                    value = 'm' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'station_h',                &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'station_h',                                    &
+                                   attribute_name = char_unit,                                     &
                                    value = 'm' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'x',                        &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'x',                                            &
+                                   attribute_name = char_unit,                                     &
                                    value = 'm' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'y',                        &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'y',                                            &
+                                   attribute_name = char_unit,                                     &
                                    value = 'm' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'E_UTM',                    &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'E_UTM',                                        &
+                                   attribute_name = char_unit,                                     &
                                    value = 'm' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'N_UTM',                    &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'N_UTM',                                        &
+                                   attribute_name = char_unit,                                     &
                                    value = 'm' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'lat',                      &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'lat',                                          &
+                                   attribute_name = char_unit,                                     &
                                    value = 'degrees_north' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'lon',                      &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'lon',                                          &
+                                   attribute_name = char_unit,                                     &
                                    value = 'degrees_east' )
 !
 !--    Long name
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'station_name',             &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'station_name',                                 &
+                                   attribute_name = char_long,                                     &
                                    value = 'station name')
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'time',                     &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'time',                                         &
+                                   attribute_name = char_long,                                     &
                                    value = 'time')
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'z',                        &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'z',                                            &
+                                   attribute_name = char_long,                                     &
                                    value = 'height above origin' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'station_h',                &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'station_h',                                    &
+                                   attribute_name = char_long,                                     &
                                    value = 'surface altitude' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'x',                        &
-                                   attribute_name = char_long,                 &
-                                   value = 'distance to origin in x-direction' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'y',                        &
-                                   attribute_name = char_long,                 &
-                                   value = 'distance to origin in y-direction' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'E_UTM',                    &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'x',                                            &
+                                   attribute_name = char_long,                                     &
+                                   value = 'distance to origin in x-direction')
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'y',                                            &
+                                   attribute_name = char_long,                                     &
+                                   value = 'distance to origin in y-direction')
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'E_UTM',                                        &
+                                   attribute_name = char_long,                                     &
                                    value = 'easting' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'N_UTM',                    &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'N_UTM',                                        &
+                                   attribute_name = char_long,                                     &
                                    value = 'northing' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'lat',                      &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'lat',                                          &
+                                   attribute_name = char_long,                                     &
                                    value = 'latitude' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'lon',                      &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'lon',                                          &
+                                   attribute_name = char_long,                                     &
                                    value = 'longitude' )
 !
 !--    Standard name
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'station_name',             &
-                                   attribute_name = char_standard,             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'station_name',                                 &
+                                   attribute_name = char_standard,                                 &
                                    value = 'platform_name')
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'time',                     &
-                                   attribute_name = char_standard,             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'time',                                         &
+                                   attribute_name = char_standard,                                 &
                                    value = 'time')
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'z',                        &
-                                   attribute_name = char_standard,             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'z',                                            &
+                                   attribute_name = char_standard,                                 &
                                    value = 'height_above_mean_sea_level' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'station_h',                &
-                                   attribute_name = char_standard,             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'station_h',                                    &
+                                   attribute_name = char_standard,                                 &
                                    value = 'surface_altitude' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'E_UTM',                    &
-                                   attribute_name = char_standard,             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'E_UTM',                                        &
+                                   attribute_name = char_standard,                                 &
                                    value = 'projection_x_coordinate' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'N_UTM',                    &
-                                   attribute_name = char_standard,             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'N_UTM',                                        &
+                                   attribute_name = char_standard,                                 &
                                    value = 'projection_y_coordinate' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'lat',                      &
-                                   attribute_name = char_standard,             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'lat',                                          &
+                                   attribute_name = char_standard,                                 &
                                    value = 'latitude' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'lon',                      &
-                                   attribute_name = char_standard,             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'lon',                                          &
+                                   attribute_name = char_standard,                                 &
                                    value = 'longitude' )
 !
 !--    Axis
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'time',                     &
-                                   attribute_name = 'axis',                    &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'time',                                         &
+                                   attribute_name = 'axis',                                        &
                                    value = 'T')
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'z',                        &
-                                   attribute_name = 'axis',                    &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'z',                                            &
+                                   attribute_name = 'axis',                                        &
                                    value = 'Z' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'x',                        &
-                                   attribute_name = 'axis',                    &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'x',                                            &
+                                   attribute_name = 'axis',                                        &
                                    value = 'X' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'y',                        &
-                                   attribute_name = 'axis',                    &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'y',                                            &
+                                   attribute_name = 'axis',                                        &
                                    value = 'Y' )
 !
 !--    Set further individual attributes for the coordinate variables.
 !--    For station name
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'station_name',             &
-                                   attribute_name = 'cf_role',                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'station_name',                                 &
+                                   attribute_name = 'cf_role',                                     &
                                    value = 'timeseries_id' )
 !
 !--    For time
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'time',                     &
-                                   attribute_name = 'calendar',                &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'time',                                         &
+                                   attribute_name = 'calendar',                                    &
                                    value = 'proleptic_gregorian' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'time',                     &
-                                   attribute_name = 'bounds',                  &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'time',                                         &
+                                   attribute_name = 'bounds',                                      &
                                    value = 'time_bounds' )
 !
 !--    For vertical reference system
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'vrs',                      &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'vrs',                                          &
+                                   attribute_name = char_long,                                     &
                                    value = 'vertical reference system' )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'vrs',                      &
-                                   attribute_name = 'system_name',             &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'vrs',                                          &
+                                   attribute_name = 'system_name',                                 &
                                    value = 'DHHN2016' )
 !
 !--    For z
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'z',                        &
-                                   attribute_name = 'positive',                &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'z',                                            &
+                                   attribute_name = 'positive',                                    &
                                    value = 'up' )
 !
 !--    For coordinate reference system
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'epsg_code',               &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'epsg_code',                                   &
                                    value = coord_ref_sys%epsg_code )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'false_easting',           &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'false_easting',                               &
                                    value = coord_ref_sys%false_easting )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'false_northing',          &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'false_northing',                              &
                                    value = coord_ref_sys%false_northing )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'grid_mapping_name',       &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'grid_mapping_name',                           &
                                    value = coord_ref_sys%grid_mapping_name )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'inverse_flattening',      &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'inverse_flattening',                          &
                                    value = coord_ref_sys%inverse_flattening )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
                                    attribute_name = 'latitude_of_projection_origin',&
                                    value = coord_ref_sys%latitude_of_projection_origin )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = char_long,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = char_long,                                     &
                                    value = coord_ref_sys%long_name )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'longitude_of_central_meridian', &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'longitude_of_central_meridian',               &
                                    value = coord_ref_sys%longitude_of_central_meridian )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'longitude_of_prime_meridian', &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'longitude_of_prime_meridian',                 &
                                    value = coord_ref_sys%longitude_of_prime_meridian )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'scale_factor_at_central_meridian', &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'scale_factor_at_central_meridian',            &
                                    value = coord_ref_sys%scale_factor_at_central_meridian )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = 'semi_major_axis',         &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = 'semi_major_axis',                             &
                                    value = coord_ref_sys%semi_major_axis )
-       return_value = dom_def_att( vmea(l)%nc_filename,                        &
-                                   variable_name = 'crs',                      &
-                                   attribute_name = char_unit,                 &
+       return_value = dom_def_att( vmea(l)%nc_filename,                                            &
+                                   variable_name = 'crs',                                          &
+                                   attribute_name = char_unit,                                     &
                                    value = coord_ref_sys%units )
 !
-!--    In case of sampled soil quantities, define further dimensions and
-!--    coordinates.
+!--    In case of sampled soil quantities, define further dimensions and coordinates.
        IF ( vmea(l)%soil_sampling )  THEN
 !
@@ -1942,8 +1807,8 @@
           ENDDO
 
-          return_value = dom_def_dim( vmea(l)%nc_filename,                     &
-                                      dimension_name = 'station_soil',         &
-                                      output_type = 'int32',                   &
-                                      bounds = (/1_iwp,vmea(l)%ns_soil_tot/),  &
+          return_value = dom_def_dim( vmea(l)%nc_filename,                                         &
+                                      dimension_name = 'station_soil',                             &
+                                      output_type = 'int32',                                       &
+                                      bounds = (/1_iwp,vmea(l)%ns_soil_tot/),                      &
                                       values_int32 = ndim )
           DEALLOCATE( ndim )
@@ -1955,8 +1820,8 @@
           ENDDO
 
-          return_value = dom_def_dim( vmea(l)%nc_filename,                     &
-                                      dimension_name = 'ntime_soil',           &
-                                      output_type = 'int32',                   &
-                                      bounds = (/1_iwp,ntimesteps/),           &
+          return_value = dom_def_dim( vmea(l)%nc_filename,                                         &
+                                      dimension_name = 'ntime_soil',                               &
+                                      output_type = 'int32',                                       &
+                                      bounds = (/1_iwp,ntimesteps/),                               &
                                       values_int32 = ndim )
           DEALLOCATE( ndim )
@@ -1964,229 +1829,229 @@
 !--       time for soil
           variable_name = 'time_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil',      &
-                                                          'ntime_soil  '/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil',                          &
+                                                          'ntime_soil  '/),                        &
                                       output_type = 'real32' )
 !
 !--       station_name for soil
           variable_name = 'station_name_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/ 'max_name_len',     &
-                                                           'station_soil' /),  &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/ 'max_name_len',                         &
+                                                           'station_soil' /),                      &
                                       output_type = 'char' )
 !
 !--       z
           variable_name = 'z_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil'/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil'/),                        &
                                       output_type = 'real32' )
 !
 !--       station_h for soil
           variable_name = 'station_h_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil'/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil'/),                        &
                                       output_type = 'real32' )
 !
 !--       x soil
           variable_name = 'x_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil'/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil'/),                        &
                                       output_type = 'real32' )
 !
 !-        y soil
           variable_name = 'y_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil'/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil'/),                        &
                                       output_type = 'real32' )
 !
 !--       E-UTM soil
           variable_name = 'E_UTM_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil'/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil'/),                        &
                                       output_type = 'real32' )
 !
 !--       N-UTM soil
           variable_name = 'N_UTM_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil'/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil'/),                        &
                                       output_type = 'real32' )
 !
 !--       latitude soil
           variable_name = 'lat_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil'/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil'/),                        &
                                       output_type = 'real32' )
 !
 !--       longitude soil
           variable_name = 'lon_soil'
-          return_value = dom_def_var( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      dimension_names = (/'station_soil'/),    &
+          return_value = dom_def_var( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      dimension_names = (/'station_soil'/),                        &
                                       output_type = 'real32' )
 !
-!--       Set attributes for the coordinate variables. Note, not all coordinates
-!--       have the same number of attributes.
+!--       Set attributes for the coordinate variables. Note, not all coordinates have the same
+!--       number of attributes.
 !--       Units
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'time_soil',             &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'time_soil',                                 &
+                                      attribute_name = char_unit,                                  &
                                       value = 'seconds since ' // origin_date_time )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'z_soil',                &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'z_soil',                                    &
+                                      attribute_name = char_unit,                                  &
                                       value = 'm' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'station_h_soil',        &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'station_h_soil',                            &
+                                      attribute_name = char_unit,                                  &
                                       value = 'm' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'x_soil',                &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'x_soil',                                    &
+                                      attribute_name = char_unit,                                  &
                                       value = 'm' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'y_soil',                &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'y_soil',                                    &
+                                      attribute_name = char_unit,                                  &
                                       value = 'm' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'E_UTM_soil',            &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'E_UTM_soil',                                &
+                                      attribute_name = char_unit,                                  &
                                       value = 'm' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'N_UTM_soil',            &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'N_UTM_soil',                                &
+                                      attribute_name = char_unit,                                  &
                                       value = 'm' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'lat_soil',              &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'lat_soil',                                  &
+                                      attribute_name = char_unit,                                  &
                                       value = 'degrees_north' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'lon_soil',              &
-                                      attribute_name = char_unit,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'lon_soil',                                  &
+                                      attribute_name = char_unit,                                  &
                                       value = 'degrees_east' )
 !
 !--       Long name
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'station_name_soil',     &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'station_name_soil',                         &
+                                      attribute_name = char_long,                                  &
                                       value = 'station name')
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'time_soil',             &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'time_soil',                                 &
+                                      attribute_name = char_long,                                  &
                                       value = 'time')
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'z_soil',                &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'z_soil',                                    &
+                                      attribute_name = char_long,                                  &
                                       value = 'height above origin' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'station_h_soil',        &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'station_h_soil',                            &
+                                      attribute_name = char_long,                                  &
                                       value = 'surface altitude' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'x_soil',                &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'x_soil',                                    &
+                                      attribute_name = char_long,                                  &
                                       value = 'distance to origin in x-direction' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'y_soil',                &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'y_soil',                                    &
+                                      attribute_name = char_long,                                  &
                                       value = 'distance to origin in y-direction' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'E_UTM_soil',            &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'E_UTM_soil',                                &
+                                      attribute_name = char_long,                                  &
                                       value = 'easting' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'N_UTM_soil',            &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'N_UTM_soil',                                &
+                                      attribute_name = char_long,                                  &
                                       value = 'northing' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'lat_soil',              &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'lat_soil',                                  &
+                                      attribute_name = char_long,                                  &
                                       value = 'latitude' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'lon_soil',              &
-                                      attribute_name = char_long,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'lon_soil',                                  &
+                                      attribute_name = char_long,                                  &
                                       value = 'longitude' )
 !
 !--       Standard name
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'station_name_soil',     &
-                                      attribute_name = char_standard,          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'station_name_soil',                         &
+                                      attribute_name = char_standard,                              &
                                       value = 'platform_name')
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'time_soil',             &
-                                      attribute_name = char_standard,          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'time_soil',                                 &
+                                      attribute_name = char_standard,                              &
                                       value = 'time')
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'z_soil',                &
-                                      attribute_name = char_standard,          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'z_soil',                                    &
+                                      attribute_name = char_standard,                              &
                                       value = 'height_above_mean_sea_level' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'station_h_soil',        &
-                                      attribute_name = char_standard,          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'station_h_soil',                            &
+                                      attribute_name = char_standard,                              &
                                       value = 'surface_altitude' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'E_UTM_soil',            &
-                                      attribute_name = char_standard,          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'E_UTM_soil',                                &
+                                      attribute_name = char_standard,                              &
                                       value = 'projection_x_coordinate' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'N_UTM_soil',            &
-                                      attribute_name = char_standard,          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'N_UTM_soil',                                &
+                                      attribute_name = char_standard,                              &
                                       value = 'projection_y_coordinate' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'lat_soil',              &
-                                      attribute_name = char_standard,          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'lat_soil',                                  &
+                                      attribute_name = char_standard,                              &
                                       value = 'latitude' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'lon_soil',              &
-                                      attribute_name = char_standard,          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'lon_soil',                                  &
+                                      attribute_name = char_standard,                              &
                                       value = 'longitude' )
 !
 !--       Axis
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'time_soil',             &
-                                      attribute_name = 'axis',                 &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'time_soil',                                 &
+                                      attribute_name = 'axis',                                     &
                                       value = 'T')
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'z_soil',                &
-                                      attribute_name = 'axis',                 &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'z_soil',                                    &
+                                      attribute_name = 'axis',                                     &
                                       value = 'Z' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'x_soil',                &
-                                      attribute_name = 'axis',                 &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'x_soil',                                    &
+                                      attribute_name = 'axis',                                     &
                                       value = 'X' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'y_soil',                &
-                                      attribute_name = 'axis',                 &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'y_soil',                                    &
+                                      attribute_name = 'axis',                                     &
                                       value = 'Y' )
 !
 !--       Set further individual attributes for the coordinate variables.
 !--       For station name soil
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'station_name_soil',     &
-                                      attribute_name = 'cf_role',              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'station_name_soil',                         &
+                                      attribute_name = 'cf_role',                                  &
                                       value = 'timeseries_id' )
 !
 !--       For time soil
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'time_soil',             &
-                                      attribute_name = 'calendar',             &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'time_soil',                                 &
+                                      attribute_name = 'calendar',                                 &
                                       value = 'proleptic_gregorian' )
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'time_soil',             &
-                                      attribute_name = 'bounds',               &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'time_soil',                                 &
+                                      attribute_name = 'bounds',                                   &
                                       value = 'time_bounds' )
 !
 !--       For z soil
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = 'z_soil',                &
-                                      attribute_name = 'positive',             &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = 'z_soil',                                    &
+                                      attribute_name = 'positive',                                 &
                                       value = 'up' )
        ENDIF
@@ -2196,60 +2061,59 @@
           variable_name = TRIM( vmea(l)%var_atts(n)%name )
 !
-!--       In order to link the correct dimension names, atmosphere and soil
-!--       variables need to be distinguished.
-          IF ( vmea(l)%soil_sampling  .AND.                                    &
+!--       In order to link the correct dimension names, atmosphere and soil variables need to be
+!--       distinguished.
+          IF ( vmea(l)%soil_sampling  .AND.                                                        &
                ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
 
-             return_value = dom_def_var( vmea(l)%nc_filename,                  &
-                                         variable_name = variable_name,        &
-                                         dimension_names = (/'station_soil',   &
-                                                             'ntime_soil  '/), &
+             return_value = dom_def_var( vmea(l)%nc_filename,                                      &
+                                         variable_name = variable_name,                            &
+                                         dimension_names = (/'station_soil',                       &
+                                                             'ntime_soil  '/),                     &
                                          output_type = 'real32' )
           ELSE
 
-             return_value = dom_def_var( vmea(l)%nc_filename,                  &
-                                         variable_name = variable_name,        &
-                                         dimension_names = (/'station',        &
-                                                             'ntime  '/),      &
+             return_value = dom_def_var( vmea(l)%nc_filename,                                      &
+                                         variable_name = variable_name,                            &
+                                         dimension_names = (/'station',                            &
+                                                             'ntime  '/),                          &
                                          output_type = 'real32' )
           ENDIF
 !
-!--       Set variable attributes. Please note, for some variables not all
-!--       attributes are defined, e.g. standard_name for the horizontal wind
-!--       components.
+!--       Set variable attributes. Please note, for some variables not all attributes are defined,
+!--       e.g. standard_name for the horizontal wind components.
           CALL vm_set_attributes( vmea(l)%var_atts(n) )
 
           IF ( vmea(l)%var_atts(n)%long_name /= 'none' )  THEN
-             return_value = dom_def_att( vmea(l)%nc_filename,                  &
-                                         variable_name = variable_name,        &
-                                         attribute_name = char_long,           &
+             return_value = dom_def_att( vmea(l)%nc_filename,                                      &
+                                         variable_name = variable_name,                            &
+                                         attribute_name = char_long,                               &
                                          value = TRIM( vmea(l)%var_atts(n)%long_name ) )
           ENDIF
           IF ( vmea(l)%var_atts(n)%standard_name /= 'none' )  THEN
-             return_value = dom_def_att( vmea(l)%nc_filename,                  &
-                                         variable_name = variable_name,        &
-                                         attribute_name = char_standard,       &
+             return_value = dom_def_att( vmea(l)%nc_filename,                                      &
+                                         variable_name = variable_name,                            &
+                                         attribute_name = char_standard,                           &
                                          value = TRIM( vmea(l)%var_atts(n)%standard_name ) )
           ENDIF
           IF ( vmea(l)%var_atts(n)%units /= 'none' )  THEN
-             return_value = dom_def_att( vmea(l)%nc_filename,                  &
-                                         variable_name = variable_name,        &
-                                         attribute_name = char_unit,           &
+             return_value = dom_def_att( vmea(l)%nc_filename,                                      &
+                                         variable_name = variable_name,                            &
+                                         attribute_name = char_unit,                               &
                                          value = TRIM( vmea(l)%var_atts(n)%units ) )
           ENDIF
 
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      attribute_name = 'grid_mapping',         &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      attribute_name = 'grid_mapping',                             &
                                       value = TRIM( vmea(l)%var_atts(n)%grid_mapping ) )
 
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      attribute_name = 'coordinates',          &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      attribute_name = 'coordinates',                              &
                                       value = TRIM( vmea(l)%var_atts(n)%coordinates ) )
 
-          return_value = dom_def_att( vmea(l)%nc_filename,                     &
-                                      variable_name = variable_name,           &
-                                      attribute_name = char_fill,              &
+          return_value = dom_def_att( vmea(l)%nc_filename,                                         &
+                                      variable_name = variable_name,                               &
+                                      attribute_name = char_fill,                                  &
                                       value = REAL( vmea(l)%var_atts(n)%fill_value, KIND=4 ) )
 
@@ -2261,13 +2125,13 @@
  END SUBROUTINE vm_init_output
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Parallel NetCDF output via data-output module.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE vm_data_output
 
-    CHARACTER(LEN=100) ::  variable_name !< name of output variable
-    CHARACTER(LEN=maximum_name_length), DIMENSION(:), ALLOCATABLE :: station_name !< string for station name, consecutively ordered
+    CHARACTER(LEN=100) ::  variable_name  !< name of output variable
+    CHARACTER(LEN=maximum_name_length), DIMENSION(:), ALLOCATABLE :: station_name  !< string for station name, consecutively ordered
 
     CHARACTER(LEN=1), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_char_target  !< target for output name arrays
@@ -2278,12 +2142,11 @@
     INTEGER(iwp)       ::  nn            !< loop index for number of characters in a name
     INTEGER            ::  return_value  !< returned status value of called function
-
     INTEGER(iwp)       ::  t_ind         !< time index
 
-    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  oro_rel                  !< relative altitude of model surface
-    REAL(wp), DIMENSION(:), POINTER               ::  output_values_1d_pointer !< pointer for 1d output array
-    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET   ::  output_values_1d_target  !< target for 1d output array
-    REAL(wp), DIMENSION(:,:), POINTER             ::  output_values_2d_pointer !< pointer for 2d output array
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_target  !< target for 2d output array
+    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  oro_rel                   !< relative altitude of model surface
+    REAL(wp), DIMENSION(:), POINTER               ::  output_values_1d_pointer  !< pointer for 1d output array
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET   ::  output_values_1d_target   !< target for 1d output array
+    REAL(wp), DIMENSION(:,:), POINTER             ::  output_values_2d_pointer  !< pointer for 2d output array
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_target   !< target for 2d output array
 
     CALL cpu_log( log_point_s(26), 'VM output', 'start' )
@@ -2301,37 +2164,36 @@
 !
 !--       Output of Easting coordinate. Before output, recalculate EUTM.
-          output_values_1d_target = init_model%origin_x                        &
-                    + REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
-                    * COS( init_model%rotation_angle * pi / 180.0_wp )         &
-                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
+          output_values_1d_target = init_model%origin_x                                            &
+                    + REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx                     &
+                    * COS( init_model%rotation_angle * pi / 180.0_wp )                             &
+                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy                     &
                     * SIN( init_model%rotation_angle * pi / 180.0_wp )
 
           output_values_1d_pointer => output_values_1d_target
 
-          return_value =                                                       &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'E_UTM',                                      &
-                                 values_realwp_1d = output_values_1d_pointer,  &
-                                 bounds_start = (/vmea(l)%start_coord_a/),     &
+          return_value =                                                                           &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'E_UTM',                                                          &
+                                 values_realwp_1d = output_values_1d_pointer,                      &
+                                 bounds_start = (/vmea(l)%start_coord_a/),                         &
                                  bounds_end   = (/vmea(l)%end_coord_a  /) )
 !
 !--       Output of Northing coordinate. Before output, recalculate NUTM.
-          output_values_1d_target = init_model%origin_y                        &
-                    - REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
-                    * SIN( init_model%rotation_angle * pi / 180.0_wp )         &
-                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
+          output_values_1d_target = init_model%origin_y                                            &
+                    - REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx                     &
+                    * SIN( init_model%rotation_angle * pi / 180.0_wp )                             &
+                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy                     &
                     * COS( init_model%rotation_angle * pi / 180.0_wp )
 
           output_values_1d_pointer => output_values_1d_target
-          return_value =                                                       &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'N_UTM',                                      &
-                                 values_realwp_1d = output_values_1d_pointer,  &
-                                 bounds_start = (/vmea(l)%start_coord_a/),     &
+          return_value =                                                                           &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'N_UTM',                                                          &
+                                 values_realwp_1d = output_values_1d_pointer,                      &
+                                 bounds_start = (/vmea(l)%start_coord_a/),                         &
                                  bounds_end   = (/vmea(l)%end_coord_a  /) )
 !
 !--       Output of relative height coordinate.
-!--       Before this is output, first define the relative orographie height
-!--       and add this to z.
+!--       Before this is output, first define the relative orographie height and add this to z.
           ALLOCATE( oro_rel(1:vmea(l)%ns) )
           DO  n = 1, vmea(l)%ns
@@ -2341,21 +2203,20 @@
           output_values_1d_target = vmea(l)%zar(1:vmea(l)%ns) + oro_rel(:)
           output_values_1d_pointer => output_values_1d_target
-          return_value =                                                       &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'z',                                          &
-                                 values_realwp_1d = output_values_1d_pointer,  &
-                                 bounds_start = (/vmea(l)%start_coord_a/),     &
+          return_value =                                                                           &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'z',                                                              &
+                                 values_realwp_1d = output_values_1d_pointer,                      &
+                                 bounds_start = (/vmea(l)%start_coord_a/),                         &
                                  bounds_end   = (/vmea(l)%end_coord_a  /) )
 !
-!--       Write surface altitude for the station. Note, since z is already
-!--       a relative observation height, station_h must be zero, in order
-!--       to obtain the observation level.
+!--       Write surface altitude for the station. Note, since z is already a relative observation
+!--       height, station_h must be zero, in order to obtain the observation level.
           output_values_1d_target = oro_rel(:)
           output_values_1d_pointer => output_values_1d_target
-          return_value =                                                       &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'station_h',                                  &
-                                 values_realwp_1d = output_values_1d_pointer,  &
-                                 bounds_start = (/vmea(l)%start_coord_a/),     &
+          return_value =                                                                           &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'station_h',                                                      &
+                                 values_realwp_1d = output_values_1d_pointer,                      &
+                                 bounds_start = (/vmea(l)%start_coord_a/),                         &
                                  bounds_end   = (/vmea(l)%end_coord_a  /) )
 
@@ -2378,9 +2239,9 @@
           output_values_2d_char_pointer => output_values_2d_char_target
 
-          return_value =                                                       &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'station_name',                               &
-                                 values_char_2d = output_values_2d_char_pointer,&
-                                 bounds_start = (/ 1,                   vmea(l)%start_coord_a /),&
+          return_value =                                                                           &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'station_name',                                                   &
+                                 values_char_2d = output_values_2d_char_pointer,                   &
+                                 bounds_start = (/ 1,                   vmea(l)%start_coord_a /),  &
                                  bounds_end   = (/ maximum_name_length, vmea(l)%end_coord_a /) )
 
@@ -2388,41 +2249,39 @@
           DEALLOCATE( output_values_2d_char_target )
 !
-!--       In case of sampled soil quantities, output also the respective
-!--       coordinate arrays.
+!--       In case of sampled soil quantities, output also the respective coordinate arrays.
           IF ( vmea(l)%soil_sampling )  THEN
              ALLOCATE( output_values_1d_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
 !
 !--          Output of Easting coordinate. Before output, recalculate EUTM.
-             output_values_1d_target = init_model%origin_x                     &
-               + REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
-               * COS( init_model%rotation_angle * pi / 180.0_wp )              &
-               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
+             output_values_1d_target = init_model%origin_x                                         &
+               + REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx                     &
+               * COS( init_model%rotation_angle * pi / 180.0_wp )                                  &
+               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy                     &
                * SIN( init_model%rotation_angle * pi / 180.0_wp )
              output_values_1d_pointer => output_values_1d_target
-             return_value =                                                    &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'E_UTM_soil',                                 &
-                                 values_realwp_1d = output_values_1d_pointer,  &
-                                 bounds_start = (/vmea(l)%start_coord_s/),     &
+             return_value =                                                                        &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'E_UTM_soil',                                                     &
+                                 values_realwp_1d = output_values_1d_pointer,                      &
+                                 bounds_start = (/vmea(l)%start_coord_s/),                         &
                                  bounds_end   = (/vmea(l)%end_coord_s  /) )
 !
 !--          Output of Northing coordinate. Before output, recalculate NUTM.
-             output_values_1d_target = init_model%origin_y                     &
-               - REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
-               * SIN( init_model%rotation_angle * pi / 180.0_wp )              &
-               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
+             output_values_1d_target = init_model%origin_y                                         &
+               - REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx                     &
+               * SIN( init_model%rotation_angle * pi / 180.0_wp )                                  &
+               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy                     &
                * COS( init_model%rotation_angle * pi / 180.0_wp )
 
              output_values_1d_pointer => output_values_1d_target
-             return_value =                                                    &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'N_UTM_soil',                                 &
-                                 values_realwp_1d = output_values_1d_pointer,  &
-                                 bounds_start = (/vmea(l)%start_coord_s/),     &
+             return_value =                                                                        &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'N_UTM_soil',                                                     &
+                                 values_realwp_1d = output_values_1d_pointer,                      &
+                                 bounds_start = (/vmea(l)%start_coord_s/),                         &
                                  bounds_end   = (/vmea(l)%end_coord_s  /) )
 !
 !--          Output of relative height coordinate.
-!--          Before this is output, first define the relative orographie height
-!--          and add this to z.
+!--          Before this is output, first define the relative orographie height and add this to z.
              ALLOCATE( oro_rel(1:vmea(l)%ns_soil) )
              DO  n = 1, vmea(l)%ns_soil
@@ -2432,21 +2291,20 @@
              output_values_1d_target = vmea(l)%depth(1:vmea(l)%ns_soil) + oro_rel(:)
              output_values_1d_pointer => output_values_1d_target
-             return_value =                                                    &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'z_soil',                                     &
-                                 values_realwp_1d = output_values_1d_pointer,  &
-                                 bounds_start = (/vmea(l)%start_coord_s/),     &
+             return_value =                                                                        &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'z_soil',                                                         &
+                                 values_realwp_1d = output_values_1d_pointer,                      &
+                                 bounds_start = (/vmea(l)%start_coord_s/),                         &
                                  bounds_end   = (/vmea(l)%end_coord_s  /) )
 !
-!--          Write surface altitude for the station. Note, since z is already
-!--          a relative observation height, station_h must be zero, in order
-!--          to obtain the observation level.
+!--          Write surface altitude for the station. Note, since z is already a relative observation
+!--          height, station_h must be zero, in order to obtain the observation level.
              output_values_1d_target = oro_rel(:)
              output_values_1d_pointer => output_values_1d_target
-             return_value =                                                    &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'station_h_soil',                             &
-                                 values_realwp_1d = output_values_1d_pointer,  &
-                                 bounds_start = (/vmea(l)%start_coord_s/),     &
+             return_value =                                                                        &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'station_h_soil',                                                 &
+                                 values_realwp_1d = output_values_1d_pointer,                      &
+                                 bounds_start = (/vmea(l)%start_coord_s/),                         &
                                  bounds_end   = (/vmea(l)%end_coord_s  /) )
 
@@ -2456,5 +2314,5 @@
 !--          Write station name
              ALLOCATE ( station_name(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
-             ALLOCATE ( output_values_2d_char_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s, &
+             ALLOCATE ( output_values_2d_char_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s,    &
                                                      1:maximum_name_length) )
 
@@ -2468,9 +2326,9 @@
              output_values_2d_char_pointer => output_values_2d_char_target
 
-             return_value =                                                    &
-                  dom_write_var( vmea(l)%nc_filename,                          &
-                                 'station_name_soil',                          &
-                                 values_char_2d = output_values_2d_char_pointer,&
-                                 bounds_start = (/ 1,                   vmea(l)%start_coord_s /),&
+             return_value =                                                                        &
+                  dom_write_var( vmea(l)%nc_filename,                                              &
+                                 'station_name_soil',                                              &
+                                 values_char_2d = output_values_2d_char_pointer,                   &
+                                 bounds_start = (/ 1,                   vmea(l)%start_coord_s /),  &
                                  bounds_end   = (/ maximum_name_length, vmea(l)%end_coord_s   /) )
 
@@ -2496,5 +2354,5 @@
 !--    Write output variables. Distinguish between atmosphere and soil variables.
        DO  n = 1, vmea(l)%nmeas
-          IF ( vmea(l)%soil_sampling  .AND.                                    &
+          IF ( vmea(l)%soil_sampling  .AND.                                                        &
             ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
 !
@@ -2505,8 +2363,8 @@
              output_values_2d_pointer => output_values_2d_target
 
-             return_value = dom_write_var( vmea(l)%nc_filename,                &
-                                           variable_name,                      &
+             return_value = dom_write_var( vmea(l)%nc_filename,                                    &
+                                           variable_name,                                          &
                                            values_realwp_2d = output_values_2d_pointer, &
-                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
+                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/),        &
                                            bounds_end   = (/vmea(l)%end_coord_s, t_ind /) )
 
@@ -2514,9 +2372,9 @@
              output_values_2d_target(t_ind,:) = vmea(l)%measured_vars_soil(:,n)
              output_values_2d_pointer => output_values_2d_target
-             return_value =                                                    &
-                      dom_write_var( vmea(l)%nc_filename,                      &
-                                     variable_name,                            &
-                                     values_realwp_2d = output_values_2d_pointer, &
-                                     bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
+             return_value =                                                                        &
+                      dom_write_var( vmea(l)%nc_filename,                                          &
+                                     variable_name,                                                &
+                                     values_realwp_2d = output_values_2d_pointer,                  &
+                                     bounds_start = (/vmea(l)%start_coord_s, t_ind/),              &
                                      bounds_end   = (/vmea(l)%end_coord_s, t_ind  /) )
              DEALLOCATE( output_values_2d_target )
@@ -2529,8 +2387,8 @@
              output_values_2d_pointer => output_values_2d_target
 
-             return_value = dom_write_var( vmea(l)%nc_filename,                &
-                                           variable_name,                      &
-                                           values_realwp_2d = output_values_2d_pointer, &
-                                           bounds_start = (/vmea(l)%start_coord_a, t_ind/), &
+             return_value = dom_write_var( vmea(l)%nc_filename,                                    &
+                                           variable_name,                                          &
+                                           values_realwp_2d = output_values_2d_pointer,            &
+                                           bounds_start = (/vmea(l)%start_coord_a, t_ind/),        &
                                            bounds_end   = (/vmea(l)%end_coord_a, t_ind/) )
 
@@ -2539,9 +2397,9 @@
              output_values_2d_target(t_ind,:) = vmea(l)%measured_vars(:,n)
              output_values_2d_pointer => output_values_2d_target
-             return_value =                                                    &
-                      dom_write_var( vmea(l)%nc_filename,                      &
-                                     variable_name,                            &
-                                     values_realwp_2d = output_values_2d_pointer, &
-                                     bounds_start = (/ vmea(l)%start_coord_a, t_ind /), &
+             return_value =                                                                        &
+                      dom_write_var( vmea(l)%nc_filename,                                          &
+                                     variable_name,                                                &
+                                     values_realwp_2d = output_values_2d_pointer,                  &
+                                     bounds_start = (/ vmea(l)%start_coord_a, t_ind /),            &
                                      bounds_end   = (/ vmea(l)%end_coord_a, t_ind /) )
 
@@ -2558,19 +2416,19 @@
 
 
-  END SUBROUTINE vm_data_output
-
-!------------------------------------------------------------------------------!
+ END SUBROUTINE vm_data_output
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Sampling of the actual quantities along the observation coordinates
-!------------------------------------------------------------------------------!
-  SUBROUTINE vm_sampling
-
-    USE radiation_model_mod,                                                   &
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE vm_sampling
+
+    USE radiation_model_mod,                                                                       &
         ONLY:  radiation
 
-    USE surface_mod,                                                           &
-        ONLY:  surf_def_h,                                                     &
-               surf_lsm_h,                                                     &
+    USE surface_mod,                                                                               &
+        ONLY:  surf_def_h,                                                                         &
+               surf_lsm_h,                                                                         &
                surf_usm_h
 
@@ -2585,10 +2443,10 @@
      INTEGER(iwp) ::  nn        !< running index over the number of chemcal species
 
-     LOGICAL ::  match_lsm !< flag indicating natural-type surface
-     LOGICAL ::  match_usm !< flag indicating urban-type surface
-
-     REAL(wp) ::  e_s      !< saturation water vapor pressure
-     REAL(wp) ::  q_s      !< saturation mixing ratio
-     REAL(wp) ::  q_wv     !< mixing ratio
+     LOGICAL ::  match_lsm  !< flag indicating natural-type surface
+     LOGICAL ::  match_usm  !< flag indicating urban-type surface
+
+     REAL(wp) ::  e_s   !< saturation water vapor pressure
+     REAL(wp) ::  q_s   !< saturation mixing ratio
+     REAL(wp) ::  q_wv  !< mixing ratio
 
      CALL cpu_log( log_point_s(27), 'VM sampling', 'start' )
@@ -2598,8 +2456,6 @@
 !
 !--     At the beginning, set _FillValues
-        IF ( ALLOCATED( vmea(l)%measured_vars      ) )                         &
-           vmea(l)%measured_vars      = vmea(l)%fillout
-        IF ( ALLOCATED( vmea(l)%measured_vars_soil ) )                         &
-           vmea(l)%measured_vars_soil = vmea(l)%fillout
+        IF ( ALLOCATED( vmea(l)%measured_vars ) ) vmea(l)%measured_vars = vmea(l)%fillout
+        IF ( ALLOCATED( vmea(l)%measured_vars_soil ) ) vmea(l)%measured_vars_soil = vmea(l)%fillout
 !
 !--     Loop over all variables measured at this site.
@@ -2624,6 +2480,5 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )     &
-                                                  - degc_to_k
+                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k ) - degc_to_k
                     ENDDO
                  ENDIF
@@ -2647,7 +2502,5 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(m,n) = ( q(k,j,i)                 &
-                                                    / ( 1.0_wp - q(k,j,i) ) )  &
-                                                  * rho_air(k)
+                       vmea(l)%measured_vars(m,n) = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) * rho_air(k)
                     ENDDO
                  ENDIF
@@ -2659,7 +2512,6 @@
 !                        j = vmea(l)%j(m)
 !                        i = vmea(l)%i(m)
-!                        vmea(l)%measured_vars(m,n) = ( q(k,j,i)                 &
-!                                                     / ( 1.0_wp - q(k,j,i) ) )  &
-!                                                   * rho_air(k)
+!                        vmea(l)%measured_vars(m,n) = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) )          &
+!                                                     * rho_air(k)
 !                     ENDDO
                  ENDIF
@@ -2672,6 +2524,6 @@
                        i = vmea(l)%i(m)
 !
-!--                    Calculate actual temperature, water vapor saturation
-!--                    pressure, and based on this  the saturation mixing ratio.
+!--                    Calculate actual temperature, water vapor saturation pressure and, based on
+!--                    this, the saturation mixing ratio.
                        e_s  = magnus( exner(k) * pt(k,j,i) )
                        q_s  = rd_d_rv * e_s / ( hyp(k) - e_s )
@@ -2711,7 +2563,7 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(m,n) = SQRT(                         &
-                                   ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
-                                   ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2   &
+                    vmea(l)%measured_vars(m,n) = SQRT(                                             &
+                                                 ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 +       &
+                                                 ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2         &
                                                      )
                  ENDDO
@@ -2723,9 +2575,8 @@
                     i = vmea(l)%i(m)
 
-                    vmea(l)%measured_vars(m,n) = 180.0_wp + 180.0_wp / pi      &
-                                               * ATAN2(                        &
-                                          0.5_wp * ( v(k,j,i) + v(k,j+1,i) ),  &
-                                          0.5_wp * ( u(k,j,i) + u(k,j,i+1) )   &
-                                                      )
+                    vmea(l)%measured_vars(m,n) = 180.0_wp + 180.0_wp / pi * ATAN2(                 &
+                                                 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ),               &
+                                                 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )                &
+                                                                                  )
                  ENDDO
 
@@ -2735,7 +2586,5 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
-                                                 ( u(k,j,i) + u(k,j,i+1) ) *   &
-                                                   pt(k,j,i)
+                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) * pt(k,j,i)
                  ENDDO
 
@@ -2745,7 +2594,5 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
-                                                 ( v(k,j,i) + v(k,j+1,i) ) *   &
-                                                   pt(k,j,i)
+                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) * pt(k,j,i)
                  ENDDO
 
@@ -2755,7 +2602,5 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
-                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
-                                                   pt(k,j,i)
+                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k-1,j,i) + w(k,j,i) ) * pt(k,j,i)
                  ENDDO
 
@@ -2766,7 +2611,5 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
-                                                    ( u(k,j,i) + u(k,j,i+1) ) *&
-                                                      q(k,j,i)
+                       vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) * q(k,j,i)
                     ENDDO
                  ENDIF
@@ -2778,7 +2621,5 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
-                                                    ( v(k,j,i) + v(k,j+1,i) ) *&
-                                                      q(k,j,i)
+                       vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) * q(k,j,i)
                     ENDDO
                  ENDIF
@@ -2790,7 +2631,5 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
-                                                    ( w(k-1,j,i) + w(k,j,i) ) *&
-                                                      q(k,j,i)
+                       vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k-1,j,i) + w(k,j,i) ) * q(k,j,i)
                     ENDDO
                  ENDIF
@@ -2801,7 +2640,6 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
-                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
-                                                 ( u(k,j,i)   + u(k,j,i+1) )
+                    vmea(l)%measured_vars(m,n) = 0.25_wp * ( w(k-1,j,i) + w(k,j,i) ) *             &
+                                                           ( u(k,j,i)   + u(k,j,i+1) )
                  ENDDO
 
@@ -2811,7 +2649,6 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
-                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
-                                                 ( v(k,j,i)   + v(k,j+1,i) )
+                    vmea(l)%measured_vars(m,n) = 0.25_wp * ( w(k-1,j,i) + w(k,j,i) ) *             &
+                                                           ( v(k,j,i)   + v(k,j+1,i) )
                  ENDDO
 
@@ -2821,15 +2658,12 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
-                                                 ( u(k,j,i)   + u(k,j,i+1) ) * &
-                                                 ( v(k,j,i)   + v(k,j+1,i) )
+                    vmea(l)%measured_vars(m,n) = 0.25_wp * ( u(k,j,i)   + u(k,j,i+1) ) *           &
+                                                           ( v(k,j,i)   + v(k,j+1,i) )
                  ENDDO
 !
-!--           Chemistry variables. List of variables may need extension.
-!--           Note, gas species in PALM are in ppm and no distinction is made
-!--           between mole-fraction and concentration quantities (all are
-!--           output in ppm so far).
-              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfno', 'mfno2',            &
-                     'mcno', 'mcno2', 'tro3' )
+!--           Chemistry variables. List of variables that may need extension. Note, gas species in
+!--           PALM are in ppm and no distinction is made between mole-fraction and concentration
+!--           quantities (all are output in ppm so far).
+              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfno', 'mfno2', 'mcno', 'mcno2', 'tro3' )
                  IF ( air_chemistry )  THEN
 !
@@ -2837,20 +2671,17 @@
 !--                 list, in order to get the internal name of the variable.
                     DO  nn = 1, UBOUND( chem_vars, 2 )
-                       IF ( TRIM( vmea(l)%var_atts(n)%name ) ==                &
+                       IF ( TRIM( vmea(l)%var_atts(n)%name ) ==                                    &
                             TRIM( chem_vars(0,nn) ) )  ind_chem = nn
                     ENDDO
 !
-!--                 Run loop over all chemical species, if the measured
-!--                 variable matches the interal name, sample the variable.
-!--                 Note, nvar as a chemistry-module variable.
+!--                 Run loop over all chemical species, if the measured variable matches the interal
+!--                 name, sample the variable. Note, nvar as a chemistry-module variable.
                     DO  nn = 1, nvar
-                       IF ( TRIM( chem_vars(1,ind_chem) ) ==                   &
-                            TRIM( chem_species(nn)%name ) )  THEN
+                       IF ( TRIM( chem_vars(1,ind_chem) ) == TRIM( chem_species(nn)%name ) )  THEN
                           DO  m = 1, vmea(l)%ns
                              k = vmea(l)%k(m)
                              j = vmea(l)%j(m)
                              i = vmea(l)%i(m)
-                             vmea(l)%measured_vars(m,n) =                      &
-                                                   chem_species(nn)%conc(k,j,i)
+                             vmea(l)%measured_vars(m,n) = chem_species(nn)%conc(k,j,i)
                           ENDDO
                        ENDIF
@@ -2861,6 +2692,6 @@
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not
-!--                 on ghost points. Hence, limit the indices.
+!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
+!--                 limit the indices.
                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                     j = MERGE( j           , nyn, j            < nyn )
@@ -2868,14 +2699,11 @@
                     i = MERGE( i           , nxr, i            < nxr )
 
-                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
-                             surf_def_h(0)%end_index(j,i)
+                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
                     ENDDO
-                    DO  mm = surf_lsm_h%start_index(j,i),                      &
-                             surf_lsm_h%end_index(j,i)
+                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm)
                     ENDDO
-                    DO  mm = surf_usm_h%start_index(j,i),                      &
-                             surf_usm_h%end_index(j,i)
+                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm)
                     ENDDO
@@ -2885,6 +2713,6 @@
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not
-!--                 on ghost points. Hence, limit the indices.
+!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
+!-                  limit the indices.
                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                     j = MERGE( j           , nyn, j            < nyn )
@@ -2892,14 +2720,11 @@
                     i = MERGE( i           , nxr, i            < nxr )
 
-                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
-                             surf_def_h(0)%end_index(j,i)
+                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
                     ENDDO
-                    DO  mm = surf_lsm_h%start_index(j,i),                      &
-                             surf_lsm_h%end_index(j,i)
+                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm)
                     ENDDO
-                    DO  mm = surf_usm_h%start_index(j,i),                      &
-                             surf_usm_h%end_index(j,i)
+                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm)
                     ENDDO
@@ -2909,6 +2734,6 @@
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not
-!--                 on ghost points. Hence, limit the indices.
+!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
+!--                 limit the indices.
                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                     j = MERGE( j           , nyn, j            < nyn )
@@ -2916,14 +2741,11 @@
                     i = MERGE( i           , nxr, i            < nxr )
 
-                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
-                             surf_def_h(0)%end_index(j,i)
+                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
                     ENDDO
-                    DO  mm = surf_lsm_h%start_index(j,i),                      &
-                             surf_lsm_h%end_index(j,i)
+                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm)
                     ENDDO
-                    DO  mm = surf_usm_h%start_index(j,i),                      &
-                             surf_usm_h%end_index(j,i)
+                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm)
                     ENDDO
@@ -2933,6 +2755,6 @@
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not
-!--                 on ghost points. Hence, limit the indices.
+!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
+!--                 limit the indices.
                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                     j = MERGE( j           , nyn, j            < nyn )
@@ -2940,14 +2762,11 @@
                     i = MERGE( i           , nxr, i            < nxr )
 
-                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
-                             surf_def_h(0)%end_index(j,i)
+                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
                     ENDDO
-                    DO  mm = surf_lsm_h%start_index(j,i),                      &
-                             surf_lsm_h%end_index(j,i)
+                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm)
                     ENDDO
-                    DO  mm = surf_usm_h%start_index(j,i),                      &
-                             surf_usm_h%end_index(j,i)
+                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm)
                     ENDDO
@@ -2957,6 +2776,6 @@
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not
-!--                 on ghost points. Hence, limit the indices.
+!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
+!--                 limit the indices.
                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                     j = MERGE( j           , nyn, j            < nyn )
@@ -2964,6 +2783,5 @@
                     i = MERGE( i           , nxr, i            < nxr )
 
-                    DO  mm = surf_lsm_h%start_index(j,i),                      &
-                             surf_lsm_h%end_index(j,i)
+                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                        vmea(l)%measured_vars(m,n) = surf_lsm_h%ghf(mm)
                     ENDDO
@@ -2973,6 +2791,6 @@
 !                  DO  m = 1, vmea(l)%ns
 ! !
-! !--                 Surface data is only available on inner subdomains, not
-! !--                 on ghost points. Hence, limit the indices.
+! !--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
+! !--                 limit the indices.
 !                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
 !                     j = MERGE( j           , nyn, j            < nyn )
@@ -2980,6 +2798,5 @@
 !                     i = MERGE( i           , nxr, i            < nxr )
 !
-!                     DO  mm = surf_lsm_h%start_index(j,i),                      &
-!                              surf_lsm_h%end_index(j,i)
+!                     DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
 !                        vmea(l)%measured_vars(m,n) = ?
 !                     ENDDO
@@ -2990,6 +2807,6 @@
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not
-!--                    on ghost points. Hence, limit the indices.
+!--                    Surface data is only available on inner subdomains, not on ghost points.
+!--                    Hence, limit the indices.
                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                        j = MERGE( j           , nyn, j            < nyn )
@@ -2997,10 +2814,8 @@
                        i = MERGE( i           , nxr, i            < nxr )
 
-                       DO  mm = surf_lsm_h%start_index(j,i),                   &
-                                surf_lsm_h%end_index(j,i)
+                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm)
                        ENDDO
-                       DO  mm = surf_usm_h%start_index(j,i),                   &
-                                surf_usm_h%end_index(j,i)
+                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm)
                        ENDDO
@@ -3012,6 +2827,6 @@
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not
-!--                    on ghost points. Hence, limit the indices.
+!--                    Surface data is only available on inner subdomains, not on ghost points.
+!--                    Hence, limit the indices.
                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                        j = MERGE( j           , nyn, j            < nyn )
@@ -3019,10 +2834,8 @@
                        i = MERGE( i           , nxr, i            < nxr )
 
-                       DO  mm = surf_lsm_h%start_index(j,i),                   &
-                                surf_lsm_h%end_index(j,i)
+                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm)
                        ENDDO
-                       DO  mm = surf_usm_h%start_index(j,i),                   &
-                                surf_usm_h%end_index(j,i)
+                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm)
                        ENDDO
@@ -3034,6 +2847,6 @@
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not
-!--                    on ghost points. Hence, limit the indices.
+!--                    Surface data is only available on inner subdomains, not on ghost points.
+!--                    Hence, limit the indices.
                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                        j = MERGE( j           , nyn, j            < nyn )
@@ -3041,10 +2854,8 @@
                        i = MERGE( i           , nxr, i            < nxr )
 
-                       DO  mm = surf_lsm_h%start_index(j,i),                   &
-                                surf_lsm_h%end_index(j,i)
+                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm)
                        ENDDO
-                       DO  mm = surf_usm_h%start_index(j,i),                   &
-                                surf_usm_h%end_index(j,i)
+                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm)
                        ENDDO
@@ -3056,6 +2867,6 @@
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not
-!--                    on ghost points. Hence, limit the indices.
+!--                    Surface data is only available on inner subdomains, not on ghost points.
+!--                    Hence, limit the indices.
                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                        j = MERGE( j           , nyn, j            < nyn )
@@ -3063,10 +2874,8 @@
                        i = MERGE( i           , nxr, i            < nxr )
 
-                       DO  mm = surf_lsm_h%start_index(j,i),                   &
-                                surf_lsm_h%end_index(j,i)
+                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm)
                        ENDDO
-                       DO  mm = surf_usm_h%start_index(j,i),                   &
-                                surf_usm_h%end_index(j,i)
+                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm)
                        ENDDO
@@ -3078,6 +2887,6 @@
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not
-!--                    on ghost points. Hence, limit the indices.
+!--                    Surface data is only available on inner subdomains, not on ghost points.
+!--                    Hence, limit the indices.
                        j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                        j = MERGE( j           , nyn, j            < nyn )
@@ -3085,10 +2894,8 @@
                        i = MERGE( i           , nxr, i            < nxr )
 
-                       DO  mm = surf_lsm_h%start_index(j,i),                   &
-                                surf_lsm_h%end_index(j,i)
+                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm)
                        ENDDO
-                       DO  mm = surf_usm_h%start_index(j,i),                   &
-                                surf_usm_h%end_index(j,i)
+                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                           vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm)
                        ENDDO
@@ -3159,6 +2966,6 @@
                           k = 0
 !
-!--                       Surface data is only available on inner subdomains,
-!--                       not on ghost points. Hence, limit the indices.
+!--                       Surface data is only available on inner subdomains, not on ghost points.
+!--                       Hence, limit the indices.
                           j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
                           j = MERGE( j           , nyn, j            < nyn )
@@ -3196,8 +3003,6 @@
                     k = vmea(l)%k_soil(m)
 
-                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
-                                surf_lsm_h%end_index(j,i)
-                    match_usm = surf_usm_h%start_index(j,i) <=                 &
-                                surf_usm_h%end_index(j,i)
+                    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
+                    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
 
                     IF ( match_lsm )  THEN
@@ -3220,6 +3025,5 @@
                     k = vmea(l)%k_soil(m)
 
-                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
-                                surf_lsm_h%end_index(j,i)
+                    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
 
                     IF ( match_lsm )  THEN
@@ -3231,4 +3035,45 @@
 
               CASE ( 'ts' ) ! surface temperature
+                 DO  m = 1, vmea(l)%ns
+!
+!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
+!--                 limit the indices.
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                    j = MERGE( j           , nyn, j            < nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                    i = MERGE( i           , nxr, i            < nxr )
+
+                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
+                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%pt_surface(mm)
+                    ENDDO
+                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
+                       vmea(l)%measured_vars(m,n) = surf_lsm_h%pt_surface(mm)
+                    ENDDO
+                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
+                       vmea(l)%measured_vars(m,n) = surf_usm_h%pt_surface(mm)
+                    ENDDO
+                 ENDDO
+
+              CASE ( 'lwp' ) ! liquid water path
+                 IF ( ASSOCIATED( ql ) )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+
+                       vmea(l)%measured_vars(m,n) = SUM( ql(nzb:nzt,j,i) * dzw(1:nzt+1) )          &
+                                                    * rho_surface
+                    ENDDO
+                 ENDIF
+
+              CASE ( 'ps' ) ! surface pressure
+                 vmea(l)%measured_vars(:,n) = surface_pressure
+
+              CASE ( 'pswrtg' ) ! platform speed above ground
+                 vmea(l)%measured_vars(:,n) = 0.0_wp
+
+              CASE ( 'pswrta' ) ! platform speed in air
+                 vmea(l)%measured_vars(:,n) = 0.0_wp
+
+              CASE ( 't_lw' ) ! water temperature
                  DO  m = 1, vmea(l)%ns
 !
@@ -3240,53 +3085,7 @@
                     i = MERGE( i           , nxr, i            < nxr )
 
-                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
-                             surf_def_h(0)%end_index(j,i)
-                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%pt_surface(mm)
-                    ENDDO
-                    DO  mm = surf_lsm_h%start_index(j,i),                      &
-                             surf_lsm_h%end_index(j,i)
-                       vmea(l)%measured_vars(m,n) = surf_lsm_h%pt_surface(mm)
-                    ENDDO
-                    DO  mm = surf_usm_h%start_index(j,i),                      &
-                             surf_usm_h%end_index(j,i)
-                       vmea(l)%measured_vars(m,n) = surf_usm_h%pt_surface(mm)
-                    ENDDO
-                 ENDDO
-
-              CASE ( 'lwp' ) ! liquid water path
-                 IF ( ASSOCIATED( ql ) )  THEN
-                    DO  m = 1, vmea(l)%ns
-                       j = vmea(l)%j(m)
-                       i = vmea(l)%i(m)
-
-                       vmea(l)%measured_vars(m,n) = SUM( ql(nzb:nzt,j,i)       &
-                                                       * dzw(1:nzt+1) )        &
-                                                  * rho_surface
-                    ENDDO
-                 ENDIF
-
-              CASE ( 'ps' ) ! surface pressure
-                 vmea(l)%measured_vars(:,n) = surface_pressure
-
-              CASE ( 'pswrtg' ) ! platform speed above ground
-                 vmea(l)%measured_vars(:,n) = 0.0_wp
-
-              CASE ( 'pswrta' ) ! platform speed in air
-                 vmea(l)%measured_vars(:,n) = 0.0_wp
-
-              CASE ( 't_lw' ) ! water temperature
-                 DO  m = 1, vmea(l)%ns
-!
-!--                 Surface data is only available on inner subdomains, not
-!--                 on ghost points. Hence, limit the indices.
-                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
-                    j = MERGE( j           , nyn, j            < nyn )
-                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
-                    i = MERGE( i           , nxr, i            < nxr )
-
-                    DO  mm = surf_lsm_h%start_index(j,i),                      &
-                             surf_lsm_h%end_index(j,i)
-                       IF ( surf_lsm_h%water_surface(m) )                      &
-                          vmea(l)%measured_vars(m,n) = t_soil_h%var_2d(nzt,m)
+                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
+                       IF ( surf_lsm_h%water_surface(m) )                                          &
+                            vmea(l)%measured_vars(m,n) = t_soil_h%var_2d(nzt,m)
                     ENDDO
 
@@ -3307,5 +3106,5 @@
      CALL cpu_log( log_point_s(27), 'VM sampling', 'stop' )
 
-  END SUBROUTINE vm_sampling
+ END SUBROUTINE vm_sampling
 
 
