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

Last change on this file since 4310 was 4226, checked in by suehring, 5 years ago

Offline nesting: data input modularized; time variable is defined relative to time_utc_init, so that input data is correctly mapped to actual model time; checks rephrased and new checks for the time dimension added; Netcdf input: routine to retrieve dimension length renamed

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