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