Home


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

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

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

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