source: palm/tags/release-6.0/SOURCE/virtual_measurement_mod.f90 @ 3862

Last change on this file since 3862 was 3473, checked in by suehring, 5 years ago

Bugfix for previous commit

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 35.9 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 3473 2018-10-30 20:50:15Z banzhafs $
27! Initial revision
28!
29! 3472 2018-10-30 20:43:50Z suehring
30!
31! Authors:
32! --------
33! @author Matthias Suehring
34!
35!
36!
37! Description:
38! ------------
39!> The module acts as an interface between 'real-world' observations and
40!> model simulations. Virtual measurements will be taken in the model at the
41!> coordinates representative for the 'real-world' measurement positions.
42!> More precisely, coordinates and measured quanties will be read from a
43!> NetCDF file which contains all required information. In the model,
44!> the same quantities (as long as all the required components are switched-on)
45!> will be sampled at the respective positions and output into an extra file,
46!> which allows for straight-forward comparison of model results with
47!> observations.
48!------------------------------------------------------------------------------!
49 MODULE virtual_measurement_mod
50
51
52    USE arrays_3d,                                                             &
53        ONLY:  q, pt, u, v, w, zu, zw
54
55    USE control_parameters,                                                    &
56        ONLY:  dz, message_string, virtual_measurement
57
58    USE cpulog,                                                                &
59        ONLY:  cpu_log, log_point
60       
61    USE grid_variables,                                                        &
62        ONLY:  dx, dy
63
64    USE indices,                                                               &
65        ONLY:  nzb, nzt, nxl, nxr, nys, nyn
66
67    USE kinds
68
69
70    IMPLICIT NONE
71
72    TYPE virt_mea
73   
74       CHARACTER(LEN=100)  ::  feature_type  !< type of the measurement
75       CHARACTER(LEN=100)  ::  site          !< name of the measurement site
76   
77       CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE ::  measured_vars_name !< name of the measured variables
78   
79       INTEGER(iwp) ::  ns    !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
80       INTEGER(iwp) ::  ntraj !< number of trajectories of a measurement
81       INTEGER(iwp) ::  nvar  !< number of measured variables
82       
83       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory or station that are no _FillValues
84       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ngp   !< number of grid points where observations for a site took place,
85                                                         !<individual for each trajectory or station that are no _FillValues
86       
87       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i  !< grid index for measurement position in x-direction
88       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j  !< grid index for measurement position in y-direction
89       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k  !< grid index for measurement position in k-direction
90           
91       LOGICAL ::  trajectory         = .FALSE. !< flag indicating that the observation is a mobile observation
92       LOGICAL ::  timseries          = .FALSE. !< flag indicating that the observation is a stationary point measurement
93       LOGICAL ::  timseries_profile  = .FALSE. !< flag indicating that the observation is a stationary profile measurement
94       
95       REAL(wp) ::  fill_eutm                         !< fill value for UTM coordinates in case of missing values
96       REAL(wp) ::  fill_nutm                         !< fill value for UTM coordinates in case of missing values
97       REAL(wp) ::  fill_zag                          !< fill value for heigth coordinates in case of missing values
98       REAL(wp) ::  fillout = -999.9                  !< fill value for output in case a observation is taken from inside a building
99       REAL(wp) ::  origin_x_obs                      !< origin of the observation in UTM coordiates in x-direction
100       REAL(wp) ::  origin_y_obs                      !< origin of the observation in UTM coordiates in y-direction
101       
102       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  xmea   !< measurement x-position in absolute UTM coordinates 
103       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ymea   !< measurement y-position in absolute UTM coordinates 
104       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zmea   !< measurement z-position in height above ground level
105       
106       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars  !< measured variables
107       
108    END TYPE virt_mea
109
110    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"                      !< dimension name for UTM coordinate easting
111    CHARACTER(LEN=11) ::  char_feature = "featureType"             !< attribute name for feature type
112    CHARACTER(LEN=10) ::  char_fillvalue = "_FillValue"            !< variable attribute name for _FillValue
113    CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
114    CHARACTER(LEN=5)  ::  char_nutm = "N_UTM"                      !< dimension name for UTM coordinate northing
115    CHARACTER(LEN=18) ::  char_numstations = "number_of_stations"  !< attribute name for number of stations
116    CHARACTER(LEN=8)  ::  char_origx = "origin_x"                  !< attribute name for station coordinate in x
117    CHARACTER(LEN=8)  ::  char_origy = "origin_y"                  !< attribute name for station coordinate in y
118    CHARACTER(LEN=4)  ::  char_site = "site"                       !< attribute name for site name
119    CHARACTER(LEN=19) ::  char_zag = "height_above_ground"         !< attribute name for height above ground variable
120    CHARACTER(LEN=10) ::  type_ts   = 'timeSeries'                 !< name of stationary point measurements
121    CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
122    CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
123   
124    CHARACTER(LEN=10), DIMENSION(1:53), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
125       (/ 'hfls      ',  & ! surface latent heat flux (W/m2)
126          'hfss      ',  & ! surface sensible heat flux (W/m2)
127          'hur       ',  & ! relative humidity (-)
128          'hus       ',  & ! specific humidity (g/kg)
129          'haa       ',  & ! absolute atmospheric humidity (kg/m3)
130          'mcpm1     ',  & ! mass concentration of PM1 (kg/m3)
131          'mcpm2p5   ',  & ! mass concentration of PM2.5 (kg/m3)
132          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
133          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
134          'mcco      ',  & ! mass concentration of CO (kg/m3)
135          'mcco2     ',  & ! mass concentration of CO2 (kg/m3)
136          'mcbcda    ',  & ! mass concentration of black carbon paritcles (kg/m3)
137          'ncaa      ',  & ! number concentation of particles (1/m3)
138          'mfco2     ',  & ! mole fraction of CO (mol/mol)
139          'mfco2     ',  & ! mole fraction of CO2 (mol/mol)
140          'mfch4     ',  & ! mole fraction of methane (mol/mol)
141          'mfnh3     ',  & ! mole fraction of amonia (mol/mol)
142          'mfno      ',  & ! mole fraction of nitrogen monoxide (mol/mol)
143          'mfno2     ',  & ! mole fraction of nitrogen dioxide (mol/mol)
144          'mfso2     ',  & ! mole fraction of sulfur dioxide (mol/mol)
145          'mfh20     ',  & ! mole fraction of water (mol/mol)
146          'plev      ',  & ! ? air pressure - hydrostaic + perturbation?
147          'rlds      ',  & ! surface downward longwave flux  (W/m2)
148          'rlus      ',  & ! surface upward longwave flux (W/m2)
149          'rsds      ',  & ! surface downward shortwave flux (W/m2)
150          'rsus      ',  & ! surface upward shortwave flux (W/m2)
151          'ta        ',  & ! air temperature (degree C)
152          't_va      ',  & ! virtual accoustic temperature (K)
153          'theta     ',  & ! potential temperature (K)
154          'tro3      ',  & ! mole fraction of ozone air (mol/mol)
155          'ts        ',  & ! scaling parameter of temperature (K)
156          'wspeed    ',  & ! ? wind speed - horizontal?
157          'wdir      ',  & ! wind direction
158          'us        ',  & ! friction velocity
159          'msoil     ',  & ! ? soil moisture - which depth? 
160          'tsoil     ',  & ! ? soil temperature - which depth?                                                               
161          'u         ',  & ! u-component
162          'ua        ',  & ! eastward wind (is there any difference to u?)
163          'uw        ',  & ! ? vertical momentum flux - total ?
164          'utheta    ',  & ! ? horizontal heat flux - total ?
165          'uv        ',  & ! upward-northward horizontal momentum flux
166          'v         ',  & ! v-component
167          'va        ',  & ! northward wind (is there any difference to v?)
168          'vw        ',  & ! ? vertical momentum flux - total ?
169          'vtheta    ',  & ! ? horizontal heat flux - total ?
170          'w         ',  & ! w-component
171          'wtheta    ',  & ! ? vertical heat flux - total ?
172          'rld       ',  & ! downward longwave radiative flux (W/m2)
173          'rlu       ',  & ! upnward longwave radiative flux (W/m2)
174          'rsd       ',  & ! downward shortwave radiative flux (W/m2)
175          'rsu       ',  & ! upward shortwave radiative flux (W/m2)
176          'rsddif    ',  & ! downward shortwave diffuse radiative flux (W/m2)
177          'rnds      '   & ! surface net downward radiative flux (W/m2)
178       /)
179                                                             
180    INTEGER(iwp) ::  id_vm                           !< NetCDF file id for virtual measurements
181    INTEGER(iwp) ::  nvm = 0                         !< number of virtual measurements
182    INTEGER(iwp) ::  observation_coverage_xy = 0     !< horizontal distance from the measurement point where observations should be taken in the surrounding
183    INTEGER(iwp) ::  observation_coverage_z  = 0     !< vertical distance from the measurement point where observations should be taken in the surrounding
184   
185    LOGICAL ::  use_virtual_measurement = .FALSE. !< Namelist parameter
186    LOGICAL ::  global_attribute = .TRUE.         !< flag indicating a global attribute
187   
188    REAL(wp) ::  vm_time_start = 0.0              !< time after virtual measurements should start
189   
190
191    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea !< virtual measurement data structure
192   
193   
194    INTERFACE vm_check_parameters
195       MODULE PROCEDURE vm_check_parameters
196    END INTERFACE vm_check_parameters
197   
198    INTERFACE vm_init
199       MODULE PROCEDURE vm_init
200    END INTERFACE vm_init
201   
202    INTERFACE vm_parin
203       MODULE PROCEDURE vm_parin
204    END INTERFACE vm_parin
205   
206    INTERFACE vm_sampling
207       MODULE PROCEDURE vm_sampling
208    END INTERFACE vm_sampling
209
210    SAVE
211
212    PRIVATE
213
214!
215!-- Public interfaces
216    PUBLIC  vm_check_parameters, vm_init, vm_parin, vm_sampling
217
218!
219!-- Public variables
220    PUBLIC  vmea, vm_time_start
221
222 CONTAINS
223
224
225!------------------------------------------------------------------------------!
226! Description:
227! ------------
228!> Check parameters for virtual measurement module
229!------------------------------------------------------------------------------!
230 SUBROUTINE vm_check_parameters
231
232    USE control_parameters,                                                    &
233        ONLY:  message_string, virtual_measurement
234 
235    USE netcdf_data_input_mod,                                                 &
236        ONLY:  input_pids_static
237       
238    IMPLICIT NONE
239   
240!
241!-- In case virtual measurements are taken, a static input file is required.
242!-- This is because UTM coordinates for the PALM domain origin are required
243!-- for correct mapping of the measurements.
244!-- ToDo: Revise this later and remove this requirement.
245    IF ( virtual_measurement  .AND.  .NOT. input_pids_static )  THEN
246       message_string = 'If virtual measurements are taken a static input ' // &
247                        'file is mandatory.'
248       CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
249    ENDIF
250 
251 END SUBROUTINE vm_check_parameters
252 
253!------------------------------------------------------------------------------!
254! Description:
255! ------------
256!> Read namelist for the virtual measurement module
257!------------------------------------------------------------------------------!
258 SUBROUTINE vm_parin
259 
260    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
261 
262    NAMELIST /virtual_measurement_parameters/  use_virtual_measurement,        &
263                                               vm_time_start
264
265    line = ' '
266
267!
268!-- Try to find stg package
269    REWIND ( 11 )
270    line = ' '
271    DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
272       READ ( 11, '(A)', END=20 )  line
273    ENDDO
274    BACKSPACE ( 11 )
275
276!
277!-- Read namelist
278    READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
279
280!
281!-- Set flag that indicates that the virtual measurement module is switched on
282    IF ( use_virtual_measurement )  virtual_measurement = .TRUE.
283   
284    GOTO 20
285
286 10 BACKSPACE( 11 )
287    READ( 11 , '(A)') line
288    CALL parin_fail_message( 'virtual_measurement_parameters', line )
289
290 20 CONTINUE
291 
292 END SUBROUTINE vm_parin
293
294
295!------------------------------------------------------------------------------!
296! Description:
297! ------------
298!> Initialize virtual measurements: read coordiante arrays and measured
299!> variables, set indicies indicating the measurement points, read further
300!> attributes, etc..
301!------------------------------------------------------------------------------!
302 SUBROUTINE vm_init
303
304    USE arrays_3d,                                                             &
305        ONLY:  zu, zw
306 
307    USE control_parameters,                                                    &
308        ONLY:  message_string
309       
310    USE grid_variables,                                                        &
311        ONLY:  ddx, ddy, dx, dy
312       
313    USE indices,                                                               &
314        ONLY:  nxl, nxr, nyn, nys
315 
316    USE netcdf_data_input_mod,                                                 &
317        ONLY:  init_model, input_file_vm,                                      &
318               netcdf_data_input_get_dimension_length,                         &
319               netcdf_data_input_att, netcdf_data_input_var
320               
321    USE surface_mod,                                                           &
322        ONLY:  get_topography_top_index_ji
323       
324    IMPLICIT NONE
325   
326    CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
327    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
328    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed
329   
330    LOGICAL ::  on_pe !< flag indicating that the respective measurement coordinate is on subdomain
331   
332    INTEGER(iwp) ::  dim_eutm  !< dimension size of UTM easting coordinate
333    INTEGER(iwp) ::  dim_nutm  !< dimension size of UTM northing coordinate
334    INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate
335    INTEGER(iwp) ::  dim_zag   !< dimension size of height coordinate
336    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
337    INTEGER(iwp) ::  icov      !< index range where observations should be taken in x-direction
338    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
339    INTEGER(iwp) ::  i_prev    !< grid index along x for UTM coordinate at previous observation time step
340    INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
341    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
342    INTEGER(iwp) ::  jcov      !< index range where observations should be taken in y-direction
343    INTEGER(iwp) ::  j_prev    !< grid index along y for UTM coordinate at previous observation time step
344    INTEGER(iwp) ::  js        !< grid index of real observation point of the respective station in y-direction
345    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
346    INTEGER(iwp) ::  kcov      !< index range where observations should be taken in z-direction
347    INTEGER(iwp) ::  ks        !< grid index of real observation point of the respective station in z-direction
348    INTEGER(iwp) ::  k_prev    !< grid index along z for UTM coordinate at previous observation time step
349    INTEGER(iwp) ::  ksurf     !< topography top index
350    INTEGER(iwp) ::  l         !< running index over all stations
351    INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
352    INTEGER(iwp) ::  ll        !< running index over all measured variables in file
353    INTEGER(iwp) ::  lll       !< running index over all allowed variables
354    INTEGER(iwp) ::  n         !< running index over trajectory coordinates
355    INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
356    INTEGER(iwp) ::  t         !< running index over number of trajectories
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!-- ALLOCATE data structure
371    ALLOCATE( vmea(1:nvm) )
372   
373!    print*, "nvm", nvm
374!
375!-- Read station coordinates and further attributes.
376!-- Note all coordinates are in UTM coordinates.
377    DO  l = 1, nvm
378!
379!--    Determine suffix which contains the ID
380       IF( l < 10 )  THEN
381          WRITE( dum, '(I1)')  l
382       ELSEIF( l < 100 )  THEN
383          WRITE( dum, '(I2)')  l
384       ELSEIF( l < 1000 )  THEN
385          WRITE( dum, '(I3)')  l
386       ELSEIF( l < 10000 )  THEN
387          WRITE( dum, '(I4)')  l
388       ELSEIF( l < 100000 )  THEN
389          WRITE( dum, '(I5)')  l
390       ENDIF
391       
392       CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx            &
393                                   // TRIM( dum ), id_vm, '', global_attribute,&
394                                   '', '' )
395       CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy            &
396                                   // TRIM( dum ), id_vm, '', global_attribute,&
397                                   '', '' )
398       CALL netcdf_data_input_att( vmea(l)%site,         char_site             &
399                                   // TRIM( dum ), id_vm, '', global_attribute,&
400                                   '', '' )
401       CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature          &
402                                   // TRIM( dum ), id_vm, '', global_attribute,&
403                                   '', '' )
404       
405!
406!---   Set logicals depending on the type of the measurement
407       IF ( INDEX( vmea(l)%feature_type, type_tspr     ) /= 0 )  THEN
408          vmea(l)%timseries_profile = .TRUE.
409       ELSEIF ( INDEX( vmea(l)%feature_type, type_ts   ) /= 0 )  THEN
410          vmea(l)%timseries         = .TRUE.
411       ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 )  THEN
412          vmea(l)%trajectory        = .TRUE.
413       ELSE
414!
415!--       Give error message
416          message_string = 'Attribue featureType = ' //                        &
417                           TRIM( vmea(l)%feature_type ) //                     &
418                           ' is not allowed.' 
419          CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 )
420       ENDIF
421!
422!--    Read string with all measured variables at this station
423       measured_variables_file = ''
424       CALL netcdf_data_input_var( measured_variables_file,                    &
425                                   char_mv // TRIM( dum ), id_vm )
426!
427!--    Count the number of measured variables which match with the variables
428!--    which are allowed to be measured in PALM. Please note, for some
429!--    NetCDF interal reasons characters end with a NULL, i.e. also empty
430!--    characters contain a NULL. Therefore, check the strings for a Null to
431!--    get the correct character length in order to compare them with the list
432!--    of allowed variables.
433       vmea(l)%nvar = 0
434       DO ll = 1, SIZE( measured_variables_file )
435          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
436               measured_variables_file(ll)(1:1) /= ' ')  THEN
437!
438!--          Obtain character length of the character
439             len_char = 1
440             DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
441                 .AND.  measured_variables_file(ll)(len_char:len_char) /= ' ' )
442                len_char = len_char + 1
443             ENDDO
444             len_char = len_char - 1
445!
446!--          Now, compare the measured variable with the list of allowed
447!--          variables.
448             DO  lll= 1, SIZE( list_allowed_variables )
449                IF ( measured_variables_file(ll)(1:len_char) ==                &
450                     TRIM( list_allowed_variables(lll) ) )  THEN
451                   vmea(l)%nvar = vmea(l)%nvar + 1
452                   measured_variables(vmea(l)%nvar) =                          &
453                                       measured_variables_file(ll)(1:len_char)
454                ENDIF
455             ENDDO
456          ENDIF
457       ENDDO
458!
459!--    Allocate array for the measured variables names for the station l.
460       ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
461
462       DO  ll = 1, vmea(l)%nvar
463          vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
464       ENDDO
465       
466!        print*, "numvars", vmea(l)%nvar, vmea(l)%measured_vars_name(1:vmea(l)%nvar)
467!
468!--    For the actual measurement ID read the UTM coordinates. Based on these,
469!--    define the index space on each subdomain where measurements should be
470!--    taken. Note, the entire coordinate arrays will not be stored on data
471!--    type as this would exceed memory requirements, particularly for
472!--    trajectory measurements. If no variable will be virtually measured,
473!--    skip the reading.
474       IF ( vmea(l)%nvar > 0 )  THEN
475!
476!--       For stationary measurements UTM coordinates are just one value and
477!--       its dimension is "station", while for mobile measurements UTM
478!--       coordinates are arrays. First, inquire dimension length for
479!--       UTM coordinates.
480          IF ( vmea(l)%trajectory )  THEN
481!
482!--          For non-stationary measurements read the number of trajectories
483             CALL netcdf_data_input_get_dimension_length( id_vm,              &
484                                                          vmea(l)%ntraj,      &
485                                                          "traj" //           &
486                                                          TRIM( dum ) )
487             CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime,   &
488                                                          "ntime" //          &
489                                                          TRIM( dum ) )
490!
491!--       For stationary measurements the dimension for UTM coordinates is 1
492          ELSE
493             vmea(l)%ntraj  = 1
494             dim_ntime = 1
495          ENDIF
496         
497!
498!-        Allocate array which defines individual time frame for each
499!--       trajectory or station
500          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
501          ALLOCATE( vmea(l)%ngp(1:vmea(l)%ntraj)   )
502!
503!--       Allocate temporary arrays for UTM and height coordinates. Note,
504!--       on file UTM coordinates might be 1D or 2D variables
505          ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
506          ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
507          ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime)  )
508!
509!--       Read _FillValue attributes
510          CALL netcdf_data_input_att( fill_eutm, char_fillvalue,               &
511                                      id_vm, '', .NOT. global_attribute, '',   &
512                                      char_eutm // TRIM( dum ) )
513          CALL netcdf_data_input_att( fill_nutm, char_fillvalue,               &
514                                      id_vm, '', .NOT. global_attribute, '',   &
515                                      char_nutm // TRIM( dum ) )
516          CALL netcdf_data_input_att( fill_zag, char_fillvalue,                &
517                                      id_vm, '', .NOT. global_attribute, '',   &
518                                      char_zag  // TRIM( dum ) )
519!
520!--       Read UTM and height coordinates coordinates for all trajectories and
521!--       times.
522          IF ( vmea(l)%trajectory )  THEN
523             CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm,  &
524                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
525             CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm,  &
526                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
527             CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), id_vm,  &
528                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
529          ELSE
530             CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm )
531             CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm )
532             CALL netcdf_data_input_var( z_ag(1,:),  char_zag  // TRIM( dum ), id_vm )
533          ENDIF
534!
535!--       Based on UTM coordinates, check if the measurement station or parts
536!--       of the trajectory is on subdomain. This case, setup grid index space
537!--       sample these quantities.
538          ns = 0
539          DO  t = 1, vmea(l)%ntraj
540!
541!--          Determine the individual time coordinate length for each station and
542!--          trajectory. This is required as several stations and trajectories
543!--          are merged into one file but they do not have the same number of
544!--          points in time, hence, missing values may occur and cannot be
545!--          processed further.
546             vmea(l)%dim_t(t) = 0
547             DO  n = 1, dim_ntime
548                IF ( e_utm(t,n) /= fill_eutm  .AND.                            &
549                     n_utm(t,n) /= fill_nutm  .AND.                            &
550                     z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
551             ENDDO
552!             
553!--          First, compute relative x- and y-coordinates with respect to the
554!--          lower-left origin of the model domain, which is the difference
555!--          betwen UTM coordinates.
556!              write(9,*) l, "before Eutm", e_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_x
557!              write(9,*) l, "before Nutm", n_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_y
558             
559             e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t))         &
560                                         - init_model%origin_x
561             n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t))         &
562                                         - init_model%origin_y
563!              write(9,*) l, "Eutm", e_utm(t,1:vmea(l)%dim_t(t))
564!              write(9,*) l, "Nutm", n_utm(t,1:vmea(l)%dim_t(t))
565!              write(9,*)
566!
567!--          Compute grid indices relative to origin and check if these are
568!--          on the subdomain. Note, virtual measurements will be taken also
569!--          at grid points surrounding the station, hence, check also for
570!--          these grid points.
571             vmea(l)%ngp(t) = 0
572             k_prev = -999
573             j_prev = -999
574             i_prev = -999
575             DO  n = 1, vmea(l)%dim_t(t)
576                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
577                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
578!
579!--             Is the observation point on subdomain?
580                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.                   &
581                          js >= nys  .AND.  js <= nyn )
582!
583!--             If the measurement is on subdomain, determine vertical index
584!--             which refers to the observation height above ground level.
585                ks = k_prev
586                IF ( on_pe )  THEN
587!                    write(9,*) " I am onpe"
588                   ksurf = get_topography_top_index_ji( js, is, 's' )
589                   ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
590                ENDIF
591!
592!--             Count the number of observation points in index space on
593!--             subdomain. Only increment if grid indices are different from
594!--             the previous one.
595                IF ( on_pe  .AND.  is /= i_prev  .AND.  js /= j_prev  .AND.    &
596                                   ks /= k_prev )  THEN
597                   ns             = ns             + 1
598                   vmea(l)%ngp(t) = vmea(l)%ngp(t) + 1
599                ENDIF
600
601!--             Store arrays for next iteration - avoid double counting
602                i_prev = is
603                j_prev = js
604                k_prev = ks
605             ENDDO
606             
607          ENDDO
608
609!
610!--       Store number of observation points on subdomain and allocate index
611!--       arrays.
612          vmea(l)%ns = ns
613         
614          ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
615          ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
616          ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
617         
618!           print*, "Num ns: ", vmea(l)%ns, "per traj", vmea(l)%ngp(:)
619!
620!--       Repeat the prior loop and save the grid indices relevant for
621!--       sampling.
622          ns = 0
623          DO  t = 1, vmea(l)%ntraj
624!
625!--          Compute grid indices relative to origin and check if these are
626!--          on the subdomain. Note, virtual measurements will be taken also
627!--          at grid points surrounding the station, hence, check also for
628!--          these grid points.
629             k_prev = -999
630             j_prev = -999
631             i_prev = -999
632             DO  n = 1, vmea(l)%dim_t(t)
633                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
634                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
635!
636!--             Is the observation point on subdomain?
637                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.                   &
638                          js >= nys  .AND.  js <= nyn )
639!
640!--             If the measurement is on subdomain, determine vertical index
641!--             which refers to the observation height above ground level.
642                ks = k_prev
643                IF ( on_pe )  THEN
644                   ksurf = get_topography_top_index_ji( js, is, 's' )
645                   ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
646                ENDIF
647!
648!--             Count the number of observation points in index space on
649!--             subdomain. Only increment if grid indices are different from
650!--             the previous one.
651                IF ( on_pe  .AND.  is /= i_prev  .AND.  js /= j_prev  .AND.    &
652                                   ks /= k_prev )  THEN
653                   ns             = ns + 1
654                   vmea(l)%i(ns)  = is
655                   vmea(l)%j(ns)  = js
656                   vmea(l)%k(ns)  = ks
657!                    write(9,*) TRIM( vmea(l)%feature_type), l, "ns", ns, "ijk", vmea(l)%i(ns), vmea(l)%j(ns), vmea(l)%k(ns)
658                ENDIF
659!
660!--             Store arrays for next iteration - avoid double counting
661                i_prev = is
662                j_prev = js
663                k_prev = ks
664             ENDDO
665             
666          ENDDO
667!
668!--       Allocate array to save the sampled values.
669!--       Todo: Is it better to allocate for all variables at a station
670!--       and store all the values before writing, or sample the variables
671!--       directly in the data output?
672          ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) )
673!
674!--       Initialize with _FillValue
675          vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) = vmea(l)%fillout
676!
677!--       Deallocate temporary coordinate arrays
678          IF ( ALLOCATED( e_utm ) )  DEALLOCATE( e_utm )
679          IF ( ALLOCATED( n_utm ) )  DEALLOCATE( n_utm )
680          IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( z_ag  )
681       ENDIF
682    ENDDO
683    flush(9)
684   
685!
686!-- Close input file for virtual measurements. Therefore, just call
687!-- the read attribute routine with the "close" option.
688    CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '',              &
689                                global_attribute, 'close', '' )
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:  pt
702
703    USE surface_mod
704   
705     IMPLICIT NONE
706     
707     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
708     
709     INTEGER(iwp) ::  l !<
710     INTEGER(iwp) ::  m !<
711     INTEGER(iwp) ::  var !<
712     
713     INTEGER(iwp) ::  mm, j, i
714
715!      write(9,*) "sampling"
716!      flush(9)
717!
718!--  Loop over all stations. For each possible variable loop over all
719!--  observation points
720     DO  l = 1, nvm
721!
722!--     Loop over all measured variables. Please note, for the moment
723!--     the same indices for scalar and velocity components are used.
724!--     ToDo: Revise this later.
725!         DO  m = 1, vmea(l)%ns
726!            j = vmea(l)%j(m)
727!            i = vmea(l)%i(m)
728!           
729!            IF ( i >= nxl  .AND.  i <= nxr  .AND.                               &
730!                 j >= nys  .AND.  j <= nyn )  THEN
731!               IF ( surf_def_h(0)%start_index(j,i) <= &
732!                    surf_def_h(0)%end_index(j,i) )  THEN
733!                 
734!                  write(9,*) "sampled"
735!                  flush(9)
736!                  mm = surf_def_h(0)%end_index(j,i)
737!               
738!                  surf_def_h(0)%pt_surface(mm) = -99.0
739!               ENDIF
740!            ENDIF
741!         ENDDO
742!         DO  var = 1, vmea(l)%nvar
743!            trimvar = TRIM( vmea(l)%measured_vars_name(var) )
744!           
745!            IF ( TRIM( trimvar ) == 'theta' )  THEN           
746!               DO  m = 1, vmea(l)%ns
747!                  vmea(l)%measured_vars(var,m) = pt(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
748!               ENDDO
749!            ENDIF           
750!            IF ( TRIM( trimvar ) == 'w' )  THEN           
751!               DO  m = 1, vmea(l)%ns
752!                  vmea(l)%measured_vars(var,m) = w(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
753!               ENDDO
754!            ENDIF
755!            IF ( TRIM( trimvar ) == 'v' )  THEN           
756!               DO  m = 1, vmea(l)%ns
757!                  vmea(l)%measured_vars(var,m) = v(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
758!               ENDDO
759!            ENDIF
760!            IF ( TRIM( trimvar ) == 'u' )  THEN           
761!               DO  m = 1, vmea(l)%ns
762!                  vmea(l)%measured_vars(var,m) = u(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
763!               ENDDO
764!            ENDIF
765!         ENDDO
766
767     ENDDO
768     
769  END SUBROUTINE vm_sampling
770 
771
772 END MODULE virtual_measurement_mod
Note: See TracBrowser for help on using the repository browser.