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

Last change on this file since 3655 was 3522, checked in by suehring, 3 years ago

Sampling of variables extended in virtual measurement module

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 48.6 KB
Line 
1!> @virtual_measurement_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: virtual_measurement_mod.f90 3522 2018-11-13 12:14:36Z knoop $
27! Sampling of variables
28!
29! 3494 2018-11-06 14:51:27Z suehring
30! Bugfixing
31!
32! 3473 2018-10-30 20:50:15Z suehring
33! Initial revision
34!
35! Authors:
36! --------
37! @author Matthias Suehring
38!
39! Description:
40! ------------
41!> The module acts as an interface between 'real-world' observations and
42!> model simulations. Virtual measurements will be taken in the model at the
43!> coordinates representative for the 'real-world' measurement positions.
44!> More precisely, coordinates and measured quanties will be read from a
45!> NetCDF file which contains all required information. In the model,
46!> the same quantities (as long as all the required components are switched-on)
47!> will be sampled at the respective positions and output into an extra file,
48!> which allows for straight-forward comparison of model results with
49!> observations.
50!>
51!> @todo list_of_allowed variables needs careful checking
52!> @todo output (binary or NetCDF) needs to be implemented
53!> @todo clean-up anything from current test modus
54!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
55!>       the (UC)2 standard
56!> @note Fluxes are not processed
57!------------------------------------------------------------------------------!
58 MODULE virtual_measurement_mod
59
60    USE arrays_3d,                                                             &
61        ONLY:  q, pt, u, v, w, zu, zw
62
63    USE chem_modules,                                                          &
64        ONLY:  nspec
65
66    USE chemistry_model_mod,                                                   &
67        ONLY:  chem_species
68       
69    USE control_parameters,                                                    &
70        ONLY:  air_chemistry, dz, humidity, neutral, message_string,           &
71               virtual_measurement
72
73    USE cpulog,                                                                &
74        ONLY:  cpu_log, log_point
75       
76    USE grid_variables,                                                        &
77        ONLY:  dx, dy
78
79    USE indices,                                                               &
80        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny
81
82    USE kinds
83
84
85    IMPLICIT NONE
86
87    TYPE virt_mea
88   
89       CHARACTER(LEN=100)  ::  feature_type  !< type of the measurement
90       CHARACTER(LEN=100)  ::  site          !< name of the measurement site
91   
92       CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE ::  measured_vars_name !< name of the measured variables
93   
94       INTEGER(iwp) ::  ns = 0 !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
95       INTEGER(iwp) ::  ntraj  !< number of trajectories of a measurement
96       INTEGER(iwp) ::  nvar   !< number of measured variables
97       
98       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory or station that are no _FillValues
99       
100       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i  !< grid index for measurement position in x-direction
101       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j  !< grid index for measurement position in y-direction
102       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k  !< grid index for measurement position in k-direction
103           
104       LOGICAL ::  trajectory         = .FALSE. !< flag indicating that the observation is a mobile observation
105       LOGICAL ::  timseries          = .FALSE. !< flag indicating that the observation is a stationary point measurement
106       LOGICAL ::  timseries_profile  = .FALSE. !< flag indicating that the observation is a stationary profile measurement
107       
108       REAL(wp) ::  fill_eutm                         !< fill value for UTM coordinates in case of missing values
109       REAL(wp) ::  fill_nutm                         !< fill value for UTM coordinates in case of missing values
110       REAL(wp) ::  fill_zag                          !< fill value for heigth coordinates in case of missing values
111       REAL(wp) ::  fillout = -999.9                  !< fill value for output in case a observation is taken from inside a building
112       REAL(wp) ::  origin_x_obs                      !< origin of the observation in UTM coordiates in x-direction
113       REAL(wp) ::  origin_y_obs                      !< origin of the observation in UTM coordiates in y-direction
114             
115       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars  !< measured variables
116       
117    END TYPE virt_mea
118
119    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"                      !< dimension name for UTM coordinate easting
120    CHARACTER(LEN=11) ::  char_feature = "featureType"             !< attribute name for feature type
121    CHARACTER(LEN=10) ::  char_fillvalue = "_FillValue"            !< variable attribute name for _FillValue
122    CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
123    CHARACTER(LEN=5)  ::  char_nutm = "N_UTM"                      !< dimension name for UTM coordinate northing
124    CHARACTER(LEN=18) ::  char_numstations = "number_of_stations"  !< attribute name for number of stations
125    CHARACTER(LEN=8)  ::  char_origx = "origin_x"                  !< attribute name for station coordinate in x
126    CHARACTER(LEN=8)  ::  char_origy = "origin_y"                  !< attribute name for station coordinate in y
127    CHARACTER(LEN=4)  ::  char_site = "site"                       !< attribute name for site name
128    CHARACTER(LEN=19) ::  char_zag = "height_above_ground"         !< attribute name for height above ground variable
129    CHARACTER(LEN=10) ::  type_ts   = 'timeSeries'                 !< name of stationary point measurements
130    CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
131    CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
132!
133!-- MS: List requires careful revision!
134    CHARACTER(LEN=10), DIMENSION(1:47), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
135       (/ 'hfls      ',  & ! surface latent heat flux (W/m2)
136          'hfss      ',  & ! surface sensible heat flux (W/m2)
137          'hur       ',  & ! relative humidity (-)
138          'hus       ',  & ! specific humidity (g/kg)
139          'haa       ',  & ! absolute atmospheric humidity (kg/m3)
140          'mcpm1     ',  & ! mass concentration of PM1 (kg/m3)
141          'mcpm2p5   ',  & ! mass concentration of PM2.5 (kg/m3)
142          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
143          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
144          'mcco      ',  & ! mass concentration of CO (kg/m3)
145          'mcco2     ',  & ! mass concentration of CO2 (kg/m3)
146          'mcbcda    ',  & ! mass concentration of black carbon paritcles (kg/m3)
147          'ncaa      ',  & ! number concentation of particles (1/m3)
148          'mfco      ',  & ! mole fraction of CO (mol/mol)
149          'mfco2     ',  & ! mole fraction of CO2 (mol/mol)
150          'mfch4     ',  & ! mole fraction of methane (mol/mol)
151          'mfnh3     ',  & ! mole fraction of amonia (mol/mol)
152          'mfno      ',  & ! mole fraction of nitrogen monoxide (mol/mol)
153          'mfno2     ',  & ! mole fraction of nitrogen dioxide (mol/mol)
154          'mfso2     ',  & ! mole fraction of sulfur dioxide (mol/mol)
155          'mfh20     ',  & ! mole fraction of water (mol/mol)
156          'plev      ',  & ! ? air pressure - hydrostaic + perturbation?
157          'rlds      ',  & ! surface downward longwave flux  (W/m2)
158          'rlus      ',  & ! surface upward longwave flux (W/m2)
159          'rsds      ',  & ! surface downward shortwave flux (W/m2)
160          'rsus      ',  & ! surface upward shortwave flux (W/m2)
161          'ta        ',  & ! air temperature (degree C)
162          't_va      ',  & ! virtual accoustic temperature (K)
163          'theta     ',  & ! potential temperature (K)
164          'tro3      ',  & ! mole fraction of ozone air (mol/mol)
165          'ts        ',  & ! scaling parameter of temperature (K)
166          'wspeed    ',  & ! ? wind speed - horizontal?
167          'wdir      ',  & ! wind direction
168          'us        ',  & ! friction velocity
169          'msoil     ',  & ! ? soil moisture - which depth? 
170          'tsoil     ',  & ! ? soil temperature - which depth?                                                               
171          'u         ',  & ! u-component
172          'ua        ',  & ! eastward wind (is there any difference to u?)
173          'v         ',  & ! v-component
174          'va        ',  & ! northward wind (is there any difference to v?)
175          'w         ',  & ! w-component
176          'rld       ',  & ! downward longwave radiative flux (W/m2)
177          'rlu       ',  & ! upnward longwave radiative flux (W/m2)
178          'rsd       ',  & ! downward shortwave radiative flux (W/m2)
179          'rsu       ',  & ! upward shortwave radiative flux (W/m2)
180          'rsddif    ',  & ! downward shortwave diffuse radiative flux (W/m2)
181          'rnds      '   & ! surface net downward radiative flux (W/m2)
182       /)
183                                                             
184    INTEGER(iwp) ::  id_vm                           !< NetCDF file id for virtual measurements
185    INTEGER(iwp) ::  nvm = 0                         !< number of virtual measurements
186    INTEGER(iwp) ::  observation_coverage_xy = 0     !< horizontal distance from the measurement point where observations should be taken in the surrounding
187    INTEGER(iwp) ::  observation_coverage_z  = 0     !< vertical distance from the measurement point where observations should be taken in the surrounding
188   
189    LOGICAL ::  use_virtual_measurement = .FALSE. !< Namelist parameter
190    LOGICAL ::  global_attribute = .TRUE.         !< flag indicating a global attribute
191   
192    REAL(wp) ::  vm_time_start = 0.0              !< time after virtual measurements should start
193   
194
195    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea !< virtual measurement data structure
196   
197   
198    INTERFACE vm_check_parameters
199       MODULE PROCEDURE vm_check_parameters
200    END INTERFACE vm_check_parameters
201   
202    INTERFACE vm_init
203       MODULE PROCEDURE vm_init
204    END INTERFACE vm_init
205   
206    INTERFACE vm_parin
207       MODULE PROCEDURE vm_parin
208    END INTERFACE vm_parin
209   
210    INTERFACE vm_sampling
211       MODULE PROCEDURE vm_sampling
212    END INTERFACE vm_sampling
213
214    SAVE
215
216    PRIVATE
217
218!
219!-- Public interfaces
220    PUBLIC  vm_check_parameters, vm_init, vm_parin, vm_sampling
221
222!
223!-- Public variables
224    PUBLIC  vmea, vm_time_start
225
226 CONTAINS
227
228
229!------------------------------------------------------------------------------!
230! Description:
231! ------------
232!> Check parameters for virtual measurement module
233!------------------------------------------------------------------------------!
234 SUBROUTINE vm_check_parameters
235
236    USE control_parameters,                                                    &
237        ONLY:  message_string, virtual_measurement
238 
239    USE netcdf_data_input_mod,                                                 &
240        ONLY:  input_pids_static
241       
242    IMPLICIT NONE
243   
244!
245!-- In case virtual measurements are taken, a static input file is required.
246!-- This is because UTM coordinates for the PALM domain origin are required
247!-- for correct mapping of the measurements.
248!-- ToDo: Revise this later and remove this requirement.
249    IF ( virtual_measurement  .AND.  .NOT. input_pids_static )  THEN
250       message_string = 'If virtual measurements are taken a static input ' // &
251                        'file is mandatory.'
252       CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
253    ENDIF
254 
255 END SUBROUTINE vm_check_parameters
256 
257!------------------------------------------------------------------------------!
258! Description:
259! ------------
260!> Read namelist for the virtual measurement module
261!------------------------------------------------------------------------------!
262 SUBROUTINE vm_parin
263 
264    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
265 
266    NAMELIST /virtual_measurement_parameters/  use_virtual_measurement,        &
267                                               vm_time_start
268
269    line = ' '
270
271!
272!-- Try to find stg package
273    REWIND ( 11 )
274    line = ' '
275    DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
276       READ ( 11, '(A)', END=20 )  line
277    ENDDO
278    BACKSPACE ( 11 )
279
280!
281!-- Read namelist
282    READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
283
284!
285!-- Set flag that indicates that the virtual measurement module is switched on
286    IF ( use_virtual_measurement )  virtual_measurement = .TRUE.
287   
288    GOTO 20
289
290 10 BACKSPACE( 11 )
291    READ( 11 , '(A)') line
292    CALL parin_fail_message( 'virtual_measurement_parameters', line )
293
294 20 CONTINUE
295 
296 END SUBROUTINE vm_parin
297
298
299!------------------------------------------------------------------------------!
300! Description:
301! ------------
302!> Initialize virtual measurements: read coordiante arrays and measured
303!> variables, set indicies indicating the measurement points, read further
304!> attributes, etc..
305!------------------------------------------------------------------------------!
306 SUBROUTINE vm_init
307
308    USE arrays_3d,                                                             &
309        ONLY:  zu, zw
310       
311    USE grid_variables,                                                        &
312        ONLY:  ddx, ddy, dx, dy
313       
314    USE indices,                                                               &
315        ONLY:  nxl, nxr, nyn, nys
316 
317    USE netcdf_data_input_mod,                                                 &
318        ONLY:  init_model, input_file_vm,                                      &
319               netcdf_data_input_get_dimension_length,                         &
320               netcdf_data_input_att, netcdf_data_input_var
321               
322    USE surface_mod,                                                           &
323        ONLY:  get_topography_top_index_ji
324       
325    IMPLICIT NONE
326   
327    CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
328    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
329    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed   
330   
331    INTEGER(iwp) ::  dim_eutm  !< dimension size of UTM easting coordinate
332    INTEGER(iwp) ::  dim_nutm  !< dimension size of UTM northing coordinate
333    INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate
334    INTEGER(iwp) ::  dim_zag   !< dimension size of height coordinate
335    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
336    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
337    INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
338    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
339    INTEGER(iwp) ::  js        !< grid index of real observation point of the respective station in y-direction
340    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
341    INTEGER(iwp) ::  kl        !< lower vertical index of surrounding grid points of an observation coordinate
342    INTEGER(iwp) ::  ks        !< grid index of real observation point of the respective station in z-direction
343    INTEGER(iwp) ::  ksurf     !< topography top index
344    INTEGER(iwp) ::  ku        !< upper vertical index of surrounding grid points of an observation coordinate
345    INTEGER(iwp) ::  l         !< running index over all stations
346    INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
347    INTEGER(iwp) ::  ll        !< running index over all measured variables in file
348    INTEGER(iwp) ::  lll       !< running index over all allowed variables
349    INTEGER(iwp) ::  n         !< running index over trajectory coordinates
350    INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
351    INTEGER(iwp) ::  t         !< running index over number of trajectories
352   
353    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
354   
355    LOGICAL ::  chem_include !< flag indicating that chemical species is considered in modelled mechanism
356    LOGICAL ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
357   
358    REAL(wp)     ::  fill_eutm !< _FillValue for coordinate array E_UTM
359    REAL(wp)     ::  fill_nutm !< _FillValue for coordinate array N_UTM
360    REAL(wp)     ::  fill_zag  !< _FillValue for height coordinate
361   
362    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm !< easting UTM coordinate, temporary variable
363    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm !< northing UTM coordinate, temporary variable,
364    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  z_ag  !< height coordinate relative to origin_z, temporary variable
365!
366!-- Obtain number of virtual measurement stations
367    CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm,   &
368                                global_attribute, 'open', '' )
369                               
370!     write(9,*) "num stationi", nvm
371!     flush(9)
372!
373!-- ALLOCATE data structure
374    ALLOCATE( vmea(1:nvm) )
375!
376!-- Allocate flag array
377    ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
378    meas_flag = 0
379 
380!
381!-- Read station coordinates and further attributes.
382!-- Note all coordinates are in UTM coordinates.
383    DO  l = 1, nvm
384!
385!--    Determine suffix which contains the ID
386       IF( l < 10 )  THEN
387          WRITE( dum, '(I1)')  l
388       ELSEIF( l < 100 )  THEN
389          WRITE( dum, '(I2)')  l
390       ELSEIF( l < 1000 )  THEN
391          WRITE( dum, '(I3)')  l
392       ELSEIF( l < 10000 )  THEN
393          WRITE( dum, '(I4)')  l
394       ELSEIF( l < 100000 )  THEN
395          WRITE( dum, '(I5)')  l
396       ENDIF
397       
398       CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx            &
399                                   // TRIM( dum ), id_vm, '', global_attribute,&
400                                   '', '' )
401       CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy            &
402                                   // TRIM( dum ), id_vm, '', global_attribute,&
403                                   '', '' )
404       CALL netcdf_data_input_att( vmea(l)%site,         char_site             &
405                                   // TRIM( dum ), id_vm, '', global_attribute,&
406                                   '', '' )
407       CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature          &
408                                   // TRIM( dum ), id_vm, '', global_attribute,&
409                                   '', '' )
410       
411!
412!---   Set logicals depending on the type of the measurement
413       IF ( INDEX( vmea(l)%feature_type, type_tspr     ) /= 0 )  THEN
414          vmea(l)%timseries_profile = .TRUE.
415       ELSEIF ( INDEX( vmea(l)%feature_type, type_ts   ) /= 0 )  THEN
416          vmea(l)%timseries         = .TRUE.
417       ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 )  THEN
418          vmea(l)%trajectory        = .TRUE.
419       ELSE
420!
421!--       Give error message
422          message_string = 'Attribue featureType = ' //                        &
423                           TRIM( vmea(l)%feature_type ) //                     &
424                           ' is not allowed.' 
425          CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 )
426       ENDIF
427!
428!--    Read string with all measured variables at this station
429       measured_variables_file = ''
430       CALL netcdf_data_input_var( measured_variables_file,                    &
431                                   char_mv // TRIM( dum ), id_vm )
432!
433!--    Count the number of measured variables which match with the variables
434!--    which are allowed to be measured in PALM. Please note, for some
435!--    NetCDF interal reasons characters end with a NULL, i.e. also empty
436!--    characters contain a NULL. Therefore, check the strings for a Null to
437!--    get the correct character length in order to compare them with the list
438!--    of allowed variables.
439       vmea(l)%nvar = 0
440       DO ll = 1, SIZE( measured_variables_file )
441          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
442               measured_variables_file(ll)(1:1) /= ' ')  THEN
443!
444!--          Obtain character length of the character
445             len_char = 1
446             DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
447                 .AND.  measured_variables_file(ll)(len_char:len_char) /= ' ' )
448                len_char = len_char + 1
449             ENDDO
450             len_char = len_char - 1
451!
452!--          Now, compare the measured variable with the list of allowed
453!--          variables.
454             DO  lll= 1, SIZE( list_allowed_variables )
455                IF ( measured_variables_file(ll)(1:len_char) ==                &
456                     TRIM( list_allowed_variables(lll) ) )  THEN
457                   vmea(l)%nvar = vmea(l)%nvar + 1
458                   measured_variables(vmea(l)%nvar) =                          &
459                                       measured_variables_file(ll)(1:len_char)
460                ENDIF
461             ENDDO
462          ENDIF
463       ENDDO
464!
465!--    Allocate array for the measured variables names for the station l.
466       ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
467
468       DO  ll = 1, vmea(l)%nvar
469          vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
470       ENDDO
471!
472!--    In case of chemistry, check if species is considered in the modelled
473!--    chemistry mechanism.
474       IF ( air_chemistry )  THEN
475          DO  ll = 1, vmea(l)%nvar
476             chem_include = .FALSE.
477             DO  n = 1, nspec
478                IF ( TRIM( vmea(l)%measured_vars_name(ll) ) ==                 &
479                     TRIM( chem_species(n)%name ) )  chem_include = .TRUE.
480             ENDDO
481             IF ( .NOT. chem_include )  THEN
482                message_string = TRIM( vmea(l)%measured_vars_name(ll) ) //     &
483                                 ' is not considered in the modelled '  //     &
484                                 'chemistry mechanism'
485                CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
486             ENDIF
487          ENDDO
488       ENDIF
489!
490!--    For the actual measurement ID read the UTM coordinates. Based on these,
491!--    define the index space on each subdomain where measurements should be
492!--    taken. Note, the entire coordinate arrays will not be stored on data
493!--    type as this would exceed memory requirements, particularly for
494!--    trajectory measurements. If no variable will be virtually measured,
495!--    skip the reading.
496       IF ( vmea(l)%nvar > 0 )  THEN
497!
498!--       For stationary measurements UTM coordinates are just one value and
499!--       its dimension is "station", while for mobile measurements UTM
500!--       coordinates are arrays. First, inquire dimension length for
501!--       UTM coordinates.
502          IF ( vmea(l)%trajectory )  THEN
503!
504!--          For non-stationary measurements read the number of trajectories
505             CALL netcdf_data_input_get_dimension_length( id_vm,              &
506                                                          vmea(l)%ntraj,      &
507                                                          "traj" //           &
508                                                          TRIM( dum ) )
509             CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime,   &
510                                                          "ntime" //          &
511                                                          TRIM( dum ) )
512!
513!--       For stationary measurements the dimension for UTM coordinates is 1
514          ELSE
515             vmea(l)%ntraj  = 1
516             dim_ntime = 1
517          ENDIF
518         
519!
520!-        Allocate array which defines individual time frame for each
521!--       trajectory or station
522          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
523!
524!--       Allocate temporary arrays for UTM and height coordinates. Note,
525!--       on file UTM coordinates might be 1D or 2D variables
526          ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
527          ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
528          ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime)  )
529!
530!--       Read _FillValue attributes
531          CALL netcdf_data_input_att( fill_eutm, char_fillvalue,               &
532                                      id_vm, '', .NOT. global_attribute, '',   &
533                                      char_eutm // TRIM( dum ) )
534          CALL netcdf_data_input_att( fill_nutm, char_fillvalue,               &
535                                      id_vm, '', .NOT. global_attribute, '',   &
536                                      char_nutm // TRIM( dum ) )
537          CALL netcdf_data_input_att( fill_zag, char_fillvalue,                &
538                                      id_vm, '', .NOT. global_attribute, '',   &
539                                      char_zag  // TRIM( dum ) )
540!
541!--       Read UTM and height coordinates coordinates for all trajectories and
542!--       times.
543          IF ( vmea(l)%trajectory )  THEN
544             CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm,  &
545                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
546             CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm,  &
547                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
548             CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), id_vm,  &
549                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
550          ELSE
551             CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm )
552             CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm )
553             CALL netcdf_data_input_var( z_ag(1,:),  char_zag  // TRIM( dum ), id_vm )
554          ENDIF
555!
556!-- For testing:
557          e_utm = e_utm - e_utm(1,1)
558          n_utm = n_utm - n_utm(1,1)
559         
560!             
561!--       First, compute relative x- and y-coordinates with respect to the
562!--       lower-left origin of the model domain, which is the difference
563!--       betwen UTM coordinates. 
564!           e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t))            &
565!                                       - init_model%origin_x
566!           n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t))            &
567!                                       - init_model%origin_y
568
569!
570!--       Based on UTM coordinates, check if the measurement station or parts
571!--       of the trajectory is on subdomain. This case, setup grid index space
572!--       sample these quantities.
573          meas_flag = 0
574          DO  t = 1, vmea(l)%ntraj
575!
576!--          Determine the individual time coordinate length for each station and
577!--          trajectory. This is required as several stations and trajectories
578!--          are merged into one file but they do not have the same number of
579!--          points in time, hence, missing values may occur and cannot be
580!--          processed further.
581             vmea(l)%dim_t(t) = 0
582             DO  n = 1, dim_ntime
583                IF ( e_utm(t,n) /= fill_eutm  .AND.                            &
584                     n_utm(t,n) /= fill_nutm  .AND.                            &
585                     z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
586             ENDDO
587
588!
589!--          Compute grid indices relative to origin and check if these are
590!--          on the subdomain. Note, virtual measurements will be taken also
591!--          at grid points surrounding the station, hence, check also for
592!--          these grid points.
593             DO  n = 1, vmea(l)%dim_t(t)
594                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
595                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
596!
597!--             Is the observation point on subdomain?
598                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.                   &
599                          js >= nys  .AND.  js <= nyn )
600!
601!--             Check if observation coordinate is on subdomain
602                IF ( on_pe )  THEN
603!
604!--                Determine vertical index which correspond to the observation
605!--                height.
606                   ksurf = get_topography_top_index_ji( js, is, 's' )
607                   ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
608!
609!--                Set mask array at the observation coordinates. Also, flag the
610!--                surrounding coordinate points, but first check whether the
611!--                surrounding coordinate points are on the subdomain.
612                   kl = MERGE( ks-1, ks, ks-1 >= nzb  .AND. ks-1 > ksurf )
613                   ku = MERGE( ks+1, ks, ks+1 <= nzt+1 )
614                 
615                   meas_flag(kl:ku,js-1:js+1,is-1:is+1) =                      &
616                               IBSET( meas_flag(kl:ku,js-1:js+1,is-1:is+1), 0 ) 
617                ENDIF
618             ENDDO
619             
620          ENDDO
621!
622!--       Based on the flag array count the number of observation points.
623          ns = 0
624          DO  is = nxl-1, nxr+1
625             DO  js = nys-1, nyn+1
626                DO  ks = nzb, nzt+1
627                   ns = ns + MERGE( 1, 0, BTEST( meas_flag(ks,js,is), 0 ) )
628                ENDDO
629             ENDDO
630          ENDDO
631!
632!--       Store number of observation points on subdomain and allocate index
633!--       arrays.
634          vmea(l)%ns = ns
635          ns = 0
636         
637          ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
638          ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
639          ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
640!
641!--       Based on the flag array store the grid indices which correspond to
642!--       the observation coordinates.
643          DO  is = nxl-1, nxr+1
644             DO  js = nys-1, nyn+1
645                DO  ks = nzb, nzt+1
646                   IF ( BTEST( meas_flag(ks,js,is), 0 ) )  THEN
647                      ns = ns + 1
648                      vmea(l)%i(ns) = is
649                      vmea(l)%j(ns) = js
650                      vmea(l)%k(ns) = ks
651                   ENDIF
652                ENDDO
653             ENDDO
654          ENDDO
655         
656!           write(9,*) l, "NS: ", ns
657!           IF ( ns > 0 )   THEN
658!              write(9,*) "i ", vmea(l)%i
659!              write(9,*) "j ", vmea(l)%j
660!              write(9,*) "k ", vmea(l)%k
661!              write(9,*)
662!           ENDIF
663!
664!--       Allocate array to save the sampled values.
665!--       Todo: Is it better to allocate for all variables at a station
666!--       and store all the values before writing, or sample the variables
667!--       directly in the data output?
668          ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) )
669!
670!--       Initialize with _FillValue
671          vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) = vmea(l)%fillout
672!
673!--       Deallocate temporary coordinate arrays
674          IF ( ALLOCATED( e_utm ) )  DEALLOCATE( e_utm )
675          IF ( ALLOCATED( n_utm ) )  DEALLOCATE( n_utm )
676          IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( z_ag  )
677       ENDIF
678    ENDDO
679!     flush(9)
680   
681!
682!-- Close input file for virtual measurements. Therefore, just call
683!-- the read attribute routine with the "close" option.
684    CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '',              &
685                                global_attribute, 'close', '' )
686!                               
687!-- Dellocate flag array
688    DEALLOCATE( meas_flag )
689       
690  END SUBROUTINE vm_init
691 
692 
693!------------------------------------------------------------------------------!
694! Description:
695! ------------
696!> Sampling of the actual quantities along the observation coordinates
697!------------------------------------------------------------------------------!
698  SUBROUTINE vm_sampling
699
700    USE arrays_3d,                                                             &
701        ONLY:  exner, pt, q, u, v, w
702
703    USE basic_constants_and_equations_mod,                                     &
704        ONLY:  pi
705   
706    USE radiation_model_mod,                                                   &
707        ONLY:  radiation 
708
709    USE surface_mod,                                                           &
710        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
711   
712     IMPLICIT NONE
713     
714     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
715     
716     INTEGER(iwp) ::  i  !< grid index in x-direction
717     INTEGER(iwp) ::  j  !< grid index in y-direction
718     INTEGER(iwp) ::  k  !< grid index in z-direction
719     INTEGER(iwp) ::  l  !< running index over the number of stations
720     INTEGER(iwp) ::  m  !< running index over all virtual observation coordinates
721     INTEGER(iwp) ::  mm !< index of surface element which corresponds to the virtual observation coordinate
722     INTEGER(iwp) ::  n  !< running index over all measured variables at a station
723     INTEGER(iwp) ::  nn !< running index over the number of chemcal species
724!
725!--  Loop over all stations. For each possible variable loop over all
726!--  observation points
727     DO  l = 1, nvm
728!
729!--     Loop over all measured variables at this station. Please note,
730!--     velocity components are interpolated onto scalar grid. 
731        DO  n = 1, vmea(l)%nvar
732       
733           SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
734           
735              CASE ( 'theta' )
736                 IF ( .NOT. neutral )  THEN
737                    DO  m = 1, vmea(l)%ns
738                       k = vmea(l)%k(m)
739                       j = vmea(l)%j(m)
740                       i = vmea(l)%i(m)
741                       vmea(l)%measured_vars(n,m) = pt(k,j,i)
742                    ENDDO
743                 ENDIF
744                 
745              CASE ( 'ta', 't_va' )
746                 IF ( .NOT. neutral )  THEN
747                    DO  m = 1, vmea(l)%ns
748                       k = vmea(l)%k(m)
749                       j = vmea(l)%j(m)
750                       i = vmea(l)%i(m)
751                       vmea(l)%measured_vars(n,m) = pt(k,j,i) * exner( k )
752                    ENDDO
753                 ENDIF
754                 
755              CASE ( 'hus', 'haa' )
756                 IF ( humidity )  THEN
757                    DO  m = 1, vmea(l)%ns
758                       k = vmea(l)%k(m)
759                       j = vmea(l)%j(m)
760                       i = vmea(l)%i(m)
761                       vmea(l)%measured_vars(n,m) = q(k,j,i)
762                    ENDDO
763                 ENDIF
764                 
765              CASE ( 'u', 'ua' )
766                 DO  m = 1, vmea(l)%ns
767                    k = vmea(l)%k(m)
768                    j = vmea(l)%j(m)
769                    i = vmea(l)%i(m)
770                    vmea(l)%measured_vars(n,m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
771                 ENDDO
772                 
773              CASE ( 'v', 'va' )
774                 DO  m = 1, vmea(l)%ns
775                    k = vmea(l)%k(m)
776                    j = vmea(l)%j(m)
777                    i = vmea(l)%i(m)
778                    vmea(l)%measured_vars(n,m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
779                 ENDDO
780                 
781              CASE ( 'w' )
782                 DO  m = 1, vmea(l)%ns
783                    k = vmea(l)%k(m)
784                    j = vmea(l)%j(m)
785                    i = vmea(l)%i(m)
786                    vmea(l)%measured_vars(n,m) = w(k,j,i) !0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
787                 ENDDO
788                 
789              CASE ( 'wspeed' )
790                 DO  m = 1, vmea(l)%ns
791                    k = vmea(l)%k(m)
792                    j = vmea(l)%j(m)
793                    i = vmea(l)%i(m)
794                    vmea(l)%measured_vars(n,m) = SQRT(                         &
795                                   ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
796                                   ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2   &
797                                                     )
798                 ENDDO
799                 
800              CASE ( 'wdir' )
801                 DO  m = 1, vmea(l)%ns
802                    k = vmea(l)%k(m)
803                    j = vmea(l)%j(m)
804                    i = vmea(l)%i(m)
805                   
806                    vmea(l)%measured_vars(n,m) = ATAN2(                        &
807                                       - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ),   &
808                                       - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )    &
809                                                      ) * 180.0_wp / pi
810                 ENDDO
811!
812!--           MS: list of variables needs extension.
813              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mcco', 'mcco2', 'mcbcda',  &
814                     'ncaa', 'mfco2', 'mfco', 'mfch4', 'mfnh3', 'mfno',        &
815                     'mfno2', 'mfso2', 'mfh20', 'tr03' )
816                 IF ( air_chemistry )  THEN                 
817                    DO  nn = 1, nspec                   
818                       IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
819                            TRIM( chem_species(nn)%name ) )  THEN                           
820                          DO  m = 1, vmea(l)%ns             
821                             k = vmea(l)%k(m)
822                             j = vmea(l)%j(m)
823                             i = vmea(l)%i(m)                   
824                             vmea(l)%measured_vars(n,m) =                      &
825                                                   chem_species(nn)%conc(k,j,i)
826                          ENDDO
827                       ENDIF
828                    ENDDO
829                 ENDIF
830                 
831              CASE ( 'us' )
832                 DO  m = 1, vmea(l)%ns
833!
834!--                 Surface data is only available on inner subdomains, not
835!--                 on ghost points. Hence, limit the indices.
836                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
837                    j = MERGE( j           , nyn, j            > nyn )
838                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
839                    i = MERGE( i           , nxr, i            > nxr )
840                   
841                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
842                             surf_def_h(0)%end_index(j,i)
843                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%us(mm)
844                    ENDDO
845                    DO  mm = surf_lsm_h%start_index(j,i),                      &
846                             surf_lsm_h%end_index(j,i)
847                       vmea(l)%measured_vars(n,m) = surf_lsm_h%us(mm)
848                    ENDDO
849                    DO  mm = surf_usm_h%start_index(j,i),                      &
850                             surf_usm_h%end_index(j,i)
851                       vmea(l)%measured_vars(n,m) = surf_usm_h%us(mm)
852                    ENDDO
853                 ENDDO
854                 
855              CASE ( 'ts' )
856                 DO  m = 1, vmea(l)%ns
857!
858!--                 Surface data is only available on inner subdomains, not
859!--                 on ghost points. Hence, limit the indices.
860                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
861                    j = MERGE( j           , nyn, j            > nyn )
862                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
863                    i = MERGE( i           , nxr, i            > nxr )
864                   
865                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
866                             surf_def_h(0)%end_index(j,i)
867                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%ts(mm)
868                    ENDDO
869                    DO  mm = surf_lsm_h%start_index(j,i),                      &
870                             surf_lsm_h%end_index(j,i)
871                       vmea(l)%measured_vars(n,m) = surf_lsm_h%ts(mm)
872                    ENDDO
873                    DO  mm = surf_usm_h%start_index(j,i),                      &
874                             surf_usm_h%end_index(j,i)
875                       vmea(l)%measured_vars(n,m) = surf_usm_h%ts(mm)
876                    ENDDO
877                 ENDDO
878                 
879              CASE ( 'hfls' )
880                 DO  m = 1, vmea(l)%ns
881!
882!--                 Surface data is only available on inner subdomains, not
883!--                 on ghost points. Hence, limit the indices.
884                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
885                    j = MERGE( j           , nyn, j            > nyn )
886                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
887                    i = MERGE( i           , nxr, i            > nxr )
888                   
889                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
890                             surf_def_h(0)%end_index(j,i)
891                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%qsws(mm)
892                    ENDDO
893                    DO  mm = surf_lsm_h%start_index(j,i),                      &
894                             surf_lsm_h%end_index(j,i)
895                       vmea(l)%measured_vars(n,m) = surf_lsm_h%qsws(mm)
896                    ENDDO
897                    DO  mm = surf_usm_h%start_index(j,i),                      &
898                             surf_usm_h%end_index(j,i)
899                       vmea(l)%measured_vars(n,m) = surf_usm_h%qsws(mm)
900                    ENDDO
901                 ENDDO
902                 
903              CASE ( 'hfss' )
904                 DO  m = 1, vmea(l)%ns
905!
906!--                 Surface data is only available on inner subdomains, not
907!--                 on ghost points. Hence, limit the indices.
908                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
909                    j = MERGE( j           , nyn, j            > nyn )
910                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
911                    i = MERGE( i           , nxr, i            > nxr )
912                   
913                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
914                             surf_def_h(0)%end_index(j,i)
915                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%shf(mm)
916                    ENDDO
917                    DO  mm = surf_lsm_h%start_index(j,i),                      &
918                             surf_lsm_h%end_index(j,i)
919                       vmea(l)%measured_vars(n,m) = surf_lsm_h%shf(mm)
920                    ENDDO
921                    DO  mm = surf_usm_h%start_index(j,i),                      &
922                             surf_usm_h%end_index(j,i)
923                       vmea(l)%measured_vars(n,m) = surf_usm_h%shf(mm)
924                    ENDDO
925                 ENDDO
926                 
927              CASE ( 'rnds' )
928                 IF ( radiation )  THEN
929                    DO  m = 1, vmea(l)%ns
930!
931!--                    Surface data is only available on inner subdomains, not
932!--                    on ghost points. Hence, limit the indices.
933                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
934                       j = MERGE( j           , nyn, j            > nyn )
935                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
936                       i = MERGE( i           , nxr, i            > nxr )
937                   
938                       DO  mm = surf_lsm_h%start_index(j,i),                   &
939                                surf_lsm_h%end_index(j,i)
940                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_net(mm)
941                       ENDDO
942                       DO  mm = surf_usm_h%start_index(j,i),                   &
943                                surf_usm_h%end_index(j,i)
944                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_net(mm)
945                       ENDDO
946                    ENDDO
947                 ENDIF
948                 
949              CASE ( 'rsus', 'rsu' )
950                 IF ( radiation )  THEN
951                    DO  m = 1, vmea(l)%ns
952!
953!--                    Surface data is only available on inner subdomains, not
954!--                    on ghost points. Hence, limit the indices.
955                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
956                       j = MERGE( j           , nyn, j            > nyn )
957                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
958                       i = MERGE( i           , nxr, i            > nxr )
959                   
960                       DO  mm = surf_lsm_h%start_index(j,i),                   &
961                                surf_lsm_h%end_index(j,i)
962                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_out(mm)
963                       ENDDO
964                       DO  mm = surf_usm_h%start_index(j,i),                   &
965                                surf_usm_h%end_index(j,i)
966                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_out(mm)
967                       ENDDO
968                    ENDDO
969                 ENDIF
970                 
971              CASE ( 'rsds', 'rsd' )
972                 IF ( radiation )  THEN
973                    DO  m = 1, vmea(l)%ns
974!
975!--                    Surface data is only available on inner subdomains, not
976!--                    on ghost points. Hence, limit the indices.
977                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
978                       j = MERGE( j           , nyn, j            > nyn )
979                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
980                       i = MERGE( i           , nxr, i            > nxr )
981                   
982                       DO  mm = surf_lsm_h%start_index(j,i),                   &
983                                surf_lsm_h%end_index(j,i)
984                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_in(mm)
985                       ENDDO
986                       DO  mm = surf_usm_h%start_index(j,i),                   &
987                                surf_usm_h%end_index(j,i)
988                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_in(mm)
989                       ENDDO
990                    ENDDO
991                 ENDIF
992                 
993              CASE ( 'rlus', 'rlu' )
994                 IF ( radiation )  THEN
995                    DO  m = 1, vmea(l)%ns
996!
997!--                    Surface data is only available on inner subdomains, not
998!--                    on ghost points. Hence, limit the indices.
999                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1000                       j = MERGE( j           , nyn, j            > nyn )
1001                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1002                       i = MERGE( i           , nxr, i            > nxr )
1003                   
1004                       DO  mm = surf_lsm_h%start_index(j,i),                   &
1005                                surf_lsm_h%end_index(j,i)
1006                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_out(mm)
1007                       ENDDO
1008                       DO  mm = surf_usm_h%start_index(j,i),                   &
1009                                surf_usm_h%end_index(j,i)
1010                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_out(mm)
1011                       ENDDO
1012                    ENDDO
1013                 ENDIF
1014                 
1015              CASE ( 'rlds', 'rld' )
1016                 IF ( radiation )  THEN
1017                    DO  m = 1, vmea(l)%ns
1018!
1019!--                    Surface data is only available on inner subdomains, not
1020!--                    on ghost points. Hence, limit the indices.
1021                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1022                       j = MERGE( j           , nyn, j            > nyn )
1023                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1024                       i = MERGE( i           , nxr, i            > nxr )
1025                   
1026                       DO  mm = surf_lsm_h%start_index(j,i),                   &
1027                                surf_lsm_h%end_index(j,i)
1028                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_in(mm)
1029                       ENDDO
1030                       DO  mm = surf_usm_h%start_index(j,i),                   &
1031                                surf_usm_h%end_index(j,i)
1032                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_in(mm)
1033                       ENDDO
1034                    ENDDO
1035                 ENDIF
1036!
1037!--           More will follow ...
1038                 
1039           END SELECT
1040
1041        ENDDO
1042
1043     ENDDO
1044     
1045  END SUBROUTINE vm_sampling
1046 
1047
1048 END MODULE virtual_measurement_mod
Note: See TracBrowser for help on using the repository browser.