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

Last change on this file since 3704 was 3704, checked in by suehring, 2 years ago

Revision of virtual-measurement module and data output enabled. Further, post-processing tool added to merge distributed virtually sampled data and to output it into NetCDF files.

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 73.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! -----------------
[3704]22! - initialization revised
23! - binary data output
24! - list of allowed variables extended
[3434]25!
26! Former revisions:
27! -----------------
[3472]28! $Id: virtual_measurement_mod.f90 3704 2019-01-29 19:51:41Z suehring $
[3522]29! Sampling of variables
30!
31! 3494 2018-11-06 14:51:27Z suehring
[3494]32! Bugfixing
33!
34! 3473 2018-10-30 20:50:15Z suehring
[3473]35! Initial revision
[3434]36!
37! Authors:
38! --------
[3522]39! @author Matthias Suehring
[3434]40!
41! Description:
42! ------------
[3471]43!> The module acts as an interface between 'real-world' observations and
44!> model simulations. Virtual measurements will be taken in the model at the
[3704]45!> coordinates representative for the 'real-world' observation coordinates.
[3471]46!> More precisely, coordinates and measured quanties will be read from a
47!> NetCDF file which contains all required information. In the model,
48!> the same quantities (as long as all the required components are switched-on)
49!> will be sampled at the respective positions and output into an extra file,
50!> which allows for straight-forward comparison of model results with
51!> observations.
[3522]52!>
53!> @todo list_of_allowed variables needs careful checking
54!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
55!>       the (UC)2 standard
56!> @note Fluxes are not processed
[3434]57!------------------------------------------------------------------------------!
[3471]58 MODULE virtual_measurement_mod
[3434]59
60    USE arrays_3d,                                                             &
61        ONLY:  q, pt, u, v, w, zu, zw
62
[3522]63    USE chem_modules,                                                          &
64        ONLY:  nspec
65
66    USE chemistry_model_mod,                                                   &
67        ONLY:  chem_species
68       
[3434]69    USE control_parameters,                                                    &
[3704]70        ONLY:  air_chemistry, dz, humidity, io_blocks, io_group, neutral,      &
71               message_string, time_since_reference_point, virtual_measurement
[3434]72
73    USE cpulog,                                                                &
74        ONLY:  cpu_log, log_point
75       
76    USE grid_variables,                                                        &
77        ONLY:  dx, dy
78
79    USE indices,                                                               &
[3704]80        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny, wall_flags_0
[3434]81
82    USE kinds
[3704]83   
84    USE netcdf_data_input_mod,                                                 &
85        ONLY:  init_model
86       
87    USE pegrid
88   
89    USE surface_mod,                                                           &
90        ONLY:  surf_lsm_h, surf_usm_h
91       
92    USE land_surface_model_mod,                                                &
93        ONLY:  nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h
94       
95    USE radiation_model_mod
96       
97    USE urban_surface_mod,                                                     &
98        ONLY:  nzb_wall, nzt_wall, t_wall_h
[3434]99
100
101    IMPLICIT NONE
[3704]102   
103    TYPE virt_general
104       INTEGER(iwp) ::  id_vm     !< NetCDF file id for virtual measurements
105       INTEGER(iwp) ::  nvm = 0   !< number of virtual measurements
106    END TYPE virt_general
[3434]107
108    TYPE virt_mea
109   
[3704]110       CHARACTER(LEN=100)  ::  feature_type      !< type of the measurement
111       CHARACTER(LEN=100)  ::  filename_original !< name of the original file
112       CHARACTER(LEN=100)  ::  site              !< name of the measurement site
[3434]113   
114       CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE ::  measured_vars_name !< name of the measured variables
115   
[3704]116       INTEGER(iwp) ::  ns = 0          !< number of observation coordinates on subdomain, for atmospheric measurements
117       INTEGER(iwp) ::  ns_tot = 0      !< total number of observation coordinates, for atmospheric measurements
118       INTEGER(iwp) ::  ntraj           !< number of trajectories of a measurement
119       INTEGER(iwp) ::  nvar            !< number of measured variables (atmosphere + soil)
[3434]120       
[3704]121       INTEGER(iwp) ::  ns_soil = 0     !< number of observation coordinates on subdomain, for soil measurements
122       INTEGER(iwp) ::  ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements
123       
[3434]124       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory or station that are no _FillValues
125       
[3704]126       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< grid index for measurement position in x-direction
127       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< grid index for measurement position in y-direction
128       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< grid index for measurement position in k-direction
129       
130       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i_soil  !< grid index for measurement position in x-direction
131       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j_soil  !< grid index for measurement position in y-direction
132       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k_soil  !< grid index for measurement position in k-direction
[3434]133           
134       LOGICAL ::  trajectory         = .FALSE. !< flag indicating that the observation is a mobile observation
135       LOGICAL ::  timseries          = .FALSE. !< flag indicating that the observation is a stationary point measurement
136       LOGICAL ::  timseries_profile  = .FALSE. !< flag indicating that the observation is a stationary profile measurement
[3704]137       LOGICAL ::  soil_sampling      = .FALSE. !< flag indicating that soil state variables were sampled
[3434]138       
[3704]139       REAL(wp) ::  fill_eutm          !< fill value for UTM coordinates in case of missing values
140       REAL(wp) ::  fill_nutm          !< fill value for UTM coordinates in case of missing values
141       REAL(wp) ::  fill_zag           !< fill value for heigth coordinates in case of missing values
142       REAL(wp) ::  fillout = -999.9   !< fill value for output in case a observation is taken from inside a building
143       REAL(wp) ::  origin_x_obs       !< origin of the observation in UTM coordiates in x-direction
144       REAL(wp) ::  origin_y_obs       !< origin of the observation in UTM coordiates in y-direction
145       
146       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  z_ag           !< measurement height above ground level
147       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth          !< measurement depth in soil
[3522]148             
[3704]149       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars       !< measured variables
150       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars_soil  !< measured variables
[3434]151       
152    END TYPE virt_mea
153
154    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"                      !< dimension name for UTM coordinate easting
155    CHARACTER(LEN=11) ::  char_feature = "featureType"             !< attribute name for feature type
[3704]156   
157    ! This need to generalized
158    CHARACTER(LEN=8)  ::  char_filename = "filename"               !< attribute name for filename
159    CHARACTER(LEN=11) ::  char_soil = "soil_sample"                !< attribute name for soil sampling indication
[3434]160    CHARACTER(LEN=10) ::  char_fillvalue = "_FillValue"            !< variable attribute name for _FillValue
161    CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
162    CHARACTER(LEN=5)  ::  char_nutm = "N_UTM"                      !< dimension name for UTM coordinate northing
163    CHARACTER(LEN=18) ::  char_numstations = "number_of_stations"  !< attribute name for number of stations
164    CHARACTER(LEN=8)  ::  char_origx = "origin_x"                  !< attribute name for station coordinate in x
165    CHARACTER(LEN=8)  ::  char_origy = "origin_y"                  !< attribute name for station coordinate in y
166    CHARACTER(LEN=4)  ::  char_site = "site"                       !< attribute name for site name
167    CHARACTER(LEN=19) ::  char_zag = "height_above_ground"         !< attribute name for height above ground variable
168    CHARACTER(LEN=10) ::  type_ts   = 'timeSeries'                 !< name of stationary point measurements
169    CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
170    CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
[3704]171   
172    CHARACTER(LEN=6), DIMENSION(1:5) ::  soil_vars       = (/                  & !< list of soil variables
173                            't_soil',                                          &
174                            'm_soil',                                          &
175                            'lwc   ',                                          &
176                            'lwcs  ',                                          &
177                            'smp   '                       /)
178                           
179    CHARACTER(LEN=10), DIMENSION(0:1,1:8) ::  chem_vars = RESHAPE( (/          &
180                                              'mcpm1     ', 'kc_PM1    ',      &
181                                              'mcpm2p5   ', 'kc_PM2.5  ',      &
182                                              'mcpm25    ', 'kc_PM25   ',      &
183                                              'mcpm10    ', 'kc_PM10   ',      &
184                                              'mfno2     ', 'kc_NO2    ',      &
185                                              'mfno      ', 'kc_NO     ',      &
186                                              'tro3      ', 'kc_O3     ',      &
187                                              'mfco      ', 'kc_CO     '       &
188                                                                   /), (/ 2, 8 /) )
[3522]189!
190!-- MS: List requires careful revision!
[3704]191    CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
[3471]192       (/ 'hfls      ',  & ! surface latent heat flux (W/m2)
193          'hfss      ',  & ! surface sensible heat flux (W/m2)
194          'hur       ',  & ! relative humidity (-)
195          'hus       ',  & ! specific humidity (g/kg)
196          'haa       ',  & ! absolute atmospheric humidity (kg/m3)
197          'mcpm1     ',  & ! mass concentration of PM1 (kg/m3)
198          'mcpm2p5   ',  & ! mass concentration of PM2.5 (kg/m3)
199          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
200          'mcco      ',  & ! mass concentration of CO (kg/m3)
201          'mcco2     ',  & ! mass concentration of CO2 (kg/m3)
202          'mcbcda    ',  & ! mass concentration of black carbon paritcles (kg/m3)
203          'ncaa      ',  & ! number concentation of particles (1/m3)
[3522]204          'mfco      ',  & ! mole fraction of CO (mol/mol)
[3471]205          'mfco2     ',  & ! mole fraction of CO2 (mol/mol)
206          'mfch4     ',  & ! mole fraction of methane (mol/mol)
207          'mfnh3     ',  & ! mole fraction of amonia (mol/mol)
208          'mfno      ',  & ! mole fraction of nitrogen monoxide (mol/mol)
209          'mfno2     ',  & ! mole fraction of nitrogen dioxide (mol/mol)
210          'mfso2     ',  & ! mole fraction of sulfur dioxide (mol/mol)
211          'mfh20     ',  & ! mole fraction of water (mol/mol)
212          'plev      ',  & ! ? air pressure - hydrostaic + perturbation?
213          'rlds      ',  & ! surface downward longwave flux  (W/m2)
214          'rlus      ',  & ! surface upward longwave flux (W/m2)
215          'rsds      ',  & ! surface downward shortwave flux (W/m2)
216          'rsus      ',  & ! surface upward shortwave flux (W/m2)
217          'ta        ',  & ! air temperature (degree C)
218          't_va      ',  & ! virtual accoustic temperature (K)
219          'theta     ',  & ! potential temperature (K)
220          'tro3      ',  & ! mole fraction of ozone air (mol/mol)
221          'ts        ',  & ! scaling parameter of temperature (K)
222          'wspeed    ',  & ! ? wind speed - horizontal?
223          'wdir      ',  & ! wind direction
224          'us        ',  & ! friction velocity
225          'msoil     ',  & ! ? soil moisture - which depth? 
226          'tsoil     ',  & ! ? soil temperature - which depth?                                                               
227          'u         ',  & ! u-component
[3704]228          'utheta    ',  & ! total eastward kinematic heat flux
[3471]229          'ua        ',  & ! eastward wind (is there any difference to u?)
230          'v         ',  & ! v-component
[3704]231          'vtheta    ',  & ! total northward kinematic heat flux
[3471]232          'va        ',  & ! northward wind (is there any difference to v?)
233          'w         ',  & ! w-component
[3704]234          'wtheta    ',  & ! total vertical kinematic heat flux
[3471]235          'rld       ',  & ! downward longwave radiative flux (W/m2)
236          'rlu       ',  & ! upnward longwave radiative flux (W/m2)
237          'rsd       ',  & ! downward shortwave radiative flux (W/m2)
238          'rsu       ',  & ! upward shortwave radiative flux (W/m2)
239          'rsddif    ',  & ! downward shortwave diffuse radiative flux (W/m2)
[3704]240          'rnds      ',  & ! surface net downward radiative flux (W/m2)
241          't_soil    ',  &
242          'm_soil    ',  &
243          'lwc       ',  &
244          'lwcs      ',  &
245          'smp       '   &
[3471]246       /)
[3704]247                                                           
[3434]248   
[3704]249    LOGICAL ::  global_attribute = .TRUE.         !< flag indicating a global attribute
250    LOGICAL ::  init = .TRUE.                     !< flag indicating initialization of data output
[3434]251    LOGICAL ::  use_virtual_measurement = .FALSE. !< Namelist parameter
252   
253    REAL(wp) ::  vm_time_start = 0.0              !< time after virtual measurements should start
254
[3704]255    TYPE( virt_general )                        ::  vmea_general !< data structure which encompass general variables
[3434]256    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea !< virtual measurement data structure
257   
258    INTERFACE vm_check_parameters
259       MODULE PROCEDURE vm_check_parameters
260    END INTERFACE vm_check_parameters
261   
[3704]262    INTERFACE vm_data_output
263       MODULE PROCEDURE vm_data_output
264    END INTERFACE vm_data_output
265   
[3434]266    INTERFACE vm_init
267       MODULE PROCEDURE vm_init
268    END INTERFACE vm_init
269   
[3704]270    INTERFACE vm_last_actions
271       MODULE PROCEDURE vm_last_actions
272    END INTERFACE vm_last_actions
273   
[3434]274    INTERFACE vm_parin
275       MODULE PROCEDURE vm_parin
276    END INTERFACE vm_parin
277   
278    INTERFACE vm_sampling
279       MODULE PROCEDURE vm_sampling
280    END INTERFACE vm_sampling
281
282    SAVE
283
284    PRIVATE
285
286!
287!-- Public interfaces
[3704]288    PUBLIC  vm_check_parameters, vm_data_output, vm_init, vm_last_actions,     &
289            vm_parin, vm_sampling
[3434]290
291!
292!-- Public variables
[3704]293    PUBLIC  vmea, vmea_general, vm_time_start
[3434]294
295 CONTAINS
296
297
298!------------------------------------------------------------------------------!
299! Description:
300! ------------
[3471]301!> Check parameters for virtual measurement module
[3434]302!------------------------------------------------------------------------------!
303 SUBROUTINE vm_check_parameters
304
305    USE control_parameters,                                                    &
306        ONLY:  message_string, virtual_measurement
307 
308    USE netcdf_data_input_mod,                                                 &
309        ONLY:  input_pids_static
310       
311    IMPLICIT NONE
312   
313!
314!-- In case virtual measurements are taken, a static input file is required.
315!-- This is because UTM coordinates for the PALM domain origin are required
316!-- for correct mapping of the measurements.
317!-- ToDo: Revise this later and remove this requirement.
318    IF ( virtual_measurement  .AND.  .NOT. input_pids_static )  THEN
[3704]319       message_string = 'If virtual measurements are taken, a static input ' //&
[3434]320                        'file is mandatory.'
321       CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
322    ENDIF
323 
324 END SUBROUTINE vm_check_parameters
325 
326!------------------------------------------------------------------------------!
327! Description:
328! ------------
[3471]329!> Read namelist for the virtual measurement module
[3434]330!------------------------------------------------------------------------------!
331 SUBROUTINE vm_parin
332 
333    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
334 
335    NAMELIST /virtual_measurement_parameters/  use_virtual_measurement,        &
336                                               vm_time_start
337
338    line = ' '
339
340!
341!-- Try to find stg package
342    REWIND ( 11 )
343    line = ' '
344    DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
345       READ ( 11, '(A)', END=20 )  line
346    ENDDO
347    BACKSPACE ( 11 )
348
349!
350!-- Read namelist
351    READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
352
353!
[3471]354!-- Set flag that indicates that the virtual measurement module is switched on
[3434]355    IF ( use_virtual_measurement )  virtual_measurement = .TRUE.
356   
357    GOTO 20
358
359 10 BACKSPACE( 11 )
360    READ( 11 , '(A)') line
361    CALL parin_fail_message( 'virtual_measurement_parameters', line )
362
363 20 CONTINUE
364 
365 END SUBROUTINE vm_parin
366
367
368!------------------------------------------------------------------------------!
369! Description:
370! ------------
371!> Initialize virtual measurements: read coordiante arrays and measured
372!> variables, set indicies indicating the measurement points, read further
373!> attributes, etc..
374!------------------------------------------------------------------------------!
375 SUBROUTINE vm_init
376
377    USE arrays_3d,                                                             &
378        ONLY:  zu, zw
379       
380    USE grid_variables,                                                        &
381        ONLY:  ddx, ddy, dx, dy
382       
383    USE indices,                                                               &
384        ONLY:  nxl, nxr, nyn, nys
385 
386    USE netcdf_data_input_mod,                                                 &
[3704]387        ONLY:  init_model, input_file_vm,                                      &
388               netcdf_data_input_get_dimension_length,                         &
[3434]389               netcdf_data_input_att, netcdf_data_input_var
390               
391    USE surface_mod,                                                           &
392        ONLY:  get_topography_top_index_ji
393       
394    IMPLICIT NONE
395   
396    CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
[3704]397    CHARACTER(LEN=5)    ::  dummy_read                !< dummy string indicate station id
[3434]398    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
[3522]399    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed   
[3434]400   
[3704]401    INTEGER(iwp) ::  dim_eutm  !< dimension size of UTM easting coordinate
402    INTEGER(iwp) ::  dim_nutm  !< dimension size of UTM northing coordinate
[3434]403    INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate
[3704]404    INTEGER(iwp) ::  dim_zag   !< dimension size of height coordinate
405    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
406    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
[3434]407    INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
[3704]408    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
[3434]409    INTEGER(iwp) ::  js        !< grid index of real observation point of the respective station in y-direction
[3704]410    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
[3522]411    INTEGER(iwp) ::  kl        !< lower vertical index of surrounding grid points of an observation coordinate
[3434]412    INTEGER(iwp) ::  ks        !< grid index of real observation point of the respective station in z-direction
413    INTEGER(iwp) ::  ksurf     !< topography top index
[3522]414    INTEGER(iwp) ::  ku        !< upper vertical index of surrounding grid points of an observation coordinate
[3434]415    INTEGER(iwp) ::  l         !< running index over all stations
416    INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
417    INTEGER(iwp) ::  ll        !< running index over all measured variables in file
418    INTEGER(iwp) ::  lll       !< running index over all allowed variables
419    INTEGER(iwp) ::  n         !< running index over trajectory coordinates
420    INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
421    INTEGER(iwp) ::  t         !< running index over number of trajectories
[3704]422    INTEGER(iwp) ::  m
[3434]423   
[3704]424    INTEGER(KIND=1)::  soil_dum
425   
426    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns_all !< dummy array used to sum-up the number of observation coordinates
427   
[3522]428    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
429   
430    LOGICAL ::  chem_include !< flag indicating that chemical species is considered in modelled mechanism
431    LOGICAL ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
432   
[3434]433    REAL(wp)     ::  fill_eutm !< _FillValue for coordinate array E_UTM
434    REAL(wp)     ::  fill_nutm !< _FillValue for coordinate array N_UTM
435    REAL(wp)     ::  fill_zag  !< _FillValue for height coordinate
436   
[3437]437    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm !< easting UTM coordinate, temporary variable
438    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm !< northing UTM coordinate, temporary variable,
439    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  z_ag  !< height coordinate relative to origin_z, temporary variable
[3434]440!
[3704]441!-- Obtain number of sites. Also, pass the 'open' string, in order to initially
442!-- open the measurement driver.
443    CALL netcdf_data_input_att( vmea_general%nvm, char_numstations,            &
444                                vmea_general%id_vm, input_file_vm,             &
[3434]445                                global_attribute, 'open', '' )
[3522]446                               
[3434]447!
[3704]448!-- Allocate data structure which encompass all required information, such as
449!-- grid points indicies, absolute UTM coordinates, the measured quantities,
450!-- etc. .
451    ALLOCATE( vmea(1:vmea_general%nvm) )
[3434]452!
[3704]453!-- Allocate flag array. This dummy array is used to identify grid points
454!-- where virtual measurements should be taken. Please note, at least one
455!-- ghost point is required, in order to include also the surrounding
456!-- grid points of the original coordinate. 
[3522]457    ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
458    meas_flag = 0
459!
[3704]460!-- Loop over all sites.
461    DO  l = 1, vmea_general%nvm
[3434]462!
[3704]463!--    Determine suffix which contains the ID, ordered according to the number
464!--    of measurements.
[3434]465       IF( l < 10 )  THEN
466          WRITE( dum, '(I1)')  l
467       ELSEIF( l < 100 )  THEN
468          WRITE( dum, '(I2)')  l
469       ELSEIF( l < 1000 )  THEN
470          WRITE( dum, '(I3)')  l
471       ELSEIF( l < 10000 )  THEN
472          WRITE( dum, '(I4)')  l
473       ELSEIF( l < 100000 )  THEN
474          WRITE( dum, '(I5)')  l
475       ENDIF
[3704]476!
477!--    Read site coordinates (UTM).
478       CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx //         &
479                                   TRIM( dum ), vmea_general%id_vm, '',        &
480                                   global_attribute, '', '' )
481       CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy //         &
482                                   TRIM( dum ), vmea_general%id_vm, '',        &
483                                   global_attribute, '', '' )
484!
485!--    Read site name                 
486       CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ),     &
487                                   vmea_general%id_vm, '', global_attribute,   &
[3434]488                                   '', '' )
[3704]489!
490!--    Read type of the measurement (trajectory, profile, timeseries).
491       CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature //       &
492                                   TRIM( dum ), vmea_general%id_vm, '',        &
493                                   global_attribute, '', '' )
494!
495!--    Read the name of the original file where observational data is stored.
496       CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // &
497                                   TRIM( dum ), vmea_general%id_vm, '',        &
498                                   global_attribute, '', '' )
499!
500!--    Read a flag which indicates that also soil quantities are take at the
501!--    respective site (is part of the virtual measurement driver). 
502       CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ),         &
503                                   vmea_general%id_vm, '', global_attribute,   &
[3434]504                                   '', '' )
505!
[3704]506!--    Set flag for soil-sampling.
507       IF ( soil_dum == 1 )  vmea(l)%soil_sampling = .TRUE.
508!
[3434]509!---   Set logicals depending on the type of the measurement
510       IF ( INDEX( vmea(l)%feature_type, type_tspr     ) /= 0 )  THEN
511          vmea(l)%timseries_profile = .TRUE.
512       ELSEIF ( INDEX( vmea(l)%feature_type, type_ts   ) /= 0 )  THEN
513          vmea(l)%timseries         = .TRUE.
514       ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 )  THEN
515          vmea(l)%trajectory        = .TRUE.
[3704]516!
517!--   Give error message in case the type matches non of the pre-defined types.
[3434]518       ELSE
519          message_string = 'Attribue featureType = ' //                        &
520                           TRIM( vmea(l)%feature_type ) //                     &
521                           ' is not allowed.' 
522          CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 )
523       ENDIF
524!
[3704]525!--    Read string with all measured variables at this site
[3434]526       measured_variables_file = ''
527       CALL netcdf_data_input_var( measured_variables_file,                    &
[3704]528                                   char_mv // TRIM( dum ), vmea_general%id_vm )
[3434]529!
[3704]530!--    Count the number of measured variables. Only count variables that match
531!--    with the allowed variables.
532!--    Please note, for some NetCDF interal reasons characters end with a NULL,
533!--    i.e. also empty characters contain a NULL. Therefore, check the strings
534!--    for a NULL to get the correct character length in order to compare
535!--    them with the list of allowed variables.
536       vmea(l)%nvar   = 0
[3434]537       DO ll = 1, SIZE( measured_variables_file )
538          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
539               measured_variables_file(ll)(1:1) /= ' ')  THEN
540!
541!--          Obtain character length of the character
542             len_char = 1
543             DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
544                 .AND.  measured_variables_file(ll)(len_char:len_char) /= ' ' )
545                len_char = len_char + 1
546             ENDDO
547             len_char = len_char - 1
548!
549!--          Now, compare the measured variable with the list of allowed
550!--          variables.
551             DO  lll= 1, SIZE( list_allowed_variables )
552                IF ( measured_variables_file(ll)(1:len_char) ==                &
553                     TRIM( list_allowed_variables(lll) ) )  THEN
554                   vmea(l)%nvar = vmea(l)%nvar + 1
555                   measured_variables(vmea(l)%nvar) =                          &
556                                       measured_variables_file(ll)(1:len_char)
557                ENDIF
558             ENDDO
559          ENDIF
560       ENDDO
561!
[3704]562!--    Allocate array for the measured variables names for the respective site.
[3434]563       ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
564
565       DO  ll = 1, vmea(l)%nvar
566          vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
567       ENDDO
568!
[3522]569!--    In case of chemistry, check if species is considered in the modelled
570!--    chemistry mechanism.
[3704]571!        IF ( air_chemistry )  THEN
572!           DO  ll = 1, vmea(l)%nvar
573!              chem_include = .FALSE.
574!              DO  n = 1, nspec
575!                 IF ( TRIM( vmea(l)%measured_vars_name(ll) ) ==                 &
576!                      TRIM( chem_species(n)%name ) )  chem_include = .TRUE.
577!              ENDDO
578! !
579! !--  Revise this. It should only check for chemistry variables and not for all!
580!              IF ( .NOT. chem_include )  THEN
581!                 message_string = TRIM( vmea(l)%measured_vars_name(ll) ) //     &
582!                                  ' is not considered in the modelled '  //     &
583!                                  'chemistry mechanism'
584!                 CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
585!              ENDIF
586!           ENDDO
587!        ENDIF
[3522]588!
[3704]589!--    Read the UTM coordinates for the actual site. Based on the coordinates,
590!--    define the grid-index space on each subdomain where virtual measurements
591!--    should be taken. Note, the entire coordinate arrays will not be stored
592!--    as this would exceed memory requirements, particularly for trajectory
593!--    measurements.
[3434]594       IF ( vmea(l)%nvar > 0 )  THEN
595!
596!--       For stationary measurements UTM coordinates are just one value and
597!--       its dimension is "station", while for mobile measurements UTM
[3704]598!--       coordinates are arrays depending on the number of trajectories and
599!--       time, according to (UC)2 standard. First, inquire dimension length
600!--       of the UTM coordinates.
[3434]601          IF ( vmea(l)%trajectory )  THEN
602!
603!--          For non-stationary measurements read the number of trajectories
[3704]604!--          and the number of time coordinates.
605             CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, &
[3434]606                                                          vmea(l)%ntraj,      &
607                                                          "traj" //           &
608                                                          TRIM( dum ) )
[3704]609             CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, &
610                                                          dim_ntime,          &
[3434]611                                                          "ntime" //          &
612                                                          TRIM( dum ) )
613!
[3704]614!--       For stationary measurements the dimension for UTM and time
615!--       coordinates is 1.
[3434]616          ELSE
617             vmea(l)%ntraj  = 1
618             dim_ntime = 1
619          ENDIF
620!
621!-        Allocate array which defines individual time frame for each
[3704]622!--       trajectory or station.
[3434]623          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
624!
625!--       Allocate temporary arrays for UTM and height coordinates. Note,
626!--       on file UTM coordinates might be 1D or 2D variables
[3437]627          ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
628          ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
629          ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime)  )
[3434]630!
[3704]631!--       Read _FillValue attributes of the coordinate dimensions.
[3434]632          CALL netcdf_data_input_att( fill_eutm, char_fillvalue,               &
[3704]633                                      vmea_general%id_vm, '',                  &
634                                      .NOT. global_attribute, '',              &
[3434]635                                      char_eutm // TRIM( dum ) )
636          CALL netcdf_data_input_att( fill_nutm, char_fillvalue,               &
[3704]637                                      vmea_general%id_vm, '',                  &
638                                      .NOT. global_attribute, '',              &
[3434]639                                      char_nutm // TRIM( dum ) )
640          CALL netcdf_data_input_att( fill_zag, char_fillvalue,                &
[3704]641                                      vmea_general%id_vm, '',                  &
642                                      .NOT. global_attribute, '',              &
[3434]643                                      char_zag  // TRIM( dum ) )
644!
645!--       Read UTM and height coordinates coordinates for all trajectories and
646!--       times.
[3437]647          IF ( vmea(l)%trajectory )  THEN
[3704]648             CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ),      &
649                                         vmea_general%id_vm,                   &
[3437]650                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
[3704]651             CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ),      &
652                                         vmea_general%id_vm,                   &
[3437]653                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
[3704]654             CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ),        &
655                                         vmea_general%id_vm,                   &
[3437]656                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
657          ELSE
[3704]658             CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), &
659                                         vmea_general%id_vm )
660             CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), &
661                                         vmea_general%id_vm )
662             CALL netcdf_data_input_var( z_ag(1,:),  char_zag  // TRIM( dum ), &
663                                         vmea_general%id_vm )
664          ENDIF         
[3434]665!
666!--       Based on UTM coordinates, check if the measurement station or parts
667!--       of the trajectory is on subdomain. This case, setup grid index space
668!--       sample these quantities.
[3522]669          meas_flag = 0
[3434]670          DO  t = 1, vmea(l)%ntraj
[3704]671!             
672!--          First, compute relative x- and y-coordinates with respect to the
673!--          lower-left origin of the model domain, which is the difference
674!--          betwen UTM coordinates. Note, if the origin is not correct, the
675!--          virtual sites will be misplaced.
676             e_utm(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
677             n_utm(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
[3434]678!
679!--          Determine the individual time coordinate length for each station and
680!--          trajectory. This is required as several stations and trajectories
681!--          are merged into one file but they do not have the same number of
682!--          points in time, hence, missing values may occur and cannot be
[3704]683!--          processed further. This is actually a work-around for the specific
684!--          (UC)2 dataset, but it won't harm in anyway.
[3434]685             vmea(l)%dim_t(t) = 0
686             DO  n = 1, dim_ntime
[3437]687                IF ( e_utm(t,n) /= fill_eutm  .AND.                            &
688                     n_utm(t,n) /= fill_nutm  .AND.                            &
689                     z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
[3434]690             ENDDO
691!
692!--          Compute grid indices relative to origin and check if these are
693!--          on the subdomain. Note, virtual measurements will be taken also
694!--          at grid points surrounding the station, hence, check also for
695!--          these grid points.
[3437]696             DO  n = 1, vmea(l)%dim_t(t)
697                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
698                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
[3434]699!
700!--             Is the observation point on subdomain?
701                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.                   &
702                          js >= nys  .AND.  js <= nyn )
703!
[3522]704!--             Check if observation coordinate is on subdomain
[3434]705                IF ( on_pe )  THEN
[3522]706!
707!--                Determine vertical index which correspond to the observation
708!--                height.
[3434]709                   ksurf = get_topography_top_index_ji( js, is, 's' )
[3437]710                   ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
[3434]711!
[3522]712!--                Set mask array at the observation coordinates. Also, flag the
713!--                surrounding coordinate points, but first check whether the
714!--                surrounding coordinate points are on the subdomain.
[3704]715                   kl = MERGE( ks-1, ks, ks-1 >= nzb  .AND. ks-1 >= ksurf )
716                   ku = MERGE( ks+1, ks, ks+1 < nzt+1 )
[3522]717                 
[3704]718                   DO  i = is-1, is+1
719                      DO  j = js-1, js+1
720                         DO  k = kl, ku
721                            meas_flag(k,j,i) = MERGE(                          &
722                                             IBSET( meas_flag(k,j,i), 0 ),     &
723                                             0,                                &
724                                             BTEST( wall_flags_0(k,j,i), 0 )   &
725                                                    )
726                         ENDDO
727                      ENDDO
728                   ENDDO
[3434]729                ENDIF
730             ENDDO
731             
732          ENDDO
733!
[3704]734!--       Based on the flag array count the number of of sampling coordinates.
735!--       Please note, sampling coordinates in atmosphere and soil may be
736!--       different, as within the soil all levels will be measured.           
737!--       Hence, count individually. Start with atmoshere.
[3522]738          ns = 0
[3704]739          DO  i = nxl-1, nxr+1
740             DO  j = nys-1, nyn+1
741                DO  k = nzb, nzt+1
742                   ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
[3522]743                ENDDO
744             ENDDO
745          ENDDO
[3704]746         
[3522]747!
[3434]748!--       Store number of observation points on subdomain and allocate index
[3704]749!--       arrays as well as array containing height information.
[3434]750          vmea(l)%ns = ns
751         
752          ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
753          ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
754          ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
[3704]755          ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) )         
[3434]756!
[3522]757!--       Based on the flag array store the grid indices which correspond to
758!--       the observation coordinates.
[3704]759          ns = 0
760          DO  i = nxl-1, nxr+1
761             DO  j = nys-1, nyn+1
762                DO  k = nzb, nzt+1
763                   IF ( BTEST( meas_flag(k,j,i), 0 ) )  THEN
[3522]764                      ns = ns + 1
[3704]765                      vmea(l)%i(ns) = i
766                      vmea(l)%j(ns) = j
767                      vmea(l)%k(ns) = k
768                      vmea(l)%z_ag(ns)  = zu(k) -                              &
769                                   zw(get_topography_top_index_ji( j, i, 's' ))
[3522]770                   ENDIF
771                ENDDO
[3434]772             ENDDO
773          ENDDO
774!
[3704]775!--       Same for the soil. Based on the flag array, count the number of
776!--       sampling coordinates in soil. Sample at all soil levels in this case.
777          IF ( vmea(l)%soil_sampling )  THEN
778             DO  i = nxl, nxr
779                DO  j = nys, nyn
780                   IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
781                      IF ( surf_lsm_h%start_index(j,i) <=                      &
782                           surf_lsm_h%end_index(j,i) )  THEN
783                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
784                                                      nzt_soil - nzb_soil + 1 
785                      ENDIF
786                      IF ( surf_usm_h%start_index(j,i) <=                      &
787                           surf_usm_h%end_index(j,i) )  THEN
788                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
789                                                      nzt_wall - nzb_wall + 1 
790                      ENDIF
791                   ENDIF
792                ENDDO
793             ENDDO
794          ENDIF         
795!
796!--       Allocate index arrays as well as array containing height information
797!--       for soil.
798          IF ( vmea(l)%soil_sampling )  THEN
799             ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) )
800             ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) )
801             ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) )
802             ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) )
803          ENDIF     
804!
805!--       For soil, store the grid indices.
806          ns = 0
807          IF ( vmea(l)%soil_sampling )  THEN
808             DO  i = nxl, nxr
809                DO  j = nys, nyn
810                   IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
811                      IF ( surf_lsm_h%start_index(j,i) <=                      &
812                           surf_lsm_h%end_index(j,i) )  THEN
813                         m = surf_lsm_h%start_index(j,i)
814                         DO  k = nzb_soil, nzt_soil
815                            ns = ns + 1
816                            vmea(l)%i_soil(ns) = i
817                            vmea(l)%j_soil(ns) = j
818                            vmea(l)%k_soil(ns) = k
819                            vmea(l)%depth(ns)  = zs(k)
820                         ENDDO
821                      ENDIF
822                     
823                      IF ( surf_usm_h%start_index(j,i) <=                      &
824                           surf_usm_h%end_index(j,i) )  THEN
825                         m = surf_usm_h%start_index(j,i)
826                         DO  k = nzb_wall, nzt_wall
827                            ns = ns + 1
828                            vmea(l)%i_soil(ns) = i
829                            vmea(l)%j_soil(ns) = j
830                            vmea(l)%k_soil(ns) = k
831                            vmea(l)%depth(ns)  = surf_usm_h%zw(k,m)
832                         ENDDO
833                      ENDIF
834                   ENDIF
835                ENDDO
836             ENDDO
837          ENDIF
838!
[3434]839!--       Allocate array to save the sampled values.
[3704]840          ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) )
841         
842          IF ( vmea(l)%soil_sampling )                                         &
843             ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,           &
844                                                  1:vmea(l)%nvar) )
[3434]845!
[3704]846!--       Initialize with _FillValues
847          vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) = vmea(l)%fillout
848          IF ( vmea(l)%soil_sampling )                                         &
849             vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nvar) =    &
850                                                                vmea(l)%fillout
[3434]851!
852!--       Deallocate temporary coordinate arrays
853          IF ( ALLOCATED( e_utm ) )  DEALLOCATE( e_utm )
854          IF ( ALLOCATED( n_utm ) )  DEALLOCATE( n_utm )
855          IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( z_ag  )
[3704]856          IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( vmea(l)%dim_t )
[3434]857       ENDIF
858    ENDDO
859!
860!-- Close input file for virtual measurements. Therefore, just call
861!-- the read attribute routine with the "close" option.
[3704]862    CALL netcdf_data_input_att( vmea_general%nvm, char_numstations,            &
863                                vmea_general%id_vm, '',                        &
[3434]864                                global_attribute, 'close', '' )
[3704]865!
866!-- Sum-up the number of observation coordiates, for atmosphere first.
867!-- This is actually only required for data output.
868    ALLOCATE( ns_all(1:vmea_general%nvm) )
869    ns_all = 0   
870#if defined( __parallel )
871    CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER,  &
872                        MPI_SUM, comm2d, ierr )
873#else
874    ns_all(:) = vmea(:)%ns
875#endif
876    vmea(:)%ns_tot = ns_all(:)
877!
878!-- Now for soil
879    ns_all = 0   
880#if defined( __parallel )
881    CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,          &
882                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
883#else
884    ns_all(:) = vmea(:)%ns_soil
885#endif
886    vmea(:)%ns_soil_tot = ns_all(:)
887   
888    DEALLOCATE( ns_all )
[3522]889!                               
890!-- Dellocate flag array
891    DEALLOCATE( meas_flag )
[3704]892!
893!-- Initialize binary data output of virtual measurements.
894!-- Open binary output file.
895    CALL check_open( 27 )
896!
897!-- Output header information.
898    CALL vm_data_output
[3522]899       
[3434]900  END SUBROUTINE vm_init
901 
902 
903!------------------------------------------------------------------------------!
904! Description:
905! ------------
[3704]906!> Binary data output.
907!------------------------------------------------------------------------------!
908  SUBROUTINE vm_data_output
909   
910     USE pegrid
911   
912     IMPLICIT NONE
913         
914     INTEGER(iwp) ::  i         !< running index over IO blocks   
915     INTEGER(iwp) ::  l         !< running index over all stations
916     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
917!
918!--  Header output on each PE
919     IF ( init )  THEN
920
921        DO  i = 0, io_blocks-1
922           IF ( i == io_group )  THEN
923              WRITE ( 27 )  'number of measurements            '
924              WRITE ( 27 )  vmea_general%nvm
925
926              DO  l = 1, vmea_general%nvm
927                 WRITE ( 27 )  'site                              '
928                 WRITE ( 27 )  vmea(l)%site
929                 WRITE ( 27 )  'file                              '
930                 WRITE ( 27 )  vmea(l)%filename_original
931                 WRITE ( 27 )  'feature_type                      '
932                 WRITE ( 27 )  vmea(l)%feature_type
933                 WRITE ( 27 )  'origin_x_obs                      '
934                 WRITE ( 27 )  vmea(l)%origin_x_obs
935                 WRITE ( 27 )  'origin_y_obs                      '
936                 WRITE ( 27 )  vmea(l)%origin_y_obs
937                 WRITE ( 27 )  'total number of observation points'                               
938                 WRITE ( 27 )  vmea(l)%ns_tot
939                 WRITE ( 27 )  'number of measured variables      '
940                 WRITE ( 27 )  vmea(l)%nvar
941                 WRITE ( 27 )  'variables                         '
942                 WRITE ( 27 )  vmea(l)%measured_vars_name(:)
943                 WRITE ( 27 )  'number of observation points      '
944                 WRITE ( 27 )  vmea(l)%ns
945                 WRITE ( 27 )  'E_UTM                             '
946                 WRITE ( 27 )  init_model%origin_x +                           &
947                        REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx
948                 WRITE ( 27 )  'N_UTM                             '
949                 WRITE ( 27 )  init_model%origin_y +                           &
950                        REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy
951                 WRITE ( 27 )  'Z_AG                              '
952                 WRITE ( 27 )  vmea(l)%z_ag(1:vmea(l)%ns)
953                 WRITE ( 27 )  'soil sampling                     '
954                 WRITE ( 27 )  MERGE( 'yes                               ',    &
955                                      'no                                ',    &
956                                      vmea(l)%soil_sampling )
957 
958                 IF ( vmea(l)%soil_sampling )  THEN                 
959                    WRITE ( 27 )  'total number of soil points       '                               
960                    WRITE ( 27 )  vmea(l)%ns_soil_tot
961                    print*, "vmea(l)%ns_soil_tot", vmea(l)%ns_soil_tot
962                    WRITE ( 27 )  'number of soil points             '
963                    WRITE ( 27 )  vmea(l)%ns_soil
964                    WRITE ( 27 )  'E_UTM soil                        '
965                    WRITE ( 27 )  init_model%origin_x +                        &
966                           REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp,   &
967                                 KIND = wp ) * dx
968                    WRITE ( 27 )  'N_UTM soil                        '
969                    WRITE ( 27 )  init_model%origin_y +                        &
970                           REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp,   &
971                                 KIND = wp ) * dy
972                    WRITE ( 27 )  'DEPTH                             '
973                    WRITE ( 27 )  vmea(l)%depth(1:vmea(l)%ns_soil)
974                 ENDIF
975              ENDDO
976
977           ENDIF
978        ENDDO
979       
980#if defined( __parallel )
981        CALL MPI_BARRIER( comm2d, ierr )
982#endif
983!
984!--     After header information is written, set control flag to .FALSE.
985        init = .FALSE.
986!
987!--  Data output at each measurement timestep on each PE
988     ELSE
989        DO  i = 0, io_blocks-1
990
991           IF ( i == io_group )  THEN
992              WRITE( 27 )  'output time                       '
993              WRITE( 27 )  time_since_reference_point
994              DO  l = 1, vmea_general%nvm
995!
996!--              Skip binary writing if no observation points are defined on PE
997                 IF ( vmea(l)%ns < 1  .AND.  vmea(l)%ns_soil < 1)  CYCLE                 
998                 DO  n = 1, vmea(l)%nvar
999                    WRITE( 27 )  vmea(l)%measured_vars_name(n)
1000                    IF ( vmea(l)%soil_sampling  .AND.                           &
1001                         ANY( TRIM( vmea(l)%measured_vars_name(n))  ==          &
1002                              soil_vars ) )  THEN                   
1003                       WRITE( 27 )  vmea(l)%measured_vars_soil(:,n)
1004                    ELSE
1005                       WRITE( 27 )  vmea(l)%measured_vars(:,n)
1006                    ENDIF
1007                 ENDDO
1008           
1009              ENDDO
1010           ENDIF
1011        ENDDO
1012#if defined( __parallel )
1013        CALL MPI_BARRIER( comm2d, ierr )
1014#endif
1015     ENDIF
1016 
1017  END SUBROUTINE vm_data_output 
1018 
1019 
1020!------------------------------------------------------------------------------!
1021! Description:
1022! ------------
1023!> Write end-of-file statement as last action.
1024!------------------------------------------------------------------------------!
1025  SUBROUTINE vm_last_actions
1026   
1027     USE pegrid
1028   
1029     IMPLICIT NONE
1030         
1031     INTEGER(iwp) ::  i         !< running index over IO blocks   
1032     INTEGER(iwp) ::  l         !< running index over all stations
1033     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
1034 
1035     DO  i = 0, io_blocks-1
1036        IF ( i == io_group )  THEN
1037           WRITE( 27 )  'EOF                               '
1038        ENDIF
1039     ENDDO
1040#if defined( __parallel )
1041        CALL MPI_BARRIER( comm2d, ierr )
1042#endif
1043!
1044!--  Close binary file
1045     CALL close_file( 27 )
1046 
1047  END SUBROUTINE vm_last_actions
1048 
1049!------------------------------------------------------------------------------!
1050! Description:
1051! ------------
[3434]1052!> Sampling of the actual quantities along the observation coordinates
1053!------------------------------------------------------------------------------!
[3471]1054  SUBROUTINE vm_sampling
[3434]1055
[3522]1056    USE arrays_3d,                                                             &
1057        ONLY:  exner, pt, q, u, v, w
[3471]1058
[3522]1059    USE basic_constants_and_equations_mod,                                     &
1060        ONLY:  pi
[3434]1061   
[3522]1062    USE radiation_model_mod,                                                   &
1063        ONLY:  radiation 
1064
1065    USE surface_mod,                                                           &
1066        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
1067   
[3434]1068     IMPLICIT NONE
1069     
[3704]1070     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
[3434]1071     
[3704]1072     INTEGER(iwp) ::  i         !< grid index in x-direction
1073     INTEGER(iwp) ::  j         !< grid index in y-direction
1074     INTEGER(iwp) ::  k         !< grid index in z-direction
1075     INTEGER(iwp) ::  ind_chem  !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming
1076     INTEGER(iwp) ::  l         !< running index over the number of stations
1077     INTEGER(iwp) ::  m         !< running index over all virtual observation coordinates
1078     INTEGER(iwp) ::  mm        !< index of surface element which corresponds to the virtual observation coordinate
1079     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
1080     INTEGER(iwp) ::  nn        !< running index over the number of chemcal species
1081     
1082     LOGICAL ::  match_lsm !< flag indicating natural-type surface
1083     LOGICAL ::  match_usm !< flag indicating urban-type surface
[3434]1084!
[3704]1085!--  Loop over all sites.
1086     DO  l = 1, vmea_general%nvm
[3434]1087!
[3704]1088!--     At the beginning, set _FillValues
1089        IF ( ALLOCATED( vmea(l)%measured_vars      ) )                         &
1090           vmea(l)%measured_vars      = vmea(l)%fillout
1091        IF ( ALLOCATED( vmea(l)%measured_vars_soil ) )                         &
1092           vmea(l)%measured_vars_soil = vmea(l)%fillout
1093!
1094!--     Loop over all variables measured at this site. 
[3522]1095        DO  n = 1, vmea(l)%nvar
1096       
1097           SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
1098           
1099              CASE ( 'theta' )
1100                 IF ( .NOT. neutral )  THEN
1101                    DO  m = 1, vmea(l)%ns
1102                       k = vmea(l)%k(m)
1103                       j = vmea(l)%j(m)
1104                       i = vmea(l)%i(m)
[3704]1105                       vmea(l)%measured_vars(m,n) = pt(k,j,i)
[3522]1106                    ENDDO
1107                 ENDIF
1108                 
[3704]1109              CASE ( 'ta' )
[3522]1110                 IF ( .NOT. neutral )  THEN
1111                    DO  m = 1, vmea(l)%ns
1112                       k = vmea(l)%k(m)
1113                       j = vmea(l)%j(m)
1114                       i = vmea(l)%i(m)
[3704]1115                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )
[3522]1116                    ENDDO
1117                 ENDIF
[3704]1118             
1119              CASE ( 't_va' )
[3522]1120                 
1121              CASE ( 'hus', 'haa' )
1122                 IF ( humidity )  THEN
1123                    DO  m = 1, vmea(l)%ns
1124                       k = vmea(l)%k(m)
1125                       j = vmea(l)%j(m)
1126                       i = vmea(l)%i(m)
[3704]1127                       vmea(l)%measured_vars(m,n) = q(k,j,i)
[3522]1128                    ENDDO
1129                 ENDIF
1130                 
1131              CASE ( 'u', 'ua' )
1132                 DO  m = 1, vmea(l)%ns
1133                    k = vmea(l)%k(m)
1134                    j = vmea(l)%j(m)
1135                    i = vmea(l)%i(m)
[3704]1136                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
[3522]1137                 ENDDO
1138                 
1139              CASE ( 'v', 'va' )
1140                 DO  m = 1, vmea(l)%ns
1141                    k = vmea(l)%k(m)
1142                    j = vmea(l)%j(m)
1143                    i = vmea(l)%i(m)
[3704]1144                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
[3522]1145                 ENDDO
1146                 
1147              CASE ( 'w' )
1148                 DO  m = 1, vmea(l)%ns
[3704]1149                    k = MAX ( 1, vmea(l)%k(m) ) 
[3522]1150                    j = vmea(l)%j(m)
1151                    i = vmea(l)%i(m)
[3704]1152                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[3522]1153                 ENDDO
1154                 
1155              CASE ( 'wspeed' )
1156                 DO  m = 1, vmea(l)%ns
1157                    k = vmea(l)%k(m)
1158                    j = vmea(l)%j(m)
1159                    i = vmea(l)%i(m)
[3704]1160                    vmea(l)%measured_vars(m,n) = SQRT(                         &
[3522]1161                                   ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
1162                                   ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2   &
1163                                                     )
1164                 ENDDO
1165                 
1166              CASE ( 'wdir' )
1167                 DO  m = 1, vmea(l)%ns
1168                    k = vmea(l)%k(m)
1169                    j = vmea(l)%j(m)
1170                    i = vmea(l)%i(m)
1171                   
[3704]1172                    vmea(l)%measured_vars(m,n) = ATAN2(                        &
[3522]1173                                       - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ),   &
1174                                       - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )    &
1175                                                      ) * 180.0_wp / pi
1176                 ENDDO
[3704]1177                 
1178              CASE ( 'utheta' )
1179                 DO  m = 1, vmea(l)%ns
1180                    k = vmea(l)%k(m)
1181                    j = vmea(l)%j(m)
1182                    i = vmea(l)%i(m)
1183                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
1184                                                 ( u(k,j,i) + u(k,j,i+1) ) *   &
1185                                                   pt(k,j,i)
1186                 ENDDO
1187                 
1188              CASE ( 'vtheta' )
1189                 DO  m = 1, vmea(l)%ns
1190                    k = vmea(l)%k(m)
1191                    j = vmea(l)%j(m)
1192                    i = vmea(l)%i(m)
1193                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
1194                                                 ( v(k,j,i) + v(k,j+1,i) ) *   &
1195                                                   pt(k,j,i)
1196                 ENDDO
1197                 
1198              CASE ( 'wtheta' )
1199                 DO  m = 1, vmea(l)%ns
1200                    k = MAX ( 1, vmea(l)%k(m) )
1201                    j = vmea(l)%j(m)
1202                    i = vmea(l)%i(m)
1203                    vmea(l)%measured_vars(m,n) = 0.5_wp *                      &
1204                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
1205                                                   pt(k,j,i)
1206                 ENDDO
1207                 
1208              CASE ( 'uw' )
1209                 DO  m = 1, vmea(l)%ns
1210                    k = MAX ( 1, vmea(l)%k(m) )
1211                    j = vmea(l)%j(m)
1212                    i = vmea(l)%i(m)
1213                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
1214                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
1215                                                 ( u(k,j,i)   + u(k,j,i+1) )
1216                 ENDDO
1217                 
1218              CASE ( 'vw' )
1219                 DO  m = 1, vmea(l)%ns
1220                    k = MAX ( 1, vmea(l)%k(m) )
1221                    j = vmea(l)%j(m)
1222                    i = vmea(l)%i(m)
1223                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
1224                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
1225                                                 ( v(k,j,i)   + v(k,j+1,i) )
1226                 ENDDO
1227                 
1228              CASE ( 'uv' )
1229                 DO  m = 1, vmea(l)%ns
1230                    k = MAX ( 1, vmea(l)%k(m) )
1231                    j = vmea(l)%j(m)
1232                    i = vmea(l)%i(m)
1233                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
1234                                                 ( u(k,j,i)   + u(k,j,i+1) ) * &
1235                                                 ( v(k,j,i)   + v(k,j+1,i) )
1236                 ENDDO
[3522]1237!
[3704]1238!--           List of variables may need extension.
1239              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2',    & 
1240                     'tro3' )                     
1241                 IF ( air_chemistry )  THEN
1242!
1243!--                 First, search for the measured variable in the chem_vars
1244!--                 list, in order to get the internal name of the variable.
1245                    DO  nn = 1, UBOUND( chem_vars, 2 )
1246                       IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
1247                            TRIM( chem_vars(0,nn) ) )  ind_chem = nn
1248                    ENDDO
1249!
1250!--                 Run loop over all chemical species, if the measured
1251!--                 variable matches the interal name, sample the variable.
[3522]1252                    DO  nn = 1, nspec                   
[3704]1253                       IF ( TRIM( chem_vars(1,ind_chem) ) ==                   &
[3522]1254                            TRIM( chem_species(nn)%name ) )  THEN                           
1255                          DO  m = 1, vmea(l)%ns             
1256                             k = vmea(l)%k(m)
1257                             j = vmea(l)%j(m)
1258                             i = vmea(l)%i(m)                   
[3704]1259                             vmea(l)%measured_vars(m,n) =                      &
[3522]1260                                                   chem_species(nn)%conc(k,j,i)
1261                          ENDDO
1262                       ENDIF
1263                    ENDDO
1264                 ENDIF
1265                 
1266              CASE ( 'us' )
1267                 DO  m = 1, vmea(l)%ns
1268!
1269!--                 Surface data is only available on inner subdomains, not
1270!--                 on ghost points. Hence, limit the indices.
1271                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1272                    j = MERGE( j           , nyn, j            > nyn )
1273                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1274                    i = MERGE( i           , nxr, i            > nxr )
1275                   
1276                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
1277                             surf_def_h(0)%end_index(j,i)
[3704]1278                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
[3522]1279                    ENDDO
1280                    DO  mm = surf_lsm_h%start_index(j,i),                      &
1281                             surf_lsm_h%end_index(j,i)
[3704]1282                       vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm)
[3522]1283                    ENDDO
1284                    DO  mm = surf_usm_h%start_index(j,i),                      &
1285                             surf_usm_h%end_index(j,i)
[3704]1286                       vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm)
[3522]1287                    ENDDO
1288                 ENDDO
1289                 
1290              CASE ( 'ts' )
1291                 DO  m = 1, vmea(l)%ns
1292!
1293!--                 Surface data is only available on inner subdomains, not
1294!--                 on ghost points. Hence, limit the indices.
1295                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1296                    j = MERGE( j           , nyn, j            > nyn )
1297                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1298                    i = MERGE( i           , nxr, i            > nxr )
1299                   
1300                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
1301                             surf_def_h(0)%end_index(j,i)
[3704]1302                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
[3522]1303                    ENDDO
1304                    DO  mm = surf_lsm_h%start_index(j,i),                      &
1305                             surf_lsm_h%end_index(j,i)
[3704]1306                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm)
[3522]1307                    ENDDO
1308                    DO  mm = surf_usm_h%start_index(j,i),                      &
1309                             surf_usm_h%end_index(j,i)
[3704]1310                       vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm)
[3522]1311                    ENDDO
1312                 ENDDO
1313                 
1314              CASE ( 'hfls' )
1315                 DO  m = 1, vmea(l)%ns
1316!
1317!--                 Surface data is only available on inner subdomains, not
1318!--                 on ghost points. Hence, limit the indices.
1319                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1320                    j = MERGE( j           , nyn, j            > nyn )
1321                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1322                    i = MERGE( i           , nxr, i            > nxr )
1323                   
1324                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
1325                             surf_def_h(0)%end_index(j,i)
[3704]1326                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
[3522]1327                    ENDDO
1328                    DO  mm = surf_lsm_h%start_index(j,i),                      &
1329                             surf_lsm_h%end_index(j,i)
[3704]1330                       vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm)
[3522]1331                    ENDDO
1332                    DO  mm = surf_usm_h%start_index(j,i),                      &
1333                             surf_usm_h%end_index(j,i)
[3704]1334                       vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm)
[3522]1335                    ENDDO
1336                 ENDDO
1337                 
1338              CASE ( 'hfss' )
1339                 DO  m = 1, vmea(l)%ns
1340!
1341!--                 Surface data is only available on inner subdomains, not
1342!--                 on ghost points. Hence, limit the indices.
1343                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1344                    j = MERGE( j           , nyn, j            > nyn )
1345                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1346                    i = MERGE( i           , nxr, i            > nxr )
1347                   
1348                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
1349                             surf_def_h(0)%end_index(j,i)
[3704]1350                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
[3522]1351                    ENDDO
1352                    DO  mm = surf_lsm_h%start_index(j,i),                      &
1353                             surf_lsm_h%end_index(j,i)
[3704]1354                       vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm)
[3522]1355                    ENDDO
1356                    DO  mm = surf_usm_h%start_index(j,i),                      &
1357                             surf_usm_h%end_index(j,i)
[3704]1358                       vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm)
[3522]1359                    ENDDO
1360                 ENDDO
1361                 
1362              CASE ( 'rnds' )
1363                 IF ( radiation )  THEN
1364                    DO  m = 1, vmea(l)%ns
1365!
1366!--                    Surface data is only available on inner subdomains, not
1367!--                    on ghost points. Hence, limit the indices.
1368                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1369                       j = MERGE( j           , nyn, j            > nyn )
1370                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1371                       i = MERGE( i           , nxr, i            > nxr )
1372                   
1373                       DO  mm = surf_lsm_h%start_index(j,i),                   &
1374                                surf_lsm_h%end_index(j,i)
[3704]1375                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm)
[3522]1376                       ENDDO
1377                       DO  mm = surf_usm_h%start_index(j,i),                   &
1378                                surf_usm_h%end_index(j,i)
[3704]1379                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm)
[3522]1380                       ENDDO
1381                    ENDDO
1382                 ENDIF
1383                 
[3704]1384              CASE ( 'rsus' )
[3522]1385                 IF ( radiation )  THEN
1386                    DO  m = 1, vmea(l)%ns
1387!
1388!--                    Surface data is only available on inner subdomains, not
1389!--                    on ghost points. Hence, limit the indices.
1390                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1391                       j = MERGE( j           , nyn, j            > nyn )
1392                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1393                       i = MERGE( i           , nxr, i            > nxr )
1394                   
1395                       DO  mm = surf_lsm_h%start_index(j,i),                   &
1396                                surf_lsm_h%end_index(j,i)
[3704]1397                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm)
[3522]1398                       ENDDO
1399                       DO  mm = surf_usm_h%start_index(j,i),                   &
1400                                surf_usm_h%end_index(j,i)
[3704]1401                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm)
[3522]1402                       ENDDO
1403                    ENDDO
1404                 ENDIF
1405                 
[3704]1406              CASE ( 'rsds' )
[3522]1407                 IF ( radiation )  THEN
1408                    DO  m = 1, vmea(l)%ns
1409!
1410!--                    Surface data is only available on inner subdomains, not
1411!--                    on ghost points. Hence, limit the indices.
1412                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1413                       j = MERGE( j           , nyn, j            > nyn )
1414                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1415                       i = MERGE( i           , nxr, i            > nxr )
1416                   
1417                       DO  mm = surf_lsm_h%start_index(j,i),                   &
1418                                surf_lsm_h%end_index(j,i)
[3704]1419                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm)
[3522]1420                       ENDDO
1421                       DO  mm = surf_usm_h%start_index(j,i),                   &
1422                                surf_usm_h%end_index(j,i)
[3704]1423                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm)
[3522]1424                       ENDDO
1425                    ENDDO
1426                 ENDIF
1427                 
[3704]1428              CASE ( 'rlus' )
[3522]1429                 IF ( radiation )  THEN
1430                    DO  m = 1, vmea(l)%ns
1431!
1432!--                    Surface data is only available on inner subdomains, not
1433!--                    on ghost points. Hence, limit the indices.
1434                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1435                       j = MERGE( j           , nyn, j            > nyn )
1436                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1437                       i = MERGE( i           , nxr, i            > nxr )
1438                   
1439                       DO  mm = surf_lsm_h%start_index(j,i),                   &
1440                                surf_lsm_h%end_index(j,i)
[3704]1441                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm)
[3522]1442                       ENDDO
1443                       DO  mm = surf_usm_h%start_index(j,i),                   &
1444                                surf_usm_h%end_index(j,i)
[3704]1445                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm)
[3522]1446                       ENDDO
1447                    ENDDO
1448                 ENDIF
1449                 
[3704]1450              CASE ( 'rlds' )
[3522]1451                 IF ( radiation )  THEN
1452                    DO  m = 1, vmea(l)%ns
1453!
1454!--                    Surface data is only available on inner subdomains, not
1455!--                    on ghost points. Hence, limit the indices.
1456                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1457                       j = MERGE( j           , nyn, j            > nyn )
1458                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1459                       i = MERGE( i           , nxr, i            > nxr )
1460                   
1461                       DO  mm = surf_lsm_h%start_index(j,i),                   &
1462                                surf_lsm_h%end_index(j,i)
[3704]1463                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm)
[3522]1464                       ENDDO
1465                       DO  mm = surf_usm_h%start_index(j,i),                   &
1466                                surf_usm_h%end_index(j,i)
[3704]1467                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm)
[3522]1468                       ENDDO
1469                    ENDDO
1470                 ENDIF
[3704]1471                 
1472              CASE ( 'rsd' )
1473                 IF ( radiation )  THEN
1474                    DO  m = 1, vmea(l)%ns
1475                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 
1476                       j = vmea(l)%j(m)
1477                       i = vmea(l)%i(m)
1478                   
1479                       vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
1480                    ENDDO
1481                 ENDIF
1482                 
1483              CASE ( 'rsu' )
1484                 IF ( radiation )  THEN
1485                    DO  m = 1, vmea(l)%ns
1486                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 
1487                       j = vmea(l)%j(m)
1488                       i = vmea(l)%i(m)
1489                   
1490                       vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
1491                    ENDDO
1492                 ENDIF
1493                 
1494              CASE ( 'rlu' )
1495                 IF ( radiation )  THEN
1496                    DO  m = 1, vmea(l)%ns
1497                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 
1498                       j = vmea(l)%j(m)
1499                       i = vmea(l)%i(m)
1500                   
1501                       vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
1502                    ENDDO
1503                 ENDIF
1504                 
1505              CASE ( 'rld' )
1506                 IF ( radiation )  THEN
1507                    DO  m = 1, vmea(l)%ns
1508                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 
1509                       j = vmea(l)%j(m)
1510                       i = vmea(l)%i(m)
1511                   
1512                       vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
1513                    ENDDO
1514                 ENDIF
1515                 
1516              CASE ( 'rsddif' )
1517                 IF ( radiation )  THEN
1518                    DO  m = 1, vmea(l)%ns
1519                       j = vmea(l)%j(m)
1520                       i = vmea(l)%i(m)
1521                   
1522                       vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
1523                    ENDDO
1524                 ENDIF
1525                 
1526              CASE ( 't_soil' )
1527                 DO  m = 1, vmea(l)%ns_soil
1528                    i = vmea(l)%i_soil(m)
1529                    j = vmea(l)%j_soil(m)
1530                    k = vmea(l)%k_soil(m)
1531                   
1532                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
1533                                surf_lsm_h%end_index(j,i)
1534                    match_usm = surf_usm_h%start_index(j,i) <=                 &
1535                                surf_usm_h%end_index(j,i)
1536                               
1537                    IF ( match_lsm )  THEN
1538                       mm = surf_lsm_h%start_index(j,i)
1539                       vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm)
1540                    ENDIF
1541                   
1542                    IF ( match_usm )  THEN
1543                       mm = surf_usm_h%start_index(j,i)
1544                       vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm)
1545                    ENDIF
1546                 ENDDO
1547                 
1548              CASE ( 'm_soil' )
1549                 DO  m = 1, vmea(l)%ns_soil
1550                    i = vmea(l)%i_soil(m)
1551                    j = vmea(l)%j_soil(m)
1552                    k = vmea(l)%k_soil(m)
1553                   
1554                    match_lsm = surf_lsm_h%start_index(j,i) <=                 &
1555                                surf_lsm_h%end_index(j,i)
1556                               
1557                    IF ( match_lsm )  THEN
1558                       mm = surf_lsm_h%start_index(j,i)
1559                       vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm)
1560                    ENDIF
1561                   
1562                 ENDDO
[3522]1563!
1564!--           More will follow ...
[3704]1565
1566!
1567!--           No match found - just set a fill value
1568              CASE DEFAULT
1569                 vmea(l)%measured_vars(:,n) = vmea(l)%fillout
[3522]1570           END SELECT
1571
[3494]1572        ENDDO
[3434]1573
1574     ENDDO
[3704]1575         
[3471]1576  END SUBROUTINE vm_sampling
[3434]1577 
1578
[3471]1579 END MODULE virtual_measurement_mod
Note: See TracBrowser for help on using the repository browser.