Home


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

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

Last change on this file since 4592 was 4536, checked in by raasch, 4 years ago

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