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

Last change on this file since 3904 was 3904, checked in by gronemeier, 5 years ago

rotate coordinates of stations by given rotation_angle

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