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

Last change on this file since 4446 was 4444, checked in by raasch, 5 years ago

bugfix: cpp-directives for serial mode added

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