Home


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

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

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

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