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

Last change on this file since 3690 was 3665, checked in by raasch, 5 years ago

dummy statements added to avoid compiler warnings about unused variables, unused variables removed, ssh-call for submitting batch jobs on remote systems modified again to avoid output of login messages on specific systems

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