Home


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

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

Last change on this file since 4694 was 4671, checked in by pavelkrc, 4 years ago

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