Home


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

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

Last change on this file since 4845 was 4843, checked in by raasch, 4 years ago

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