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

Last change on this file since 4356 was 4346, checked in by motisi, 5 years ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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