source: palm/trunk/SOURCE/virtual_measurement_mod.f90 @ 4643

Last change on this file since 4643 was 4642, checked in by suehring, 4 years ago

Do not set attribute bounds for time variable, as it refers to time_bounds which is not defined for non-aggregated quantities (according to data standard)

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 154.4 KB
Line 
1!> @virtual_measurement_mod.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: virtual_measurement_mod.f90 4642 2020-08-13 15:47:33Z pavelkrc $
27! Do not set attribute bounds for time variable, as it refers to time_bounds which is not defined
28! for non-aggregated quantities (according to data standard)
29!
30! 4641 2020-08-13 09:57:07Z suehring
31! - To be in agreement with (UC)2 data standard do not list the measured variables in attribute
32!   data_content but simply set 'airmeteo'
33! - Bugfix in setting long_name attribute for variable t_va and for global attribute creation_time
34!
35! 4536 2020-05-17 17:24:13Z raasch
36! bugfix: preprocessor directive adjusted
37!
38! 4504 2020-04-20 12:11:24Z raasch
39! file re-formatted to follow the PALM coding standard
40!
41! 4481 2020-03-31 18:55:54Z maronga
42! bugfix: cpp-directives for serial mode added
43!
44! 4438 2020-03-03 20:49:28Z suehring
45! Add cpu-log points
46!
47! 4422 2020-02-24 22:45:13Z suehring
48! Missing trim()
49!
50! 4408 2020-02-14 10:04:39Z gronemeier
51! - Output of character string station_name after DOM has been enabled to
52!   output character variables
53! - Bugfix, missing coupling_char statement when opening the input file
54!
55! 4408 2020-02-14 10:04:39Z gronemeier
56! write fill_value attribute
57!
58! 4406 2020-02-13 20:06:29Z knoop
59! Bugix: removed oro_rel wrong loop bounds and removed unnecessary restart method
60!
61! 4400 2020-02-10 20:32:41Z suehring
62! Revision of the module:
63! - revised input from NetCDF setup file
64! - parallel NetCDF output via data-output module ( Tobias Gronemeier )
65! - variable attributes added
66! - further variables defined
67!
68! 4346 2019-12-18 11:55:56Z motisi
69! Introduction of wall_flags_total_0, which currently sets bits based on static
70! topography information used in wall_flags_static_0
71!
72! 4329 2019-12-10 15:46:36Z motisi
73! Renamed wall_flags_0 to wall_flags_static_0
74!
75! 4226 2019-09-10 17:03:24Z suehring
76! Netcdf input routine for dimension length renamed
77!
78! 4182 2019-08-22 15:20:23Z scharf
79! Corrected "Former revisions" section
80!
81! 4168 2019-08-16 13:50:17Z suehring
82! Replace function get_topography_top_index by topo_top_ind
83!
84! 3988 2019-05-22 11:32:37Z kanani
85! Add variables to enable steering of output interval for virtual measurements
86!
87! 3913 2019-04-17 15:12:28Z gronemeier
88! Bugfix: rotate positions of measurements before writing them into file
89!
90! 3910 2019-04-17 11:46:56Z suehring
91! Bugfix in rotation of UTM coordinates
92!
93! 3904 2019-04-16 18:22:51Z gronemeier
94! Rotate coordinates of stations by given rotation_angle
95!
96! 3876 2019-04-08 18:41:49Z knoop
97! Remove print statement
98!
99! 3854 2019-04-02 16:59:33Z suehring
100! renamed nvar to nmeas, replaced USE chem_modules by USE chem_gasphase_mod and
101! nspec by nvar
102!
103! 3766 2019-02-26 16:23:41Z raasch
104! unused variables removed
105!
106! 3718 2019-02-06 11:08:28Z suehring
107! Adjust variable name connections between UC2 and chemistry variables
108!
109! 3717 2019-02-05 17:21:16Z suehring
110! Additional check + error numbers adjusted
111!
112! 3706 2019-01-29 20:02:26Z suehring
113! unused variables removed
114!
115! 3705 2019-01-29 19:56:39Z suehring
116! - initialization revised
117! - binary data output
118! - list of allowed variables extended
119!
120! 3704 2019-01-29 19:51:41Z suehring
121! Sampling of variables
122!
123! 3473 2018-10-30 20:50:15Z suehring
124! Initial revision
125!
126! Authors:
127! --------
128! @author Matthias Suehring
129! @author Tobias Gronemeier
130!
131! Description:
132! ------------
133!> The module acts as an interface between 'real-world' observations and model simulations.
134!> Virtual measurements will be taken in the model at the coordinates representative for the
135!> 'real-world' observation coordinates. More precisely, coordinates and measured quanties will be
136!> read from a NetCDF file which contains all required information. In the model, the same
137!> quantities (as long as all the required components are switched-on) will be sampled at the
138!> respective positions and output into an extra file, which allows for straight-forward comparison
139!> of model results with observations.
140!--------------------------------------------------------------------------------------------------!
141 MODULE virtual_measurement_mod
142
143    USE arrays_3d,                                                                                 &
144        ONLY:  dzw,                                                                                &
145               exner,                                                                              &
146               hyp,                                                                                &
147               q,                                                                                  &
148               ql,                                                                                 &
149               pt,                                                                                 &
150               rho_air,                                                                            &
151               u,                                                                                  &
152               v,                                                                                  &
153               w,                                                                                  &
154               zu,                                                                                 &
155               zw
156
157    USE basic_constants_and_equations_mod,                                                         &
158        ONLY:  convert_utm_to_geographic,                                                          &
159               degc_to_k,                                                                          &
160               magnus,                                                                             &
161               pi,                                                                                 &
162               rd_d_rv
163
164    USE chem_gasphase_mod,                                                                         &
165        ONLY:  nvar
166
167    USE chem_modules,                                                                              &
168        ONLY:  chem_species
169
170    USE control_parameters,                                                                        &
171        ONLY:  air_chemistry,                                                                      &
172               coupling_char,                                                                      &
173               dz,                                                                                 &
174               end_time,                                                                           &
175               humidity,                                                                           &
176               message_string,                                                                     &
177               neutral,                                                                            &
178               origin_date_time,                                                                   &
179               rho_surface,                                                                        &
180               surface_pressure,                                                                   &
181               time_since_reference_point,                                                         &
182               virtual_measurement
183
184    USE cpulog,                                                                                    &
185        ONLY:  cpu_log,                                                                            &
186               log_point_s
187
188    USE data_output_module
189
190    USE grid_variables,                                                                            &
191        ONLY:  ddx,                                                                                &
192               ddy,                                                                                &
193               dx,                                                                                 &
194               dy
195
196    USE indices,                                                                                   &
197        ONLY:  nbgp,                                                                               &
198               nzb,                                                                                &
199               nzt,                                                                                &
200               nxl,                                                                                &
201               nxlg,                                                                               &
202               nxr,                                                                                &
203               nxrg,                                                                               &
204               nys,                                                                                &
205               nysg,                                                                               &
206               nyn,                                                                                &
207               nyng,                                                                               &
208               topo_top_ind,                                                                       &
209               wall_flags_total_0
210
211    USE kinds
212
213    USE netcdf_data_input_mod,                                                                     &
214        ONLY:  close_input_file,                                                                   &
215               coord_ref_sys,                                                                      &
216               crs_list,                                                                           &
217               get_attribute,                                                                      &
218               get_dimension_length,                                                               &
219               get_variable,                                                                       &
220               init_model,                                                                         &
221               input_file_atts,                                                                    &
222               input_file_vm,                                                                      &
223               input_pids_static,                                                                  &
224               input_pids_vm,                                                                      &
225               inquire_fill_value,                                                                 &
226               open_read_file,                                                                     &
227               pids_id
228
229    USE pegrid
230
231    USE surface_mod,                                                                               &
232        ONLY:  surf_lsm_h,                                                                         &
233               surf_usm_h
234
235    USE land_surface_model_mod,                                                                    &
236        ONLY:  m_soil_h,                                                                           &
237               nzb_soil,                                                                           &
238               nzt_soil,                                                                           &
239               t_soil_h,                                                                           &
240               zs
241
242    USE radiation_model_mod,                                                                       &
243        ONLY:  rad_lw_in,                                                                          &
244               rad_lw_out,                                                                         &
245               rad_sw_in,                                                                          &
246               rad_sw_in_diff,                                                                     &
247               rad_sw_out,                                                                         &
248               radiation_scheme
249
250    USE urban_surface_mod,                                                                         &
251        ONLY:  nzb_wall,                                                                           &
252               nzt_wall,                                                                           &
253               t_wall_h
254
255
256    IMPLICIT NONE
257
258    TYPE virt_general
259       INTEGER(iwp) ::  nvm = 0  !< number of virtual measurements
260    END TYPE virt_general
261
262    TYPE virt_var_atts
263       CHARACTER(LEN=100) ::  coordinates           !< defined longname of the variable
264       CHARACTER(LEN=100) ::  grid_mapping          !< defined longname of the variable
265       CHARACTER(LEN=100) ::  long_name             !< defined longname of the variable
266       CHARACTER(LEN=100) ::  name                  !< variable name
267       CHARACTER(LEN=100) ::  standard_name         !< defined standard name of the variable
268       CHARACTER(LEN=100) ::  units                 !< unit of the output variable
269
270       REAL(wp)           ::  fill_value = -9999.0  !< _FillValue attribute
271    END TYPE virt_var_atts
272
273    TYPE virt_mea
274       CHARACTER(LEN=100)  ::  feature_type                      !< type of the real-world measurement
275       CHARACTER(LEN=100)  ::  feature_type_out = 'timeSeries'   !< type of the virtual measurement
276                                                                 !< (all will be timeSeries, even trajectories)
277       CHARACTER(LEN=100)  ::  nc_filename                       !< name of the NetCDF output file for the station
278       CHARACTER(LEN=100)  ::  site                              !< name of the measurement site
279
280       CHARACTER(LEN=1000) ::  data_content = REPEAT(' ', 1000)  !< string of measured variables (data output only)
281
282       INTEGER(iwp) ::  end_coord_a     = 0  !< end coordinate in NetCDF file for local atmosphere observations
283       INTEGER(iwp) ::  end_coord_s     = 0  !< end coordinate in NetCDF file for local soil observations
284       INTEGER(iwp) ::  file_time_index = 0  !< time index in NetCDF output file
285       INTEGER(iwp) ::  ns              = 0  !< number of observation coordinates on subdomain, for atmospheric measurements
286       INTEGER(iwp) ::  ns_tot          = 0  !< total number of observation coordinates, for atmospheric measurements
287       INTEGER(iwp) ::  n_tr_st              !< number of trajectories / station of a measurement
288       INTEGER(iwp) ::  nmeas                !< number of measured variables (atmosphere + soil)
289       INTEGER(iwp) ::  ns_soil         = 0  !< number of observation coordinates on subdomain, for soil measurements
290       INTEGER(iwp) ::  ns_soil_tot     = 0  !< total number of observation coordinates, for soil measurements
291       INTEGER(iwp) ::  start_coord_a   = 0  !< start coordinate in NetCDF file for local atmosphere observations
292       INTEGER(iwp) ::  start_coord_s   = 0  !< start coordinate in NetCDF file for local soil observations
293
294       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t  !< number observations individual for each trajectory
295                                                          !< or station that are no _FillValues
296
297       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< grid index for measurement position in x-direction
298       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< grid index for measurement position in y-direction
299       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< grid index for measurement position in k-direction
300
301       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i_soil  !< grid index for measurement position in x-direction
302       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j_soil  !< grid index for measurement position in y-direction
303       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k_soil  !< grid index for measurement position in k-direction
304
305       LOGICAL ::  soil_sampling      = .FALSE.  !< flag indicating that soil state variables were sampled
306       LOGICAL ::  trajectory         = .FALSE.  !< flag indicating that the observation is a mobile observation
307       LOGICAL ::  timseries          = .FALSE.  !< flag indicating that the observation is a stationary point measurement
308       LOGICAL ::  timseries_profile  = .FALSE.  !< flag indicating that the observation is a stationary profile measurement
309
310       REAL(wp) ::  fill_eutm          !< fill value for UTM coordinates in case of missing values
311       REAL(wp) ::  fill_nutm          !< fill value for UTM coordinates in case of missing values
312       REAL(wp) ::  fill_zar           !< fill value for heigth coordinates in case of missing values
313       REAL(wp) ::  fillout = -9999.0  !< fill value for output in case an observation is taken e.g. from inside a building
314       REAL(wp) ::  origin_x_obs       !< origin of the observation in UTM coordiates in x-direction
315       REAL(wp) ::  origin_y_obs       !< origin of the observation in UTM coordiates in y-direction
316
317       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth         !< measurement depth in soil
318       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zar           !< measurement height above ground level
319
320       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars       !< measured variables
321       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars_soil  !< measured variables
322
323       TYPE( virt_var_atts ), DIMENSION(:), ALLOCATABLE ::  var_atts !< variable attributes
324    END TYPE virt_mea
325
326    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"            !< dimension name for UTM coordinate easting
327    CHARACTER(LEN=11) ::  char_feature = "featureType"   !< attribute name for feature type
328
329    ! This need to be generalized
330    CHARACTER(LEN=10) ::  char_fill = '_FillValue'                 !< attribute name for fill value
331    CHARACTER(LEN=9)  ::  char_long = 'long_name'                  !< attribute name for long_name
332    CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
333    CHARACTER(LEN=5)  ::  char_nutm = "N_UTM"                      !< dimension name for UTM coordinate northing
334    CHARACTER(LEN=18) ::  char_numstations = "number_of_stations"  !< attribute name for number of stations
335    CHARACTER(LEN=8)  ::  char_origx = "origin_x"                  !< attribute name for station coordinate in x
336    CHARACTER(LEN=8)  ::  char_origy = "origin_y"                  !< attribute name for station coordinate in y
337    CHARACTER(LEN=4)  ::  char_site = "site"                       !< attribute name for site name
338    CHARACTER(LEN=11) ::  char_soil = "soil_sample"                !< attribute name for soil sampling indication
339    CHARACTER(LEN=13) ::  char_standard = 'standard_name'          !< attribute name for standard_name
340    CHARACTER(LEN=9)  ::  char_station_h = "station_h"             !< variable name indicating height of the site
341    CHARACTER(LEN=5)  ::  char_unit = 'units'                      !< attribute name for standard_name
342    CHARACTER(LEN=1)  ::  char_zar = "z"                           !< attribute name indicating height above reference level
343    CHARACTER(LEN=10) ::  type_ts   = 'timeSeries'                 !< name of stationary point measurements
344    CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
345    CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
346
347    CHARACTER(LEN=6), DIMENSION(1:5) ::  soil_vars = (/ 't_soil', & !< list of soil variables
348                                                        'm_soil',                                  &
349                                                        'lwc   ',                                  &
350                                                        'lwcs  ',                                  &
351                                                        'smp   ' /)
352
353    CHARACTER(LEN=10), DIMENSION(0:1,1:8) ::  chem_vars = RESHAPE( (/ 'mcpm1     ', 'PM1       ',  &
354                                                                      'mcpm2p5   ', 'PM2.5     ',  &
355                                                                      'mcpm10    ', 'PM10      ',  &
356                                                                      'mfno2     ', 'NO2       ',  &
357                                                                      'mfno      ', 'NO        ',  &
358                                                                      'mcno2     ', 'NO2       ',  &
359                                                                      'mcno      ', 'NO        ',  &
360                                                                      'tro3      ', 'O3        '   &
361                                                                    /), (/ 2, 8 /) )
362
363    INTEGER(iwp) ::  maximum_name_length = 32  !< maximum name length of station names
364    INTEGER(iwp) ::  ntimesteps                !< number of timesteps defined in NetCDF output file
365    INTEGER(iwp) ::  off_pr              = 1   !< number of neighboring grid points (in each direction) where virtual profile
366                                               !< measurements shall be taken, in addition to the given coordinates in the driver
367    INTEGER(iwp) ::  off_ts              = 1   !< number of neighboring grid points (in each direction) where virtual timeseries
368                                               !< measurements shall be taken, in addition to the given coordinates in the driver
369    INTEGER(iwp) ::  off_tr              = 1   !< number of neighboring grid points (in each direction) where virtual trajectory
370                                               !< measurements shall be taken, in addition to the given coordinates in the driver
371    LOGICAL ::  global_attribute          = .TRUE.   !< flag indicating a global attribute
372    LOGICAL ::  initial_write_coordinates = .FALSE.  !< flag indicating a global attribute
373    LOGICAL ::  use_virtual_measurement   = .FALSE.  !< Namelist parameter
374
375    REAL(wp) ::  dt_virtual_measurement   = 0.0_wp  !< sampling interval
376    REAL(wp) ::  time_virtual_measurement = 0.0_wp  !< time since last sampling
377    REAL(wp) ::  vm_time_start            = 0.0     !< time after which sampling shall start
378
379    TYPE( virt_general )                        ::  vmea_general  !< data structure which encompasses global variables
380    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea          !< data structure containing station-specific variables
381
382    INTERFACE vm_check_parameters
383       MODULE PROCEDURE vm_check_parameters
384    END INTERFACE vm_check_parameters
385
386    INTERFACE vm_data_output
387       MODULE PROCEDURE vm_data_output
388    END INTERFACE vm_data_output
389
390    INTERFACE vm_init
391       MODULE PROCEDURE vm_init
392    END INTERFACE vm_init
393
394    INTERFACE vm_init_output
395       MODULE PROCEDURE vm_init_output
396    END INTERFACE vm_init_output
397
398    INTERFACE vm_parin
399       MODULE PROCEDURE vm_parin
400    END INTERFACE vm_parin
401
402    INTERFACE vm_sampling
403       MODULE PROCEDURE vm_sampling
404    END INTERFACE vm_sampling
405
406    SAVE
407
408    PRIVATE
409
410!
411!-- Public interfaces
412    PUBLIC  vm_check_parameters,                                                                   &
413            vm_data_output,                                                                        &
414            vm_init,                                                                               &
415            vm_init_output,                                                                        &
416            vm_parin,                                                                              &
417            vm_sampling
418
419!
420!-- Public variables
421    PUBLIC  dt_virtual_measurement,                                                                &
422            time_virtual_measurement,                                                              &
423            vmea,                                                                                  &
424            vmea_general,                                                                          &
425            vm_time_start
426
427 CONTAINS
428
429
430!--------------------------------------------------------------------------------------------------!
431! Description:
432! ------------
433!> Check parameters for virtual measurement module
434!--------------------------------------------------------------------------------------------------!
435 SUBROUTINE vm_check_parameters
436
437    IF ( .NOT. virtual_measurement )  RETURN
438!
439!-- Virtual measurements require a setup file.
440    IF ( .NOT. input_pids_vm )  THEN
441       message_string = 'If virtual measurements are taken, a setup input ' //                     &
442                        'file for the site locations is mandatory.'
443       CALL message( 'vm_check_parameters', 'PA0533', 1, 2, 0, 6, 0 )
444    ENDIF
445!
446!-- In case virtual measurements are taken, a static input file is required.
447!-- This is because UTM coordinates for the PALM domain origin are required for correct mapping of
448!-- the measurements.
449!-- ToDo: Revise this later and remove this requirement.
450    IF ( .NOT. input_pids_static )  THEN
451       message_string = 'If virtual measurements are taken, a static input file is mandatory.'
452       CALL message( 'vm_check_parameters', 'PA0534', 1, 2, 0, 6, 0 )
453    ENDIF
454
455#if !defined( __netcdf4_parallel )
456!
457!-- In case of non-parallel NetCDF the virtual measurement output is not
458!-- working. This is only designed for parallel NetCDF.
459    message_string = 'If virtual measurements are taken, parallel NetCDF is required.'
460    CALL message( 'vm_check_parameters', 'PA0708', 1, 2, 0, 6, 0 )
461#endif
462!
463!-- Check if the given number of neighboring grid points do not exceed the number
464!-- of ghost points.
465    IF ( off_pr > nbgp - 1  .OR.  off_ts > nbgp - 1  .OR.  off_tr > nbgp - 1 )  THEN
466       WRITE(message_string,*)                                                                     &
467                        'If virtual measurements are taken, the number ' //                        &
468                        'of surrounding grid points must not be larger ' //                        &
469                        'than the number of ghost points - 1, which is: ', nbgp - 1
470       CALL message( 'vm_check_parameters', 'PA0705', 1, 2, 0, 6, 0 )
471    ENDIF
472
473    IF ( dt_virtual_measurement <= 0.0 )  THEN
474       message_string = 'dt_virtual_measurement must be > 0.0'
475       CALL message( 'check_parameters', 'PA0706', 1, 2, 0, 6, 0 )
476    ENDIF
477
478 END SUBROUTINE vm_check_parameters
479
480!--------------------------------------------------------------------------------------------------!
481! Description:
482! ------------
483!> Subroutine defines variable attributes according to UC2 standard. Note, later  this list can be
484!> moved to the data-output module where it can be re-used also for other output.
485!--------------------------------------------------------------------------------------------------!
486 SUBROUTINE vm_set_attributes( output_variable )
487
488    TYPE( virt_var_atts ), INTENT(INOUT) ::  output_variable !< data structure with attributes that need to be set
489
490    output_variable%long_name     = 'none'
491    output_variable%standard_name = 'none'
492    output_variable%units         = 'none'
493    output_variable%coordinates   = 'lon lat E_UTM N_UTM x y z time station_name'
494    output_variable%grid_mapping  = 'crs'
495
496    SELECT CASE ( TRIM( output_variable%name ) )
497
498       CASE ( 'u' )
499          output_variable%long_name     = 'u wind component'
500          output_variable%units         = 'm s-1'
501
502       CASE ( 'ua' )
503          output_variable%long_name     = 'eastward wind'
504          output_variable%standard_name = 'eastward_wind'
505          output_variable%units         = 'm s-1'
506
507       CASE ( 'v' )
508          output_variable%long_name     = 'v wind component'
509          output_variable%units         = 'm s-1'
510
511       CASE ( 'va' )
512          output_variable%long_name     = 'northward wind'
513          output_variable%standard_name = 'northward_wind'
514          output_variable%units         = 'm s-1'
515
516       CASE ( 'w' )
517          output_variable%long_name     = 'w wind component'
518          output_variable%standard_name = 'upward_air_velocity'
519          output_variable%units         = 'm s-1'
520
521       CASE ( 'wspeed' )
522          output_variable%long_name     = 'wind speed'
523          output_variable%standard_name = 'wind_speed'
524          output_variable%units         = 'm s-1'
525
526       CASE ( 'wdir' )
527          output_variable%long_name     = 'wind from direction'
528          output_variable%standard_name = 'wind_from_direction'
529          output_variable%units         = 'degrees'
530
531       CASE ( 'theta' )
532          output_variable%long_name     = 'air potential temperature'
533          output_variable%standard_name = 'air_potential_temperature'
534          output_variable%units         = 'K'
535
536       CASE ( 'utheta' )
537          output_variable%long_name     = 'eastward kinematic sensible heat flux in air'
538          output_variable%units         = 'K m s-1'
539
540       CASE ( 'vtheta' )
541          output_variable%long_name     = 'northward kinematic sensible heat flux in air'
542          output_variable%units         = 'K m s-1'
543
544       CASE ( 'wtheta' )
545          output_variable%long_name     = 'upward kinematic sensible heat flux in air'
546          output_variable%units         = 'K m s-1'
547
548       CASE ( 'ta' )
549          output_variable%long_name     = 'air temperature'
550          output_variable%standard_name = 'air_temperature'
551          output_variable%units         = 'degree_C'
552
553       CASE ( 't_va' )
554          output_variable%long_name     = 'virtual acoustic temperature'
555          output_variable%units         = 'K'
556
557       CASE ( 'haa' )
558          output_variable%long_name     = 'absolute atmospheric humidity'
559          output_variable%units         = 'kg m-3'
560
561       CASE ( 'hus' )
562          output_variable%long_name     = 'specific humidity'
563          output_variable%standard_name = 'specific_humidity'
564          output_variable%units         = 'kg kg-1'
565
566       CASE ( 'hur' )
567          output_variable%long_name     = 'relative humidity'
568          output_variable%standard_name = 'relative_humidity'
569          output_variable%units         = '1'
570
571       CASE ( 'rlu' )
572          output_variable%long_name     = 'upwelling longwave flux in air'
573          output_variable%standard_name = 'upwelling_longwave_flux_in_air'
574          output_variable%units         = 'W m-2'
575
576       CASE ( 'rlus' )
577          output_variable%long_name     = 'surface upwelling longwave flux in air'
578          output_variable%standard_name = 'surface_upwelling_longwave_flux_in_air'
579          output_variable%units         = 'W m-2'
580
581       CASE ( 'rld' )
582          output_variable%long_name     = 'downwelling longwave flux in air'
583          output_variable%standard_name = 'downwelling_longwave_flux_in_air'
584          output_variable%units         = 'W m-2'
585
586       CASE ( 'rsddif' )
587          output_variable%long_name     = 'diffuse downwelling shortwave flux in air'
588          output_variable%standard_name = 'diffuse_downwelling_shortwave_flux_in_air'
589          output_variable%units         = 'W m-2'
590
591       CASE ( 'rsd' )
592          output_variable%long_name     = 'downwelling shortwave flux in air'
593          output_variable%standard_name = 'downwelling_shortwave_flux_in_air'
594          output_variable%units         = 'W m-2'
595
596       CASE ( 'rnds' )
597          output_variable%long_name     = 'surface net downward radiative flux'
598          output_variable%standard_name = 'surface_net_downward_radiative_flux'
599          output_variable%units         = 'W m-2'
600
601       CASE ( 'rsu' )
602          output_variable%long_name     = 'upwelling shortwave flux in air'
603          output_variable%standard_name = 'upwelling_shortwave_flux_in_air'
604          output_variable%units         = 'W m-2'
605
606       CASE ( 'rsus' )
607          output_variable%long_name     = 'surface upwelling shortwave flux in air'
608          output_variable%standard_name = 'surface_upwelling_shortwave_flux_in_air'
609          output_variable%units         = 'W m-2'
610
611       CASE ( 'rsds' )
612          output_variable%long_name     = 'surface downwelling shortwave flux in air'
613          output_variable%standard_name = 'surface_downwelling_shortwave_flux_in_air'
614          output_variable%units         = 'W m-2'
615
616       CASE ( 'hfss' )
617          output_variable%long_name     = 'surface upward sensible heat flux'
618          output_variable%standard_name = 'surface_upward_sensible_heat_flux'
619          output_variable%units         = 'W m-2'
620
621       CASE ( 'hfls' )
622          output_variable%long_name     = 'surface upward latent heat flux'
623          output_variable%standard_name = 'surface_upward_latent_heat_flux'
624          output_variable%units         = 'W m-2'
625
626       CASE ( 'ts' )
627          output_variable%long_name     = 'surface temperature'
628          output_variable%standard_name = 'surface_temperature'
629          output_variable%units         = 'K'
630
631       CASE ( 'thetas' )
632          output_variable%long_name     = 'surface layer temperature scale'
633          output_variable%units         = 'K'
634
635       CASE ( 'us' )
636          output_variable%long_name     = 'friction velocity'
637          output_variable%units         = 'm s-1'
638
639       CASE ( 'uw' )
640          output_variable%long_name     = 'upward eastward kinematic momentum flux in air'
641          output_variable%units         = 'm2 s-2'
642
643       CASE ( 'vw' )
644          output_variable%long_name     = 'upward northward kinematic momentum flux in air'
645          output_variable%units         = 'm2 s-2'
646
647       CASE ( 'uv' )
648          output_variable%long_name     = 'eastward northward kinematic momentum flux in air'
649          output_variable%units         = 'm2 s-2'
650
651       CASE ( 'plev' )
652          output_variable%long_name     = 'air pressure'
653          output_variable%standard_name = 'air_pressure'
654          output_variable%units         = 'Pa'
655
656       CASE ( 'm_soil' )
657          output_variable%long_name     = 'soil moisture volumetric'
658          output_variable%units         = 'm3 m-3'
659
660       CASE ( 't_soil' )
661          output_variable%long_name     = 'soil temperature'
662          output_variable%standard_name = 'soil_temperature'
663          output_variable%units         = 'degree_C'
664
665       CASE ( 'hfdg' )
666          output_variable%long_name     = 'downward heat flux at ground level in soil'
667          output_variable%standard_name = 'downward_heat_flux_at_ground_level_in_soil'
668          output_variable%units         = 'W m-2'
669
670       CASE ( 'hfds' )
671          output_variable%long_name     = 'downward heat flux in soil'
672          output_variable%standard_name = 'downward_heat_flux_in_soil'
673          output_variable%units         = 'W m-2'
674
675       CASE ( 'hfla' )
676          output_variable%long_name     = 'upward latent heat flux in air'
677          output_variable%standard_name = 'upward_latent_heat_flux_in_air'
678          output_variable%units         = 'W m-2'
679
680       CASE ( 'hfsa' )
681          output_variable%long_name     = 'upward latent heat flux in air'
682          output_variable%standard_name = 'upward_sensible_heat_flux_in_air'
683          output_variable%units         = 'W m-2'
684
685       CASE ( 'jno2' )
686          output_variable%long_name     = 'photolysis rate of nitrogen dioxide'
687          output_variable%standard_name = 'photolysis_rate_of_nitrogen_dioxide'
688          output_variable%units         = 's-1'
689
690       CASE ( 'lwcs' )
691          output_variable%long_name     = 'liquid water content of soil layer'
692          output_variable%standard_name = 'liquid_water_content_of_soil_layer'
693          output_variable%units         = 'kg m-2'
694
695       CASE ( 'lwp' )
696          output_variable%long_name     = 'liquid water path'
697          output_variable%standard_name = 'atmosphere_mass_content_of_cloud_liquid_water'
698          output_variable%units         = 'kg m-2'
699
700       CASE ( 'ps' )
701          output_variable%long_name     = 'surface air pressure'
702          output_variable%standard_name = 'surface_air_pressure'
703          output_variable%units         = 'hPa'
704
705       CASE ( 'pswrtg' )
706          output_variable%long_name     = 'platform speed wrt ground'
707          output_variable%standard_name = 'platform_speed_wrt_ground'
708          output_variable%units         = 'm s-1'
709
710       CASE ( 'pswrta' )
711          output_variable%long_name     = 'platform speed wrt air'
712          output_variable%standard_name = 'platform_speed_wrt_air'
713          output_variable%units         = 'm s-1'
714
715       CASE ( 'pwv' )
716          output_variable%long_name     = 'water vapor partial pressure in air'
717          output_variable%standard_name = 'water_vapor_partial_pressure_in_air'
718          output_variable%units         = 'hPa'
719
720       CASE ( 'ssdu' )
721          output_variable%long_name     = 'duration of sunshine'
722          output_variable%standard_name = 'duration_of_sunshine'
723          output_variable%units         = 's'
724
725       CASE ( 't_lw' )
726          output_variable%long_name     = 'land water temperature'
727          output_variable%units         = 'degree_C'
728
729       CASE ( 'tb' )
730          output_variable%long_name     = 'brightness temperature'
731          output_variable%standard_name = 'brightness_temperature'
732          output_variable%units         = 'K'
733
734       CASE ( 'uqv' )
735          output_variable%long_name     = 'eastward kinematic latent heat flux in air'
736          output_variable%units         = 'g kg-1 m s-1'
737
738       CASE ( 'vqv' )
739          output_variable%long_name     = 'northward kinematic latent heat flux in air'
740          output_variable%units         = 'g kg-1 m s-1'
741
742       CASE ( 'wqv' )
743          output_variable%long_name     = 'upward kinematic latent heat flux in air'
744          output_variable%units         = 'g kg-1 m s-1'
745
746       CASE ( 'zcb' )
747          output_variable%long_name     = 'cloud base altitude'
748          output_variable%standard_name = 'cloud_base_altitude'
749          output_variable%units         = 'm'
750
751       CASE ( 'zmla' )
752          output_variable%long_name     = 'atmosphere boundary layer thickness'
753          output_variable%standard_name = 'atmosphere_boundary_layer_thickness'
754          output_variable%units         = 'm'
755
756       CASE ( 'mcpm1' )
757          output_variable%long_name     = 'mass concentration of pm1 ambient aerosol particles in air'
758          output_variable%standard_name = 'mass_concentration_of_pm1_ambient_aerosol_particles_in_air'
759          output_variable%units         = 'kg m-3'
760
761       CASE ( 'mcpm10' )
762          output_variable%long_name     = 'mass concentration of pm10 ambient aerosol particles in air'
763          output_variable%standard_name = 'mass_concentration_of_pm10_ambient_aerosol_particles_in_air'
764          output_variable%units         = 'kg m-3'
765
766       CASE ( 'mcpm2p5' )
767          output_variable%long_name     = 'mass concentration of pm2p5 ambient aerosol particles in air'
768          output_variable%standard_name = 'mass_concentration_of_pm2p5_ambient_aerosol_particles_in_air'
769          output_variable%units         = 'kg m-3'
770
771       CASE ( 'mfno', 'mcno'  )
772          output_variable%long_name     = 'mole fraction of nitrogen monoxide in air'
773          output_variable%standard_name = 'mole_fraction_of_nitrogen_monoxide_in_air'
774          output_variable%units         = 'ppm' !'mol mol-1'
775
776       CASE ( 'mfno2', 'mcno2'  )
777          output_variable%long_name     = 'mole fraction of nitrogen dioxide in air'
778          output_variable%standard_name = 'mole_fraction_of_nitrogen_dioxide_in_air'
779          output_variable%units         = 'ppm' !'mol mol-1'
780
781       CASE ( 'ncaa' )
782          output_variable%long_name     = 'number concentration of ambient aerosol particles in air'
783          output_variable%standard_name = 'number_concentration_of_ambient_aerosol_particles_in_air'
784          output_variable%units         = 'm-3' !'mol mol-1'
785
786       CASE ( 'tro3'  )
787          output_variable%long_name     = 'mole fraction of ozone in air'
788          output_variable%standard_name = 'mole_fraction_of_ozone_in_air'
789          output_variable%units         = 'ppm' !'mol mol-1'
790
791       CASE DEFAULT
792
793    END SELECT
794
795 END SUBROUTINE vm_set_attributes
796
797
798!--------------------------------------------------------------------------------------------------!
799! Description:
800! ------------
801!> Read namelist for the virtual measurement module
802!--------------------------------------------------------------------------------------------------!
803 SUBROUTINE vm_parin
804
805    CHARACTER(LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
806
807    NAMELIST /virtual_measurement_parameters/  dt_virtual_measurement,                             &
808                                               off_ts,                                             &
809                                               off_pr,                                             &
810                                               off_tr,                                             &
811                                               use_virtual_measurement,                            &
812                                               vm_time_start
813
814    line = ' '
815!
816!-- Try to find stg package
817    REWIND ( 11 )
818    line = ' '
819    DO  WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
820       READ ( 11, '(A)', END=20 )  line
821    ENDDO
822    BACKSPACE ( 11 )
823
824!
825!-- Read namelist
826    READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
827
828!
829!-- Set flag that indicates that the virtual measurement module is switched on
830    IF ( use_virtual_measurement )  virtual_measurement = .TRUE.
831    GOTO 20
832
833 10 BACKSPACE( 11 )
834    READ( 11 , '(A)') line
835    CALL parin_fail_message( 'virtual_measurement_parameters', line )
836
837 20 CONTINUE
838
839 END SUBROUTINE vm_parin
840
841
842!--------------------------------------------------------------------------------------------------!
843! Description:
844! ------------
845!> Initialize virtual measurements: read coordiante arrays and measured variables, set indicies
846!> indicating the measurement points, read further attributes, etc..
847!--------------------------------------------------------------------------------------------------!
848 SUBROUTINE vm_init
849
850    CHARACTER(LEN=5)                  ::  dum                           !< dummy string indicating station id
851    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables_file = ''  !< array with all measured variables read from NetCDF
852    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables      = ''  !< dummy array with all measured variables that are allowed
853
854    INTEGER(iwp) ::  dim_ntime  !< dimension size of time coordinate
855    INTEGER(iwp) ::  i          !< grid index of virtual observation point in x-direction
856    INTEGER(iwp) ::  is         !< grid index of real observation point of the respective station in x-direction
857    INTEGER(iwp) ::  j          !< grid index of observation point in x-direction
858    INTEGER(iwp) ::  js         !< grid index of real observation point of the respective station in y-direction
859    INTEGER(iwp) ::  k          !< grid index of observation point in x-direction
860    INTEGER(iwp) ::  kl         !< lower vertical index of surrounding grid points of an observation coordinate
861    INTEGER(iwp) ::  ks         !< grid index of real observation point of the respective station in z-direction
862    INTEGER(iwp) ::  ksurf      !< topography top index
863    INTEGER(iwp) ::  ku         !< upper vertical index of surrounding grid points of an observation coordinate
864    INTEGER(iwp) ::  l          !< running index over all stations
865    INTEGER(iwp) ::  len_char   !< character length of single measured variables without Null character
866    INTEGER(iwp) ::  ll         !< running index over all measured variables in file
867    INTEGER(iwp) ::  m          !< running index for surface elements
868    INTEGER(iwp) ::  n          !< running index over trajectory coordinates
869    INTEGER(iwp) ::  nofill     !< dummy for nofill return value (not used)
870    INTEGER(iwp) ::  ns         !< counter variable for number of observation points on subdomain
871    INTEGER(iwp) ::  off        !< number of surrounding grid points to be sampled
872    INTEGER(iwp) ::  t          !< running index over number of trajectories
873
874    INTEGER(KIND=1)                             ::  soil_dum  !< dummy variable to input a soil flag
875
876    INTEGER(iwp), DIMENSION(:), ALLOCATABLE     ::  ns_all  !< dummy array used to sum-up the number of observation coordinates
877
878#if defined( __netcdf4_parallel )
879    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_atmos  !< number of observation points for each station on each mpi rank
880    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_soil   !< number of observation points for each station on each mpi rank
881#endif
882
883    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag  !< mask array indicating measurement positions
884
885    LOGICAL  ::  on_pe  !< flag indicating that the respective measurement coordinate is on subdomain
886
887    REAL(wp) ::  fill_eutm !< _FillValue for coordinate array E_UTM
888    REAL(wp) ::  fill_nutm !< _FillValue for coordinate array N_UTM
889    REAL(wp) ::  fill_zar  !< _FillValue for height coordinate
890
891    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm      !< easting UTM coordinate, temporary variable
892    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm_tmp  !< EUTM coordinate before rotation
893    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm      !< northing UTM coordinate, temporary variable
894    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm_tmp  !< NUTM coordinate before rotation
895    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  station_h  !< station height above reference
896    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zar        !< observation height above reference
897#if defined( __netcdf )
898!
899!-- Open the input file.
900    CALL open_read_file( TRIM( input_file_vm ) // TRIM( coupling_char ), pids_id )
901!
902!-- Obtain number of sites.
903    CALL get_attribute( pids_id, char_numstations, vmea_general%nvm, global_attribute )
904!
905!-- Allocate data structure which encompasses all required information, such as  grid points indicies,
906!-- absolute UTM coordinates, the measured quantities, etc. .
907    ALLOCATE( vmea(1:vmea_general%nvm) )
908!
909!-- Allocate flag array. This dummy array is used to identify grid points where virtual measurements
910!-- should be taken. Please note, in order to include also the surrounding grid points of the
911!-- original coordinate, ghost points are required.
912    ALLOCATE( meas_flag(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
913    meas_flag = 0
914!
915!-- Loop over all sites in the setup file.
916    DO  l = 1, vmea_general%nvm
917!
918!--    Determine suffix which contains the ID, ordered according to the number of measurements.
919       IF( l < 10 )  THEN
920          WRITE( dum, '(I1)')  l
921       ELSEIF( l < 100 )  THEN
922          WRITE( dum, '(I2)')  l
923       ELSEIF( l < 1000 )  THEN
924          WRITE( dum, '(I3)')  l
925       ELSEIF( l < 10000 )  THEN
926          WRITE( dum, '(I4)')  l
927       ELSEIF( l < 100000 )  THEN
928          WRITE( dum, '(I5)')  l
929       ENDIF
930!
931!--    Read the origin site coordinates (UTM).
932       CALL get_attribute( pids_id, char_origx // TRIM( dum ), vmea(l)%origin_x_obs, global_attribute )
933       CALL get_attribute( pids_id, char_origy // TRIM( dum ), vmea(l)%origin_y_obs, global_attribute )
934!
935!--    Read site name.
936       CALL get_attribute( pids_id, char_site // TRIM( dum ), vmea(l)%site, global_attribute )
937!
938!--    Read a flag which indicates that also soil quantities are take at the respective site
939!--    (is part of the virtual measurement driver).
940       CALL get_attribute( pids_id, char_soil // TRIM( dum ), soil_dum, global_attribute )
941!
942!--    Set flag indicating soil-sampling.
943       IF ( soil_dum == 1 )  vmea(l)%soil_sampling = .TRUE.
944!
945!--    Read type of the measurement (trajectory, profile, timeseries).
946       CALL get_attribute( pids_id, char_feature // TRIM( dum ), vmea(l)%feature_type, global_attribute )
947!
948!---   Set logicals depending on the type of the measurement
949       IF ( INDEX( vmea(l)%feature_type, type_tspr     ) /= 0 )  THEN
950          vmea(l)%timseries_profile = .TRUE.
951       ELSEIF ( INDEX( vmea(l)%feature_type, type_ts   ) /= 0 )  THEN
952          vmea(l)%timseries         = .TRUE.
953       ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 )  THEN
954          vmea(l)%trajectory        = .TRUE.
955!
956!--    Give error message in case the type matches non of the pre-defined types.
957       ELSE
958          message_string = 'Attribue featureType = ' // TRIM( vmea(l)%feature_type ) // ' is not allowed.'
959          CALL message( 'vm_init', 'PA0535', 1, 2, 0, 6, 0 )
960       ENDIF
961!
962!--    Read string with all measured variables at this site.
963       measured_variables_file = ''
964       CALL get_variable( pids_id, char_mv // TRIM( dum ), measured_variables_file )
965!
966!--    Count the number of measured variables.
967!--    Please note, for some NetCDF interal reasons, characters end with a NULL, i.e. also empty
968!--    characters contain a NULL. Therefore, check the strings for a NULL to get the correct
969!--    character length in order to compare them with the list of allowed variables.
970       vmea(l)%nmeas  = 1
971       DO  ll = 1, SIZE( measured_variables_file )
972          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.                                  &
973               measured_variables_file(ll)(1:1) /= ' ')  THEN
974!
975!--          Obtain character length of the character
976             len_char = 1
977             DO  WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)  .AND.          &
978                 measured_variables_file(ll)(len_char:len_char) /= ' ' )
979                len_char = len_char + 1
980             ENDDO
981             len_char = len_char - 1
982
983             measured_variables(vmea(l)%nmeas) = measured_variables_file(ll)(1:len_char)
984             vmea(l)%nmeas = vmea(l)%nmeas + 1
985
986          ENDIF
987       ENDDO
988       vmea(l)%nmeas = vmea(l)%nmeas - 1
989!
990!--    Allocate data-type array for the measured variables names and attributes at the respective
991!--    site.
992       ALLOCATE( vmea(l)%var_atts(1:vmea(l)%nmeas) )
993!
994!--    Store the variable names in a data structure, which assigns further attributes to this name.
995!--    Further, for data output reasons, create a string of output variables, which will be written
996!--    into the attribute data_content.
997       DO  ll = 1, vmea(l)%nmeas
998          vmea(l)%var_atts(ll)%name = TRIM( measured_variables(ll) )
999
1000!           vmea(l)%data_content = TRIM( vmea(l)%data_content ) // " " //                            &
1001!                                  TRIM( vmea(l)%var_atts(ll)%name )
1002       ENDDO
1003!
1004!--    Read all the UTM coordinates for the site. Based on the coordinates, define the grid-index
1005!--    space on each subdomain where virtual measurements should be taken. Note, the entire
1006!--    coordinate array (on the entire model domain) won't be stored as this would exceed memory
1007!--    requirements, particularly for trajectories.
1008       IF ( vmea(l)%nmeas > 0 )  THEN
1009!
1010!--       For stationary measurements UTM coordinates are just one value and its dimension is
1011!--       "station", while for mobile measurements UTM coordinates are arrays depending on the
1012!--       number of trajectories and time, according to (UC)2 standard. First, inquire dimension
1013!--       length of the UTM coordinates.
1014          IF ( vmea(l)%trajectory )  THEN
1015!
1016!--          For non-stationary measurements read the number of trajectories and the number of time
1017!--          coordinates.
1018             CALL get_dimension_length( pids_id, vmea(l)%n_tr_st, "traj" // TRIM( dum ) )
1019             CALL get_dimension_length( pids_id, dim_ntime, "ntime" // TRIM( dum ) )
1020!
1021!--       For stationary measurements the dimension for UTM is station and for the time-coordinate
1022!--       it is one.
1023          ELSE
1024             CALL get_dimension_length( pids_id, vmea(l)%n_tr_st, "station" // TRIM( dum ) )
1025             dim_ntime = 1
1026          ENDIF
1027!
1028!-        Allocate array which defines individual time/space frame for each trajectory or station.
1029          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%n_tr_st) )
1030!
1031!--       Allocate temporary arrays for UTM and height coordinates. Note, on file UTM coordinates
1032!--       might be 1D or 2D variables
1033          ALLOCATE( e_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
1034          ALLOCATE( n_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
1035          ALLOCATE( station_h(1:vmea(l)%n_tr_st,1:dim_ntime)   )
1036          ALLOCATE( zar(1:vmea(l)%n_tr_st,1:dim_ntime)         )
1037          e_utm     = 0.0_wp
1038          n_utm     = 0.0_wp
1039          station_h = 0.0_wp
1040          zar       = 0.0_wp
1041
1042          ALLOCATE( e_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
1043          ALLOCATE( n_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
1044!
1045!--       Read UTM and height coordinates for all trajectories and times. Note, in case
1046!--       these obtain any missing values, replace them with default _FillValues.
1047          CALL inquire_fill_value( pids_id, char_eutm // TRIM( dum ), nofill, fill_eutm )
1048          CALL inquire_fill_value( pids_id, char_nutm // TRIM( dum ), nofill, fill_nutm )
1049          CALL inquire_fill_value( pids_id, char_zar // TRIM( dum ), nofill, fill_zar )
1050!
1051!--       Further line is just to avoid compiler warnings. nofill might be used in future.
1052          IF ( nofill == 0  .OR.  nofill /= 0 )  CONTINUE
1053!
1054!--       Read observation coordinates. Please note, for trajectories the observation height is
1055!--       stored directly in z, while for timeSeries it is stored in z - station_h, according to
1056!--       UC2-standard.
1057          IF ( vmea(l)%trajectory )  THEN
1058             CALL get_variable( pids_id, char_eutm // TRIM( dum ), e_utm, 0, dim_ntime-1, 0,       &
1059                                vmea(l)%n_tr_st-1 )
1060             CALL get_variable( pids_id, char_nutm // TRIM( dum ), n_utm, 0, dim_ntime-1, 0,       &
1061                                vmea(l)%n_tr_st-1 )
1062             CALL get_variable( pids_id, char_zar // TRIM( dum ), zar, 0, dim_ntime-1, 0,          &
1063                                vmea(l)%n_tr_st-1 )
1064          ELSE
1065             CALL get_variable( pids_id, char_eutm // TRIM( dum ), e_utm(:,1) )
1066             CALL get_variable( pids_id, char_nutm // TRIM( dum ), n_utm(:,1) )
1067             CALL get_variable( pids_id, char_station_h // TRIM( dum ), station_h(:,1) )
1068             CALL get_variable( pids_id, char_zar // TRIM( dum ), zar(:,1) )
1069          ENDIF
1070
1071          e_utm = MERGE( e_utm, vmea(l)%fillout, e_utm /= fill_eutm )
1072          n_utm = MERGE( n_utm, vmea(l)%fillout, n_utm /= fill_nutm )
1073          zar   = MERGE( zar,   vmea(l)%fillout, zar   /= fill_zar  )
1074!
1075!--       Compute observation height above ground.
1076          zar  = zar - station_h
1077!
1078!--       Based on UTM coordinates, check if the measurement station or parts of the trajectory are
1079!--       on subdomain. This case, setup grid index space sample these quantities.
1080          meas_flag = 0
1081          DO  t = 1, vmea(l)%n_tr_st
1082!
1083!--          First, compute relative x- and y-coordinates with respect to the lower-left origin of
1084!--          the model domain, which is the difference between UTM coordinates. Note, if the origin
1085!--          is not correct, the virtual sites will be misplaced. Further, in case of an rotated
1086!--          model domain, the UTM coordinates must also be rotated.
1087             e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
1088             n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
1089             e_utm(t,1:dim_ntime) = COS( init_model%rotation_angle * pi / 180.0_wp )               &
1090                                    * e_utm_tmp(t,1:dim_ntime)                                     &
1091                                  - SIN( init_model%rotation_angle * pi / 180.0_wp )               &
1092                                    * n_utm_tmp(t,1:dim_ntime)
1093             n_utm(t,1:dim_ntime) = SIN( init_model%rotation_angle * pi / 180.0_wp )               &
1094                                    * e_utm_tmp(t,1:dim_ntime)                                     &
1095                                  + COS( init_model%rotation_angle * pi / 180.0_wp )               &
1096                                    * n_utm_tmp(t,1:dim_ntime)
1097!
1098!--          Determine the individual time coordinate length for each station and trajectory. This
1099!--          is required as several stations and trajectories are merged into one file but they do
1100!--          not have the same number of points in time, hence, missing values may occur and cannot
1101!--          be processed further. This is actually a work-around for the specific (UC)2 dataset,
1102!--          but it won't harm anyway.
1103             vmea(l)%dim_t(t) = 0
1104             DO  n = 1, dim_ntime
1105                IF ( e_utm(t,n) /= fill_eutm  .AND.  n_utm(t,n) /= fill_nutm  .AND.                &
1106                     zar(t,n)   /= fill_zar )  vmea(l)%dim_t(t) = n
1107             ENDDO
1108!
1109!--          Compute grid indices relative to origin and check if these are on the subdomain. Note,
1110!--          virtual measurements will be taken also at grid points surrounding the station, hence,
1111!--          check also for these grid points. The number of surrounding grid points is set
1112!--          according to the featureType.
1113             IF ( vmea(l)%timseries_profile )  THEN
1114                off = off_pr
1115             ELSEIF ( vmea(l)%timseries     )  THEN
1116                off = off_ts
1117             ELSEIF ( vmea(l)%trajectory    )  THEN
1118                off = off_tr
1119             ENDIF
1120
1121             DO  n = 1, vmea(l)%dim_t(t)
1122                 is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
1123                 js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
1124!
1125!--             Is the observation point on subdomain?
1126                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.  js >= nys  .AND.  js <= nyn )
1127!
1128!--             Check if observation coordinate is on subdomain.
1129                IF ( on_pe )  THEN
1130!
1131!--                Determine vertical index which corresponds to the observation height.
1132                   ksurf = topo_top_ind(js,is,0)
1133                   ks = MINLOC( ABS( zu - zw(ksurf) - zar(t,n) ), DIM = 1 ) - 1
1134!
1135!--                Set mask array at the observation coordinates. Also, flag the surrounding
1136!--                coordinate points, but first check whether the surrounding coordinate points are
1137!--                on the subdomain.
1138                   kl = MERGE( ks-off, ksurf, ks-off >= nzb  .AND. ks-off >= ksurf )
1139                   ku = MERGE( ks+off, nzt,   ks+off < nzt+1 )
1140
1141                   DO  i = is-off, is+off
1142                      DO  j = js-off, js+off
1143                         DO  k = kl, ku
1144                            meas_flag(k,j,i) = MERGE( IBSET( meas_flag(k,j,i), 0 ), 0,             &
1145                                               BTEST( wall_flags_total_0(k,j,i), 0 ) )
1146                         ENDDO
1147                      ENDDO
1148                   ENDDO
1149                ENDIF
1150             ENDDO
1151
1152          ENDDO
1153!
1154!--       Based on the flag array, count the number of sampling coordinates. Please note, sampling
1155!--       coordinates in atmosphere and soil may be different, as within the soil all levels will be
1156!--       measured. Hence, count individually. Start with atmoshere.
1157          ns = 0
1158          DO  i = nxl-off, nxr+off
1159             DO  j = nys-off, nyn+off
1160                DO  k = nzb, nzt+1
1161                   ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
1162                ENDDO
1163             ENDDO
1164          ENDDO
1165
1166!
1167!--       Store number of observation points on subdomain and allocate index arrays as well as array
1168!--       containing height information.
1169          vmea(l)%ns = ns
1170
1171          ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
1172          ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
1173          ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
1174          ALLOCATE( vmea(l)%zar(1:vmea(l)%ns) )
1175!
1176!--       Based on the flag array store the grid indices which correspond to the observation
1177!--       coordinates.
1178          ns = 0
1179          DO  i = nxl-off, nxr+off
1180             DO  j = nys-off, nyn+off
1181                DO  k = nzb, nzt+1
1182                   IF ( BTEST( meas_flag(k,j,i), 0 ) )  THEN
1183                      ns = ns + 1
1184                      vmea(l)%i(ns) = i
1185                      vmea(l)%j(ns) = j
1186                      vmea(l)%k(ns) = k
1187                      vmea(l)%zar(ns)  = zu(k) - zw(topo_top_ind(j,i,0))
1188                   ENDIF
1189                ENDDO
1190             ENDDO
1191          ENDDO
1192!
1193!--       Same for the soil. Based on the flag array, count the number of sampling coordinates in
1194!--       soil. Sample at all soil levels in this case. Please note, soil variables can only be
1195!--       sampled on subdomains, not on ghost layers.
1196          IF ( vmea(l)%soil_sampling )  THEN
1197             DO  i = nxl, nxr
1198                DO  j = nys, nyn
1199                   IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
1200                      IF ( surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) )  THEN
1201                         vmea(l)%ns_soil = vmea(l)%ns_soil + nzt_soil - nzb_soil + 1
1202                      ENDIF
1203                      IF ( surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) )  THEN
1204                         vmea(l)%ns_soil = vmea(l)%ns_soil + nzt_wall - nzb_wall + 1
1205                      ENDIF
1206                   ENDIF
1207                ENDDO
1208             ENDDO
1209          ENDIF
1210!
1211!--       Allocate index arrays as well as array containing height information for soil.
1212          IF ( vmea(l)%soil_sampling )  THEN
1213             ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) )
1214             ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) )
1215             ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) )
1216             ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil)  )
1217          ENDIF
1218!
1219!--       For soil, store the grid indices.
1220          ns = 0
1221          IF ( vmea(l)%soil_sampling )  THEN
1222             DO  i = nxl, nxr
1223                DO  j = nys, nyn
1224                   IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
1225                      IF ( surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) )  THEN
1226                         m = surf_lsm_h%start_index(j,i)
1227                         DO  k = nzb_soil, nzt_soil
1228                            ns = ns + 1
1229                            vmea(l)%i_soil(ns) = i
1230                            vmea(l)%j_soil(ns) = j
1231                            vmea(l)%k_soil(ns) = k
1232                            vmea(l)%depth(ns)  = - zs(k)
1233                         ENDDO
1234                      ENDIF
1235
1236                      IF ( surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) )  THEN
1237                         m = surf_usm_h%start_index(j,i)
1238                         DO  k = nzb_wall, nzt_wall
1239                            ns = ns + 1
1240                            vmea(l)%i_soil(ns) = i
1241                            vmea(l)%j_soil(ns) = j
1242                            vmea(l)%k_soil(ns) = k
1243                            vmea(l)%depth(ns)  = - surf_usm_h%zw(k,m)
1244                         ENDDO
1245                      ENDIF
1246                   ENDIF
1247                ENDDO
1248             ENDDO
1249          ENDIF
1250!
1251!--       Allocate array to save the sampled values.
1252          ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) )
1253
1254          IF ( vmea(l)%soil_sampling )                                                             &
1255             ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil, 1:vmea(l)%nmeas) )
1256!
1257!--       Initialize with _FillValues
1258          vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) = vmea(l)%fillout
1259          IF ( vmea(l)%soil_sampling )                                                             &
1260             vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nmeas) = vmea(l)%fillout
1261!
1262!--       Deallocate temporary coordinate arrays
1263          IF ( ALLOCATED( e_utm )     )  DEALLOCATE( e_utm )
1264          IF ( ALLOCATED( n_utm )     )  DEALLOCATE( n_utm )
1265          IF ( ALLOCATED( e_utm_tmp ) )  DEALLOCATE( e_utm_tmp )
1266          IF ( ALLOCATED( n_utm_tmp ) )  DEALLOCATE( n_utm_tmp )
1267          IF ( ALLOCATED( n_utm )     )  DEALLOCATE( n_utm )
1268          IF ( ALLOCATED( zar  )      )  DEALLOCATE( vmea(l)%dim_t )
1269          IF ( ALLOCATED( zar  )      )  DEALLOCATE( zar  )
1270          IF ( ALLOCATED( station_h ) )  DEALLOCATE( station_h )
1271
1272       ENDIF
1273    ENDDO
1274!
1275!-- Dellocate flag array
1276    DEALLOCATE( meas_flag )
1277!
1278!-- Close input file for virtual measurements.
1279    CALL close_input_file( pids_id )
1280!
1281!-- Sum-up the number of observation coordiates, for atmosphere first.
1282!-- This is actually only required for data output.
1283    ALLOCATE( ns_all(1:vmea_general%nvm) )
1284    ns_all = 0
1285#if defined( __parallel )
1286    CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm,                                   &
1287                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
1288#else
1289    ns_all(:) = vmea(:)%ns
1290#endif
1291    vmea(:)%ns_tot = ns_all(:)
1292!
1293!-- Now for soil
1294    ns_all = 0
1295#if defined( __parallel )
1296    CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,                              &
1297                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
1298#else
1299    ns_all(:) = vmea(:)%ns_soil
1300#endif
1301    vmea(:)%ns_soil_tot = ns_all(:)
1302
1303    DEALLOCATE( ns_all )
1304!
1305!-- In case of parallel NetCDF the start coordinate for each mpi rank needs to be defined, so that
1306!-- each processor knows where to write the data.
1307#if defined( __netcdf4_parallel )
1308    ALLOCATE( ns_atmos(0:numprocs-1,1:vmea_general%nvm) )
1309    ALLOCATE( ns_soil(0:numprocs-1,1:vmea_general%nvm)  )
1310    ns_atmos = 0
1311    ns_soil  = 0
1312
1313    DO  l = 1, vmea_general%nvm
1314       ns_atmos(myid,l) = vmea(l)%ns
1315       ns_soil(myid,l)  = vmea(l)%ns_soil
1316    ENDDO
1317
1318#if defined( __parallel )
1319    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_atmos, numprocs * vmea_general%nvm,                       &
1320                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
1321    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_soil, numprocs * vmea_general%nvm,                        &
1322                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
1323#else
1324    ns_atmos(0,:) = vmea(:)%ns
1325    ns_soil(0,:)  = vmea(:)%ns_soil
1326#endif
1327
1328!
1329!-- Determine the start coordinate in NetCDF file for the local arrays. Note, start coordinates are
1330!-- initialized with zero for sake of simplicity in summation. However, in NetCDF the start
1331!-- coordinates must be >= 1, so that a one needs to be added at the end.
1332    DO  l = 1, vmea_general%nvm
1333       DO  n  = 0, myid - 1
1334          vmea(l)%start_coord_a = vmea(l)%start_coord_a + ns_atmos(n,l)
1335          vmea(l)%start_coord_s = vmea(l)%start_coord_s + ns_soil(n,l)
1336       ENDDO
1337!
1338!--    Start coordinate in NetCDF starts always at one not at 0.
1339       vmea(l)%start_coord_a = vmea(l)%start_coord_a + 1
1340       vmea(l)%start_coord_s = vmea(l)%start_coord_s + 1
1341!
1342!--    Determine the local end coordinate
1343       vmea(l)%end_coord_a = vmea(l)%start_coord_a + vmea(l)%ns - 1
1344       vmea(l)%end_coord_s = vmea(l)%start_coord_s + vmea(l)%ns_soil - 1
1345    ENDDO
1346
1347    DEALLOCATE( ns_atmos )
1348    DEALLOCATE( ns_soil  )
1349
1350#endif
1351
1352#endif
1353
1354 END SUBROUTINE vm_init
1355
1356
1357!--------------------------------------------------------------------------------------------------!
1358! Description:
1359! ------------
1360!> Initialize output using data-output module
1361!--------------------------------------------------------------------------------------------------!
1362 SUBROUTINE vm_init_output
1363
1364    CHARACTER(LEN=100) ::  variable_name  !< name of output variable
1365
1366    INTEGER(iwp) ::  l             !< loop index
1367    INTEGER(iwp) ::  n             !< loop index
1368    INTEGER      ::  return_value  !< returned status value of called function
1369
1370    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ndim !< dummy to write dimension
1371
1372    REAL(wp) ::  dum_lat  !< transformed geographical coordinate (latitude)
1373    REAL(wp) ::  dum_lon  !< transformed geographical coordinate (longitude)
1374
1375!
1376!-- Determine the number of output timesteps.
1377    ntimesteps = CEILING( ( end_time - MAX( vm_time_start, time_since_reference_point )            &
1378                          ) / dt_virtual_measurement )
1379!
1380!-- Create directory where output files will be stored.
1381    CALL local_system( 'mkdir -p VM_OUTPUT' // TRIM( coupling_char ) )
1382!
1383!-- Loop over all sites.
1384    DO  l = 1, vmea_general%nvm
1385!
1386!--    Skip if no observations will be taken for this site.
1387       IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
1388!
1389!--    Define output file.
1390       WRITE( vmea(l)%nc_filename, '(A,I4.4)' ) 'VM_OUTPUT' // TRIM( coupling_char ) // '/' //     &
1391              'site', l
1392
1393       return_value = dom_def_file( vmea(l)%nc_filename, 'netcdf4-parallel' )
1394!
1395!--    Define global attributes.
1396!--    Before, transform UTM into geographical coordinates.
1397       CALL convert_utm_to_geographic( crs_list, vmea(l)%origin_x_obs, vmea(l)%origin_y_obs,       &
1398                                       dum_lon, dum_lat )
1399
1400       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'site',                   &
1401                                   value = TRIM( vmea(l)%site ) )
1402       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'title',                  &
1403                                   value = 'Virtual measurement output')
1404       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'source',                 &
1405                                   value = 'PALM-4U')
1406       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'institution',            &
1407                                   value = input_file_atts%institution )
1408       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'acronym',                &
1409                                   value = input_file_atts%acronym )
1410       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'author',                 &
1411                                   value = input_file_atts%author )
1412       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'contact_person',         &
1413                                   value = input_file_atts%author )
1414       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'iop',                    &
1415                                   value = input_file_atts%campaign )
1416       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'campaign',               &
1417                                   value = 'PALM-4U' )
1418       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_time ',           &
1419                                   value = origin_date_time)
1420       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'location',               &
1421                                   value = input_file_atts%location )
1422       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_x',               &
1423                                   value = vmea(l)%origin_x_obs )
1424       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_y',               &
1425                                   value = vmea(l)%origin_y_obs )
1426       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_lon',             &
1427                                   value = dum_lon )
1428       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_lat',             &
1429                                   value = dum_lat )
1430       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_z', value = 0.0 )
1431       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'rotation_angle',         &
1432                                   value = input_file_atts%rotation_angle )
1433       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'featureType',            &
1434                                   value = TRIM( vmea(l)%feature_type_out ) )
1435       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'data_content',           &
1436                                   value = TRIM( vmea(l)%data_content ) )
1437       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'creation_time',          &
1438                                   value = input_file_atts%creation_time )
1439       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'version', value = 1 ) !input_file_atts%version
1440       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'Conventions',            &
1441                                   value = input_file_atts%conventions )
1442       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'dependencies',           &
1443                                   value = input_file_atts%dependencies )
1444       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'history',                &
1445                                   value = input_file_atts%history )
1446       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'references',             &
1447                                   value = input_file_atts%references )
1448       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'comment',                &
1449                                   value = input_file_atts%comment )
1450       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'keywords',               &
1451                                   value = input_file_atts%keywords )
1452       return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'licence',                &
1453                                   value = '[UC]2 Open Licence; see [UC]2 ' //                     &
1454                                           'data policy available at ' //                          &
1455                                           'www.uc2-program.org/uc2_data_policy.pdf' )
1456!
1457!--    Define dimensions.
1458!--    station
1459       ALLOCATE( ndim(1:vmea(l)%ns_tot) )
1460       DO  n = 1, vmea(l)%ns_tot
1461          ndim(n) = n
1462       ENDDO
1463       return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'station',                &
1464                                   output_type = 'int32', bounds = (/1_iwp, vmea(l)%ns_tot/),      &
1465                                   values_int32 = ndim )
1466       DEALLOCATE( ndim )
1467!
1468!--    ntime
1469       ALLOCATE( ndim(1:ntimesteps) )
1470       DO  n = 1, ntimesteps
1471          ndim(n) = n
1472       ENDDO
1473
1474       return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'ntime',                  &
1475                                   output_type = 'int32', bounds = (/1_iwp, ntimesteps/),          &
1476                                   values_int32 = ndim )
1477       DEALLOCATE( ndim )
1478!
1479!--    nv
1480       ALLOCATE( ndim(1:2) )
1481       DO  n = 1, 2
1482          ndim(n) = n
1483       ENDDO
1484
1485       return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'nv',                     &
1486                                   output_type = 'int32', bounds = (/1_iwp, 2_iwp/),               &
1487                                   values_int32 = ndim )
1488       DEALLOCATE( ndim )
1489!
1490!--    maximum name length
1491       ALLOCATE( ndim(1:maximum_name_length) )
1492       DO  n = 1, maximum_name_length
1493          ndim(n) = n
1494       ENDDO
1495
1496       return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'max_name_len',           &
1497                                   output_type = 'int32',                                          &
1498                                   bounds = (/1_iwp, maximum_name_length /), values_int32 = ndim )
1499       DEALLOCATE( ndim )
1500!
1501!--    Define coordinate variables.
1502!--    time
1503       variable_name = 'time'
1504       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1505                                   dimension_names = (/ 'station  ', 'ntime    '/),                &
1506                                   output_type = 'real32' )
1507!
1508!--    station_name
1509       variable_name = 'station_name'
1510       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1511                                   dimension_names = (/ 'max_name_len', 'station     ' /),         &
1512                                   output_type = 'char' )
1513!
1514!--    vrs (vertical reference system)
1515       variable_name = 'vrs'
1516       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1517                                   dimension_names = (/ 'station' /), output_type = 'int8' )
1518!
1519!--    crs (coordinate reference system)
1520       variable_name = 'crs'
1521       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1522                                   dimension_names = (/ 'station' /), output_type = 'int8' )
1523!
1524!--    z
1525       variable_name = 'z'
1526       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1527                                   dimension_names = (/'station'/), output_type = 'real32' )
1528!
1529!--    station_h
1530       variable_name = 'station_h'
1531       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1532                                   dimension_names = (/'station'/), output_type = 'real32' )
1533!
1534!--    x
1535       variable_name = 'x'
1536       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1537                                   dimension_names = (/'station'/), output_type = 'real32' )
1538!
1539!--    y
1540       variable_name = 'y'
1541       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1542                                   dimension_names = (/'station'/), output_type = 'real32' )
1543!
1544!--    E-UTM
1545       variable_name = 'E_UTM'
1546       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1547                                   dimension_names = (/'station'/), output_type = 'real32' )
1548!
1549!--    N-UTM
1550       variable_name = 'N_UTM'
1551       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1552                                   dimension_names = (/'station'/), output_type = 'real32' )
1553!
1554!--    latitude
1555       variable_name = 'lat'
1556       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1557                                   dimension_names = (/'station'/), output_type = 'real32' )
1558!
1559!--    longitude
1560       variable_name = 'lon'
1561       return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,             &
1562                                   dimension_names = (/'station'/), output_type = 'real32' )
1563!
1564!--    Set attributes for the coordinate variables. Note, not all coordinates have the same number
1565!--    of attributes.
1566!--    Units
1567       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time',                    &
1568                                   attribute_name = char_unit, value = 'seconds since ' //         &
1569                                   origin_date_time )
1570       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z',                       &
1571                                   attribute_name = char_unit, value = 'm' )
1572       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h',               &
1573                                   attribute_name = char_unit, value = 'm' )
1574       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x',                       &
1575                                   attribute_name = char_unit, value = 'm' )
1576       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y',                       &
1577                                   attribute_name = char_unit, value = 'm' )
1578       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM',                   &
1579                                   attribute_name = char_unit, value = 'm' )
1580       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM',                   &
1581                                   attribute_name = char_unit, value = 'm' )
1582       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat',                     &
1583                                   attribute_name = char_unit, value = 'degrees_north' )
1584       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon',                     &
1585                                   attribute_name = char_unit, value = 'degrees_east' )
1586!
1587!--    Long name
1588       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name',            &
1589                                   attribute_name = char_long, value = 'station name')
1590       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time',                    &
1591                                   attribute_name = char_long, value = 'time')
1592       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z',                       &
1593                                   attribute_name = char_long, value = 'height above origin' )
1594       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h',               &
1595                                   attribute_name = char_long, value = 'surface altitude' )
1596       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x',                       &
1597                                   attribute_name = char_long,                                     &
1598                                   value = 'distance to origin in x-direction')
1599       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y',                       &
1600                                   attribute_name = char_long,                                     &
1601                                   value = 'distance to origin in y-direction')
1602       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM',                   &
1603                                   attribute_name = char_long, value = 'easting' )
1604       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM',                   &
1605                                   attribute_name = char_long, value = 'northing' )
1606       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat',                     &
1607                                   attribute_name = char_long, value = 'latitude' )
1608       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon',                     &
1609                                   attribute_name = char_long, value = 'longitude' )
1610!
1611!--    Standard name
1612       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name',            &
1613                                   attribute_name = char_standard, value = 'platform_name')
1614       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time',                    &
1615                                   attribute_name = char_standard, value = 'time')
1616       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z',                       &
1617                                   attribute_name = char_standard,                                 &
1618                                   value = 'height_above_mean_sea_level' )
1619       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h',               &
1620                                   attribute_name = char_standard, value = 'surface_altitude' )
1621       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM',                   &
1622                                   attribute_name = char_standard,                                 &
1623                                   value = 'projection_x_coordinate' )
1624       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM',                   &
1625                                   attribute_name = char_standard,                                 &
1626                                   value = 'projection_y_coordinate' )
1627       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat',                     &
1628                                   attribute_name = char_standard, value = 'latitude' )
1629       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon',                     &
1630                                   attribute_name = char_standard, value = 'longitude' )
1631!
1632!--    Axis
1633       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time',                    &
1634                                   attribute_name = 'axis', value = 'T')
1635       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z',                       &
1636                                   attribute_name = 'axis', value = 'Z' )
1637       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x',                       &
1638                                   attribute_name = 'axis', value = 'X' )
1639       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y',                       &
1640                                   attribute_name = 'axis', value = 'Y' )
1641!
1642!--    Set further individual attributes for the coordinate variables.
1643!--    For station name
1644       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name',            &
1645                                   attribute_name = 'cf_role', value = 'timeseries_id' )
1646!
1647!--    For time
1648       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time',                    &
1649                                   attribute_name = 'calendar', value = 'proleptic_gregorian' )
1650!        return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time',                    &
1651!                                    attribute_name = 'bounds', value = 'time_bounds' )
1652!
1653!--    For vertical reference system
1654       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'vrs',                     &
1655                                   attribute_name = char_long, value = 'vertical reference system' )
1656       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'vrs',                     &
1657                                   attribute_name = 'system_name', value = 'DHHN2016' )
1658!
1659!--    For z
1660       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z',                       &
1661                                   attribute_name = 'positive', value = 'up' )
1662!
1663!--    For coordinate reference system
1664       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1665                                   attribute_name = 'epsg_code', value = coord_ref_sys%epsg_code )
1666       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1667                                   attribute_name = 'false_easting',                               &
1668                                   value = coord_ref_sys%false_easting )
1669       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1670                                   attribute_name = 'false_northing',                              &
1671                                   value = coord_ref_sys%false_northing )
1672       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1673                                   attribute_name = 'grid_mapping_name',                           &
1674                                   value = coord_ref_sys%grid_mapping_name )
1675       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1676                                   attribute_name = 'inverse_flattening',                          &
1677                                   value = coord_ref_sys%inverse_flattening )
1678       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1679                                   attribute_name = 'latitude_of_projection_origin',&
1680                                   value = coord_ref_sys%latitude_of_projection_origin )
1681       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1682                                   attribute_name = char_long, value = coord_ref_sys%long_name )
1683       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1684                                   attribute_name = 'longitude_of_central_meridian',               &
1685                                   value = coord_ref_sys%longitude_of_central_meridian )
1686       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1687                                   attribute_name = 'longitude_of_prime_meridian',                 &
1688                                   value = coord_ref_sys%longitude_of_prime_meridian )
1689       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1690                                   attribute_name = 'scale_factor_at_central_meridian',            &
1691                                   value = coord_ref_sys%scale_factor_at_central_meridian )
1692       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1693                                   attribute_name = 'semi_major_axis',                             &
1694                                   value = coord_ref_sys%semi_major_axis )
1695       return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs',                     &
1696                                   attribute_name = char_unit, value = coord_ref_sys%units )
1697!
1698!--    In case of sampled soil quantities, define further dimensions and coordinates.
1699       IF ( vmea(l)%soil_sampling )  THEN
1700!
1701!--       station for soil
1702          ALLOCATE( ndim(1:vmea(l)%ns_soil_tot) )
1703          DO  n = 1, vmea(l)%ns_soil_tot
1704             ndim(n) = n
1705          ENDDO
1706
1707          return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'station_soil',        &
1708                                      output_type = 'int32',                                       &
1709                                      bounds = (/1_iwp,vmea(l)%ns_soil_tot/), values_int32 = ndim )
1710          DEALLOCATE( ndim )
1711!
1712!--       ntime for soil
1713          ALLOCATE( ndim(1:ntimesteps) )
1714          DO  n = 1, ntimesteps
1715             ndim(n) = n
1716          ENDDO
1717
1718          return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'ntime_soil',          &
1719                                      output_type = 'int32', bounds = (/1_iwp,ntimesteps/),        &
1720                                      values_int32 = ndim )
1721          DEALLOCATE( ndim )
1722!
1723!--       time for soil
1724          variable_name = 'time_soil'
1725          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1726                                      dimension_names = (/'station_soil', 'ntime_soil  '/),        &
1727                                      output_type = 'real32' )
1728!
1729!--       station_name for soil
1730          variable_name = 'station_name_soil'
1731          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1732                                      dimension_names = (/ 'max_name_len', 'station_soil' /),      &
1733                                      output_type = 'char' )
1734!
1735!--       z
1736          variable_name = 'z_soil'
1737          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1738                                      dimension_names = (/'station_soil'/), output_type = 'real32' )
1739!
1740!--       station_h for soil
1741          variable_name = 'station_h_soil'
1742          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1743                                      dimension_names = (/'station_soil'/), output_type = 'real32' )
1744!
1745!--       x soil
1746          variable_name = 'x_soil'
1747          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1748                                      dimension_names = (/'station_soil'/), output_type = 'real32' )
1749!
1750!-        y soil
1751          variable_name = 'y_soil'
1752          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1753                                      dimension_names = (/'station_soil'/), output_type = 'real32' )
1754!
1755!--       E-UTM soil
1756          variable_name = 'E_UTM_soil'
1757          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1758                                      dimension_names = (/'station_soil'/), output_type = 'real32' )
1759!
1760!--       N-UTM soil
1761          variable_name = 'N_UTM_soil'
1762          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1763                                      dimension_names = (/'station_soil'/), output_type = 'real32' )
1764!
1765!--       latitude soil
1766          variable_name = 'lat_soil'
1767          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1768                                      dimension_names = (/'station_soil'/), output_type = 'real32' )
1769!
1770!--       longitude soil
1771          variable_name = 'lon_soil'
1772          return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,          &
1773                                      dimension_names = (/'station_soil'/), output_type = 'real32' )
1774!
1775!--       Set attributes for the coordinate variables. Note, not all coordinates have the same
1776!--       number of attributes.
1777!--       Units
1778          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil',            &
1779                                      attribute_name = char_unit, value = 'seconds since ' //      &
1780                                      origin_date_time )
1781          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil',               &
1782                                      attribute_name = char_unit, value = 'm' )
1783          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h_soil',       &
1784                                      attribute_name = char_unit, value = 'm' )
1785          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x_soil',               &
1786                                      attribute_name = char_unit, value = 'm' )
1787          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y_soil',               &
1788                                      attribute_name = char_unit, value = 'm' )
1789          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM_soil',           &
1790                                      attribute_name = char_unit, value = 'm' )
1791          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM_soil',           &
1792                                      attribute_name = char_unit, value = 'm' )
1793          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat_soil',             &
1794                                      attribute_name = char_unit, value = 'degrees_north' )
1795          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon_soil',             &
1796                                      attribute_name = char_unit, value = 'degrees_east' )
1797!
1798!--       Long name
1799          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name_soil',    &
1800                                      attribute_name = char_long, value = 'station name')
1801          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil',            &
1802                                      attribute_name = char_long, value = 'time')
1803          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil',               &
1804                                      attribute_name = char_long, value = 'height above origin' )
1805          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h_soil',       &
1806                                      attribute_name = char_long, value = 'surface altitude' )
1807          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x_soil',               &
1808                                      attribute_name = char_long,                                  &
1809                                      value = 'distance to origin in x-direction' )
1810          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y_soil',               &
1811                                      attribute_name = char_long,                                  &
1812                                      value = 'distance to origin in y-direction' )
1813          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM_soil',           &
1814                                      attribute_name = char_long, value = 'easting' )
1815          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM_soil',           &
1816                                      attribute_name = char_long, value = 'northing' )
1817          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat_soil',             &
1818                                      attribute_name = char_long, value = 'latitude' )
1819          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon_soil',             &
1820                                      attribute_name = char_long, value = 'longitude' )
1821!
1822!--       Standard name
1823          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name_soil',    &
1824                                      attribute_name = char_standard, value = 'platform_name')
1825          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil',            &
1826                                      attribute_name = char_standard, value = 'time')
1827          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil',               &
1828                                      attribute_name = char_standard,                              &
1829                                      value = 'height_above_mean_sea_level' )
1830          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h_soil',       &
1831                                      attribute_name = char_standard, value = 'surface_altitude' )
1832          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM_soil',           &
1833                                      attribute_name = char_standard,                              &
1834                                      value = 'projection_x_coordinate' )
1835          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM_soil',           &
1836                                      attribute_name = char_standard,                              &
1837                                      value = 'projection_y_coordinate' )
1838          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat_soil',             &
1839                                      attribute_name = char_standard, value = 'latitude' )
1840          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon_soil',             &
1841                                      attribute_name = char_standard, value = 'longitude' )
1842!
1843!--       Axis
1844          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil',            &
1845                                      attribute_name = 'axis', value = 'T')
1846          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil',               &
1847                                      attribute_name = 'axis', value = 'Z' )
1848          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x_soil',               &
1849                                      attribute_name = 'axis', value = 'X' )
1850          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y_soil',               &
1851                                      attribute_name = 'axis', value = 'Y' )
1852!
1853!--       Set further individual attributes for the coordinate variables.
1854!--       For station name soil
1855          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name_soil',    &
1856                                      attribute_name = 'cf_role', value = 'timeseries_id' )
1857!
1858!--       For time soil
1859          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil',            &
1860                                      attribute_name = 'calendar', value = 'proleptic_gregorian' )
1861!           return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil',            &
1862!                                       attribute_name = 'bounds', value = 'time_bounds' )
1863!
1864!--       For z soil
1865          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil',               &
1866                                      attribute_name = 'positive', value = 'up' )
1867       ENDIF
1868!
1869!--    Define variables that shall be sampled.
1870       DO  n = 1, vmea(l)%nmeas
1871          variable_name = TRIM( vmea(l)%var_atts(n)%name )
1872!
1873!--       In order to link the correct dimension names, atmosphere and soil variables need to be
1874!--       distinguished.
1875          IF ( vmea(l)%soil_sampling  .AND.                                                        &
1876               ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
1877
1878             return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,       &
1879                                         dimension_names = (/'station_soil', 'ntime_soil  '/),     &
1880                                         output_type = 'real32' )
1881          ELSE
1882
1883             return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name,       &
1884                                         dimension_names = (/'station', 'ntime  '/),               &
1885                                         output_type = 'real32' )
1886          ENDIF
1887!
1888!--       Set variable attributes. Please note, for some variables not all attributes are defined,
1889!--       e.g. standard_name for the horizontal wind components.
1890          CALL vm_set_attributes( vmea(l)%var_atts(n) )
1891
1892          IF ( vmea(l)%var_atts(n)%long_name /= 'none' )  THEN
1893             return_value = dom_def_att( vmea(l)%nc_filename,  variable_name = variable_name,      &
1894                                         attribute_name = char_long,                               &
1895                                         value = TRIM( vmea(l)%var_atts(n)%long_name ) )
1896          ENDIF
1897          IF ( vmea(l)%var_atts(n)%standard_name /= 'none' )  THEN
1898             return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name,       &
1899                                         attribute_name = char_standard,                           &
1900                                         value = TRIM( vmea(l)%var_atts(n)%standard_name ) )
1901          ENDIF
1902          IF ( vmea(l)%var_atts(n)%units /= 'none' )  THEN
1903             return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name,       &
1904                                         attribute_name = char_unit,                               &
1905                                         value = TRIM( vmea(l)%var_atts(n)%units ) )
1906          ENDIF
1907
1908          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name,          &
1909                                      attribute_name = 'grid_mapping',                             &
1910                                      value = TRIM( vmea(l)%var_atts(n)%grid_mapping ) )
1911
1912          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name,          &
1913                                      attribute_name = 'coordinates',                              &
1914                                      value = TRIM( vmea(l)%var_atts(n)%coordinates ) )
1915
1916          return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name,          &
1917                                      attribute_name = char_fill,                                  &
1918                                      value = REAL( vmea(l)%var_atts(n)%fill_value, KIND=4 ) )
1919
1920       ENDDO  ! loop over variables per site
1921
1922    ENDDO  ! loop over sites
1923
1924
1925 END SUBROUTINE vm_init_output
1926
1927!--------------------------------------------------------------------------------------------------!
1928! Description:
1929! ------------
1930!> Parallel NetCDF output via data-output module.
1931!--------------------------------------------------------------------------------------------------!
1932 SUBROUTINE vm_data_output
1933
1934    CHARACTER(LEN=100) ::  variable_name  !< name of output variable
1935    CHARACTER(LEN=maximum_name_length), DIMENSION(:), ALLOCATABLE :: station_name  !< string for station name, consecutively ordered
1936
1937    CHARACTER(LEN=1), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_char_target  !< target for output name arrays
1938    CHARACTER(LEN=1), DIMENSION(:,:), POINTER             ::  output_values_2d_char_pointer !< pointer for output name arrays
1939
1940    INTEGER(iwp)       ::  l             !< loop index for the number of sites
1941    INTEGER(iwp)       ::  n             !< loop index for observation points
1942    INTEGER(iwp)       ::  nn            !< loop index for number of characters in a name
1943    INTEGER            ::  return_value  !< returned status value of called function
1944    INTEGER(iwp)       ::  t_ind         !< time index
1945
1946    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  dum_lat                   !< transformed geographical coordinate (latitude)
1947    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  dum_lon                   !< transformed geographical coordinate (longitude)
1948    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  oro_rel                   !< relative altitude of model surface
1949    REAL(wp), DIMENSION(:), POINTER               ::  output_values_1d_pointer  !< pointer for 1d output array
1950    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET   ::  output_values_1d_target   !< target for 1d output array
1951    REAL(wp), DIMENSION(:,:), POINTER             ::  output_values_2d_pointer  !< pointer for 2d output array
1952    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_target   !< target for 2d output array
1953
1954    CALL cpu_log( log_point_s(26), 'VM output', 'start' )
1955!
1956!-- At the first call of this routine write the spatial coordinates.
1957    IF ( .NOT. initial_write_coordinates )  THEN
1958!
1959!--    Write spatial coordinates.
1960       DO  l = 1, vmea_general%nvm
1961!
1962!--       Skip if no observations were taken.
1963          IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
1964
1965          ALLOCATE( output_values_1d_target(vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
1966!
1967!--       Output of Easting coordinate. Before output, recalculate EUTM.
1968          output_values_1d_target = init_model%origin_x                                            &
1969                    + REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx                     &
1970                    * COS( init_model%rotation_angle * pi / 180.0_wp )                             &
1971                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy                     &
1972                    * SIN( init_model%rotation_angle * pi / 180.0_wp )
1973
1974          output_values_1d_pointer => output_values_1d_target
1975
1976          return_value = dom_write_var( vmea(l)%nc_filename, 'E_UTM',                              &
1977                                        values_realwp_1d = output_values_1d_pointer,               &
1978                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
1979                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
1980!
1981!--       Output of Northing coordinate. Before output, recalculate NUTM.
1982          output_values_1d_target = init_model%origin_y                                            &
1983                    - REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx                     &
1984                    * SIN( init_model%rotation_angle * pi / 180.0_wp )                             &
1985                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy                     &
1986                    * COS( init_model%rotation_angle * pi / 180.0_wp )
1987
1988          output_values_1d_pointer => output_values_1d_target
1989          return_value = dom_write_var( vmea(l)%nc_filename, 'N_UTM',                              &
1990                                        values_realwp_1d = output_values_1d_pointer,               &
1991                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
1992                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
1993!
1994!--       Output of longitude and latitude coordinate. Before output, convert it.
1995          ALLOCATE( dum_lat(1:vmea(l)%ns) )
1996          ALLOCATE( dum_lon(1:vmea(l)%ns) )
1997
1998          DO  n = 1, vmea(l)%ns
1999             CALL convert_utm_to_geographic( crs_list,                                             &
2000                                             init_model%origin_x                                   &
2001                                           + REAL( vmea(l)%i(n) + 0.5_wp, KIND = wp ) * dx         &
2002                                           * COS( init_model%rotation_angle * pi / 180.0_wp )      &
2003                                           + REAL( vmea(l)%j(n) + 0.5_wp, KIND = wp ) * dy         &
2004                                           * SIN( init_model%rotation_angle * pi / 180.0_wp ),     &
2005                                             init_model%origin_y                                   &
2006                                           - REAL( vmea(l)%i(n) + 0.5_wp, KIND = wp ) * dx         &
2007                                           * SIN( init_model%rotation_angle * pi / 180.0_wp )      &
2008                                           + REAL( vmea(l)%j(n) + 0.5_wp, KIND = wp ) * dy         &
2009                                           * COS( init_model%rotation_angle * pi / 180.0_wp ),     &
2010                                             dum_lon(n), dum_lat(n) )
2011          ENDDO
2012
2013          output_values_1d_target = dum_lat
2014          output_values_1d_pointer => output_values_1d_target
2015          return_value = dom_write_var( vmea(l)%nc_filename, 'lat',                                &
2016                                        values_realwp_1d = output_values_1d_pointer,               &
2017                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
2018                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
2019
2020          output_values_1d_target = dum_lon
2021          output_values_1d_pointer => output_values_1d_target
2022          return_value = dom_write_var( vmea(l)%nc_filename, 'lon',                                &
2023                                        values_realwp_1d = output_values_1d_pointer,               &
2024                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
2025                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
2026          DEALLOCATE( dum_lat )
2027          DEALLOCATE( dum_lon )
2028!
2029!--       Output of relative height coordinate.
2030!--       Before this is output, first define the relative orographie height and add this to z.
2031          ALLOCATE( oro_rel(1:vmea(l)%ns) )
2032          DO  n = 1, vmea(l)%ns
2033             oro_rel(n) = zw(topo_top_ind(vmea(l)%j(n),vmea(l)%i(n),3))
2034          ENDDO
2035
2036          output_values_1d_target = vmea(l)%zar(1:vmea(l)%ns) + oro_rel(:)
2037          output_values_1d_pointer => output_values_1d_target
2038          return_value = dom_write_var( vmea(l)%nc_filename, 'z',                                  &
2039                                        values_realwp_1d = output_values_1d_pointer,               &
2040                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
2041                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
2042!
2043!--       Write surface altitude for the station. Note, since z is already a relative observation
2044!--       height, station_h must be zero, in order to obtain the observation level.
2045          output_values_1d_target = oro_rel(:)
2046          output_values_1d_pointer => output_values_1d_target
2047          return_value = dom_write_var( vmea(l)%nc_filename, 'station_h',                          &
2048                                        values_realwp_1d = output_values_1d_pointer,               &
2049                                        bounds_start = (/vmea(l)%start_coord_a/),                  &
2050                                        bounds_end   = (/vmea(l)%end_coord_a  /) )
2051
2052          DEALLOCATE( oro_rel )
2053          DEALLOCATE( output_values_1d_target )
2054!
2055!--       Write station name
2056          ALLOCATE ( station_name(vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
2057          ALLOCATE ( output_values_2d_char_target(vmea(l)%start_coord_a:vmea(l)%end_coord_a, &
2058                                                  1:maximum_name_length) )
2059
2060          DO  n = vmea(l)%start_coord_a, vmea(l)%end_coord_a
2061             station_name(n) = REPEAT( ' ', maximum_name_length )
2062             WRITE( station_name(n), '(A,I10.10)') "station", n
2063             DO  nn = 1, maximum_name_length
2064                output_values_2d_char_target(n,nn) = station_name(n)(nn:nn)
2065             ENDDO
2066          ENDDO
2067
2068          output_values_2d_char_pointer => output_values_2d_char_target
2069
2070          return_value = dom_write_var( vmea(l)%nc_filename, 'station_name',                       &
2071                                        values_char_2d = output_values_2d_char_pointer,            &
2072                                        bounds_start = (/ 1, vmea(l)%start_coord_a /),             &
2073                                        bounds_end   = (/ maximum_name_length,                     &
2074                                        vmea(l)%end_coord_a /) )
2075
2076          DEALLOCATE( station_name )
2077          DEALLOCATE( output_values_2d_char_target )
2078!
2079!--       In case of sampled soil quantities, output also the respective coordinate arrays.
2080          IF ( vmea(l)%soil_sampling )  THEN
2081             ALLOCATE( output_values_1d_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
2082!
2083!--          Output of Easting coordinate. Before output, recalculate EUTM.
2084             output_values_1d_target = init_model%origin_x                                         &
2085               + REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx                     &
2086               * COS( init_model%rotation_angle * pi / 180.0_wp )                                  &
2087               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy                     &
2088               * SIN( init_model%rotation_angle * pi / 180.0_wp )
2089             output_values_1d_pointer => output_values_1d_target
2090             return_value = dom_write_var( vmea(l)%nc_filename, 'E_UTM_soil',                      &
2091                                           values_realwp_1d = output_values_1d_pointer,            &
2092                                           bounds_start = (/vmea(l)%start_coord_s/),               &
2093                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
2094!
2095!--          Output of Northing coordinate. Before output, recalculate NUTM.
2096             output_values_1d_target = init_model%origin_y                                         &
2097               - REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx                &
2098               * SIN( init_model%rotation_angle * pi / 180.0_wp )                                  &
2099               + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy                &
2100               * COS( init_model%rotation_angle * pi / 180.0_wp )
2101
2102             output_values_1d_pointer => output_values_1d_target
2103             return_value = dom_write_var( vmea(l)%nc_filename, 'N_UTM_soil',                      &
2104                                           values_realwp_1d = output_values_1d_pointer,            &
2105                                           bounds_start = (/vmea(l)%start_coord_s/),               &
2106                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
2107!
2108!--          Output of longitude and latitude coordinate. Before output, convert it.
2109             ALLOCATE( dum_lat(1:vmea(l)%ns_soil) )
2110             ALLOCATE( dum_lon(1:vmea(l)%ns_soil) )
2111
2112             DO  n = 1, vmea(l)%ns_soil
2113                CALL convert_utm_to_geographic( crs_list,                                          &
2114                                                init_model%origin_x                                &
2115                                              + REAL( vmea(l)%i_soil(n) + 0.5_wp, KIND = wp ) * dx &
2116                                              * COS( init_model%rotation_angle * pi / 180.0_wp )   &
2117                                              + REAL( vmea(l)%j_soil(n) + 0.5_wp, KIND = wp ) * dy &
2118                                              * SIN( init_model%rotation_angle * pi / 180.0_wp ),  &
2119                                                init_model%origin_y                                &
2120                                              - REAL( vmea(l)%i_soil(n) + 0.5_wp, KIND = wp ) * dx &
2121                                              * SIN( init_model%rotation_angle * pi / 180.0_wp )   &
2122                                              + REAL( vmea(l)%j_soil(n) + 0.5_wp, KIND = wp ) * dy &
2123                                              * COS( init_model%rotation_angle * pi / 180.0_wp ),  &
2124                                                dum_lon(n), dum_lat(n) )
2125             ENDDO
2126
2127             output_values_1d_target = dum_lat
2128             output_values_1d_pointer => output_values_1d_target
2129             return_value = dom_write_var( vmea(l)%nc_filename, 'lat_soil',                        &
2130                                           values_realwp_1d = output_values_1d_pointer,            &
2131                                           bounds_start = (/vmea(l)%start_coord_s/),               &
2132                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
2133
2134             output_values_1d_target = dum_lon
2135             output_values_1d_pointer => output_values_1d_target
2136             return_value = dom_write_var( vmea(l)%nc_filename, 'lon_soil',                        &
2137                                           values_realwp_1d = output_values_1d_pointer,            &
2138                                           bounds_start = (/vmea(l)%start_coord_s/),               &
2139                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
2140             DEALLOCATE( dum_lat )
2141             DEALLOCATE( dum_lon )
2142!
2143!--          Output of relative height coordinate.
2144!--          Before this is output, first define the relative orographie height and add this to z.
2145             ALLOCATE( oro_rel(1:vmea(l)%ns_soil) )
2146             DO  n = 1, vmea(l)%ns_soil
2147                oro_rel(n) = zw(topo_top_ind(vmea(l)%j_soil(n),vmea(l)%i_soil(n),3))
2148             ENDDO
2149
2150             output_values_1d_target = vmea(l)%depth(1:vmea(l)%ns_soil) + oro_rel(:)
2151             output_values_1d_pointer => output_values_1d_target
2152             return_value = dom_write_var( vmea(l)%nc_filename, 'z_soil',                          &
2153                                           values_realwp_1d = output_values_1d_pointer,            &
2154                                           bounds_start = (/vmea(l)%start_coord_s/),               &
2155                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
2156!
2157!--          Write surface altitude for the station. Note, since z is already a relative observation
2158!--          height, station_h must be zero, in order to obtain the observation level.
2159             output_values_1d_target = oro_rel(:)
2160             output_values_1d_pointer => output_values_1d_target
2161             return_value = dom_write_var( vmea(l)%nc_filename, 'station_h_soil',                  &
2162                                           values_realwp_1d = output_values_1d_pointer,            &
2163                                           bounds_start = (/vmea(l)%start_coord_s/),               &
2164                                           bounds_end   = (/vmea(l)%end_coord_s  /) )
2165
2166             DEALLOCATE( oro_rel )
2167             DEALLOCATE( output_values_1d_target )
2168!
2169!--          Write station name
2170             ALLOCATE ( station_name(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
2171             ALLOCATE ( output_values_2d_char_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s,    &
2172                                                     1:maximum_name_length) )
2173
2174             DO  n = vmea(l)%start_coord_s, vmea(l)%end_coord_s
2175                station_name(n) = REPEAT( ' ', maximum_name_length )
2176                WRITE( station_name(n), '(A,I10.10)') "station", n
2177                DO  nn = 1, maximum_name_length
2178                   output_values_2d_char_target(n,nn) = station_name(n)(nn:nn)
2179                ENDDO
2180             ENDDO
2181             output_values_2d_char_pointer => output_values_2d_char_target
2182
2183             return_value = dom_write_var( vmea(l)%nc_filename, 'station_name_soil',               &
2184                                           values_char_2d = output_values_2d_char_pointer,         &
2185                                           bounds_start = (/ 1, vmea(l)%start_coord_s /),          &
2186                                           bounds_end   = (/ maximum_name_length,                  &
2187                                           vmea(l)%end_coord_s   /) )
2188
2189             DEALLOCATE( station_name )
2190             DEALLOCATE( output_values_2d_char_target )
2191
2192          ENDIF
2193
2194       ENDDO  ! loop over sites
2195
2196       initial_write_coordinates = .TRUE.
2197    ENDIF
2198!
2199!-- Loop over all sites.
2200    DO  l = 1, vmea_general%nvm
2201!
2202!--    Skip if no observations were taken.
2203       IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
2204!
2205!--    Determine time index in file.
2206       t_ind = vmea(l)%file_time_index + 1
2207!
2208!--    Write output variables. Distinguish between atmosphere and soil variables.
2209       DO  n = 1, vmea(l)%nmeas
2210          IF ( vmea(l)%soil_sampling  .AND.                                                        &
2211            ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
2212!
2213!--          Write time coordinate to file
2214             variable_name = 'time_soil'
2215             ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
2216             output_values_2d_target(t_ind,:) = time_since_reference_point
2217             output_values_2d_pointer => output_values_2d_target
2218
2219             return_value = dom_write_var( vmea(l)%nc_filename, variable_name,                     &
2220                                           values_realwp_2d = output_values_2d_pointer,            &
2221                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/),        &
2222                                           bounds_end   = (/vmea(l)%end_coord_s, t_ind /) )
2223
2224             variable_name = TRIM( vmea(l)%var_atts(n)%name )
2225             output_values_2d_target(t_ind,:) = vmea(l)%measured_vars_soil(:,n)
2226             output_values_2d_pointer => output_values_2d_target
2227             return_value = dom_write_var( vmea(l)%nc_filename, variable_name,                     &
2228                                           values_realwp_2d = output_values_2d_pointer,            &
2229                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/),        &
2230                                           bounds_end   = (/vmea(l)%end_coord_s, t_ind  /) )
2231             DEALLOCATE( output_values_2d_target )
2232          ELSE
2233!
2234!--          Write time coordinate to file
2235             variable_name = 'time'
2236             ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
2237             output_values_2d_target(t_ind,:) = time_since_reference_point
2238             output_values_2d_pointer => output_values_2d_target
2239
2240             return_value = dom_write_var( vmea(l)%nc_filename, variable_name,                     &
2241                                           values_realwp_2d = output_values_2d_pointer,            &
2242                                           bounds_start = (/vmea(l)%start_coord_a, t_ind/),        &
2243                                           bounds_end   = (/vmea(l)%end_coord_a, t_ind/) )
2244
2245             variable_name = TRIM( vmea(l)%var_atts(n)%name )
2246
2247             output_values_2d_target(t_ind,:) = vmea(l)%measured_vars(:,n)
2248             output_values_2d_pointer => output_values_2d_target
2249             return_value = dom_write_var( vmea(l)%nc_filename, variable_name,                     &
2250                                           values_realwp_2d = output_values_2d_pointer,            &
2251                                           bounds_start = (/ vmea(l)%start_coord_a, t_ind /),      &
2252                                           bounds_end   = (/ vmea(l)%end_coord_a, t_ind /) )
2253
2254             DEALLOCATE( output_values_2d_target )
2255          ENDIF
2256       ENDDO
2257!
2258!--    Update number of written time indices
2259       vmea(l)%file_time_index = t_ind
2260
2261    ENDDO  ! loop over sites
2262
2263    CALL cpu_log( log_point_s(26), 'VM output', 'stop' )
2264
2265
2266 END SUBROUTINE vm_data_output
2267
2268!--------------------------------------------------------------------------------------------------!
2269! Description:
2270! ------------
2271!> Sampling of the actual quantities along the observation coordinates
2272!--------------------------------------------------------------------------------------------------!
2273 SUBROUTINE vm_sampling
2274
2275    USE radiation_model_mod,                                                                       &
2276        ONLY:  radiation
2277
2278    USE surface_mod,                                                                               &
2279        ONLY:  surf_def_h,                                                                         &
2280               surf_lsm_h,                                                                         &
2281               surf_usm_h
2282
2283     INTEGER(iwp) ::  i         !< grid index in x-direction
2284     INTEGER(iwp) ::  j         !< grid index in y-direction
2285     INTEGER(iwp) ::  k         !< grid index in z-direction
2286     INTEGER(iwp) ::  ind_chem  !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming
2287     INTEGER(iwp) ::  l         !< running index over the number of stations
2288     INTEGER(iwp) ::  m         !< running index over all virtual observation coordinates
2289     INTEGER(iwp) ::  mm        !< index of surface element which corresponds to the virtual observation coordinate
2290     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
2291     INTEGER(iwp) ::  nn        !< running index over the number of chemcal species
2292
2293     LOGICAL ::  match_lsm  !< flag indicating natural-type surface
2294     LOGICAL ::  match_usm  !< flag indicating urban-type surface
2295
2296     REAL(wp) ::  e_s   !< saturation water vapor pressure
2297     REAL(wp) ::  q_s   !< saturation mixing ratio
2298     REAL(wp) ::  q_wv  !< mixing ratio
2299
2300     CALL cpu_log( log_point_s(27), 'VM sampling', 'start' )
2301!
2302!--  Loop over all sites.
2303     DO  l = 1, vmea_general%nvm
2304!
2305!--     At the beginning, set _FillValues
2306        IF ( ALLOCATED( vmea(l)%measured_vars ) ) vmea(l)%measured_vars = vmea(l)%fillout
2307        IF ( ALLOCATED( vmea(l)%measured_vars_soil ) ) vmea(l)%measured_vars_soil = vmea(l)%fillout
2308!
2309!--     Loop over all variables measured at this site.
2310        DO  n = 1, vmea(l)%nmeas
2311
2312           SELECT CASE ( TRIM( vmea(l)%var_atts(n)%name ) )
2313
2314              CASE ( 'theta' ) ! potential temperature
2315                 IF ( .NOT. neutral )  THEN
2316                    DO  m = 1, vmea(l)%ns
2317                       k = vmea(l)%k(m)
2318                       j = vmea(l)%j(m)
2319                       i = vmea(l)%i(m)
2320                       vmea(l)%measured_vars(m,n) = pt(k,j,i)
2321                    ENDDO
2322                 ENDIF
2323
2324              CASE ( 'ta' ) ! absolute temperature
2325                 IF ( .NOT. neutral )  THEN
2326                    DO  m = 1, vmea(l)%ns
2327                       k = vmea(l)%k(m)
2328                       j = vmea(l)%j(m)
2329                       i = vmea(l)%i(m)
2330                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k ) - degc_to_k
2331                    ENDDO
2332                 ENDIF
2333
2334              CASE ( 't_va' )
2335
2336              CASE ( 'hus' ) ! mixing ratio
2337                 IF ( humidity )  THEN
2338                    DO  m = 1, vmea(l)%ns
2339                       k = vmea(l)%k(m)
2340                       j = vmea(l)%j(m)
2341                       i = vmea(l)%i(m)
2342                       vmea(l)%measured_vars(m,n) = q(k,j,i)
2343                    ENDDO
2344                 ENDIF
2345
2346              CASE ( 'haa' ) ! absolute humidity
2347                 IF ( humidity )  THEN
2348                    DO  m = 1, vmea(l)%ns
2349                       k = vmea(l)%k(m)
2350                       j = vmea(l)%j(m)
2351                       i = vmea(l)%i(m)
2352                       vmea(l)%measured_vars(m,n) = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) * rho_air(k)
2353                    ENDDO
2354                 ENDIF
2355
2356              CASE ( 'pwv' ) ! water vapor partial pressure
2357                 IF ( humidity )  THEN
2358!                     DO  m = 1, vmea(l)%ns
2359!                        k = vmea(l)%k(m)
2360!                        j = vmea(l)%j(m)
2361!                        i = vmea(l)%i(m)
2362!                        vmea(l)%measured_vars(m,n) = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) )          &
2363!                                                     * rho_air(k)
2364!                     ENDDO
2365                 ENDIF
2366
2367              CASE ( 'hur' ) ! relative humidity
2368                 IF ( humidity )  THEN
2369                    DO  m = 1, vmea(l)%ns
2370                       k = vmea(l)%k(m)
2371                       j = vmea(l)%j(m)
2372                       i = vmea(l)%i(m)
2373!
2374!--                    Calculate actual temperature, water vapor saturation pressure and, based on
2375!--                    this, the saturation mixing ratio.
2376                       e_s  = magnus( exner(k) * pt(k,j,i) )
2377                       q_s  = rd_d_rv * e_s / ( hyp(k) - e_s )
2378                       q_wv = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) * rho_air(k)
2379
2380                       vmea(l)%measured_vars(m,n) = q_wv / ( q_s + 1E-10_wp )
2381                    ENDDO
2382                 ENDIF
2383
2384              CASE ( 'u', 'ua' ) ! u-component
2385                 DO  m = 1, vmea(l)%ns
2386                    k = vmea(l)%k(m)
2387                    j = vmea(l)%j(m)
2388                    i = vmea(l)%i(m)
2389                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
2390                 ENDDO
2391
2392              CASE ( 'v', 'va' ) ! v-component
2393                 DO  m = 1, vmea(l)%ns
2394                    k = vmea(l)%k(m)
2395                    j = vmea(l)%j(m)
2396                    i = vmea(l)%i(m)
2397                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
2398                 ENDDO
2399
2400              CASE ( 'w' ) ! w-component
2401                 DO  m = 1, vmea(l)%ns
2402                    k = MAX ( 1, vmea(l)%k(m) )
2403                    j = vmea(l)%j(m)
2404                    i = vmea(l)%i(m)
2405                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
2406                 ENDDO
2407
2408              CASE ( 'wspeed' ) ! horizontal wind speed
2409                 DO  m = 1, vmea(l)%ns
2410                    k = vmea(l)%k(m)
2411                    j = vmea(l)%j(m)
2412                    i = vmea(l)%i(m)
2413                    vmea(l)%measured_vars(m,n) = SQRT(   ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 &
2414                                                       + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 &
2415                                                     )
2416                 ENDDO
2417
2418              CASE ( 'wdir' ) ! wind direction
2419                 DO  m = 1, vmea(l)%ns
2420                    k = vmea(l)%k(m)
2421                    j = vmea(l)%j(m)
2422                    i = vmea(l)%i(m)
2423
2424                    vmea(l)%measured_vars(m,n) = 180.0_wp + 180.0_wp / pi * ATAN2(                 &
2425                                                               0.5_wp * ( v(k,j,i) + v(k,j+1,i) ), &
2426                                                               0.5_wp * ( u(k,j,i) + u(k,j,i+1) )  &
2427                                                                                 )
2428                 ENDDO
2429
2430              CASE ( 'utheta' )
2431                 DO  m = 1, vmea(l)%ns
2432                    k = vmea(l)%k(m)
2433                    j = vmea(l)%j(m)
2434                    i = vmea(l)%i(m)
2435                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) * pt(k,j,i)
2436                 ENDDO
2437
2438              CASE ( 'vtheta' )
2439                 DO  m = 1, vmea(l)%ns
2440                    k = vmea(l)%k(m)
2441                    j = vmea(l)%j(m)
2442                    i = vmea(l)%i(m)
2443                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) * pt(k,j,i)
2444                 ENDDO
2445
2446              CASE ( 'wtheta' )
2447                 DO  m = 1, vmea(l)%ns
2448                    k = MAX ( 1, vmea(l)%k(m) )
2449                    j = vmea(l)%j(m)
2450                    i = vmea(l)%i(m)
2451                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k-1,j,i) + w(k,j,i) ) * pt(k,j,i)
2452                 ENDDO
2453
2454              CASE ( 'uqv' )
2455                 IF ( humidity )  THEN
2456                    DO  m = 1, vmea(l)%ns
2457                       k = vmea(l)%k(m)
2458                       j = vmea(l)%j(m)
2459                       i = vmea(l)%i(m)
2460                       vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) * q(k,j,i)
2461                    ENDDO
2462                 ENDIF
2463
2464              CASE ( 'vqv' )
2465                 IF ( humidity )  THEN
2466                    DO  m = 1, vmea(l)%ns
2467                       k = vmea(l)%k(m)
2468                       j = vmea(l)%j(m)
2469                       i = vmea(l)%i(m)
2470                       vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) * q(k,j,i)
2471                    ENDDO
2472                 ENDIF
2473
2474              CASE ( 'wqv' )
2475                 IF ( humidity )  THEN
2476                    DO  m = 1, vmea(l)%ns
2477                       k = MAX ( 1, vmea(l)%k(m) )
2478                       j = vmea(l)%j(m)
2479                       i = vmea(l)%i(m)
2480                       vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k-1,j,i) + w(k,j,i) ) * q(k,j,i)
2481                    ENDDO
2482                 ENDIF
2483
2484              CASE ( 'uw' )
2485                 DO  m = 1, vmea(l)%ns
2486                    k = MAX ( 1, vmea(l)%k(m) )
2487                    j = vmea(l)%j(m)
2488                    i = vmea(l)%i(m)
2489                    vmea(l)%measured_vars(m,n) = 0.25_wp * ( w(k-1,j,i) + w(k,j,i) ) *             &
2490                                                           ( u(k,j,i)   + u(k,j,i+1) )
2491                 ENDDO
2492
2493              CASE ( 'vw' )
2494                 DO  m = 1, vmea(l)%ns
2495                    k = MAX ( 1, vmea(l)%k(m) )
2496                    j = vmea(l)%j(m)
2497                    i = vmea(l)%i(m)
2498                    vmea(l)%measured_vars(m,n) = 0.25_wp * ( w(k-1,j,i) + w(k,j,i) ) *             &
2499                                                           ( v(k,j,i)   + v(k,j+1,i) )
2500                 ENDDO
2501
2502              CASE ( 'uv' )
2503                 DO  m = 1, vmea(l)%ns
2504                    k = vmea(l)%k(m)
2505                    j = vmea(l)%j(m)
2506                    i = vmea(l)%i(m)
2507                    vmea(l)%measured_vars(m,n) = 0.25_wp * ( u(k,j,i)   + u(k,j,i+1) ) *           &
2508                                                           ( v(k,j,i)   + v(k,j+1,i) )
2509                 ENDDO
2510!
2511!--           Chemistry variables. List of variables that may need extension. Note, gas species in
2512!--           PALM are in ppm and no distinction is made between mole-fraction and concentration
2513!--           quantities (all are output in ppm so far).
2514              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfno', 'mfno2', 'mcno', 'mcno2', 'tro3' )
2515                 IF ( air_chemistry )  THEN
2516!
2517!--                 First, search for the measured variable in the chem_vars
2518!--                 list, in order to get the internal name of the variable.
2519                    DO  nn = 1, UBOUND( chem_vars, 2 )
2520                       IF ( TRIM( vmea(l)%var_atts(n)%name ) ==                                    &
2521                            TRIM( chem_vars(0,nn) ) )  ind_chem = nn
2522                    ENDDO
2523!
2524!--                 Run loop over all chemical species, if the measured variable matches the interal
2525!--                 name, sample the variable. Note, nvar as a chemistry-module variable.
2526                    DO  nn = 1, nvar
2527                       IF ( TRIM( chem_vars(1,ind_chem) ) == TRIM( chem_species(nn)%name ) )  THEN
2528                          DO  m = 1, vmea(l)%ns
2529                             k = vmea(l)%k(m)
2530                             j = vmea(l)%j(m)
2531                             i = vmea(l)%i(m)
2532                             vmea(l)%measured_vars(m,n) = chem_species(nn)%conc(k,j,i)
2533                          ENDDO
2534                       ENDIF
2535                    ENDDO
2536                 ENDIF
2537
2538              CASE ( 'us' ) ! friction velocity
2539                 DO  m = 1, vmea(l)%ns
2540!
2541!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
2542!--                 limit the indices.
2543                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2544                    j = MERGE( j           , nyn, j            < nyn )
2545                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2546                    i = MERGE( i           , nxr, i            < nxr )
2547
2548                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
2549                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
2550                    ENDDO
2551                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2552                       vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm)
2553                    ENDDO
2554                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2555                       vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm)
2556                    ENDDO
2557                 ENDDO
2558
2559              CASE ( 'thetas' ) ! scaling parameter temperature
2560                 DO  m = 1, vmea(l)%ns
2561!
2562!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
2563!-                  limit the indices.
2564                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2565                    j = MERGE( j           , nyn, j            < nyn )
2566                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2567                    i = MERGE( i           , nxr, i            < nxr )
2568
2569                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
2570                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
2571                    ENDDO
2572                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2573                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm)
2574                    ENDDO
2575                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2576                       vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm)
2577                    ENDDO
2578                 ENDDO
2579
2580              CASE ( 'hfls' ) ! surface latent heat flux
2581                 DO  m = 1, vmea(l)%ns
2582!
2583!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
2584!--                 limit the indices.
2585                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2586                    j = MERGE( j           , nyn, j            < nyn )
2587                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2588                    i = MERGE( i           , nxr, i            < nxr )
2589
2590                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
2591                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
2592                    ENDDO
2593                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2594                       vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm)
2595                    ENDDO
2596                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2597                       vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm)
2598                    ENDDO
2599                 ENDDO
2600
2601              CASE ( 'hfss' ) ! surface sensible heat flux
2602                 DO  m = 1, vmea(l)%ns
2603!
2604!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
2605!--                 limit the indices.
2606                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2607                    j = MERGE( j           , nyn, j            < nyn )
2608                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2609                    i = MERGE( i           , nxr, i            < nxr )
2610
2611                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
2612                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
2613                    ENDDO
2614                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2615                       vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm)
2616                    ENDDO
2617                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2618                       vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm)
2619                    ENDDO
2620                 ENDDO
2621
2622              CASE ( 'hfdg' ) ! ground heat flux
2623                 DO  m = 1, vmea(l)%ns
2624!
2625!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
2626!--                 limit the indices.
2627                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2628                    j = MERGE( j           , nyn, j            < nyn )
2629                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2630                    i = MERGE( i           , nxr, i            < nxr )
2631
2632                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2633                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ghf(mm)
2634                    ENDDO
2635                 ENDDO
2636
2637              CASE ( 'lwcs' )  ! liquid water of soil layer
2638!                  DO  m = 1, vmea(l)%ns
2639! !
2640! !--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
2641! !--                 limit the indices.
2642!                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2643!                     j = MERGE( j           , nyn, j            < nyn )
2644!                     i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2645!                     i = MERGE( i           , nxr, i            < nxr )
2646!
2647!                     DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2648!                        vmea(l)%measured_vars(m,n) = ?
2649!                     ENDDO
2650!                  ENDDO
2651
2652              CASE ( 'rnds' ) ! surface net radiation
2653                 IF ( radiation )  THEN
2654                    DO  m = 1, vmea(l)%ns
2655!
2656!--                    Surface data is only available on inner subdomains, not on ghost points.
2657!--                    Hence, limit the indices.
2658                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2659                       j = MERGE( j           , nyn, j            < nyn )
2660                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2661                       i = MERGE( i           , nxr, i            < nxr )
2662
2663                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2664                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm)
2665                       ENDDO
2666                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2667                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm)
2668                       ENDDO
2669                    ENDDO
2670                 ENDIF
2671
2672              CASE ( 'rsus' ) ! surface shortwave out
2673                 IF ( radiation )  THEN
2674                    DO  m = 1, vmea(l)%ns
2675!
2676!--                    Surface data is only available on inner subdomains, not on ghost points.
2677!--                    Hence, limit the indices.
2678                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2679                       j = MERGE( j           , nyn, j            < nyn )
2680                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2681                       i = MERGE( i           , nxr, i            < nxr )
2682
2683                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2684                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm)
2685                       ENDDO
2686                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2687                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm)
2688                       ENDDO
2689                    ENDDO
2690                 ENDIF
2691
2692              CASE ( 'rsds' ) ! surface shortwave in
2693                 IF ( radiation )  THEN
2694                    DO  m = 1, vmea(l)%ns
2695!
2696!--                    Surface data is only available on inner subdomains, not on ghost points.
2697!--                    Hence, limit the indices.
2698                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2699                       j = MERGE( j           , nyn, j            < nyn )
2700                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2701                       i = MERGE( i           , nxr, i            < nxr )
2702
2703                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2704                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm)
2705                       ENDDO
2706                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2707                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm)
2708                       ENDDO
2709                    ENDDO
2710                 ENDIF
2711
2712              CASE ( 'rlus' ) ! surface longwave out
2713                 IF ( radiation )  THEN
2714                    DO  m = 1, vmea(l)%ns
2715!
2716!--                    Surface data is only available on inner subdomains, not on ghost points.
2717!--                    Hence, limit the indices.
2718                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2719                       j = MERGE( j           , nyn, j            < nyn )
2720                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2721                       i = MERGE( i           , nxr, i            < nxr )
2722
2723                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2724                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm)
2725                       ENDDO
2726                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2727                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm)
2728                       ENDDO
2729                    ENDDO
2730                 ENDIF
2731
2732              CASE ( 'rlds' ) ! surface longwave in
2733                 IF ( radiation )  THEN
2734                    DO  m = 1, vmea(l)%ns
2735!
2736!--                    Surface data is only available on inner subdomains, not on ghost points.
2737!--                    Hence, limit the indices.
2738                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2739                       j = MERGE( j           , nyn, j            < nyn )
2740                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2741                       i = MERGE( i           , nxr, i            < nxr )
2742
2743                       DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2744                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm)
2745                       ENDDO
2746                       DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2747                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm)
2748                       ENDDO
2749                    ENDDO
2750                 ENDIF
2751
2752              CASE ( 'rsd' ) ! shortwave in
2753                 IF ( radiation )  THEN
2754                    IF ( radiation_scheme /= 'rrtmg' )  THEN
2755                       DO  m = 1, vmea(l)%ns
2756                          k = 0
2757                          j = vmea(l)%j(m)
2758                          i = vmea(l)%i(m)
2759                          vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
2760                       ENDDO
2761                    ELSE
2762                       DO  m = 1, vmea(l)%ns
2763                          k = vmea(l)%k(m)
2764                          j = vmea(l)%j(m)
2765                          i = vmea(l)%i(m)
2766                          vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
2767                       ENDDO
2768                    ENDIF
2769                 ENDIF
2770
2771              CASE ( 'rsu' ) ! shortwave out
2772                 IF ( radiation )  THEN
2773                    IF ( radiation_scheme /= 'rrtmg' )  THEN
2774                       DO  m = 1, vmea(l)%ns
2775                          k = 0
2776                          j = vmea(l)%j(m)
2777                          i = vmea(l)%i(m)
2778                          vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
2779                       ENDDO
2780                    ELSE
2781                       DO  m = 1, vmea(l)%ns
2782                          k = vmea(l)%k(m)
2783                          j = vmea(l)%j(m)
2784                          i = vmea(l)%i(m)
2785                          vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
2786                       ENDDO
2787                    ENDIF
2788                 ENDIF
2789
2790              CASE ( 'rlu' ) ! longwave out
2791                 IF ( radiation )  THEN
2792                    IF ( radiation_scheme /= 'rrtmg' )  THEN
2793                       DO  m = 1, vmea(l)%ns
2794                          k = 0
2795                          j = vmea(l)%j(m)
2796                          i = vmea(l)%i(m)
2797                          vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
2798                       ENDDO
2799                    ELSE
2800                       DO  m = 1, vmea(l)%ns
2801                          k = vmea(l)%k(m)
2802                          j = vmea(l)%j(m)
2803                          i = vmea(l)%i(m)
2804                          vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
2805                       ENDDO
2806                    ENDIF
2807                 ENDIF
2808
2809              CASE ( 'rld' ) ! longwave in
2810                 IF ( radiation )  THEN
2811                    IF ( radiation_scheme /= 'rrtmg' )  THEN
2812                       DO  m = 1, vmea(l)%ns
2813                          k = 0
2814!
2815!--                       Surface data is only available on inner subdomains, not on ghost points.
2816!--                       Hence, limit the indices.
2817                          j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2818                          j = MERGE( j           , nyn, j            < nyn )
2819                          i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2820                          i = MERGE( i           , nxr, i            < nxr )
2821
2822                          vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
2823                       ENDDO
2824                    ELSE
2825                       DO  m = 1, vmea(l)%ns
2826                          k = vmea(l)%k(m)
2827                          j = vmea(l)%j(m)
2828                          i = vmea(l)%i(m)
2829                          vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
2830                       ENDDO
2831                    ENDIF
2832                 ENDIF
2833
2834              CASE ( 'rsddif' ) ! shortwave in, diffuse part
2835                 IF ( radiation )  THEN
2836                    DO  m = 1, vmea(l)%ns
2837                       j = vmea(l)%j(m)
2838                       i = vmea(l)%i(m)
2839
2840                       vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
2841                    ENDDO
2842                 ENDIF
2843
2844              CASE ( 't_soil' ) ! soil and wall temperature
2845                 DO  m = 1, vmea(l)%ns_soil
2846                    j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
2847                    j = MERGE( j                , nyn, j                 < nyn )
2848                    i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
2849                    i = MERGE( i                , nxr, i                 < nxr )
2850                    k = vmea(l)%k_soil(m)
2851
2852                    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
2853                    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
2854
2855                    IF ( match_lsm )  THEN
2856                       mm = surf_lsm_h%start_index(j,i)
2857                       vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm)
2858                    ENDIF
2859
2860                    IF ( match_usm )  THEN
2861                       mm = surf_usm_h%start_index(j,i)
2862                       vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm)
2863                    ENDIF
2864                 ENDDO
2865
2866              CASE ( 'm_soil' ) ! soil moisture
2867                 DO  m = 1, vmea(l)%ns_soil
2868                    j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
2869                    j = MERGE( j                , nyn, j                 < nyn )
2870                    i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
2871                    i = MERGE( i                , nxr, i                 < nxr )
2872                    k = vmea(l)%k_soil(m)
2873
2874                    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
2875
2876                    IF ( match_lsm )  THEN
2877                       mm = surf_lsm_h%start_index(j,i)
2878                       vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm)
2879                    ENDIF
2880
2881                 ENDDO
2882
2883              CASE ( 'ts' ) ! surface temperature
2884                 DO  m = 1, vmea(l)%ns
2885!
2886!--                 Surface data is only available on inner subdomains, not on ghost points. Hence,
2887!--                 limit the indices.
2888                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2889                    j = MERGE( j           , nyn, j            < nyn )
2890                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2891                    i = MERGE( i           , nxr, i            < nxr )
2892
2893                    DO  mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
2894                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%pt_surface(mm)
2895                    ENDDO
2896                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2897                       vmea(l)%measured_vars(m,n) = surf_lsm_h%pt_surface(mm)
2898                    ENDDO
2899                    DO  mm = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
2900                       vmea(l)%measured_vars(m,n) = surf_usm_h%pt_surface(mm)
2901                    ENDDO
2902                 ENDDO
2903
2904              CASE ( 'lwp' ) ! liquid water path
2905                 IF ( ASSOCIATED( ql ) )  THEN
2906                    DO  m = 1, vmea(l)%ns
2907                       j = vmea(l)%j(m)
2908                       i = vmea(l)%i(m)
2909
2910                       vmea(l)%measured_vars(m,n) = SUM( ql(nzb:nzt,j,i) * dzw(1:nzt+1) )          &
2911                                                    * rho_surface
2912                    ENDDO
2913                 ENDIF
2914
2915              CASE ( 'ps' ) ! surface pressure
2916                 vmea(l)%measured_vars(:,n) = surface_pressure
2917
2918              CASE ( 'pswrtg' ) ! platform speed above ground
2919                 vmea(l)%measured_vars(:,n) = 0.0_wp
2920
2921              CASE ( 'pswrta' ) ! platform speed in air
2922                 vmea(l)%measured_vars(:,n) = 0.0_wp
2923
2924              CASE ( 't_lw' ) ! water temperature
2925                 DO  m = 1, vmea(l)%ns
2926!
2927!--                 Surface data is only available on inner subdomains, not
2928!--                 on ghost points. Hence, limit the indices.
2929                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2930                    j = MERGE( j           , nyn, j            < nyn )
2931                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2932                    i = MERGE( i           , nxr, i            < nxr )
2933
2934                    DO  mm = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
2935                       IF ( surf_lsm_h%water_surface(m) )                                          &
2936                            vmea(l)%measured_vars(m,n) = t_soil_h%var_2d(nzt,m)
2937                    ENDDO
2938
2939                 ENDDO
2940!
2941!--           More will follow ...
2942              CASE ( 'ncaa' )
2943!
2944!--           No match found - just set a fill value
2945              CASE DEFAULT
2946                 vmea(l)%measured_vars(:,n) = vmea(l)%fillout
2947           END SELECT
2948
2949        ENDDO
2950
2951     ENDDO
2952
2953     CALL cpu_log( log_point_s(27), 'VM sampling', 'stop' )
2954
2955 END SUBROUTINE vm_sampling
2956
2957
2958 END MODULE virtual_measurement_mod
Note: See TracBrowser for help on using the repository browser.