Home


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

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

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

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