Home


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

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

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

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