Home


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

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

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

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