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 4329 2019-12-10 15:46:36Z suehring $ |
---|
27 | ! Renamed wall_flags_0 to wall_flags_static_0 |
---|
28 | ! |
---|
29 | ! 4226 2019-09-10 17:03:24Z suehring |
---|
30 | ! Netcdf input routine for dimension length renamed |
---|
31 | ! |
---|
32 | ! 4182 2019-08-22 15:20:23Z scharf |
---|
33 | ! Corrected "Former revisions" section |
---|
34 | ! |
---|
35 | ! 4168 2019-08-16 13:50:17Z suehring |
---|
36 | ! Replace function get_topography_top_index by topo_top_ind |
---|
37 | ! |
---|
38 | ! 3988 2019-05-22 11:32:37Z kanani |
---|
39 | ! Add variables to enable steering of output interval for virtual measurements |
---|
40 | ! |
---|
41 | ! 3913 2019-04-17 15:12:28Z gronemeier |
---|
42 | ! Bugfix: rotate positions of measurements before writing them into file |
---|
43 | ! |
---|
44 | ! 3910 2019-04-17 11:46:56Z suehring |
---|
45 | ! Bugfix in rotation of UTM coordinates |
---|
46 | ! |
---|
47 | ! 3904 2019-04-16 18:22:51Z gronemeier |
---|
48 | ! Rotate coordinates of stations by given rotation_angle |
---|
49 | ! |
---|
50 | ! 3876 2019-04-08 18:41:49Z knoop |
---|
51 | ! Remove print statement |
---|
52 | ! |
---|
53 | ! 3854 2019-04-02 16:59:33Z suehring |
---|
54 | ! renamed nvar to nmeas, replaced USE chem_modules by USE chem_gasphase_mod and |
---|
55 | ! nspec by nvar |
---|
56 | ! |
---|
57 | ! 3766 2019-02-26 16:23:41Z raasch |
---|
58 | ! unused variables removed |
---|
59 | ! |
---|
60 | ! 3718 2019-02-06 11:08:28Z suehring |
---|
61 | ! Adjust variable name connections between UC2 and chemistry variables |
---|
62 | ! |
---|
63 | ! 3717 2019-02-05 17:21:16Z suehring |
---|
64 | ! Additional check + error numbers adjusted |
---|
65 | ! |
---|
66 | ! 3706 2019-01-29 20:02:26Z suehring |
---|
67 | ! unused variables removed |
---|
68 | ! |
---|
69 | ! 3705 2019-01-29 19:56:39Z suehring |
---|
70 | ! - initialization revised |
---|
71 | ! - binary data output |
---|
72 | ! - list of allowed variables extended |
---|
73 | ! |
---|
74 | ! 3704 2019-01-29 19:51:41Z suehring |
---|
75 | ! Sampling of variables |
---|
76 | ! |
---|
77 | ! 3473 2018-10-30 20:50:15Z suehring |
---|
78 | ! Initial revision |
---|
79 | ! |
---|
80 | ! Authors: |
---|
81 | ! -------- |
---|
82 | ! @author Matthias Suehring |
---|
83 | ! |
---|
84 | ! Description: |
---|
85 | ! ------------ |
---|
86 | !> The module acts as an interface between 'real-world' observations and |
---|
87 | !> model simulations. Virtual measurements will be taken in the model at the |
---|
88 | !> coordinates representative for the 'real-world' observation coordinates. |
---|
89 | !> More precisely, coordinates and measured quanties will be read from a |
---|
90 | !> NetCDF file which contains all required information. In the model, |
---|
91 | !> the same quantities (as long as all the required components are switched-on) |
---|
92 | !> will be sampled at the respective positions and output into an extra file, |
---|
93 | !> which allows for straight-forward comparison of model results with |
---|
94 | !> observations. |
---|
95 | !> |
---|
96 | !> @todo list_of_allowed variables needs careful checking |
---|
97 | !> @todo Check if sign of surface fluxes for heat, radiation, etc., follows |
---|
98 | !> the (UC)2 standard |
---|
99 | !> @note Fluxes are not processed |
---|
100 | !------------------------------------------------------------------------------! |
---|
101 | MODULE virtual_measurement_mod |
---|
102 | |
---|
103 | USE arrays_3d, & |
---|
104 | ONLY: q, pt, u, v, w, zu, zw |
---|
105 | |
---|
106 | USE basic_constants_and_equations_mod, & |
---|
107 | ONLY: pi |
---|
108 | |
---|
109 | USE chem_gasphase_mod, & |
---|
110 | ONLY: nvar |
---|
111 | |
---|
112 | USE chem_modules, & |
---|
113 | ONLY: chem_species |
---|
114 | |
---|
115 | USE control_parameters, & |
---|
116 | ONLY: air_chemistry, dz, humidity, io_blocks, io_group, neutral, & |
---|
117 | message_string, time_since_reference_point, virtual_measurement |
---|
118 | |
---|
119 | USE cpulog, & |
---|
120 | ONLY: cpu_log, log_point |
---|
121 | |
---|
122 | USE grid_variables, & |
---|
123 | ONLY: dx, dy |
---|
124 | |
---|
125 | USE indices, & |
---|
126 | ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny, topo_top_ind, wall_flags_static_0 |
---|
127 | |
---|
128 | USE kinds |
---|
129 | |
---|
130 | USE netcdf_data_input_mod, & |
---|
131 | ONLY: init_model |
---|
132 | |
---|
133 | USE pegrid |
---|
134 | |
---|
135 | USE surface_mod, & |
---|
136 | ONLY: surf_lsm_h, surf_usm_h |
---|
137 | |
---|
138 | USE land_surface_model_mod, & |
---|
139 | ONLY: nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h |
---|
140 | |
---|
141 | USE radiation_model_mod |
---|
142 | |
---|
143 | USE urban_surface_mod, & |
---|
144 | ONLY: nzb_wall, nzt_wall, t_wall_h |
---|
145 | |
---|
146 | |
---|
147 | IMPLICIT NONE |
---|
148 | |
---|
149 | TYPE virt_general |
---|
150 | INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements |
---|
151 | INTEGER(iwp) :: nvm = 0 !< number of virtual measurements |
---|
152 | END TYPE virt_general |
---|
153 | |
---|
154 | TYPE virt_mea |
---|
155 | |
---|
156 | CHARACTER(LEN=100) :: feature_type !< type of the measurement |
---|
157 | CHARACTER(LEN=100) :: filename_original !< name of the original file |
---|
158 | CHARACTER(LEN=100) :: site !< name of the measurement site |
---|
159 | |
---|
160 | CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables |
---|
161 | |
---|
162 | INTEGER(iwp) :: ns = 0 !< number of observation coordinates on subdomain, for atmospheric measurements |
---|
163 | INTEGER(iwp) :: ns_tot = 0 !< total number of observation coordinates, for atmospheric measurements |
---|
164 | INTEGER(iwp) :: ntraj !< number of trajectories of a measurement |
---|
165 | INTEGER(iwp) :: nmeas !< number of measured variables (atmosphere + soil) |
---|
166 | |
---|
167 | INTEGER(iwp) :: ns_soil = 0 !< number of observation coordinates on subdomain, for soil measurements |
---|
168 | INTEGER(iwp) :: ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements |
---|
169 | |
---|
170 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues |
---|
171 | |
---|
172 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction |
---|
173 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction |
---|
174 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction |
---|
175 | |
---|
176 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_soil !< grid index for measurement position in x-direction |
---|
177 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_soil !< grid index for measurement position in y-direction |
---|
178 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_soil !< grid index for measurement position in k-direction |
---|
179 | |
---|
180 | LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation |
---|
181 | LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement |
---|
182 | LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement |
---|
183 | LOGICAL :: soil_sampling = .FALSE. !< flag indicating that soil state variables were sampled |
---|
184 | |
---|
185 | REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values |
---|
186 | REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values |
---|
187 | REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values |
---|
188 | REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building |
---|
189 | REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction |
---|
190 | REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction |
---|
191 | |
---|
192 | REAL(wp), DIMENSION(:), ALLOCATABLE :: z_ag !< measurement height above ground level |
---|
193 | REAL(wp), DIMENSION(:), ALLOCATABLE :: depth !< measurement depth in soil |
---|
194 | |
---|
195 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables |
---|
196 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars_soil !< measured variables |
---|
197 | |
---|
198 | END TYPE virt_mea |
---|
199 | |
---|
200 | CHARACTER(LEN=5) :: char_eutm = "E_UTM" !< dimension name for UTM coordinate easting |
---|
201 | CHARACTER(LEN=11) :: char_feature = "featureType" !< attribute name for feature type |
---|
202 | |
---|
203 | ! This need to generalized |
---|
204 | CHARACTER(LEN=8) :: char_filename = "filename" !< attribute name for filename |
---|
205 | CHARACTER(LEN=11) :: char_soil = "soil_sample" !< attribute name for soil sampling indication |
---|
206 | CHARACTER(LEN=10) :: char_fillvalue = "_FillValue" !< variable attribute name for _FillValue |
---|
207 | CHARACTER(LEN=18) :: char_mv = "measured_variables" !< variable name for the array with the measured variable names |
---|
208 | CHARACTER(LEN=5) :: char_nutm = "N_UTM" !< dimension name for UTM coordinate northing |
---|
209 | CHARACTER(LEN=18) :: char_numstations = "number_of_stations" !< attribute name for number of stations |
---|
210 | CHARACTER(LEN=8) :: char_origx = "origin_x" !< attribute name for station coordinate in x |
---|
211 | CHARACTER(LEN=8) :: char_origy = "origin_y" !< attribute name for station coordinate in y |
---|
212 | CHARACTER(LEN=4) :: char_site = "site" !< attribute name for site name |
---|
213 | CHARACTER(LEN=19) :: char_zag = "height_above_ground" !< attribute name for height above ground variable |
---|
214 | CHARACTER(LEN=10) :: type_ts = 'timeSeries' !< name of stationary point measurements |
---|
215 | CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements |
---|
216 | CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements |
---|
217 | |
---|
218 | CHARACTER(LEN=6), DIMENSION(1:5) :: soil_vars = (/ & !< list of soil variables |
---|
219 | 't_soil', & |
---|
220 | 'm_soil', & |
---|
221 | 'lwc ', & |
---|
222 | 'lwcs ', & |
---|
223 | 'smp ' /) |
---|
224 | |
---|
225 | CHARACTER(LEN=10), DIMENSION(0:1,1:8) :: chem_vars = RESHAPE( (/ & |
---|
226 | 'mcpm1 ', 'PM1 ', & |
---|
227 | 'mcpm2p5 ', 'PM2.5 ', & |
---|
228 | 'mcpm25 ', 'PM25 ', & |
---|
229 | 'mcpm10 ', 'PM10 ', & |
---|
230 | 'mfno2 ', 'NO2 ', & |
---|
231 | 'mfno ', 'NO ', & |
---|
232 | 'tro3 ', 'O3 ', & |
---|
233 | 'mfco ', 'CO ' & |
---|
234 | /), (/ 2, 8 /) ) |
---|
235 | ! |
---|
236 | !-- MS: List requires careful revision! |
---|
237 | CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM |
---|
238 | (/ 'hfls ', & ! surface latent heat flux (W/m2) |
---|
239 | 'hfss ', & ! surface sensible heat flux (W/m2) |
---|
240 | 'hur ', & ! relative humidity (-) |
---|
241 | 'hus ', & ! specific humidity (g/kg) |
---|
242 | 'haa ', & ! absolute atmospheric humidity (kg/m3) |
---|
243 | 'mcpm1 ', & ! mass concentration of PM1 (kg/m3) |
---|
244 | 'mcpm2p5 ', & ! mass concentration of PM2.5 (kg/m3) |
---|
245 | 'mcpm10 ', & ! mass concentration of PM10 (kg/m3) |
---|
246 | 'mcco ', & ! mass concentration of CO (kg/m3) |
---|
247 | 'mcco2 ', & ! mass concentration of CO2 (kg/m3) |
---|
248 | 'mcbcda ', & ! mass concentration of black carbon paritcles (kg/m3) |
---|
249 | 'ncaa ', & ! number concentation of particles (1/m3) |
---|
250 | 'mfco ', & ! mole fraction of CO (mol/mol) |
---|
251 | 'mfco2 ', & ! mole fraction of CO2 (mol/mol) |
---|
252 | 'mfch4 ', & ! mole fraction of methane (mol/mol) |
---|
253 | 'mfnh3 ', & ! mole fraction of amonia (mol/mol) |
---|
254 | 'mfno ', & ! mole fraction of nitrogen monoxide (mol/mol) |
---|
255 | 'mfno2 ', & ! mole fraction of nitrogen dioxide (mol/mol) |
---|
256 | 'mfso2 ', & ! mole fraction of sulfur dioxide (mol/mol) |
---|
257 | 'mfh20 ', & ! mole fraction of water (mol/mol) |
---|
258 | 'plev ', & ! ? air pressure - hydrostaic + perturbation? |
---|
259 | 'rlds ', & ! surface downward longwave flux (W/m2) |
---|
260 | 'rlus ', & ! surface upward longwave flux (W/m2) |
---|
261 | 'rsds ', & ! surface downward shortwave flux (W/m2) |
---|
262 | 'rsus ', & ! surface upward shortwave flux (W/m2) |
---|
263 | 'ta ', & ! air temperature (degree C) |
---|
264 | 't_va ', & ! virtual accoustic temperature (K) |
---|
265 | 'theta ', & ! potential temperature (K) |
---|
266 | 'tro3 ', & ! mole fraction of ozone air (mol/mol) |
---|
267 | 'ts ', & ! scaling parameter of temperature (K) |
---|
268 | 'wspeed ', & ! ? wind speed - horizontal? |
---|
269 | 'wdir ', & ! wind direction |
---|
270 | 'us ', & ! friction velocity |
---|
271 | 'msoil ', & ! ? soil moisture - which depth? |
---|
272 | 'tsoil ', & ! ? soil temperature - which depth? |
---|
273 | 'u ', & ! u-component |
---|
274 | 'utheta ', & ! total eastward kinematic heat flux |
---|
275 | 'ua ', & ! eastward wind (is there any difference to u?) |
---|
276 | 'v ', & ! v-component |
---|
277 | 'vtheta ', & ! total northward kinematic heat flux |
---|
278 | 'va ', & ! northward wind (is there any difference to v?) |
---|
279 | 'w ', & ! w-component |
---|
280 | 'wtheta ', & ! total vertical kinematic heat flux |
---|
281 | 'rld ', & ! downward longwave radiative flux (W/m2) |
---|
282 | 'rlu ', & ! upnward longwave radiative flux (W/m2) |
---|
283 | 'rsd ', & ! downward shortwave radiative flux (W/m2) |
---|
284 | 'rsu ', & ! upward shortwave radiative flux (W/m2) |
---|
285 | 'rsddif ', & ! downward shortwave diffuse radiative flux (W/m2) |
---|
286 | 'rnds ', & ! surface net downward radiative flux (W/m2) |
---|
287 | 't_soil ', & |
---|
288 | 'm_soil ', & |
---|
289 | 'lwc ', & |
---|
290 | 'lwcs ', & |
---|
291 | 'smp ' & |
---|
292 | /) |
---|
293 | |
---|
294 | |
---|
295 | LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute |
---|
296 | LOGICAL :: init = .TRUE. !< flag indicating initialization of data output |
---|
297 | LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter |
---|
298 | |
---|
299 | REAL(wp) :: dt_virtual_measurement = 0.0_wp !< namelist parameter |
---|
300 | REAL(wp) :: time_virtual_measurement = 0.0_wp !< time since last 3d output |
---|
301 | REAL(wp) :: vm_time_start = 0.0 !< time after virtual measurements should start |
---|
302 | |
---|
303 | TYPE( virt_general ) :: vmea_general !< data structure which encompass general variables |
---|
304 | TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< virtual measurement data structure |
---|
305 | |
---|
306 | INTERFACE vm_check_parameters |
---|
307 | MODULE PROCEDURE vm_check_parameters |
---|
308 | END INTERFACE vm_check_parameters |
---|
309 | |
---|
310 | INTERFACE vm_data_output |
---|
311 | MODULE PROCEDURE vm_data_output |
---|
312 | END INTERFACE vm_data_output |
---|
313 | |
---|
314 | INTERFACE vm_init |
---|
315 | MODULE PROCEDURE vm_init |
---|
316 | END INTERFACE vm_init |
---|
317 | |
---|
318 | INTERFACE vm_last_actions |
---|
319 | MODULE PROCEDURE vm_last_actions |
---|
320 | END INTERFACE vm_last_actions |
---|
321 | |
---|
322 | INTERFACE vm_parin |
---|
323 | MODULE PROCEDURE vm_parin |
---|
324 | END INTERFACE vm_parin |
---|
325 | |
---|
326 | INTERFACE vm_sampling |
---|
327 | MODULE PROCEDURE vm_sampling |
---|
328 | END INTERFACE vm_sampling |
---|
329 | |
---|
330 | SAVE |
---|
331 | |
---|
332 | PRIVATE |
---|
333 | |
---|
334 | ! |
---|
335 | !-- Public interfaces |
---|
336 | PUBLIC vm_check_parameters, vm_data_output, vm_init, vm_last_actions, & |
---|
337 | vm_parin, vm_sampling |
---|
338 | |
---|
339 | ! |
---|
340 | !-- Public variables |
---|
341 | PUBLIC dt_virtual_measurement, time_virtual_measurement, vmea, vmea_general, vm_time_start |
---|
342 | |
---|
343 | CONTAINS |
---|
344 | |
---|
345 | |
---|
346 | !------------------------------------------------------------------------------! |
---|
347 | ! Description: |
---|
348 | ! ------------ |
---|
349 | !> Check parameters for virtual measurement module |
---|
350 | !------------------------------------------------------------------------------! |
---|
351 | SUBROUTINE vm_check_parameters |
---|
352 | |
---|
353 | USE control_parameters, & |
---|
354 | ONLY: message_string, virtual_measurement |
---|
355 | |
---|
356 | USE netcdf_data_input_mod, & |
---|
357 | ONLY: input_pids_static, input_pids_vm |
---|
358 | |
---|
359 | IMPLICIT NONE |
---|
360 | |
---|
361 | ! |
---|
362 | !-- Virtual measurements require a setup file. |
---|
363 | IF ( virtual_measurement .AND. .NOT. input_pids_vm ) THEN |
---|
364 | message_string = 'If virtual measurements are taken, a setup input ' // & |
---|
365 | 'file for the site locations is mandatory.' |
---|
366 | CALL message( 'vm_check_parameters', 'PA0533', 1, 2, 0, 6, 0 ) |
---|
367 | ENDIF |
---|
368 | ! |
---|
369 | !-- In case virtual measurements are taken, a static input file is required. |
---|
370 | !-- This is because UTM coordinates for the PALM domain origin are required |
---|
371 | !-- for correct mapping of the measurements. |
---|
372 | !-- ToDo: Revise this later and remove this requirement. |
---|
373 | IF ( virtual_measurement .AND. .NOT. input_pids_static ) THEN |
---|
374 | message_string = 'If virtual measurements are taken, a static input ' //& |
---|
375 | 'file is mandatory.' |
---|
376 | CALL message( 'vm_check_parameters', 'PA0534', 1, 2, 0, 6, 0 ) |
---|
377 | ENDIF |
---|
378 | |
---|
379 | END SUBROUTINE vm_check_parameters |
---|
380 | |
---|
381 | !------------------------------------------------------------------------------! |
---|
382 | ! Description: |
---|
383 | ! ------------ |
---|
384 | !> Read namelist for the virtual measurement module |
---|
385 | !------------------------------------------------------------------------------! |
---|
386 | SUBROUTINE vm_parin |
---|
387 | |
---|
388 | CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file |
---|
389 | |
---|
390 | NAMELIST /virtual_measurement_parameters/ dt_virtual_measurement, & |
---|
391 | use_virtual_measurement, & |
---|
392 | vm_time_start |
---|
393 | |
---|
394 | line = ' ' |
---|
395 | |
---|
396 | ! |
---|
397 | !-- Try to find stg package |
---|
398 | REWIND ( 11 ) |
---|
399 | line = ' ' |
---|
400 | DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 ) |
---|
401 | READ ( 11, '(A)', END=20 ) line |
---|
402 | ENDDO |
---|
403 | BACKSPACE ( 11 ) |
---|
404 | |
---|
405 | ! |
---|
406 | !-- Read namelist |
---|
407 | READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 ) |
---|
408 | |
---|
409 | ! |
---|
410 | !-- Set flag that indicates that the virtual measurement module is switched on |
---|
411 | IF ( use_virtual_measurement ) virtual_measurement = .TRUE. |
---|
412 | |
---|
413 | GOTO 20 |
---|
414 | |
---|
415 | 10 BACKSPACE( 11 ) |
---|
416 | READ( 11 , '(A)') line |
---|
417 | CALL parin_fail_message( 'virtual_measurement_parameters', line ) |
---|
418 | |
---|
419 | 20 CONTINUE |
---|
420 | |
---|
421 | END SUBROUTINE vm_parin |
---|
422 | |
---|
423 | |
---|
424 | !------------------------------------------------------------------------------! |
---|
425 | ! Description: |
---|
426 | ! ------------ |
---|
427 | !> Initialize virtual measurements: read coordiante arrays and measured |
---|
428 | !> variables, set indicies indicating the measurement points, read further |
---|
429 | !> attributes, etc.. |
---|
430 | !------------------------------------------------------------------------------! |
---|
431 | SUBROUTINE vm_init |
---|
432 | |
---|
433 | USE arrays_3d, & |
---|
434 | ONLY: zu, zw |
---|
435 | |
---|
436 | USE grid_variables, & |
---|
437 | ONLY: ddx, ddy, dx, dy |
---|
438 | |
---|
439 | USE indices, & |
---|
440 | ONLY: nxl, nxr, nyn, nys |
---|
441 | |
---|
442 | USE netcdf_data_input_mod, & |
---|
443 | ONLY: get_dimension_length, & |
---|
444 | init_model, & |
---|
445 | input_file_vm, & |
---|
446 | netcdf_data_input_att, & |
---|
447 | netcdf_data_input_var |
---|
448 | |
---|
449 | IMPLICIT NONE |
---|
450 | |
---|
451 | CHARACTER(LEN=5) :: dum !< dummy string indicate station id |
---|
452 | CHARACTER(LEN=10), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF |
---|
453 | CHARACTER(LEN=10), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed |
---|
454 | |
---|
455 | INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate |
---|
456 | INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction |
---|
457 | INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction |
---|
458 | INTEGER(iwp) :: j !< grid index of observation point in x-direction |
---|
459 | INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction |
---|
460 | INTEGER(iwp) :: k !< grid index of observation point in x-direction |
---|
461 | INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate |
---|
462 | INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction |
---|
463 | INTEGER(iwp) :: ksurf !< topography top index |
---|
464 | INTEGER(iwp) :: ku !< upper vertical index of surrounding grid points of an observation coordinate |
---|
465 | INTEGER(iwp) :: l !< running index over all stations |
---|
466 | INTEGER(iwp) :: len_char !< character length of single measured variables without Null character |
---|
467 | INTEGER(iwp) :: ll !< running index over all measured variables in file |
---|
468 | INTEGER(iwp) :: lll !< running index over all allowed variables |
---|
469 | INTEGER(iwp) :: n !< running index over trajectory coordinates |
---|
470 | INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain |
---|
471 | INTEGER(iwp) :: t !< running index over number of trajectories |
---|
472 | INTEGER(iwp) :: m |
---|
473 | |
---|
474 | INTEGER(KIND=1):: soil_dum |
---|
475 | |
---|
476 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ns_all !< dummy array used to sum-up the number of observation coordinates |
---|
477 | |
---|
478 | INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: meas_flag !< mask array indicating measurement positions |
---|
479 | |
---|
480 | ! LOGICAL :: chem_include !< flag indicating that chemical species is considered in modelled mechanism |
---|
481 | LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain |
---|
482 | |
---|
483 | REAL(wp) :: fill_eutm !< _FillValue for coordinate array E_UTM |
---|
484 | REAL(wp) :: fill_nutm !< _FillValue for coordinate array N_UTM |
---|
485 | REAL(wp) :: fill_zag !< _FillValue for height coordinate |
---|
486 | |
---|
487 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm !< easting UTM coordinate, temporary variable |
---|
488 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm !< northing UTM coordinate, temporary variable |
---|
489 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm_tmp !< EUTM coordinate before rotation |
---|
490 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm_tmp !< NUTM coordinate before rotation |
---|
491 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z_ag !< height coordinate relative to origin_z, temporary variable |
---|
492 | ! |
---|
493 | !-- Obtain number of sites. Also, pass the 'open' string, in order to initially |
---|
494 | !-- open the measurement driver. |
---|
495 | CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, & |
---|
496 | vmea_general%id_vm, input_file_vm, & |
---|
497 | global_attribute, 'open', '' ) |
---|
498 | ! |
---|
499 | !-- Allocate data structure which encompass all required information, such as |
---|
500 | !-- grid points indicies, absolute UTM coordinates, the measured quantities, |
---|
501 | !-- etc. . |
---|
502 | ALLOCATE( vmea(1:vmea_general%nvm) ) |
---|
503 | ! |
---|
504 | !-- Allocate flag array. This dummy array is used to identify grid points |
---|
505 | !-- where virtual measurements should be taken. Please note, at least one |
---|
506 | !-- ghost point is required, in order to include also the surrounding |
---|
507 | !-- grid points of the original coordinate. |
---|
508 | ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) |
---|
509 | meas_flag = 0 |
---|
510 | ! |
---|
511 | !-- Loop over all sites. |
---|
512 | DO l = 1, vmea_general%nvm |
---|
513 | ! |
---|
514 | !-- Determine suffix which contains the ID, ordered according to the number |
---|
515 | !-- of measurements. |
---|
516 | IF( l < 10 ) THEN |
---|
517 | WRITE( dum, '(I1)') l |
---|
518 | ELSEIF( l < 100 ) THEN |
---|
519 | WRITE( dum, '(I2)') l |
---|
520 | ELSEIF( l < 1000 ) THEN |
---|
521 | WRITE( dum, '(I3)') l |
---|
522 | ELSEIF( l < 10000 ) THEN |
---|
523 | WRITE( dum, '(I4)') l |
---|
524 | ELSEIF( l < 100000 ) THEN |
---|
525 | WRITE( dum, '(I5)') l |
---|
526 | ENDIF |
---|
527 | ! |
---|
528 | !-- Read site coordinates (UTM). |
---|
529 | CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx // & |
---|
530 | TRIM( dum ), vmea_general%id_vm, '', & |
---|
531 | global_attribute, '', '' ) |
---|
532 | CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy // & |
---|
533 | TRIM( dum ), vmea_general%id_vm, '', & |
---|
534 | global_attribute, '', '' ) |
---|
535 | ! |
---|
536 | !-- Read site name |
---|
537 | CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ), & |
---|
538 | vmea_general%id_vm, '', global_attribute, & |
---|
539 | '', '' ) |
---|
540 | ! |
---|
541 | !-- Read type of the measurement (trajectory, profile, timeseries). |
---|
542 | CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature // & |
---|
543 | TRIM( dum ), vmea_general%id_vm, '', & |
---|
544 | global_attribute, '', '' ) |
---|
545 | ! |
---|
546 | !-- Read the name of the original file where observational data is stored. |
---|
547 | CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // & |
---|
548 | TRIM( dum ), vmea_general%id_vm, '', & |
---|
549 | global_attribute, '', '' ) |
---|
550 | ! |
---|
551 | !-- Read a flag which indicates that also soil quantities are take at the |
---|
552 | !-- respective site (is part of the virtual measurement driver). |
---|
553 | CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ), & |
---|
554 | vmea_general%id_vm, '', global_attribute, & |
---|
555 | '', '' ) |
---|
556 | ! |
---|
557 | !-- Set flag for soil-sampling. |
---|
558 | IF ( soil_dum == 1 ) vmea(l)%soil_sampling = .TRUE. |
---|
559 | ! |
---|
560 | !--- Set logicals depending on the type of the measurement |
---|
561 | IF ( INDEX( vmea(l)%feature_type, type_tspr ) /= 0 ) THEN |
---|
562 | vmea(l)%timseries_profile = .TRUE. |
---|
563 | ELSEIF ( INDEX( vmea(l)%feature_type, type_ts ) /= 0 ) THEN |
---|
564 | vmea(l)%timseries = .TRUE. |
---|
565 | ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN |
---|
566 | vmea(l)%trajectory = .TRUE. |
---|
567 | ! |
---|
568 | !-- Give error message in case the type matches non of the pre-defined types. |
---|
569 | ELSE |
---|
570 | message_string = 'Attribue featureType = ' // & |
---|
571 | TRIM( vmea(l)%feature_type ) // & |
---|
572 | ' is not allowed.' |
---|
573 | CALL message( 'vm_init', 'PA0535', 1, 2, 0, 6, 0 ) |
---|
574 | ENDIF |
---|
575 | ! |
---|
576 | !-- Read string with all measured variables at this site |
---|
577 | measured_variables_file = '' |
---|
578 | CALL netcdf_data_input_var( measured_variables_file, & |
---|
579 | char_mv // TRIM( dum ), vmea_general%id_vm ) |
---|
580 | ! |
---|
581 | !-- Count the number of measured variables. Only count variables that match |
---|
582 | !-- with the allowed variables. |
---|
583 | !-- Please note, for some NetCDF interal reasons characters end with a NULL, |
---|
584 | !-- i.e. also empty characters contain a NULL. Therefore, check the strings |
---|
585 | !-- for a NULL to get the correct character length in order to compare |
---|
586 | !-- them with the list of allowed variables. |
---|
587 | vmea(l)%nmeas = 0 |
---|
588 | DO ll = 1, SIZE( measured_variables_file ) |
---|
589 | IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. & |
---|
590 | measured_variables_file(ll)(1:1) /= ' ') THEN |
---|
591 | ! |
---|
592 | !-- Obtain character length of the character |
---|
593 | len_char = 1 |
---|
594 | DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)& |
---|
595 | .AND. measured_variables_file(ll)(len_char:len_char) /= ' ' ) |
---|
596 | len_char = len_char + 1 |
---|
597 | ENDDO |
---|
598 | len_char = len_char - 1 |
---|
599 | ! |
---|
600 | !-- Now, compare the measured variable with the list of allowed |
---|
601 | !-- variables. |
---|
602 | DO lll= 1, SIZE( list_allowed_variables ) |
---|
603 | IF ( measured_variables_file(ll)(1:len_char) == & |
---|
604 | TRIM( list_allowed_variables(lll) ) ) THEN |
---|
605 | vmea(l)%nmeas = vmea(l)%nmeas + 1 |
---|
606 | measured_variables(vmea(l)%nmeas) = & |
---|
607 | measured_variables_file(ll)(1:len_char) |
---|
608 | ENDIF |
---|
609 | ENDDO |
---|
610 | ENDIF |
---|
611 | ENDDO |
---|
612 | |
---|
613 | |
---|
614 | ! |
---|
615 | !-- Allocate array for the measured variables names for the respective site. |
---|
616 | ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nmeas) ) |
---|
617 | |
---|
618 | DO ll = 1, vmea(l)%nmeas |
---|
619 | vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) ) |
---|
620 | ENDDO |
---|
621 | ! |
---|
622 | !-- In case of chemistry, check if species is considered in the modelled |
---|
623 | !-- chemistry mechanism. |
---|
624 | ! IF ( air_chemistry ) THEN |
---|
625 | ! DO ll = 1, vmea(l)%nmeas |
---|
626 | ! chem_include = .FALSE. |
---|
627 | ! DO n = 1, nvar |
---|
628 | ! IF ( TRIM( vmea(l)%measured_vars_name(ll) ) == & |
---|
629 | ! TRIM( chem_species(n)%name ) ) chem_include = .TRUE. |
---|
630 | ! ENDDO |
---|
631 | ! ! |
---|
632 | ! !-- Revise this. It should only check for chemistry variables and not for all! |
---|
633 | ! IF ( .NOT. chem_include ) THEN |
---|
634 | ! message_string = TRIM( vmea(l)%measured_vars_name(ll) ) // & |
---|
635 | ! ' is not considered in the modelled ' // & |
---|
636 | ! 'chemistry mechanism' |
---|
637 | ! CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 ) |
---|
638 | ! ENDIF |
---|
639 | ! ENDDO |
---|
640 | ! ENDIF |
---|
641 | ! |
---|
642 | !-- Read the UTM coordinates for the actual site. Based on the coordinates, |
---|
643 | !-- define the grid-index space on each subdomain where virtual measurements |
---|
644 | !-- should be taken. Note, the entire coordinate arrays will not be stored |
---|
645 | !-- as this would exceed memory requirements, particularly for trajectory |
---|
646 | !-- measurements. |
---|
647 | IF ( vmea(l)%nmeas > 0 ) THEN |
---|
648 | ! |
---|
649 | !-- For stationary measurements UTM coordinates are just one value and |
---|
650 | !-- its dimension is "station", while for mobile measurements UTM |
---|
651 | !-- coordinates are arrays depending on the number of trajectories and |
---|
652 | !-- time, according to (UC)2 standard. First, inquire dimension length |
---|
653 | !-- of the UTM coordinates. |
---|
654 | IF ( vmea(l)%trajectory ) THEN |
---|
655 | ! |
---|
656 | !-- For non-stationary measurements read the number of trajectories |
---|
657 | !-- and the number of time coordinates. |
---|
658 | CALL get_dimension_length( vmea_general%id_vm, & |
---|
659 | vmea(l)%ntraj, & |
---|
660 | "traj" // & |
---|
661 | TRIM( dum ) ) |
---|
662 | CALL get_dimension_length( vmea_general%id_vm, & |
---|
663 | dim_ntime, & |
---|
664 | "ntime" // & |
---|
665 | TRIM( dum ) ) |
---|
666 | ! |
---|
667 | !-- For stationary measurements the dimension for UTM and time |
---|
668 | !-- coordinates is 1. |
---|
669 | ELSE |
---|
670 | vmea(l)%ntraj = 1 |
---|
671 | dim_ntime = 1 |
---|
672 | ENDIF |
---|
673 | ! |
---|
674 | !- Allocate array which defines individual time frame for each |
---|
675 | !-- trajectory or station. |
---|
676 | ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) ) |
---|
677 | ! |
---|
678 | !-- Allocate temporary arrays for UTM and height coordinates. Note, |
---|
679 | !-- on file UTM coordinates might be 1D or 2D variables |
---|
680 | ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) ) |
---|
681 | ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) ) |
---|
682 | ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime) ) |
---|
683 | |
---|
684 | ALLOCATE( e_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) ) |
---|
685 | ALLOCATE( n_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) ) |
---|
686 | ! |
---|
687 | !-- Read _FillValue attributes of the coordinate dimensions. |
---|
688 | CALL netcdf_data_input_att( fill_eutm, char_fillvalue, & |
---|
689 | vmea_general%id_vm, '', & |
---|
690 | .NOT. global_attribute, '', & |
---|
691 | char_eutm // TRIM( dum ) ) |
---|
692 | CALL netcdf_data_input_att( fill_nutm, char_fillvalue, & |
---|
693 | vmea_general%id_vm, '', & |
---|
694 | .NOT. global_attribute, '', & |
---|
695 | char_nutm // TRIM( dum ) ) |
---|
696 | CALL netcdf_data_input_att( fill_zag, char_fillvalue, & |
---|
697 | vmea_general%id_vm, '', & |
---|
698 | .NOT. global_attribute, '', & |
---|
699 | char_zag // TRIM( dum ) ) |
---|
700 | ! |
---|
701 | !-- Read UTM and height coordinates coordinates for all trajectories and |
---|
702 | !-- times. |
---|
703 | IF ( vmea(l)%trajectory ) THEN |
---|
704 | CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), & |
---|
705 | vmea_general%id_vm, & |
---|
706 | 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) |
---|
707 | CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), & |
---|
708 | vmea_general%id_vm, & |
---|
709 | 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) |
---|
710 | CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), & |
---|
711 | vmea_general%id_vm, & |
---|
712 | 0, dim_ntime-1, 0, vmea(l)%ntraj-1 ) |
---|
713 | ELSE |
---|
714 | CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), & |
---|
715 | vmea_general%id_vm ) |
---|
716 | CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), & |
---|
717 | vmea_general%id_vm ) |
---|
718 | CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), & |
---|
719 | vmea_general%id_vm ) |
---|
720 | ENDIF |
---|
721 | ! |
---|
722 | !-- Based on UTM coordinates, check if the measurement station or parts |
---|
723 | !-- of the trajectory is on subdomain. This case, setup grid index space |
---|
724 | !-- sample these quantities. |
---|
725 | meas_flag = 0 |
---|
726 | DO t = 1, vmea(l)%ntraj |
---|
727 | ! |
---|
728 | !-- First, compute relative x- and y-coordinates with respect to the |
---|
729 | !-- lower-left origin of the model domain, which is the difference |
---|
730 | !-- between UTM coordinates. Note, if the origin is not correct, the |
---|
731 | !-- virtual sites will be misplaced. |
---|
732 | e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x |
---|
733 | n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y |
---|
734 | e_utm(t,1:dim_ntime) = COS( init_model%rotation_angle * pi / 180.0_wp ) & |
---|
735 | * e_utm_tmp(t,1:dim_ntime) & |
---|
736 | - SIN( init_model%rotation_angle * pi / 180.0_wp ) & |
---|
737 | * n_utm_tmp(t,1:dim_ntime) |
---|
738 | n_utm(t,1:dim_ntime) = SIN( init_model%rotation_angle * pi / 180.0_wp ) & |
---|
739 | * e_utm_tmp(t,1:dim_ntime) & |
---|
740 | + COS( init_model%rotation_angle * pi / 180.0_wp ) & |
---|
741 | * n_utm_tmp(t,1:dim_ntime) |
---|
742 | ! |
---|
743 | !-- Determine the individual time coordinate length for each station and |
---|
744 | !-- trajectory. This is required as several stations and trajectories |
---|
745 | !-- are merged into one file but they do not have the same number of |
---|
746 | !-- points in time, hence, missing values may occur and cannot be |
---|
747 | !-- processed further. This is actually a work-around for the specific |
---|
748 | !-- (UC)2 dataset, but it won't harm in anyway. |
---|
749 | vmea(l)%dim_t(t) = 0 |
---|
750 | DO n = 1, dim_ntime |
---|
751 | IF ( e_utm(t,n) /= fill_eutm .AND. & |
---|
752 | n_utm(t,n) /= fill_nutm .AND. & |
---|
753 | z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n |
---|
754 | ENDDO |
---|
755 | ! |
---|
756 | !-- Compute grid indices relative to origin and check if these are |
---|
757 | !-- on the subdomain. Note, virtual measurements will be taken also |
---|
758 | !-- at grid points surrounding the station, hence, check also for |
---|
759 | !-- these grid points. |
---|
760 | DO n = 1, vmea(l)%dim_t(t) |
---|
761 | is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp ) |
---|
762 | js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp ) |
---|
763 | ! |
---|
764 | !-- Is the observation point on subdomain? |
---|
765 | on_pe = ( is >= nxl .AND. is <= nxr .AND. & |
---|
766 | js >= nys .AND. js <= nyn ) |
---|
767 | ! |
---|
768 | !-- Check if observation coordinate is on subdomain |
---|
769 | IF ( on_pe ) THEN |
---|
770 | ! |
---|
771 | !-- Determine vertical index which correspond to the observation |
---|
772 | !-- height. |
---|
773 | ksurf = topo_top_ind(js,is,0) |
---|
774 | ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1 |
---|
775 | ! |
---|
776 | !-- Set mask array at the observation coordinates. Also, flag the |
---|
777 | !-- surrounding coordinate points, but first check whether the |
---|
778 | !-- surrounding coordinate points are on the subdomain. |
---|
779 | kl = MERGE( ks-1, ks, ks-1 >= nzb .AND. ks-1 >= ksurf ) |
---|
780 | ku = MERGE( ks+1, ks, ks+1 < nzt+1 ) |
---|
781 | |
---|
782 | DO i = is-1, is+1 |
---|
783 | DO j = js-1, js+1 |
---|
784 | DO k = kl, ku |
---|
785 | meas_flag(k,j,i) = MERGE( & |
---|
786 | IBSET( meas_flag(k,j,i), 0 ), & |
---|
787 | 0, & |
---|
788 | BTEST( wall_flags_static_0(k,j,i), 0 ) & |
---|
789 | ) |
---|
790 | ENDDO |
---|
791 | ENDDO |
---|
792 | ENDDO |
---|
793 | ENDIF |
---|
794 | ENDDO |
---|
795 | |
---|
796 | ENDDO |
---|
797 | ! |
---|
798 | !-- Based on the flag array count the number of of sampling coordinates. |
---|
799 | !-- Please note, sampling coordinates in atmosphere and soil may be |
---|
800 | !-- different, as within the soil all levels will be measured. |
---|
801 | !-- Hence, count individually. Start with atmoshere. |
---|
802 | ns = 0 |
---|
803 | DO i = nxl-1, nxr+1 |
---|
804 | DO j = nys-1, nyn+1 |
---|
805 | DO k = nzb, nzt+1 |
---|
806 | ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) ) |
---|
807 | ENDDO |
---|
808 | ENDDO |
---|
809 | ENDDO |
---|
810 | |
---|
811 | ! |
---|
812 | !-- Store number of observation points on subdomain and allocate index |
---|
813 | !-- arrays as well as array containing height information. |
---|
814 | vmea(l)%ns = ns |
---|
815 | |
---|
816 | ALLOCATE( vmea(l)%i(1:vmea(l)%ns) ) |
---|
817 | ALLOCATE( vmea(l)%j(1:vmea(l)%ns) ) |
---|
818 | ALLOCATE( vmea(l)%k(1:vmea(l)%ns) ) |
---|
819 | ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) ) |
---|
820 | ! |
---|
821 | !-- Based on the flag array store the grid indices which correspond to |
---|
822 | !-- the observation coordinates. |
---|
823 | ns = 0 |
---|
824 | DO i = nxl-1, nxr+1 |
---|
825 | DO j = nys-1, nyn+1 |
---|
826 | DO k = nzb, nzt+1 |
---|
827 | IF ( BTEST( meas_flag(k,j,i), 0 ) ) THEN |
---|
828 | ns = ns + 1 |
---|
829 | vmea(l)%i(ns) = i |
---|
830 | vmea(l)%j(ns) = j |
---|
831 | vmea(l)%k(ns) = k |
---|
832 | vmea(l)%z_ag(ns) = zu(k) - zw(topo_top_ind(j,i,0)) |
---|
833 | ENDIF |
---|
834 | ENDDO |
---|
835 | ENDDO |
---|
836 | ENDDO |
---|
837 | ! |
---|
838 | !-- Same for the soil. Based on the flag array, count the number of |
---|
839 | !-- sampling coordinates in soil. Sample at all soil levels in this case. |
---|
840 | IF ( vmea(l)%soil_sampling ) THEN |
---|
841 | DO i = nxl, nxr |
---|
842 | DO j = nys, nyn |
---|
843 | IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN |
---|
844 | IF ( surf_lsm_h%start_index(j,i) <= & |
---|
845 | surf_lsm_h%end_index(j,i) ) THEN |
---|
846 | vmea(l)%ns_soil = vmea(l)%ns_soil + & |
---|
847 | nzt_soil - nzb_soil + 1 |
---|
848 | ENDIF |
---|
849 | IF ( surf_usm_h%start_index(j,i) <= & |
---|
850 | surf_usm_h%end_index(j,i) ) THEN |
---|
851 | vmea(l)%ns_soil = vmea(l)%ns_soil + & |
---|
852 | nzt_wall - nzb_wall + 1 |
---|
853 | ENDIF |
---|
854 | ENDIF |
---|
855 | ENDDO |
---|
856 | ENDDO |
---|
857 | ENDIF |
---|
858 | ! |
---|
859 | !-- Allocate index arrays as well as array containing height information |
---|
860 | !-- for soil. |
---|
861 | IF ( vmea(l)%soil_sampling ) THEN |
---|
862 | ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) ) |
---|
863 | ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) ) |
---|
864 | ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) ) |
---|
865 | ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) ) |
---|
866 | ENDIF |
---|
867 | ! |
---|
868 | !-- For soil, store the grid indices. |
---|
869 | ns = 0 |
---|
870 | IF ( vmea(l)%soil_sampling ) THEN |
---|
871 | DO i = nxl, nxr |
---|
872 | DO j = nys, nyn |
---|
873 | IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN |
---|
874 | IF ( surf_lsm_h%start_index(j,i) <= & |
---|
875 | surf_lsm_h%end_index(j,i) ) THEN |
---|
876 | m = surf_lsm_h%start_index(j,i) |
---|
877 | DO k = nzb_soil, nzt_soil |
---|
878 | ns = ns + 1 |
---|
879 | vmea(l)%i_soil(ns) = i |
---|
880 | vmea(l)%j_soil(ns) = j |
---|
881 | vmea(l)%k_soil(ns) = k |
---|
882 | vmea(l)%depth(ns) = zs(k) |
---|
883 | ENDDO |
---|
884 | ENDIF |
---|
885 | |
---|
886 | IF ( surf_usm_h%start_index(j,i) <= & |
---|
887 | surf_usm_h%end_index(j,i) ) THEN |
---|
888 | m = surf_usm_h%start_index(j,i) |
---|
889 | DO k = nzb_wall, nzt_wall |
---|
890 | ns = ns + 1 |
---|
891 | vmea(l)%i_soil(ns) = i |
---|
892 | vmea(l)%j_soil(ns) = j |
---|
893 | vmea(l)%k_soil(ns) = k |
---|
894 | vmea(l)%depth(ns) = surf_usm_h%zw(k,m) |
---|
895 | ENDDO |
---|
896 | ENDIF |
---|
897 | ENDIF |
---|
898 | ENDDO |
---|
899 | ENDDO |
---|
900 | ENDIF |
---|
901 | ! |
---|
902 | !-- Allocate array to save the sampled values. |
---|
903 | ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) ) |
---|
904 | |
---|
905 | IF ( vmea(l)%soil_sampling ) & |
---|
906 | ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil, & |
---|
907 | 1:vmea(l)%nmeas) ) |
---|
908 | ! |
---|
909 | !-- Initialize with _FillValues |
---|
910 | vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) = vmea(l)%fillout |
---|
911 | IF ( vmea(l)%soil_sampling ) & |
---|
912 | vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nmeas) = & |
---|
913 | vmea(l)%fillout |
---|
914 | ! |
---|
915 | !-- Deallocate temporary coordinate arrays |
---|
916 | IF ( ALLOCATED( e_utm ) ) DEALLOCATE( e_utm ) |
---|
917 | IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm ) |
---|
918 | IF ( ALLOCATED( e_utm_tmp ) ) DEALLOCATE( e_utm_tmp ) |
---|
919 | IF ( ALLOCATED( n_utm_tmp ) ) DEALLOCATE( n_utm_tmp ) |
---|
920 | IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm ) |
---|
921 | IF ( ALLOCATED( z_ag ) ) DEALLOCATE( z_ag ) |
---|
922 | IF ( ALLOCATED( z_ag ) ) DEALLOCATE( vmea(l)%dim_t ) |
---|
923 | ENDIF |
---|
924 | ENDDO |
---|
925 | ! |
---|
926 | !-- Close input file for virtual measurements. Therefore, just call |
---|
927 | !-- the read attribute routine with the "close" option. |
---|
928 | CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, & |
---|
929 | vmea_general%id_vm, '', & |
---|
930 | global_attribute, 'close', '' ) |
---|
931 | ! |
---|
932 | !-- Sum-up the number of observation coordiates, for atmosphere first. |
---|
933 | !-- This is actually only required for data output. |
---|
934 | ALLOCATE( ns_all(1:vmea_general%nvm) ) |
---|
935 | ns_all = 0 |
---|
936 | #if defined( __parallel ) |
---|
937 | CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER, & |
---|
938 | MPI_SUM, comm2d, ierr ) |
---|
939 | #else |
---|
940 | ns_all(:) = vmea(:)%ns |
---|
941 | #endif |
---|
942 | vmea(:)%ns_tot = ns_all(:) |
---|
943 | ! |
---|
944 | !-- Now for soil |
---|
945 | ns_all = 0 |
---|
946 | #if defined( __parallel ) |
---|
947 | CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm, & |
---|
948 | MPI_INTEGER, MPI_SUM, comm2d, ierr ) |
---|
949 | #else |
---|
950 | ns_all(:) = vmea(:)%ns_soil |
---|
951 | #endif |
---|
952 | vmea(:)%ns_soil_tot = ns_all(:) |
---|
953 | |
---|
954 | DEALLOCATE( ns_all ) |
---|
955 | ! |
---|
956 | !-- Dellocate flag array |
---|
957 | DEALLOCATE( meas_flag ) |
---|
958 | ! |
---|
959 | !-- Initialize binary data output of virtual measurements. |
---|
960 | !-- Open binary output file. |
---|
961 | CALL check_open( 27 ) |
---|
962 | ! |
---|
963 | !-- Output header information. |
---|
964 | CALL vm_data_output |
---|
965 | |
---|
966 | END SUBROUTINE vm_init |
---|
967 | |
---|
968 | |
---|
969 | !------------------------------------------------------------------------------! |
---|
970 | ! Description: |
---|
971 | ! ------------ |
---|
972 | !> Binary data output. |
---|
973 | !------------------------------------------------------------------------------! |
---|
974 | SUBROUTINE vm_data_output |
---|
975 | |
---|
976 | USE pegrid |
---|
977 | |
---|
978 | IMPLICIT NONE |
---|
979 | |
---|
980 | INTEGER(iwp) :: i !< running index over IO blocks |
---|
981 | INTEGER(iwp) :: l !< running index over all stations |
---|
982 | INTEGER(iwp) :: n !< running index over all measured variables at a station |
---|
983 | ! |
---|
984 | !-- Header output on each PE |
---|
985 | IF ( init ) THEN |
---|
986 | |
---|
987 | DO i = 0, io_blocks-1 |
---|
988 | IF ( i == io_group ) THEN |
---|
989 | WRITE ( 27 ) 'number of measurements ' |
---|
990 | WRITE ( 27 ) vmea_general%nvm |
---|
991 | |
---|
992 | DO l = 1, vmea_general%nvm |
---|
993 | WRITE ( 27 ) 'site ' |
---|
994 | WRITE ( 27 ) vmea(l)%site |
---|
995 | WRITE ( 27 ) 'file ' |
---|
996 | WRITE ( 27 ) vmea(l)%filename_original |
---|
997 | WRITE ( 27 ) 'feature_type ' |
---|
998 | WRITE ( 27 ) vmea(l)%feature_type |
---|
999 | WRITE ( 27 ) 'origin_x_obs ' |
---|
1000 | WRITE ( 27 ) vmea(l)%origin_x_obs |
---|
1001 | WRITE ( 27 ) 'origin_y_obs ' |
---|
1002 | WRITE ( 27 ) vmea(l)%origin_y_obs |
---|
1003 | WRITE ( 27 ) 'total number of observation points' |
---|
1004 | WRITE ( 27 ) vmea(l)%ns_tot |
---|
1005 | WRITE ( 27 ) 'number of measured variables ' |
---|
1006 | WRITE ( 27 ) vmea(l)%nmeas |
---|
1007 | WRITE ( 27 ) 'variables ' |
---|
1008 | WRITE ( 27 ) vmea(l)%measured_vars_name(:) |
---|
1009 | WRITE ( 27 ) 'number of observation points ' |
---|
1010 | WRITE ( 27 ) vmea(l)%ns |
---|
1011 | WRITE ( 27 ) 'E_UTM ' |
---|
1012 | WRITE ( 27 ) init_model%origin_x & |
---|
1013 | + REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx & |
---|
1014 | * COS( init_model%rotation_angle * pi / 180.0_wp ) & |
---|
1015 | + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy & |
---|
1016 | * SIN( init_model%rotation_angle * pi / 180.0_wp ) |
---|
1017 | WRITE ( 27 ) 'N_UTM ' |
---|
1018 | WRITE ( 27 ) init_model%origin_y & |
---|
1019 | - REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx & |
---|
1020 | * SIN( init_model%rotation_angle * pi / 180.0_wp ) & |
---|
1021 | + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy & |
---|
1022 | * COS( init_model%rotation_angle * pi / 180.0_wp ) |
---|
1023 | WRITE ( 27 ) 'Z_AG ' |
---|
1024 | WRITE ( 27 ) vmea(l)%z_ag(1:vmea(l)%ns) |
---|
1025 | WRITE ( 27 ) 'soil sampling ' |
---|
1026 | WRITE ( 27 ) MERGE( 'yes ', & |
---|
1027 | 'no ', & |
---|
1028 | vmea(l)%soil_sampling ) |
---|
1029 | |
---|
1030 | IF ( vmea(l)%soil_sampling ) THEN |
---|
1031 | WRITE ( 27 ) 'total number of soil points ' |
---|
1032 | WRITE ( 27 ) vmea(l)%ns_soil_tot |
---|
1033 | WRITE ( 27 ) 'number of soil points ' |
---|
1034 | WRITE ( 27 ) vmea(l)%ns_soil |
---|
1035 | WRITE ( 27 ) 'E_UTM soil ' |
---|
1036 | WRITE ( 27 ) init_model%origin_x & |
---|
1037 | + REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, & |
---|
1038 | KIND = wp ) * dx & |
---|
1039 | * COS( init_model%rotation_angle * pi / 180.0_wp ) & |
---|
1040 | + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, & |
---|
1041 | KIND = wp ) * dy & |
---|
1042 | * SIN( init_model%rotation_angle * pi / 180.0_wp ) |
---|
1043 | WRITE ( 27 ) 'N_UTM soil ' |
---|
1044 | WRITE ( 27 ) init_model%origin_y & |
---|
1045 | - REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, & |
---|
1046 | KIND = wp ) * dx & |
---|
1047 | * SIN( init_model%rotation_angle * pi / 180.0_wp ) & |
---|
1048 | + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, & |
---|
1049 | KIND = wp ) * dy & |
---|
1050 | * COS( init_model%rotation_angle * pi / 180.0_wp ) |
---|
1051 | WRITE ( 27 ) 'DEPTH ' |
---|
1052 | WRITE ( 27 ) vmea(l)%depth(1:vmea(l)%ns_soil) |
---|
1053 | ENDIF |
---|
1054 | ENDDO |
---|
1055 | |
---|
1056 | ENDIF |
---|
1057 | ENDDO |
---|
1058 | |
---|
1059 | #if defined( __parallel ) |
---|
1060 | CALL MPI_BARRIER( comm2d, ierr ) |
---|
1061 | #endif |
---|
1062 | ! |
---|
1063 | !-- After header information is written, set control flag to .FALSE. |
---|
1064 | init = .FALSE. |
---|
1065 | ! |
---|
1066 | !-- Data output at each measurement timestep on each PE |
---|
1067 | ELSE |
---|
1068 | DO i = 0, io_blocks-1 |
---|
1069 | |
---|
1070 | IF ( i == io_group ) THEN |
---|
1071 | WRITE( 27 ) 'output time ' |
---|
1072 | WRITE( 27 ) time_since_reference_point |
---|
1073 | DO l = 1, vmea_general%nvm |
---|
1074 | ! |
---|
1075 | !-- Skip binary writing if no observation points are defined on PE |
---|
1076 | IF ( vmea(l)%ns < 1 .AND. vmea(l)%ns_soil < 1) CYCLE |
---|
1077 | DO n = 1, vmea(l)%nmeas |
---|
1078 | WRITE( 27 ) vmea(l)%measured_vars_name(n) |
---|
1079 | IF ( vmea(l)%soil_sampling .AND. & |
---|
1080 | ANY( TRIM( vmea(l)%measured_vars_name(n)) == & |
---|
1081 | soil_vars ) ) THEN |
---|
1082 | WRITE( 27 ) vmea(l)%measured_vars_soil(:,n) |
---|
1083 | ELSE |
---|
1084 | WRITE( 27 ) vmea(l)%measured_vars(:,n) |
---|
1085 | ENDIF |
---|
1086 | ENDDO |
---|
1087 | |
---|
1088 | ENDDO |
---|
1089 | ENDIF |
---|
1090 | ENDDO |
---|
1091 | #if defined( __parallel ) |
---|
1092 | CALL MPI_BARRIER( comm2d, ierr ) |
---|
1093 | #endif |
---|
1094 | ENDIF |
---|
1095 | |
---|
1096 | END SUBROUTINE vm_data_output |
---|
1097 | |
---|
1098 | |
---|
1099 | !------------------------------------------------------------------------------! |
---|
1100 | ! Description: |
---|
1101 | ! ------------ |
---|
1102 | !> Write end-of-file statement as last action. |
---|
1103 | !------------------------------------------------------------------------------! |
---|
1104 | SUBROUTINE vm_last_actions |
---|
1105 | |
---|
1106 | USE pegrid |
---|
1107 | |
---|
1108 | IMPLICIT NONE |
---|
1109 | |
---|
1110 | INTEGER(iwp) :: i !< running index over IO blocks |
---|
1111 | |
---|
1112 | DO i = 0, io_blocks-1 |
---|
1113 | IF ( i == io_group ) THEN |
---|
1114 | WRITE( 27 ) 'EOF ' |
---|
1115 | ENDIF |
---|
1116 | ENDDO |
---|
1117 | #if defined( __parallel ) |
---|
1118 | CALL MPI_BARRIER( comm2d, ierr ) |
---|
1119 | #endif |
---|
1120 | ! |
---|
1121 | !-- Close binary file |
---|
1122 | CALL close_file( 27 ) |
---|
1123 | |
---|
1124 | END SUBROUTINE vm_last_actions |
---|
1125 | |
---|
1126 | !------------------------------------------------------------------------------! |
---|
1127 | ! Description: |
---|
1128 | ! ------------ |
---|
1129 | !> Sampling of the actual quantities along the observation coordinates |
---|
1130 | !------------------------------------------------------------------------------! |
---|
1131 | SUBROUTINE vm_sampling |
---|
1132 | |
---|
1133 | USE arrays_3d, & |
---|
1134 | ONLY: exner, pt, q, u, v, w |
---|
1135 | |
---|
1136 | USE radiation_model_mod, & |
---|
1137 | ONLY: radiation |
---|
1138 | |
---|
1139 | USE surface_mod, & |
---|
1140 | ONLY: surf_def_h, surf_lsm_h, surf_usm_h |
---|
1141 | |
---|
1142 | IMPLICIT NONE |
---|
1143 | |
---|
1144 | INTEGER(iwp) :: i !< grid index in x-direction |
---|
1145 | INTEGER(iwp) :: j !< grid index in y-direction |
---|
1146 | INTEGER(iwp) :: k !< grid index in z-direction |
---|
1147 | INTEGER(iwp) :: ind_chem !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming |
---|
1148 | INTEGER(iwp) :: l !< running index over the number of stations |
---|
1149 | INTEGER(iwp) :: m !< running index over all virtual observation coordinates |
---|
1150 | INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate |
---|
1151 | INTEGER(iwp) :: n !< running index over all measured variables at a station |
---|
1152 | INTEGER(iwp) :: nn !< running index over the number of chemcal species |
---|
1153 | |
---|
1154 | LOGICAL :: match_lsm !< flag indicating natural-type surface |
---|
1155 | LOGICAL :: match_usm !< flag indicating urban-type surface |
---|
1156 | ! |
---|
1157 | !-- Loop over all sites. |
---|
1158 | DO l = 1, vmea_general%nvm |
---|
1159 | ! |
---|
1160 | !-- At the beginning, set _FillValues |
---|
1161 | IF ( ALLOCATED( vmea(l)%measured_vars ) ) & |
---|
1162 | vmea(l)%measured_vars = vmea(l)%fillout |
---|
1163 | IF ( ALLOCATED( vmea(l)%measured_vars_soil ) ) & |
---|
1164 | vmea(l)%measured_vars_soil = vmea(l)%fillout |
---|
1165 | ! |
---|
1166 | !-- Loop over all variables measured at this site. |
---|
1167 | DO n = 1, vmea(l)%nmeas |
---|
1168 | |
---|
1169 | SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) ) |
---|
1170 | |
---|
1171 | CASE ( 'theta' ) |
---|
1172 | IF ( .NOT. neutral ) THEN |
---|
1173 | DO m = 1, vmea(l)%ns |
---|
1174 | k = vmea(l)%k(m) |
---|
1175 | j = vmea(l)%j(m) |
---|
1176 | i = vmea(l)%i(m) |
---|
1177 | vmea(l)%measured_vars(m,n) = pt(k,j,i) |
---|
1178 | ENDDO |
---|
1179 | ENDIF |
---|
1180 | |
---|
1181 | CASE ( 'ta' ) |
---|
1182 | IF ( .NOT. neutral ) THEN |
---|
1183 | DO m = 1, vmea(l)%ns |
---|
1184 | k = vmea(l)%k(m) |
---|
1185 | j = vmea(l)%j(m) |
---|
1186 | i = vmea(l)%i(m) |
---|
1187 | vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k ) |
---|
1188 | ENDDO |
---|
1189 | ENDIF |
---|
1190 | |
---|
1191 | CASE ( 't_va' ) |
---|
1192 | |
---|
1193 | CASE ( 'hus', 'haa' ) |
---|
1194 | IF ( humidity ) THEN |
---|
1195 | DO m = 1, vmea(l)%ns |
---|
1196 | k = vmea(l)%k(m) |
---|
1197 | j = vmea(l)%j(m) |
---|
1198 | i = vmea(l)%i(m) |
---|
1199 | vmea(l)%measured_vars(m,n) = q(k,j,i) |
---|
1200 | ENDDO |
---|
1201 | ENDIF |
---|
1202 | |
---|
1203 | CASE ( 'u', 'ua' ) |
---|
1204 | DO m = 1, vmea(l)%ns |
---|
1205 | k = vmea(l)%k(m) |
---|
1206 | j = vmea(l)%j(m) |
---|
1207 | i = vmea(l)%i(m) |
---|
1208 | vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) |
---|
1209 | ENDDO |
---|
1210 | |
---|
1211 | CASE ( 'v', 'va' ) |
---|
1212 | DO m = 1, vmea(l)%ns |
---|
1213 | k = vmea(l)%k(m) |
---|
1214 | j = vmea(l)%j(m) |
---|
1215 | i = vmea(l)%i(m) |
---|
1216 | vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) |
---|
1217 | ENDDO |
---|
1218 | |
---|
1219 | CASE ( 'w' ) |
---|
1220 | DO m = 1, vmea(l)%ns |
---|
1221 | k = MAX ( 1, vmea(l)%k(m) ) |
---|
1222 | j = vmea(l)%j(m) |
---|
1223 | i = vmea(l)%i(m) |
---|
1224 | vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) |
---|
1225 | ENDDO |
---|
1226 | |
---|
1227 | CASE ( 'wspeed' ) |
---|
1228 | DO m = 1, vmea(l)%ns |
---|
1229 | k = vmea(l)%k(m) |
---|
1230 | j = vmea(l)%j(m) |
---|
1231 | i = vmea(l)%i(m) |
---|
1232 | vmea(l)%measured_vars(m,n) = SQRT( & |
---|
1233 | ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + & |
---|
1234 | ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 & |
---|
1235 | ) |
---|
1236 | ENDDO |
---|
1237 | |
---|
1238 | CASE ( 'wdir' ) |
---|
1239 | DO m = 1, vmea(l)%ns |
---|
1240 | k = vmea(l)%k(m) |
---|
1241 | j = vmea(l)%j(m) |
---|
1242 | i = vmea(l)%i(m) |
---|
1243 | |
---|
1244 | vmea(l)%measured_vars(m,n) = ATAN2( & |
---|
1245 | - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ), & |
---|
1246 | - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) & |
---|
1247 | ) * 180.0_wp / pi |
---|
1248 | ENDDO |
---|
1249 | |
---|
1250 | CASE ( 'utheta' ) |
---|
1251 | DO m = 1, vmea(l)%ns |
---|
1252 | k = vmea(l)%k(m) |
---|
1253 | j = vmea(l)%j(m) |
---|
1254 | i = vmea(l)%i(m) |
---|
1255 | vmea(l)%measured_vars(m,n) = 0.5_wp * & |
---|
1256 | ( u(k,j,i) + u(k,j,i+1) ) * & |
---|
1257 | pt(k,j,i) |
---|
1258 | ENDDO |
---|
1259 | |
---|
1260 | CASE ( 'vtheta' ) |
---|
1261 | DO m = 1, vmea(l)%ns |
---|
1262 | k = vmea(l)%k(m) |
---|
1263 | j = vmea(l)%j(m) |
---|
1264 | i = vmea(l)%i(m) |
---|
1265 | vmea(l)%measured_vars(m,n) = 0.5_wp * & |
---|
1266 | ( v(k,j,i) + v(k,j+1,i) ) * & |
---|
1267 | pt(k,j,i) |
---|
1268 | ENDDO |
---|
1269 | |
---|
1270 | CASE ( 'wtheta' ) |
---|
1271 | DO m = 1, vmea(l)%ns |
---|
1272 | k = MAX ( 1, vmea(l)%k(m) ) |
---|
1273 | j = vmea(l)%j(m) |
---|
1274 | i = vmea(l)%i(m) |
---|
1275 | vmea(l)%measured_vars(m,n) = 0.5_wp * & |
---|
1276 | ( w(k-1,j,i) + w(k,j,i) ) * & |
---|
1277 | pt(k,j,i) |
---|
1278 | ENDDO |
---|
1279 | |
---|
1280 | CASE ( 'uw' ) |
---|
1281 | DO m = 1, vmea(l)%ns |
---|
1282 | k = MAX ( 1, vmea(l)%k(m) ) |
---|
1283 | j = vmea(l)%j(m) |
---|
1284 | i = vmea(l)%i(m) |
---|
1285 | vmea(l)%measured_vars(m,n) = 0.25_wp * & |
---|
1286 | ( w(k-1,j,i) + w(k,j,i) ) * & |
---|
1287 | ( u(k,j,i) + u(k,j,i+1) ) |
---|
1288 | ENDDO |
---|
1289 | |
---|
1290 | CASE ( 'vw' ) |
---|
1291 | DO m = 1, vmea(l)%ns |
---|
1292 | k = MAX ( 1, vmea(l)%k(m) ) |
---|
1293 | j = vmea(l)%j(m) |
---|
1294 | i = vmea(l)%i(m) |
---|
1295 | vmea(l)%measured_vars(m,n) = 0.25_wp * & |
---|
1296 | ( w(k-1,j,i) + w(k,j,i) ) * & |
---|
1297 | ( v(k,j,i) + v(k,j+1,i) ) |
---|
1298 | ENDDO |
---|
1299 | |
---|
1300 | CASE ( 'uv' ) |
---|
1301 | DO m = 1, vmea(l)%ns |
---|
1302 | k = MAX ( 1, vmea(l)%k(m) ) |
---|
1303 | j = vmea(l)%j(m) |
---|
1304 | i = vmea(l)%i(m) |
---|
1305 | vmea(l)%measured_vars(m,n) = 0.25_wp * & |
---|
1306 | ( u(k,j,i) + u(k,j,i+1) ) * & |
---|
1307 | ( v(k,j,i) + v(k,j+1,i) ) |
---|
1308 | ENDDO |
---|
1309 | ! |
---|
1310 | !-- List of variables may need extension. |
---|
1311 | CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2', & |
---|
1312 | 'tro3' ) |
---|
1313 | IF ( air_chemistry ) THEN |
---|
1314 | ! |
---|
1315 | !-- First, search for the measured variable in the chem_vars |
---|
1316 | !-- list, in order to get the internal name of the variable. |
---|
1317 | DO nn = 1, UBOUND( chem_vars, 2 ) |
---|
1318 | IF ( TRIM( vmea(l)%measured_vars_name(m) ) == & |
---|
1319 | TRIM( chem_vars(0,nn) ) ) ind_chem = nn |
---|
1320 | ENDDO |
---|
1321 | ! |
---|
1322 | !-- Run loop over all chemical species, if the measured |
---|
1323 | !-- variable matches the interal name, sample the variable. |
---|
1324 | DO nn = 1, nvar |
---|
1325 | IF ( TRIM( chem_vars(1,ind_chem) ) == & |
---|
1326 | TRIM( chem_species(nn)%name ) ) THEN |
---|
1327 | DO m = 1, vmea(l)%ns |
---|
1328 | k = vmea(l)%k(m) |
---|
1329 | j = vmea(l)%j(m) |
---|
1330 | i = vmea(l)%i(m) |
---|
1331 | vmea(l)%measured_vars(m,n) = & |
---|
1332 | chem_species(nn)%conc(k,j,i) |
---|
1333 | ENDDO |
---|
1334 | ENDIF |
---|
1335 | ENDDO |
---|
1336 | ENDIF |
---|
1337 | |
---|
1338 | CASE ( 'us' ) |
---|
1339 | DO m = 1, vmea(l)%ns |
---|
1340 | ! |
---|
1341 | !-- Surface data is only available on inner subdomains, not |
---|
1342 | !-- on ghost points. Hence, limit the indices. |
---|
1343 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1344 | j = MERGE( j , nyn, j > nyn ) |
---|
1345 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1346 | i = MERGE( i , nxr, i > nxr ) |
---|
1347 | |
---|
1348 | DO mm = surf_def_h(0)%start_index(j,i), & |
---|
1349 | surf_def_h(0)%end_index(j,i) |
---|
1350 | vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm) |
---|
1351 | ENDDO |
---|
1352 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1353 | surf_lsm_h%end_index(j,i) |
---|
1354 | vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm) |
---|
1355 | ENDDO |
---|
1356 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1357 | surf_usm_h%end_index(j,i) |
---|
1358 | vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm) |
---|
1359 | ENDDO |
---|
1360 | ENDDO |
---|
1361 | |
---|
1362 | CASE ( 'ts' ) |
---|
1363 | DO m = 1, vmea(l)%ns |
---|
1364 | ! |
---|
1365 | !-- Surface data is only available on inner subdomains, not |
---|
1366 | !-- on ghost points. Hence, limit the indices. |
---|
1367 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1368 | j = MERGE( j , nyn, j > nyn ) |
---|
1369 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1370 | i = MERGE( i , nxr, i > nxr ) |
---|
1371 | |
---|
1372 | DO mm = surf_def_h(0)%start_index(j,i), & |
---|
1373 | surf_def_h(0)%end_index(j,i) |
---|
1374 | vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm) |
---|
1375 | ENDDO |
---|
1376 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1377 | surf_lsm_h%end_index(j,i) |
---|
1378 | vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm) |
---|
1379 | ENDDO |
---|
1380 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1381 | surf_usm_h%end_index(j,i) |
---|
1382 | vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm) |
---|
1383 | ENDDO |
---|
1384 | ENDDO |
---|
1385 | |
---|
1386 | CASE ( 'hfls' ) |
---|
1387 | DO m = 1, vmea(l)%ns |
---|
1388 | ! |
---|
1389 | !-- Surface data is only available on inner subdomains, not |
---|
1390 | !-- on ghost points. Hence, limit the indices. |
---|
1391 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1392 | j = MERGE( j , nyn, j > nyn ) |
---|
1393 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1394 | i = MERGE( i , nxr, i > nxr ) |
---|
1395 | |
---|
1396 | DO mm = surf_def_h(0)%start_index(j,i), & |
---|
1397 | surf_def_h(0)%end_index(j,i) |
---|
1398 | vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm) |
---|
1399 | ENDDO |
---|
1400 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1401 | surf_lsm_h%end_index(j,i) |
---|
1402 | vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm) |
---|
1403 | ENDDO |
---|
1404 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1405 | surf_usm_h%end_index(j,i) |
---|
1406 | vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm) |
---|
1407 | ENDDO |
---|
1408 | ENDDO |
---|
1409 | |
---|
1410 | CASE ( 'hfss' ) |
---|
1411 | DO m = 1, vmea(l)%ns |
---|
1412 | ! |
---|
1413 | !-- Surface data is only available on inner subdomains, not |
---|
1414 | !-- on ghost points. Hence, limit the indices. |
---|
1415 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1416 | j = MERGE( j , nyn, j > nyn ) |
---|
1417 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1418 | i = MERGE( i , nxr, i > nxr ) |
---|
1419 | |
---|
1420 | DO mm = surf_def_h(0)%start_index(j,i), & |
---|
1421 | surf_def_h(0)%end_index(j,i) |
---|
1422 | vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm) |
---|
1423 | ENDDO |
---|
1424 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1425 | surf_lsm_h%end_index(j,i) |
---|
1426 | vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm) |
---|
1427 | ENDDO |
---|
1428 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1429 | surf_usm_h%end_index(j,i) |
---|
1430 | vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm) |
---|
1431 | ENDDO |
---|
1432 | ENDDO |
---|
1433 | |
---|
1434 | CASE ( 'rnds' ) |
---|
1435 | IF ( radiation ) THEN |
---|
1436 | DO m = 1, vmea(l)%ns |
---|
1437 | ! |
---|
1438 | !-- Surface data is only available on inner subdomains, not |
---|
1439 | !-- on ghost points. Hence, limit the indices. |
---|
1440 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1441 | j = MERGE( j , nyn, j > nyn ) |
---|
1442 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1443 | i = MERGE( i , nxr, i > nxr ) |
---|
1444 | |
---|
1445 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1446 | surf_lsm_h%end_index(j,i) |
---|
1447 | vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm) |
---|
1448 | ENDDO |
---|
1449 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1450 | surf_usm_h%end_index(j,i) |
---|
1451 | vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm) |
---|
1452 | ENDDO |
---|
1453 | ENDDO |
---|
1454 | ENDIF |
---|
1455 | |
---|
1456 | CASE ( 'rsus' ) |
---|
1457 | IF ( radiation ) THEN |
---|
1458 | DO m = 1, vmea(l)%ns |
---|
1459 | ! |
---|
1460 | !-- Surface data is only available on inner subdomains, not |
---|
1461 | !-- on ghost points. Hence, limit the indices. |
---|
1462 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1463 | j = MERGE( j , nyn, j > nyn ) |
---|
1464 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1465 | i = MERGE( i , nxr, i > nxr ) |
---|
1466 | |
---|
1467 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1468 | surf_lsm_h%end_index(j,i) |
---|
1469 | vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm) |
---|
1470 | ENDDO |
---|
1471 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1472 | surf_usm_h%end_index(j,i) |
---|
1473 | vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm) |
---|
1474 | ENDDO |
---|
1475 | ENDDO |
---|
1476 | ENDIF |
---|
1477 | |
---|
1478 | CASE ( 'rsds' ) |
---|
1479 | IF ( radiation ) THEN |
---|
1480 | DO m = 1, vmea(l)%ns |
---|
1481 | ! |
---|
1482 | !-- Surface data is only available on inner subdomains, not |
---|
1483 | !-- on ghost points. Hence, limit the indices. |
---|
1484 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1485 | j = MERGE( j , nyn, j > nyn ) |
---|
1486 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1487 | i = MERGE( i , nxr, i > nxr ) |
---|
1488 | |
---|
1489 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1490 | surf_lsm_h%end_index(j,i) |
---|
1491 | vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm) |
---|
1492 | ENDDO |
---|
1493 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1494 | surf_usm_h%end_index(j,i) |
---|
1495 | vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm) |
---|
1496 | ENDDO |
---|
1497 | ENDDO |
---|
1498 | ENDIF |
---|
1499 | |
---|
1500 | CASE ( 'rlus' ) |
---|
1501 | IF ( radiation ) THEN |
---|
1502 | DO m = 1, vmea(l)%ns |
---|
1503 | ! |
---|
1504 | !-- Surface data is only available on inner subdomains, not |
---|
1505 | !-- on ghost points. Hence, limit the indices. |
---|
1506 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1507 | j = MERGE( j , nyn, j > nyn ) |
---|
1508 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1509 | i = MERGE( i , nxr, i > nxr ) |
---|
1510 | |
---|
1511 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1512 | surf_lsm_h%end_index(j,i) |
---|
1513 | vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm) |
---|
1514 | ENDDO |
---|
1515 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1516 | surf_usm_h%end_index(j,i) |
---|
1517 | vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm) |
---|
1518 | ENDDO |
---|
1519 | ENDDO |
---|
1520 | ENDIF |
---|
1521 | |
---|
1522 | CASE ( 'rlds' ) |
---|
1523 | IF ( radiation ) THEN |
---|
1524 | DO m = 1, vmea(l)%ns |
---|
1525 | ! |
---|
1526 | !-- Surface data is only available on inner subdomains, not |
---|
1527 | !-- on ghost points. Hence, limit the indices. |
---|
1528 | j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys ) |
---|
1529 | j = MERGE( j , nyn, j > nyn ) |
---|
1530 | i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl ) |
---|
1531 | i = MERGE( i , nxr, i > nxr ) |
---|
1532 | |
---|
1533 | DO mm = surf_lsm_h%start_index(j,i), & |
---|
1534 | surf_lsm_h%end_index(j,i) |
---|
1535 | vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm) |
---|
1536 | ENDDO |
---|
1537 | DO mm = surf_usm_h%start_index(j,i), & |
---|
1538 | surf_usm_h%end_index(j,i) |
---|
1539 | vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm) |
---|
1540 | ENDDO |
---|
1541 | ENDDO |
---|
1542 | ENDIF |
---|
1543 | |
---|
1544 | CASE ( 'rsd' ) |
---|
1545 | IF ( radiation ) THEN |
---|
1546 | DO m = 1, vmea(l)%ns |
---|
1547 | k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) |
---|
1548 | j = vmea(l)%j(m) |
---|
1549 | i = vmea(l)%i(m) |
---|
1550 | |
---|
1551 | vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i) |
---|
1552 | ENDDO |
---|
1553 | ENDIF |
---|
1554 | |
---|
1555 | CASE ( 'rsu' ) |
---|
1556 | IF ( radiation ) THEN |
---|
1557 | DO m = 1, vmea(l)%ns |
---|
1558 | k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) |
---|
1559 | j = vmea(l)%j(m) |
---|
1560 | i = vmea(l)%i(m) |
---|
1561 | |
---|
1562 | vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i) |
---|
1563 | ENDDO |
---|
1564 | ENDIF |
---|
1565 | |
---|
1566 | CASE ( 'rlu' ) |
---|
1567 | IF ( radiation ) THEN |
---|
1568 | DO m = 1, vmea(l)%ns |
---|
1569 | k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) |
---|
1570 | j = vmea(l)%j(m) |
---|
1571 | i = vmea(l)%i(m) |
---|
1572 | |
---|
1573 | vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i) |
---|
1574 | ENDDO |
---|
1575 | ENDIF |
---|
1576 | |
---|
1577 | CASE ( 'rld' ) |
---|
1578 | IF ( radiation ) THEN |
---|
1579 | DO m = 1, vmea(l)%ns |
---|
1580 | k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) |
---|
1581 | j = vmea(l)%j(m) |
---|
1582 | i = vmea(l)%i(m) |
---|
1583 | |
---|
1584 | vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i) |
---|
1585 | ENDDO |
---|
1586 | ENDIF |
---|
1587 | |
---|
1588 | CASE ( 'rsddif' ) |
---|
1589 | IF ( radiation ) THEN |
---|
1590 | DO m = 1, vmea(l)%ns |
---|
1591 | j = vmea(l)%j(m) |
---|
1592 | i = vmea(l)%i(m) |
---|
1593 | |
---|
1594 | vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i) |
---|
1595 | ENDDO |
---|
1596 | ENDIF |
---|
1597 | |
---|
1598 | CASE ( 't_soil' ) |
---|
1599 | DO m = 1, vmea(l)%ns_soil |
---|
1600 | i = vmea(l)%i_soil(m) |
---|
1601 | j = vmea(l)%j_soil(m) |
---|
1602 | k = vmea(l)%k_soil(m) |
---|
1603 | |
---|
1604 | match_lsm = surf_lsm_h%start_index(j,i) <= & |
---|
1605 | surf_lsm_h%end_index(j,i) |
---|
1606 | match_usm = surf_usm_h%start_index(j,i) <= & |
---|
1607 | surf_usm_h%end_index(j,i) |
---|
1608 | |
---|
1609 | IF ( match_lsm ) THEN |
---|
1610 | mm = surf_lsm_h%start_index(j,i) |
---|
1611 | vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm) |
---|
1612 | ENDIF |
---|
1613 | |
---|
1614 | IF ( match_usm ) THEN |
---|
1615 | mm = surf_usm_h%start_index(j,i) |
---|
1616 | vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm) |
---|
1617 | ENDIF |
---|
1618 | ENDDO |
---|
1619 | |
---|
1620 | CASE ( 'm_soil' ) |
---|
1621 | DO m = 1, vmea(l)%ns_soil |
---|
1622 | i = vmea(l)%i_soil(m) |
---|
1623 | j = vmea(l)%j_soil(m) |
---|
1624 | k = vmea(l)%k_soil(m) |
---|
1625 | |
---|
1626 | match_lsm = surf_lsm_h%start_index(j,i) <= & |
---|
1627 | surf_lsm_h%end_index(j,i) |
---|
1628 | |
---|
1629 | IF ( match_lsm ) THEN |
---|
1630 | mm = surf_lsm_h%start_index(j,i) |
---|
1631 | vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm) |
---|
1632 | ENDIF |
---|
1633 | |
---|
1634 | ENDDO |
---|
1635 | ! |
---|
1636 | !-- More will follow ... |
---|
1637 | |
---|
1638 | ! |
---|
1639 | !-- No match found - just set a fill value |
---|
1640 | CASE DEFAULT |
---|
1641 | vmea(l)%measured_vars(:,n) = vmea(l)%fillout |
---|
1642 | END SELECT |
---|
1643 | |
---|
1644 | ENDDO |
---|
1645 | |
---|
1646 | ENDDO |
---|
1647 | |
---|
1648 | END SUBROUTINE vm_sampling |
---|
1649 | |
---|
1650 | |
---|
1651 | END MODULE virtual_measurement_mod |
---|