Home


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

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

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

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