Home


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

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

Last change on this file since 4679 was 4671, checked in by pavelkrc, 4 years ago

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