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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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