Home


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

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

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

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