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

Last change on this file since 4406 was 4406, checked in by knoop, 20 months ago

Bugix: removed oro_rel wrong loop bounds and removed unnecessary restart method

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