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

Last change on this file since 3988 was 3988, checked in by kanani, 2 years ago

enable steering of output interval for virtual measurements

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