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

Last change on this file since 4438 was 4438, checked in by suehring, 19 months ago

Synthetic turbulence: performance optimizations - random numbers only defined and computed locally, option to compute velocity seeds locally without need of global communication; paralell random number generator: new routine to initialize 1D random number arrays; virtual measurements: CPU-log points added

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