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

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

Bugfix, do not set boundary conditions for pt in neutral runs; enable checks + additional checks for virtual measurements

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