Home


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

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

Last change on this file since 4901 was 4888, checked in by suehring, 4 years ago

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