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

Last change on this file since 3705 was 3705, checked in by suehring, 2 years ago

last commit documented

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