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

Last change on this file since 3876 was 3876, checked in by knoop, 2 years ago

Moved "photolysis_scheme", "chem_species" and "phot_frequen" to chem_modules

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