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

Last change on this file since 4408 was 4408, checked in by gronemeier, 3 years ago

write fill_value attribute in virtual-measurements module; enable character-array output in data-output module

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 164.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 4408 2020-02-14 10:04:39Z gronemeier $
27! write fill_value attribute
28!
29! 4406 2020-02-13 20:06:29Z knoop
30! Bugix: removed oro_rel wrong loop bounds and removed unnecessary restart method
31!
32! 4400 2020-02-10 20:32:41Z suehring
33! Revision of the module:
34! - revised input from NetCDF setup file
35! - parallel NetCDF output via data-output module ( Tobias Gronemeier )
36! - variable attributes added
37! - further variables defined
38!
39! 4346 2019-12-18 11:55:56Z motisi
40! Introduction of wall_flags_total_0, which currently sets bits based on static
41! topography information used in wall_flags_static_0
42!
43! 4329 2019-12-10 15:46:36Z motisi
44! Renamed wall_flags_0 to wall_flags_static_0
45!
46! 4226 2019-09-10 17:03:24Z suehring
47! Netcdf input routine for dimension length renamed
48!
49! 4182 2019-08-22 15:20:23Z scharf
50! Corrected "Former revisions" section
51!
52! 4168 2019-08-16 13:50:17Z suehring
53! Replace function get_topography_top_index by topo_top_ind
54!
55! 3988 2019-05-22 11:32:37Z kanani
56! Add variables to enable steering of output interval for virtual measurements
57!
58! 3913 2019-04-17 15:12:28Z gronemeier
59! Bugfix: rotate positions of measurements before writing them into file
60!
61! 3910 2019-04-17 11:46:56Z suehring
62! Bugfix in rotation of UTM coordinates
63!
64! 3904 2019-04-16 18:22:51Z gronemeier
65! Rotate coordinates of stations by given rotation_angle
66!
67! 3876 2019-04-08 18:41:49Z knoop
68! Remove print statement
69!
70! 3854 2019-04-02 16:59:33Z suehring
71! renamed nvar to nmeas, replaced USE chem_modules by USE chem_gasphase_mod and
72! nspec by nvar
73!
74! 3766 2019-02-26 16:23:41Z raasch
75! unused variables removed
76!
77! 3718 2019-02-06 11:08:28Z suehring
78! Adjust variable name connections between UC2 and chemistry variables
79!
80! 3717 2019-02-05 17:21:16Z suehring
81! Additional check + error numbers adjusted
82!
83! 3706 2019-01-29 20:02:26Z suehring
84! unused variables removed
85!
86! 3705 2019-01-29 19:56:39Z suehring
87! - initialization revised
88! - binary data output
89! - list of allowed variables extended
90!
91! 3704 2019-01-29 19:51:41Z suehring
92! Sampling of variables
93!
94! 3473 2018-10-30 20:50:15Z suehring
95! Initial revision
96!
97! Authors:
98! --------
99! @author Matthias Suehring
100! @author Tobias Gronemeier
101!
102! Description:
103! ------------
104!> The module acts as an interface between 'real-world' observations and
105!> model simulations. Virtual measurements will be taken in the model at the
106!> coordinates representative for the 'real-world' observation coordinates.
107!> More precisely, coordinates and measured quanties will be read from a
108!> NetCDF file which contains all required information. In the model,
109!> the same quantities (as long as all the required components are switched-on)
110!> will be sampled at the respective positions and output into an extra file,
111!> which allows for straight-forward comparison of model results with
112!> observations.
113!>
114!> @todo Check why there is an error when _FillValue attributes are added via
115!>       dom.
116!> @todo Output of character variable station_name (dom hasn't this feature
117!>       yet implemented).
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, 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, 32_iwp/),                 &
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. DOM needs to be enabled to define CHARACTER variables.
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
2263    INTEGER(iwp)       ::  l             !< loop index
2264    INTEGER(iwp)       ::  n             !< loop index
2265    INTEGER            ::  return_value  !< returned status value of called function
2266
2267    INTEGER(iwp)       ::  t_ind         !< time index
2268
2269    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  oro_rel                  !< relative altitude of model surface
2270    REAL(wp), DIMENSION(:), POINTER               ::  output_values_1d_pointer !< pointer for 1d output array
2271    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET   ::  output_values_1d_target  !< target for 1d output array
2272    REAL(wp), DIMENSION(:,:), POINTER             ::  output_values_2d_pointer !< pointer for 2d output array
2273    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_target  !< target for 2d output array
2274
2275!
2276!-- At the first call of this routine write the spatial coordinates.
2277    IF ( .NOT. initial_write_coordinates )  THEN
2278!
2279!--    Write spatial coordinates.
2280       DO  l = 1, vmea_general%nvm
2281!
2282!--       Skip if no observations were taken.
2283          IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
2284
2285          ALLOCATE( output_values_1d_target(vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
2286!
2287!--       Output of Easting coordinate. Before output, recalculate EUTM.
2288          output_values_1d_target = init_model%origin_x                        &
2289                    + REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
2290                    * COS( init_model%rotation_angle * pi / 180.0_wp )         &
2291                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
2292                    * SIN( init_model%rotation_angle * pi / 180.0_wp )
2293
2294          output_values_1d_pointer => output_values_1d_target
2295
2296          return_value =                                                       &
2297                  dom_write_var( vmea(l)%nc_filename,                          &
2298                                 'E_UTM',                                      &
2299                                 values_realwp_1d = output_values_1d_pointer,  &
2300                                 bounds_start = (/vmea(l)%start_coord_a/),     &
2301                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
2302!
2303!--       Output of Northing coordinate. Before output, recalculate NUTM.
2304          output_values_1d_target = init_model%origin_y                        &
2305                    - REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
2306                    * SIN( init_model%rotation_angle * pi / 180.0_wp )         &
2307                    + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
2308                    * COS( init_model%rotation_angle * pi / 180.0_wp )
2309
2310          output_values_1d_pointer => output_values_1d_target
2311          return_value =                                                       &
2312                  dom_write_var( vmea(l)%nc_filename,                          &
2313                                 'N_UTM',                                      &
2314                                 values_realwp_1d = output_values_1d_pointer,  &
2315                                 bounds_start = (/vmea(l)%start_coord_a/),     &
2316                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
2317!
2318!--       Output of relative height coordinate.
2319!--       Before this is output, first define the relative orographie height
2320!--       and add this to z.
2321          ALLOCATE( oro_rel(1:vmea(l)%ns) )
2322          DO  n = 1, vmea(l)%ns
2323             oro_rel(n) = zw(topo_top_ind(vmea(l)%j(n),vmea(l)%i(n),3))
2324          ENDDO
2325
2326          output_values_1d_target = vmea(l)%zar(1:vmea(l)%ns) + oro_rel(:)
2327          output_values_1d_pointer => output_values_1d_target
2328          return_value =                                                       &
2329                  dom_write_var( vmea(l)%nc_filename,                          &
2330                                 'z',                                          &
2331                                 values_realwp_1d = output_values_1d_pointer,  &
2332                                 bounds_start = (/vmea(l)%start_coord_a/),     &
2333                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
2334!
2335!--       Write surface altitude for the station. Note, since z is already
2336!--       a relative observation height, station_h must be zero, in order
2337!--       to obtain the observation level.
2338          output_values_1d_target = oro_rel(:)
2339          output_values_1d_pointer => output_values_1d_target
2340          return_value =                                                       &
2341                  dom_write_var( vmea(l)%nc_filename,                          &
2342                                 'station_h',                                  &
2343                                 values_realwp_1d = output_values_1d_pointer,  &
2344                                 bounds_start = (/vmea(l)%start_coord_a/),     &
2345                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
2346
2347          DEALLOCATE( oro_rel )
2348          DEALLOCATE( output_values_1d_target )
2349!
2350!--       In case of sampled soil quantities, output also the respective
2351!--       coordinate arrays.
2352          IF ( vmea(l)%soil_sampling )  THEN
2353             ALLOCATE( output_values_1d_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
2354!
2355!--          Output of Easting coordinate. Before output, recalculate EUTM.
2356             output_values_1d_target = init_model%origin_x                     &
2357               + REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
2358               * COS( init_model%rotation_angle * pi / 180.0_wp )              &
2359               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
2360               * SIN( init_model%rotation_angle * pi / 180.0_wp )
2361             output_values_1d_pointer => output_values_1d_target
2362             return_value =                                                    &
2363                  dom_write_var( vmea(l)%nc_filename,                          &
2364                                 'E_UTM_soil',                                 &
2365                                 values_realwp_1d = output_values_1d_pointer,  &
2366                                 bounds_start = (/vmea(l)%start_coord_s/),     &
2367                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
2368!
2369!--          Output of Northing coordinate. Before output, recalculate NUTM.
2370             output_values_1d_target = init_model%origin_y                     &
2371               - REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
2372               * SIN( init_model%rotation_angle * pi / 180.0_wp )              &
2373               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
2374               * COS( init_model%rotation_angle * pi / 180.0_wp )
2375
2376             output_values_1d_pointer => output_values_1d_target
2377             return_value =                                                    &
2378                  dom_write_var( vmea(l)%nc_filename,                          &
2379                                 'N_UTM_soil',                                 &
2380                                 values_realwp_1d = output_values_1d_pointer,  &
2381                                 bounds_start = (/vmea(l)%start_coord_s/),     &
2382                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
2383!
2384!--          Output of relative height coordinate.
2385!--          Before this is output, first define the relative orographie height
2386!--          and add this to z.
2387             ALLOCATE( oro_rel(1:vmea(l)%ns_soil) )
2388             DO  n = 1, vmea(l)%ns_soil
2389                oro_rel(n) = zw(topo_top_ind(vmea(l)%j_soil(n),vmea(l)%i_soil(n),3))
2390             ENDDO
2391
2392             output_values_1d_target = vmea(l)%depth(1:vmea(l)%ns_soil) + oro_rel(:)
2393             output_values_1d_pointer => output_values_1d_target
2394             return_value =                                                    &
2395                  dom_write_var( vmea(l)%nc_filename,                          &
2396                                 'z_soil',                                     &
2397                                 values_realwp_1d = output_values_1d_pointer,  &
2398                                 bounds_start = (/vmea(l)%start_coord_s/),     &
2399                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
2400!
2401!--          Write surface altitude for the station. Note, since z is already
2402!--          a relative observation height, station_h must be zero, in order
2403!--          to obtain the observation level.
2404             output_values_1d_target = oro_rel(:)
2405             output_values_1d_pointer => output_values_1d_target
2406             return_value =                                                    &
2407                  dom_write_var( vmea(l)%nc_filename,                          &
2408                                 'station_h_soil',                             &
2409                                 values_realwp_1d = output_values_1d_pointer,  &
2410                                 bounds_start = (/vmea(l)%start_coord_s/),     &
2411                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
2412
2413             DEALLOCATE( oro_rel )
2414             DEALLOCATE( output_values_1d_target )
2415!
2416!--          Write the stations name
2417
2418          ENDIF
2419
2420       ENDDO  ! loop over sites
2421
2422       initial_write_coordinates = .TRUE.
2423    ENDIF
2424!
2425!-- Loop over all sites.
2426    DO  l = 1, vmea_general%nvm
2427!
2428!--    Skip if no observations were taken.
2429       IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
2430!
2431!--    Determine time index in file.
2432       t_ind = vmea(l)%file_time_index + 1
2433!
2434!--    Write output variables. Distinguish between atmosphere and soil variables.
2435       DO  n = 1, vmea(l)%nmeas
2436          IF ( vmea(l)%soil_sampling  .AND.                                    &
2437            ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
2438!
2439!--          Write time coordinate to file
2440             variable_name = 'time_soil'
2441             ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
2442             output_values_2d_target(t_ind,:) = time_since_reference_point
2443             output_values_2d_pointer => output_values_2d_target
2444
2445             return_value = dom_write_var( vmea(l)%nc_filename,                &
2446                                           variable_name,                      &
2447                                           values_realwp_2d = output_values_2d_pointer, &
2448                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
2449                                           bounds_end   = (/vmea(l)%end_coord_s, t_ind /) )
2450
2451             variable_name = TRIM( vmea(l)%var_atts(n)%name )
2452             output_values_2d_target(t_ind,:) = vmea(l)%measured_vars_soil(:,n)
2453             output_values_2d_pointer => output_values_2d_target
2454             return_value =                                                    &
2455                      dom_write_var( vmea(l)%nc_filename,                      &
2456                                     variable_name,                            &
2457                                     values_realwp_2d = output_values_2d_pointer, &
2458                                     bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
2459                                     bounds_end   = (/vmea(l)%end_coord_s, t_ind  /) )
2460             DEALLOCATE( output_values_2d_target )
2461          ELSE
2462!
2463!--          Write time coordinate to file
2464             variable_name = 'time'
2465             ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
2466             output_values_2d_target(t_ind,:) = time_since_reference_point
2467             output_values_2d_pointer => output_values_2d_target
2468
2469             return_value = dom_write_var( vmea(l)%nc_filename,                &
2470                                           variable_name,                      &
2471                                           values_realwp_2d = output_values_2d_pointer, &
2472                                           bounds_start = (/vmea(l)%start_coord_a, t_ind/), &
2473                                           bounds_end   = (/vmea(l)%end_coord_a, t_ind/) )
2474
2475             variable_name = TRIM( vmea(l)%var_atts(n)%name )
2476
2477             output_values_2d_target(t_ind,:) = vmea(l)%measured_vars(:,n)
2478             output_values_2d_pointer => output_values_2d_target
2479             return_value =                                                    &
2480                      dom_write_var( vmea(l)%nc_filename,                      &
2481                                     variable_name,                            &
2482                                     values_realwp_2d = output_values_2d_pointer, &
2483                                     bounds_start = (/ vmea(l)%start_coord_a, t_ind /), &
2484                                     bounds_end   = (/ vmea(l)%end_coord_a, t_ind /) )
2485
2486             DEALLOCATE( output_values_2d_target )
2487          ENDIF
2488       ENDDO
2489!
2490!--    Update number of written time indices
2491       vmea(l)%file_time_index = t_ind
2492
2493    ENDDO  ! loop over sites
2494
2495
2496  END SUBROUTINE vm_data_output
2497
2498!------------------------------------------------------------------------------!
2499! Description:
2500! ------------
2501!> Sampling of the actual quantities along the observation coordinates
2502!------------------------------------------------------------------------------!
2503  SUBROUTINE vm_sampling
2504
2505    USE radiation_model_mod,                                                   &
2506        ONLY:  radiation
2507
2508    USE surface_mod,                                                           &
2509        ONLY:  surf_def_h,                                                     &
2510               surf_lsm_h,                                                     &
2511               surf_usm_h
2512
2513     INTEGER(iwp) ::  i         !< grid index in x-direction
2514     INTEGER(iwp) ::  j         !< grid index in y-direction
2515     INTEGER(iwp) ::  k         !< grid index in z-direction
2516     INTEGER(iwp) ::  ind_chem  !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming
2517     INTEGER(iwp) ::  l         !< running index over the number of stations
2518     INTEGER(iwp) ::  m         !< running index over all virtual observation coordinates
2519     INTEGER(iwp) ::  mm        !< index of surface element which corresponds to the virtual observation coordinate
2520     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
2521     INTEGER(iwp) ::  nn        !< running index over the number of chemcal species
2522
2523     LOGICAL ::  match_lsm !< flag indicating natural-type surface
2524     LOGICAL ::  match_usm !< flag indicating urban-type surface
2525
2526     REAL(wp) ::  e_s      !< saturation water vapor pressure
2527     REAL(wp) ::  q_s      !< saturation mixing ratio
2528     REAL(wp) ::  q_wv     !< mixing ratio
2529!
2530!--  Loop over all sites.
2531     DO  l = 1, vmea_general%nvm
2532!
2533!--     At the beginning, set _FillValues
2534        IF ( ALLOCATED( vmea(l)%measured_vars      ) )                         &
2535           vmea(l)%measured_vars      = vmea(l)%fillout
2536        IF ( ALLOCATED( vmea(l)%measured_vars_soil ) )                         &
2537           vmea(l)%measured_vars_soil = vmea(l)%fillout
2538!
2539!--     Loop over all variables measured at this site.
2540        DO  n = 1, vmea(l)%nmeas
2541
2542           SELECT CASE ( TRIM( vmea(l)%var_atts(n)%name ) )
2543
2544              CASE ( 'theta' ) ! potential temperature
2545                 IF ( .NOT. neutral )  THEN
2546                    DO  m = 1, vmea(l)%ns
2547                       k = vmea(l)%k(m)
2548                       j = vmea(l)%j(m)
2549                       i = vmea(l)%i(m)
2550                       vmea(l)%measured_vars(m,n) = pt(k,j,i)
2551                    ENDDO
2552                 ENDIF
2553
2554              CASE ( 'ta' ) ! absolute temperature
2555                 IF ( .NOT. neutral )  THEN
2556                    DO  m = 1, vmea(l)%ns
2557                       k = vmea(l)%k(m)
2558                       j = vmea(l)%j(m)
2559                       i = vmea(l)%i(m)
2560                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )     &
2561                                                  - degc_to_k
2562                    ENDDO
2563                 ENDIF
2564
2565              CASE ( 't_va' )
2566
2567              CASE ( 'hus' ) ! mixing ratio
2568                 IF ( humidity )  THEN
2569                    DO  m = 1, vmea(l)%ns
2570                       k = vmea(l)%k(m)
2571                       j = vmea(l)%j(m)
2572                       i = vmea(l)%i(m)
2573                       vmea(l)%measured_vars(m,n) = q(k,j,i)
2574                    ENDDO
2575                 ENDIF
2576
2577              CASE ( 'haa' ) ! absolute humidity
2578                 IF ( humidity )  THEN
2579                    DO  m = 1, vmea(l)%ns
2580                       k = vmea(l)%k(m)
2581                       j = vmea(l)%j(m)
2582                       i = vmea(l)%i(m)
2583                       vmea(l)%measured_vars(m,n) = ( q(k,j,i)                 &
2584                                                    / ( 1.0_wp - q(k,j,i) ) )  &
2585                                                  * rho_air(k)
2586                    ENDDO
2587                 ENDIF
2588
2589              CASE ( 'pwv' ) ! water vapor partial pressure
2590                 IF ( humidity )  THEN
2591!                     DO  m = 1, vmea(l)%ns
2592!                        k = vmea(l)%k(m)
2593!                        j = vmea(l)%j(m)
2594!                        i = vmea(l)%i(m)
2595!                        vmea(l)%measured_vars(m,n) = ( q(k,j,i)                 &
2596!                                                     / ( 1.0_wp - q(k,j,i) ) )  &
2597!                                                   * rho_air(k)
2598!                     ENDDO
2599                 ENDIF
2600
2601              CASE ( 'hur' ) ! relative humidity
2602                 IF ( humidity )  THEN
2603                    DO  m = 1, vmea(l)%ns
2604                       k = vmea(l)%k(m)
2605                       j = vmea(l)%j(m)
2606                       i = vmea(l)%i(m)
2607!
2608!--                    Calculate actual temperature, water vapor saturation
2609!--                    pressure, and based on this  the saturation mixing ratio.
2610                       e_s  = magnus( exner(k) * pt(k,j,i) )
2611                       q_s  = rd_d_rv * e_s / ( hyp(k) - e_s )
2612                       q_wv = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) * rho_air(k)
2613
2614                       vmea(l)%measured_vars(m,n) = q_wv / ( q_s + 1E-10_wp )
2615                    ENDDO
2616                 ENDIF
2617
2618              CASE ( 'u', 'ua' ) ! u-component
2619                 DO  m = 1, vmea(l)%ns
2620                    k = vmea(l)%k(m)
2621                    j = vmea(l)%j(m)
2622                    i = vmea(l)%i(m)
2623                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
2624                 ENDDO
2625
2626              CASE ( 'v', 'va' ) ! v-component
2627                 DO  m = 1, vmea(l)%ns
2628                    k = vmea(l)%k(m)
2629                    j = vmea(l)%j(m)
2630                    i = vmea(l)%i(m)
2631                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
2632                 ENDDO
2633
2634              CASE ( 'w' ) ! w-component
2635                 DO  m = 1, vmea(l)%ns
2636                    k = MAX ( 1, vmea(l)%k(m) )
2637                    j = vmea(l)%j(m)
2638                    i = vmea(l)%i(m)
2639                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
2640                 ENDDO
2641
2642              CASE ( 'wspeed' ) ! horizontal wind speed
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) = SQRT(                         &
2648                                   ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
2649                                   ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2   &
2650                                                     )
2651                 ENDDO
2652
2653              CASE ( 'wdir' ) ! wind direction
2654                 DO  m = 1, vmea(l)%ns
2655                    k = vmea(l)%k(m)
2656                    j = vmea(l)%j(m)
2657                    i = vmea(l)%i(m)
2658
2659                    vmea(l)%measured_vars(m,n) = 180.0_wp + 180.0_wp / pi      &
2660                                               * ATAN2(                        &
2661                                          0.5_wp * ( v(k,j,i) + v(k,j+1,i) ),  &
2662                                          0.5_wp * ( u(k,j,i) + u(k,j,i+1) )   &
2663                                                      )
2664                 ENDDO
2665
2666              CASE ( 'utheta' )
2667                 DO  m = 1, vmea(l)%ns
2668                    k = vmea(l)%k(m)
2669                    j = vmea(l)%j(m)
2670                    i = vmea(l)%i(m)
2671                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
2672                                                 ( u(k,j,i) + u(k,j,i+1) ) *   &
2673                                                   pt(k,j,i)
2674                 ENDDO
2675
2676              CASE ( 'vtheta' )
2677                 DO  m = 1, vmea(l)%ns
2678                    k = vmea(l)%k(m)
2679                    j = vmea(l)%j(m)
2680                    i = vmea(l)%i(m)
2681                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
2682                                                 ( v(k,j,i) + v(k,j+1,i) ) *   &
2683                                                   pt(k,j,i)
2684                 ENDDO
2685
2686              CASE ( 'wtheta' )
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 *                      &
2692                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
2693                                                   pt(k,j,i)
2694                 ENDDO
2695
2696              CASE ( 'uqv' )
2697                 IF ( humidity )  THEN
2698                    DO  m = 1, vmea(l)%ns
2699                       k = vmea(l)%k(m)
2700                       j = vmea(l)%j(m)
2701                       i = vmea(l)%i(m)
2702                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
2703                                                    ( u(k,j,i) + u(k,j,i+1) ) *&
2704                                                      q(k,j,i)
2705                    ENDDO
2706                 ENDIF
2707
2708              CASE ( 'vqv' )
2709                 IF ( humidity )  THEN
2710                    DO  m = 1, vmea(l)%ns
2711                       k = vmea(l)%k(m)
2712                       j = vmea(l)%j(m)
2713                       i = vmea(l)%i(m)
2714                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
2715                                                    ( v(k,j,i) + v(k,j+1,i) ) *&
2716                                                      q(k,j,i)
2717                    ENDDO
2718                 ENDIF
2719
2720              CASE ( 'wqv' )
2721                 IF ( humidity )  THEN
2722                    DO  m = 1, vmea(l)%ns
2723                       k = MAX ( 1, vmea(l)%k(m) )
2724                       j = vmea(l)%j(m)
2725                       i = vmea(l)%i(m)
2726                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
2727                                                    ( w(k-1,j,i) + w(k,j,i) ) *&
2728                                                      q(k,j,i)
2729                    ENDDO
2730                 ENDIF
2731
2732              CASE ( 'uw' )
2733                 DO  m = 1, vmea(l)%ns
2734                    k = MAX ( 1, vmea(l)%k(m) )
2735                    j = vmea(l)%j(m)
2736                    i = vmea(l)%i(m)
2737                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
2738                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
2739                                                 ( u(k,j,i)   + u(k,j,i+1) )
2740                 ENDDO
2741
2742              CASE ( 'vw' )
2743                 DO  m = 1, vmea(l)%ns
2744                    k = MAX ( 1, vmea(l)%k(m) )
2745                    j = vmea(l)%j(m)
2746                    i = vmea(l)%i(m)
2747                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
2748                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
2749                                                 ( v(k,j,i)   + v(k,j+1,i) )
2750                 ENDDO
2751
2752              CASE ( 'uv' )
2753                 DO  m = 1, vmea(l)%ns
2754                    k = vmea(l)%k(m)
2755                    j = vmea(l)%j(m)
2756                    i = vmea(l)%i(m)
2757                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
2758                                                 ( u(k,j,i)   + u(k,j,i+1) ) * &
2759                                                 ( v(k,j,i)   + v(k,j+1,i) )
2760                 ENDDO
2761!
2762!--           Chemistry variables. List of variables may need extension.
2763!--           Note, gas species in PALM are in ppm and no distinction is made
2764!--           between mole-fraction and concentration quantities (all are
2765!--           output in ppm so far).
2766              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfno', 'mfno2',            &
2767                     'mcno', 'mcno2', 'tro3' )
2768                 IF ( air_chemistry )  THEN
2769!
2770!--                 First, search for the measured variable in the chem_vars
2771!--                 list, in order to get the internal name of the variable.
2772                    DO  nn = 1, UBOUND( chem_vars, 2 )
2773                       IF ( TRIM( vmea(l)%var_atts(n)%name ) ==                &
2774                            TRIM( chem_vars(0,nn) ) )  ind_chem = nn
2775                    ENDDO
2776!
2777!--                 Run loop over all chemical species, if the measured
2778!--                 variable matches the interal name, sample the variable.
2779!--                 Note, nvar as a chemistry-module variable.
2780                    DO  nn = 1, nvar
2781                       IF ( TRIM( chem_vars(1,ind_chem) ) ==                   &
2782                            TRIM( chem_species(nn)%name ) )  THEN
2783                          DO  m = 1, vmea(l)%ns
2784                             k = vmea(l)%k(m)
2785                             j = vmea(l)%j(m)
2786                             i = vmea(l)%i(m)
2787                             vmea(l)%measured_vars(m,n) =                      &
2788                                                   chem_species(nn)%conc(k,j,i)
2789                          ENDDO
2790                       ENDIF
2791                    ENDDO
2792                 ENDIF
2793
2794              CASE ( 'us' ) ! friction velocity
2795                 DO  m = 1, vmea(l)%ns
2796!
2797!--                 Surface data is only available on inner subdomains, not
2798!--                 on ghost points. Hence, limit the indices.
2799                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2800                    j = MERGE( j           , nyn, j            < nyn )
2801                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2802                    i = MERGE( i           , nxr, i            < nxr )
2803
2804                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
2805                             surf_def_h(0)%end_index(j,i)
2806                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
2807                    ENDDO
2808                    DO  mm = surf_lsm_h%start_index(j,i),                      &
2809                             surf_lsm_h%end_index(j,i)
2810                       vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm)
2811                    ENDDO
2812                    DO  mm = surf_usm_h%start_index(j,i),                      &
2813                             surf_usm_h%end_index(j,i)
2814                       vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm)
2815                    ENDDO
2816                 ENDDO
2817
2818              CASE ( 'thetas' ) ! scaling parameter temperature
2819                 DO  m = 1, vmea(l)%ns
2820!
2821!--                 Surface data is only available on inner subdomains, not
2822!--                 on ghost points. Hence, limit the indices.
2823                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2824                    j = MERGE( j           , nyn, j            < nyn )
2825                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2826                    i = MERGE( i           , nxr, i            < nxr )
2827
2828                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
2829                             surf_def_h(0)%end_index(j,i)
2830                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
2831                    ENDDO
2832                    DO  mm = surf_lsm_h%start_index(j,i),                      &
2833                             surf_lsm_h%end_index(j,i)
2834                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm)
2835                    ENDDO
2836                    DO  mm = surf_usm_h%start_index(j,i),                      &
2837                             surf_usm_h%end_index(j,i)
2838                       vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm)
2839                    ENDDO
2840                 ENDDO
2841
2842              CASE ( 'hfls' ) ! surface latent heat flux
2843                 DO  m = 1, vmea(l)%ns
2844!
2845!--                 Surface data is only available on inner subdomains, not
2846!--                 on ghost points. Hence, limit the indices.
2847                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2848                    j = MERGE( j           , nyn, j            < nyn )
2849                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2850                    i = MERGE( i           , nxr, i            < nxr )
2851
2852                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
2853                             surf_def_h(0)%end_index(j,i)
2854                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
2855                    ENDDO
2856                    DO  mm = surf_lsm_h%start_index(j,i),                      &
2857                             surf_lsm_h%end_index(j,i)
2858                       vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm)
2859                    ENDDO
2860                    DO  mm = surf_usm_h%start_index(j,i),                      &
2861                             surf_usm_h%end_index(j,i)
2862                       vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm)
2863                    ENDDO
2864                 ENDDO
2865
2866              CASE ( 'hfss' ) ! surface sensible heat flux
2867                 DO  m = 1, vmea(l)%ns
2868!
2869!--                 Surface data is only available on inner subdomains, not
2870!--                 on ghost points. Hence, limit the indices.
2871                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2872                    j = MERGE( j           , nyn, j            < nyn )
2873                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2874                    i = MERGE( i           , nxr, i            < nxr )
2875
2876                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
2877                             surf_def_h(0)%end_index(j,i)
2878                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
2879                    ENDDO
2880                    DO  mm = surf_lsm_h%start_index(j,i),                      &
2881                             surf_lsm_h%end_index(j,i)
2882                       vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm)
2883                    ENDDO
2884                    DO  mm = surf_usm_h%start_index(j,i),                      &
2885                             surf_usm_h%end_index(j,i)
2886                       vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm)
2887                    ENDDO
2888                 ENDDO
2889
2890              CASE ( 'hfdg' ) ! ground heat flux
2891                 DO  m = 1, vmea(l)%ns
2892!
2893!--                 Surface data is only available on inner subdomains, not
2894!--                 on ghost points. Hence, limit the indices.
2895                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2896                    j = MERGE( j           , nyn, j            < nyn )
2897                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2898                    i = MERGE( i           , nxr, i            < nxr )
2899
2900                    DO  mm = surf_lsm_h%start_index(j,i),                      &
2901                             surf_lsm_h%end_index(j,i)
2902                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ghf(mm)
2903                    ENDDO
2904                 ENDDO
2905
2906              CASE ( 'lwcs' )  ! liquid water of soil layer
2907!                  DO  m = 1, vmea(l)%ns
2908! !
2909! !--                 Surface data is only available on inner subdomains, not
2910! !--                 on ghost points. Hence, limit the indices.
2911!                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2912!                     j = MERGE( j           , nyn, j            < nyn )
2913!                     i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2914!                     i = MERGE( i           , nxr, i            < nxr )
2915!
2916!                     DO  mm = surf_lsm_h%start_index(j,i),                      &
2917!                              surf_lsm_h%end_index(j,i)
2918!                        vmea(l)%measured_vars(m,n) = ?
2919!                     ENDDO
2920!                  ENDDO
2921
2922              CASE ( 'rnds' ) ! surface net radiation
2923                 IF ( radiation )  THEN
2924                    DO  m = 1, vmea(l)%ns
2925!
2926!--                    Surface data is only available on inner subdomains, not
2927!--                    on ghost points. Hence, limit the indices.
2928                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2929                       j = MERGE( j           , nyn, j            < nyn )
2930                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2931                       i = MERGE( i           , nxr, i            < nxr )
2932
2933                       DO  mm = surf_lsm_h%start_index(j,i),                   &
2934                                surf_lsm_h%end_index(j,i)
2935                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm)
2936                       ENDDO
2937                       DO  mm = surf_usm_h%start_index(j,i),                   &
2938                                surf_usm_h%end_index(j,i)
2939                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm)
2940                       ENDDO
2941                    ENDDO
2942                 ENDIF
2943
2944              CASE ( 'rsus' ) ! surface shortwave out
2945                 IF ( radiation )  THEN
2946                    DO  m = 1, vmea(l)%ns
2947!
2948!--                    Surface data is only available on inner subdomains, not
2949!--                    on ghost points. Hence, limit the indices.
2950                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2951                       j = MERGE( j           , nyn, j            < nyn )
2952                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2953                       i = MERGE( i           , nxr, i            < nxr )
2954
2955                       DO  mm = surf_lsm_h%start_index(j,i),                   &
2956                                surf_lsm_h%end_index(j,i)
2957                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm)
2958                       ENDDO
2959                       DO  mm = surf_usm_h%start_index(j,i),                   &
2960                                surf_usm_h%end_index(j,i)
2961                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm)
2962                       ENDDO
2963                    ENDDO
2964                 ENDIF
2965
2966              CASE ( 'rsds' ) ! surface shortwave in
2967                 IF ( radiation )  THEN
2968                    DO  m = 1, vmea(l)%ns
2969!
2970!--                    Surface data is only available on inner subdomains, not
2971!--                    on ghost points. Hence, limit the indices.
2972                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2973                       j = MERGE( j           , nyn, j            < nyn )
2974                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2975                       i = MERGE( i           , nxr, i            < nxr )
2976
2977                       DO  mm = surf_lsm_h%start_index(j,i),                   &
2978                                surf_lsm_h%end_index(j,i)
2979                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm)
2980                       ENDDO
2981                       DO  mm = surf_usm_h%start_index(j,i),                   &
2982                                surf_usm_h%end_index(j,i)
2983                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm)
2984                       ENDDO
2985                    ENDDO
2986                 ENDIF
2987
2988              CASE ( 'rlus' ) ! surface longwave out
2989                 IF ( radiation )  THEN
2990                    DO  m = 1, vmea(l)%ns
2991!
2992!--                    Surface data is only available on inner subdomains, not
2993!--                    on ghost points. Hence, limit the indices.
2994                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
2995                       j = MERGE( j           , nyn, j            < nyn )
2996                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
2997                       i = MERGE( i           , nxr, i            < nxr )
2998
2999                       DO  mm = surf_lsm_h%start_index(j,i),                   &
3000                                surf_lsm_h%end_index(j,i)
3001                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm)
3002                       ENDDO
3003                       DO  mm = surf_usm_h%start_index(j,i),                   &
3004                                surf_usm_h%end_index(j,i)
3005                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm)
3006                       ENDDO
3007                    ENDDO
3008                 ENDIF
3009
3010              CASE ( 'rlds' ) ! surface longwave in
3011                 IF ( radiation )  THEN
3012                    DO  m = 1, vmea(l)%ns
3013!
3014!--                    Surface data is only available on inner subdomains, not
3015!--                    on ghost points. Hence, limit the indices.
3016                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
3017                       j = MERGE( j           , nyn, j            < nyn )
3018                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
3019                       i = MERGE( i           , nxr, i            < nxr )
3020
3021                       DO  mm = surf_lsm_h%start_index(j,i),                   &
3022                                surf_lsm_h%end_index(j,i)
3023                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm)
3024                       ENDDO
3025                       DO  mm = surf_usm_h%start_index(j,i),                   &
3026                                surf_usm_h%end_index(j,i)
3027                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm)
3028                       ENDDO
3029                    ENDDO
3030                 ENDIF
3031
3032              CASE ( 'rsd' ) ! shortwave in
3033                 IF ( radiation )  THEN
3034                    IF ( radiation_scheme /= 'rrtmg' )  THEN
3035                       DO  m = 1, vmea(l)%ns
3036                          k = 0
3037                          j = vmea(l)%j(m)
3038                          i = vmea(l)%i(m)
3039                          vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
3040                       ENDDO
3041                    ELSE
3042                       DO  m = 1, vmea(l)%ns
3043                          k = vmea(l)%k(m)
3044                          j = vmea(l)%j(m)
3045                          i = vmea(l)%i(m)
3046                          vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
3047                       ENDDO
3048                    ENDIF
3049                 ENDIF
3050
3051              CASE ( 'rsu' ) ! shortwave out
3052                 IF ( radiation )  THEN
3053                    IF ( radiation_scheme /= 'rrtmg' )  THEN
3054                       DO  m = 1, vmea(l)%ns
3055                          k = 0
3056                          j = vmea(l)%j(m)
3057                          i = vmea(l)%i(m)
3058                          vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
3059                       ENDDO
3060                    ELSE
3061                       DO  m = 1, vmea(l)%ns
3062                          k = vmea(l)%k(m)
3063                          j = vmea(l)%j(m)
3064                          i = vmea(l)%i(m)
3065                          vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
3066                       ENDDO
3067                    ENDIF
3068                 ENDIF
3069
3070              CASE ( 'rlu' ) ! longwave out
3071                 IF ( radiation )  THEN
3072                    IF ( radiation_scheme /= 'rrtmg' )  THEN
3073                       DO  m = 1, vmea(l)%ns
3074                          k = 0
3075                          j = vmea(l)%j(m)
3076                          i = vmea(l)%i(m)
3077                          vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
3078                       ENDDO
3079                    ELSE
3080                       DO  m = 1, vmea(l)%ns
3081                          k = vmea(l)%k(m)
3082                          j = vmea(l)%j(m)
3083                          i = vmea(l)%i(m)
3084                          vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
3085                       ENDDO
3086                    ENDIF
3087                 ENDIF
3088
3089              CASE ( 'rld' ) ! longwave in
3090                 IF ( radiation )  THEN
3091                    IF ( radiation_scheme /= 'rrtmg' )  THEN
3092                       DO  m = 1, vmea(l)%ns
3093                          k = 0
3094!
3095!--                       Surface data is only available on inner subdomains,
3096!--                       not on ghost points. Hence, limit the indices.
3097                          j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
3098                          j = MERGE( j           , nyn, j            < nyn )
3099                          i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
3100                          i = MERGE( i           , nxr, i            < nxr )
3101
3102                          vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
3103                       ENDDO
3104                    ELSE
3105                       DO  m = 1, vmea(l)%ns
3106                          k = vmea(l)%k(m)
3107                          j = vmea(l)%j(m)
3108                          i = vmea(l)%i(m)
3109                          vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
3110                       ENDDO
3111                    ENDIF
3112                 ENDIF
3113
3114              CASE ( 'rsddif' ) ! shortwave in, diffuse part
3115                 IF ( radiation )  THEN
3116                    DO  m = 1, vmea(l)%ns
3117                       j = vmea(l)%j(m)
3118                       i = vmea(l)%i(m)
3119
3120                       vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
3121                    ENDDO
3122                 ENDIF
3123
3124              CASE ( 't_soil' ) ! soil and wall temperature
3125                 DO  m = 1, vmea(l)%ns_soil
3126                    j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
3127                    j = MERGE( j                , nyn, j                 < nyn )
3128                    i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
3129                    i = MERGE( i                , nxr, i                 < nxr )
3130                    k = vmea(l)%k_soil(m)
3131
3132                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
3133                                surf_lsm_h%end_index(j,i)
3134                    match_usm = surf_usm_h%start_index(j,i) <=                 &
3135                                surf_usm_h%end_index(j,i)
3136
3137                    IF ( match_lsm )  THEN
3138                       mm = surf_lsm_h%start_index(j,i)
3139                       vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm)
3140                    ENDIF
3141
3142                    IF ( match_usm )  THEN
3143                       mm = surf_usm_h%start_index(j,i)
3144                       vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm)
3145                    ENDIF
3146                 ENDDO
3147
3148              CASE ( 'm_soil' ) ! soil moisture
3149                 DO  m = 1, vmea(l)%ns_soil
3150                    j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
3151                    j = MERGE( j                , nyn, j                 < nyn )
3152                    i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
3153                    i = MERGE( i                , nxr, i                 < nxr )
3154                    k = vmea(l)%k_soil(m)
3155
3156                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
3157                                surf_lsm_h%end_index(j,i)
3158
3159                    IF ( match_lsm )  THEN
3160                       mm = surf_lsm_h%start_index(j,i)
3161                       vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm)
3162                    ENDIF
3163
3164                 ENDDO
3165
3166              CASE ( 'ts' ) ! surface temperature
3167                 DO  m = 1, vmea(l)%ns
3168!
3169!--                 Surface data is only available on inner subdomains, not
3170!--                 on ghost points. Hence, limit the indices.
3171                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
3172                    j = MERGE( j           , nyn, j            < nyn )
3173                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
3174                    i = MERGE( i           , nxr, i            < nxr )
3175
3176                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
3177                             surf_def_h(0)%end_index(j,i)
3178                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%pt_surface(mm)
3179                    ENDDO
3180                    DO  mm = surf_lsm_h%start_index(j,i),                      &
3181                             surf_lsm_h%end_index(j,i)
3182                       vmea(l)%measured_vars(m,n) = surf_lsm_h%pt_surface(mm)
3183                    ENDDO
3184                    DO  mm = surf_usm_h%start_index(j,i),                      &
3185                             surf_usm_h%end_index(j,i)
3186                       vmea(l)%measured_vars(m,n) = surf_usm_h%pt_surface(mm)
3187                    ENDDO
3188                 ENDDO
3189
3190              CASE ( 'lwp' ) ! liquid water path
3191                 IF ( ASSOCIATED( ql ) )  THEN
3192                    DO  m = 1, vmea(l)%ns
3193                       j = vmea(l)%j(m)
3194                       i = vmea(l)%i(m)
3195
3196                       vmea(l)%measured_vars(m,n) = SUM( ql(nzb:nzt,j,i)       &
3197                                                       * dzw(1:nzt+1) )        &
3198                                                  * rho_surface
3199                    ENDDO
3200                 ENDIF
3201
3202              CASE ( 'ps' ) ! surface pressure
3203                 vmea(l)%measured_vars(:,n) = surface_pressure
3204
3205              CASE ( 'pswrtg' ) ! platform speed above ground
3206                 vmea(l)%measured_vars(:,n) = 0.0_wp
3207
3208              CASE ( 'pswrta' ) ! platform speed in air
3209                 vmea(l)%measured_vars(:,n) = 0.0_wp
3210
3211              CASE ( 't_lw' ) ! water temperature
3212                 DO  m = 1, vmea(l)%ns
3213!
3214!--                 Surface data is only available on inner subdomains, not
3215!--                 on ghost points. Hence, limit the indices.
3216                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
3217                    j = MERGE( j           , nyn, j            < nyn )
3218                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
3219                    i = MERGE( i           , nxr, i            < nxr )
3220
3221                    DO  mm = surf_lsm_h%start_index(j,i),                      &
3222                             surf_lsm_h%end_index(j,i)
3223                       IF ( surf_lsm_h%water_surface(m) )                      &
3224                          vmea(l)%measured_vars(m,n) = t_soil_h%var_2d(nzt,m)
3225                    ENDDO
3226
3227                 ENDDO
3228!
3229!--           More will follow ...
3230
3231!
3232!--           No match found - just set a fill value
3233              CASE DEFAULT
3234                 vmea(l)%measured_vars(:,n) = vmea(l)%fillout
3235           END SELECT
3236
3237        ENDDO
3238
3239     ENDDO
3240
3241  END SUBROUTINE vm_sampling
3242
3243
3244 END MODULE virtual_measurement_mod
Note: See TracBrowser for help on using the repository browser.