source: palm/trunk/SOURCE/virtual_measurement_mod.f90

Last change on this file was 4888, checked in by suehring, 5 months ago

Timeseries output of ghf and r_a as sum of all horizontally upward-facing surfaces in LSM+USM simulations; virtual measurements: Reduce number of additional grid point output in z for timeseries data to a reasonable number

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