Home


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

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

Last change on this file since 4842 was 4828, checked in by Giersch, 3 years ago

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