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

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

module for virtual measurements added (in a preliminary state); new public routines to input NetCDF data directly from modules

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