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

Last change on this file since 4329 was 4329, checked in by motisi, 18 months ago

Renamed wall_flags_0 to wall_flags_static_0

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