Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

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

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