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

Last change on this file since 3833 was 3833, checked in by forkel, 2 years ago

removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead

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