Home


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

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

Last change on this file since 4788 was 4764, checked in by suehring, 4 years ago

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