Home

Context Navigation

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

Last change on this file since 3500 was 3494, checked in by suehring, 6 years ago
Surface output revised and some bugs are fixed + new post-processing tool to convert binary surface output to Paraview readable VTK files
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 35.3 KB

Line
1	!> @virtual_measurement_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: virtual_measurement_mod.f90 3494 2018-11-06 14:51:27Z gronemeier $
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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |