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

Last change on this file since 4421 was 4421, checked in by suehring, 3 years ago

Output of character station_name (required for UC2 standard); bugfix - missing coupling_char for opening NetCDF file

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