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

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

typo and print statement removed

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