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

Last change on this file since 3913 was 3913, checked in by gronemeier, 2 years ago

Bugfix: rotate positions of measurements before writing them into file (virtual_measurement_mod)

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