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

Last change on this file since 3787 was 3766, checked in by raasch, 5 years ago

unused_variables removed, bugfix in im_define_netcdf_grid argument list, trim added to avoid truncation compiler warnings, save attribute added to local targets to avoid outlive pointer target warning, first argument removed from module_interface_rrd_*, file module_interface reformatted with respect to coding standards, bugfix in surface_data_output_rrd_local (variable k removed)

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