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

Last change on this file since 4400 was 4400, checked in by suehring, 20 months ago

Revision of the virtual-measurement module: data input from NetCDF file; removed binary output - instead parallel NetCDF output using the new data-output module; variable attributes added; further variables added that can be sampled, file connections added; Functions for coordinate transformation moved to basic_constants_and_equations; netcdf_data_input_mod: unused routines netcdf_data_input_att and netcdf_data_input_var removed; new routines to inquire fill values added; Preprocessing script (palm_cvd) to setup virtual-measurement input files provided; postprocessor combine_virtual_measurements deleted

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