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

Last change on this file since 4434 was 4422, checked in by suehring, 5 years ago

missing trim() in open statement

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