Home


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

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

Last change on this file since 4535 was 4504, checked in by raasch, 5 years ago

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