Home

Context Navigation

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

Last change on this file since 3579 was 3522, checked in by suehring, 6 years ago
Sampling of variables extended in virtual measurement module
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 48.6 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 3522 2018-11-13 12:14:36Z suehring $
27	! Sampling of variables
28	!
29	! 3494 2018-11-06 14:51:27Z suehring
30	! Bugfixing
31	!
32	! 3473 2018-10-30 20:50:15Z suehring
33	! Initial revision
34	!
35	! Authors:
36	! --------
37	! @author Matthias Suehring
38	!
39	! Description:
40	! ------------
41	!> The module acts as an interface between 'real-world' observations and
42	!> model simulations. Virtual measurements will be taken in the model at the
43	!> coordinates representative for the 'real-world' measurement positions.
44	!> More precisely, coordinates and measured quanties will be read from a
45	!> NetCDF file which contains all required information. In the model,
46	!> the same quantities (as long as all the required components are switched-on)
47	!> will be sampled at the respective positions and output into an extra file,
48	!> which allows for straight-forward comparison of model results with
49	!> observations.
50	!>
51	!> @todo list_of_allowed variables needs careful checking
52	!> @todo output (binary or NetCDF) needs to be implemented
53	!> @todo clean-up anything from current test modus
54	!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
55	!> the (UC)2 standard
56	!> @note Fluxes are not processed
57	!------------------------------------------------------------------------------!
58	MODULE virtual_measurement_mod
59
60	USE arrays_3d, &
61	ONLY: q, pt, u, v, w, zu, zw
62
63	USE chem_modules, &
64	ONLY: nspec
65
66	USE chemistry_model_mod, &
67	ONLY: chem_species
68
69	USE control_parameters, &
70	ONLY: air_chemistry, dz, humidity, neutral, message_string, &
71	virtual_measurement
72
73	USE cpulog, &
74	ONLY: cpu_log, log_point
75
76	USE grid_variables, &
77	ONLY: dx, dy
78
79	USE indices, &
80	ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny
81
82	USE kinds
83
84
85	IMPLICIT NONE
86
87	TYPE virt_mea
88
89	CHARACTER(LEN=100) :: feature_type !< type of the measurement
90	CHARACTER(LEN=100) :: site !< name of the measurement site
91
92	CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables
93
94	INTEGER(iwp) :: ns = 0 !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
95	INTEGER(iwp) :: ntraj !< number of trajectories of a measurement
96	INTEGER(iwp) :: nvar !< number of measured variables
97
98	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues
99
100	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction
101	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction
102	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction
103
104	LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation
105	LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement
106	LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement
107
108	REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values
109	REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values
110	REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values
111	REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building
112	REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction
113	REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction
114
115	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables
116
117	END TYPE virt_mea
118
119	CHARACTER(LEN=5) :: char_eutm = "E_UTM" !< dimension name for UTM coordinate easting
120	CHARACTER(LEN=11) :: char_feature = "featureType" !< attribute name for feature type
121	CHARACTER(LEN=10) :: char_fillvalue = "_FillValue" !< variable attribute name for _FillValue
122	CHARACTER(LEN=18) :: char_mv = "measured_variables" !< variable name for the array with the measured variable names
123	CHARACTER(LEN=5) :: char_nutm = "N_UTM" !< dimension name for UTM coordinate northing
124	CHARACTER(LEN=18) :: char_numstations = "number_of_stations" !< attribute name for number of stations
125	CHARACTER(LEN=8) :: char_origx = "origin_x" !< attribute name for station coordinate in x
126	CHARACTER(LEN=8) :: char_origy = "origin_y" !< attribute name for station coordinate in y
127	CHARACTER(LEN=4) :: char_site = "site" !< attribute name for site name
128	CHARACTER(LEN=19) :: char_zag = "height_above_ground" !< attribute name for height above ground variable
129	CHARACTER(LEN=10) :: type_ts = 'timeSeries' !< name of stationary point measurements
130	CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements
131	CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements
132	!
133	!-- MS: List requires careful revision!
134	CHARACTER(LEN=10), DIMENSION(1:47), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM
135	(/ 'hfls ', & ! surface latent heat flux (W/m2)
136	'hfss ', & ! surface sensible heat flux (W/m2)
137	'hur ', & ! relative humidity (-)
138	'hus ', & ! specific humidity (g/kg)
139	'haa ', & ! absolute atmospheric humidity (kg/m3)
140	'mcpm1 ', & ! mass concentration of PM1 (kg/m3)
141	'mcpm2p5 ', & ! mass concentration of PM2.5 (kg/m3)
142	'mcpm10 ', & ! mass concentration of PM10 (kg/m3)
143	'mcpm10 ', & ! mass concentration of PM10 (kg/m3)
144	'mcco ', & ! mass concentration of CO (kg/m3)
145	'mcco2 ', & ! mass concentration of CO2 (kg/m3)
146	'mcbcda ', & ! mass concentration of black carbon paritcles (kg/m3)
147	'ncaa ', & ! number concentation of particles (1/m3)
148	'mfco ', & ! mole fraction of CO (mol/mol)
149	'mfco2 ', & ! mole fraction of CO2 (mol/mol)
150	'mfch4 ', & ! mole fraction of methane (mol/mol)
151	'mfnh3 ', & ! mole fraction of amonia (mol/mol)
152	'mfno ', & ! mole fraction of nitrogen monoxide (mol/mol)
153	'mfno2 ', & ! mole fraction of nitrogen dioxide (mol/mol)
154	'mfso2 ', & ! mole fraction of sulfur dioxide (mol/mol)
155	'mfh20 ', & ! mole fraction of water (mol/mol)
156	'plev ', & ! ? air pressure - hydrostaic + perturbation?
157	'rlds ', & ! surface downward longwave flux (W/m2)
158	'rlus ', & ! surface upward longwave flux (W/m2)
159	'rsds ', & ! surface downward shortwave flux (W/m2)
160	'rsus ', & ! surface upward shortwave flux (W/m2)
161	'ta ', & ! air temperature (degree C)
162	't_va ', & ! virtual accoustic temperature (K)
163	'theta ', & ! potential temperature (K)
164	'tro3 ', & ! mole fraction of ozone air (mol/mol)
165	'ts ', & ! scaling parameter of temperature (K)
166	'wspeed ', & ! ? wind speed - horizontal?
167	'wdir ', & ! wind direction
168	'us ', & ! friction velocity
169	'msoil ', & ! ? soil moisture - which depth?
170	'tsoil ', & ! ? soil temperature - which depth?
171	'u ', & ! u-component
172	'ua ', & ! eastward wind (is there any difference to u?)
173	'v ', & ! v-component
174	'va ', & ! northward wind (is there any difference to v?)
175	'w ', & ! w-component
176	'rld ', & ! downward longwave radiative flux (W/m2)
177	'rlu ', & ! upnward longwave radiative flux (W/m2)
178	'rsd ', & ! downward shortwave radiative flux (W/m2)
179	'rsu ', & ! upward shortwave radiative flux (W/m2)
180	'rsddif ', & ! downward shortwave diffuse radiative flux (W/m2)
181	'rnds ' & ! surface net downward radiative flux (W/m2)
182	/)
183
184	INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements
185	INTEGER(iwp) :: nvm = 0 !< number of virtual measurements
186	INTEGER(iwp) :: observation_coverage_xy = 0 !< horizontal distance from the measurement point where observations should be taken in the surrounding
187	INTEGER(iwp) :: observation_coverage_z = 0 !< vertical distance from the measurement point where observations should be taken in the surrounding
188
189	LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter
190	LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute
191
192	REAL(wp) :: vm_time_start = 0.0 !< time after virtual measurements should start
193
194
195	TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< virtual measurement data structure
196
197
198	INTERFACE vm_check_parameters
199	MODULE PROCEDURE vm_check_parameters
200	END INTERFACE vm_check_parameters
201
202	INTERFACE vm_init
203	MODULE PROCEDURE vm_init
204	END INTERFACE vm_init
205
206	INTERFACE vm_parin
207	MODULE PROCEDURE vm_parin
208	END INTERFACE vm_parin
209
210	INTERFACE vm_sampling
211	MODULE PROCEDURE vm_sampling
212	END INTERFACE vm_sampling
213
214	SAVE
215
216	PRIVATE
217
218	!
219	!-- Public interfaces
220	PUBLIC vm_check_parameters, vm_init, vm_parin, vm_sampling
221
222	!
223	!-- Public variables
224	PUBLIC vmea, vm_time_start
225
226	CONTAINS
227
228
229	!------------------------------------------------------------------------------!
230	! Description:
231	! ------------
232	!> Check parameters for virtual measurement module
233	!------------------------------------------------------------------------------!
234	SUBROUTINE vm_check_parameters
235
236	USE control_parameters, &
237	ONLY: message_string, virtual_measurement
238
239	USE netcdf_data_input_mod, &
240	ONLY: input_pids_static
241
242	IMPLICIT NONE
243
244	!
245	!-- In case virtual measurements are taken, a static input file is required.
246	!-- This is because UTM coordinates for the PALM domain origin are required
247	!-- for correct mapping of the measurements.
248	!-- ToDo: Revise this later and remove this requirement.
249	IF ( virtual_measurement .AND. .NOT. input_pids_static ) THEN
250	message_string = 'If virtual measurements are taken a static input ' // &
251	'file is mandatory.'
252	CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
253	ENDIF
254
255	END SUBROUTINE vm_check_parameters
256
257	!------------------------------------------------------------------------------!
258	! Description:
259	! ------------
260	!> Read namelist for the virtual measurement module
261	!------------------------------------------------------------------------------!
262	SUBROUTINE vm_parin
263
264	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
265
266	NAMELIST /virtual_measurement_parameters/ use_virtual_measurement, &
267	vm_time_start
268
269	line = ' '
270
271	!
272	!-- Try to find stg package
273	REWIND ( 11 )
274	line = ' '
275	DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
276	READ ( 11, '(A)', END=20 ) line
277	ENDDO
278	BACKSPACE ( 11 )
279
280	!
281	!-- Read namelist
282	READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
283
284	!
285	!-- Set flag that indicates that the virtual measurement module is switched on
286	IF ( use_virtual_measurement ) virtual_measurement = .TRUE.
287
288	GOTO 20
289
290	10 BACKSPACE( 11 )
291	READ( 11 , '(A)') line
292	CALL parin_fail_message( 'virtual_measurement_parameters', line )
293
294	20 CONTINUE
295
296	END SUBROUTINE vm_parin
297
298
299	!------------------------------------------------------------------------------!
300	! Description:
301	! ------------
302	!> Initialize virtual measurements: read coordiante arrays and measured
303	!> variables, set indicies indicating the measurement points, read further
304	!> attributes, etc..
305	!------------------------------------------------------------------------------!
306	SUBROUTINE vm_init
307
308	USE arrays_3d, &
309	ONLY: zu, zw
310
311	USE grid_variables, &
312	ONLY: ddx, ddy, dx, dy
313
314	USE indices, &
315	ONLY: nxl, nxr, nyn, nys
316
317	USE netcdf_data_input_mod, &
318	ONLY: init_model, input_file_vm, &
319	netcdf_data_input_get_dimension_length, &
320	netcdf_data_input_att, netcdf_data_input_var
321
322	USE surface_mod, &
323	ONLY: get_topography_top_index_ji
324
325	IMPLICIT NONE
326
327	CHARACTER(LEN=5) :: dum !< dummy string indicate station id
328	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF
329	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed
330
331	INTEGER(iwp) :: dim_eutm !< dimension size of UTM easting coordinate
332	INTEGER(iwp) :: dim_nutm !< dimension size of UTM northing coordinate
333	INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate
334	INTEGER(iwp) :: dim_zag !< dimension size of height coordinate
335	INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction
336	INTEGER(iwp) :: ii !< running index over all coordinate points of a measurement
337	INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction
338	INTEGER(iwp) :: j !< grid index of observation point in x-direction
339	INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction
340	INTEGER(iwp) :: k !< grid index of observation point in x-direction
341	INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate
342	INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction
343	INTEGER(iwp) :: ksurf !< topography top index
344	INTEGER(iwp) :: ku !< upper vertical index of surrounding grid points of an observation coordinate
345	INTEGER(iwp) :: l !< running index over all stations
346	INTEGER(iwp) :: len_char !< character length of single measured variables without Null character
347	INTEGER(iwp) :: ll !< running index over all measured variables in file
348	INTEGER(iwp) :: lll !< running index over all allowed variables
349	INTEGER(iwp) :: n !< running index over trajectory coordinates
350	INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain
351	INTEGER(iwp) :: t !< running index over number of trajectories
352
353	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: meas_flag !< mask array indicating measurement positions
354
355	LOGICAL :: chem_include !< flag indicating that chemical species is considered in modelled mechanism
356	LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain
357
358	REAL(wp) :: fill_eutm !< _FillValue for coordinate array E_UTM
359	REAL(wp) :: fill_nutm !< _FillValue for coordinate array N_UTM
360	REAL(wp) :: fill_zag !< _FillValue for height coordinate
361
362	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm !< easting UTM coordinate, temporary variable
363	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm !< northing UTM coordinate, temporary variable,
364	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z_ag !< height coordinate relative to origin_z, temporary variable
365	!
366	!-- Obtain number of virtual measurement stations
367	CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm, &
368	global_attribute, 'open', '' )
369
370	! write(9,*) "num stationi", nvm
371	! flush(9)
372	!
373	!-- ALLOCATE data structure
374	ALLOCATE( vmea(1:nvm) )
375	!
376	!-- Allocate flag array
377	ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
378	meas_flag = 0
379
380	!
381	!-- Read station coordinates and further attributes.
382	!-- Note all coordinates are in UTM coordinates.
383	DO l = 1, nvm
384	!
385	!-- Determine suffix which contains the ID
386	IF( l < 10 ) THEN
387	WRITE( dum, '(I1)') l
388	ELSEIF( l < 100 ) THEN
389	WRITE( dum, '(I2)') l
390	ELSEIF( l < 1000 ) THEN
391	WRITE( dum, '(I3)') l
392	ELSEIF( l < 10000 ) THEN
393	WRITE( dum, '(I4)') l
394	ELSEIF( l < 100000 ) THEN
395	WRITE( dum, '(I5)') l
396	ENDIF
397
398	CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx &
399	// TRIM( dum ), id_vm, '', global_attribute,&
400	'', '' )
401	CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy &
402	// TRIM( dum ), id_vm, '', global_attribute,&
403	'', '' )
404	CALL netcdf_data_input_att( vmea(l)%site, char_site &
405	// TRIM( dum ), id_vm, '', global_attribute,&
406	'', '' )
407	CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature &
408	// TRIM( dum ), id_vm, '', global_attribute,&
409	'', '' )
410
411	!
412	!--- Set logicals depending on the type of the measurement
413	IF ( INDEX( vmea(l)%feature_type, type_tspr ) /= 0 ) THEN
414	vmea(l)%timseries_profile = .TRUE.
415	ELSEIF ( INDEX( vmea(l)%feature_type, type_ts ) /= 0 ) THEN
416	vmea(l)%timseries = .TRUE.
417	ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN
418	vmea(l)%trajectory = .TRUE.
419	ELSE
420	!
421	!-- Give error message
422	message_string = 'Attribue featureType = ' // &
423	TRIM( vmea(l)%feature_type ) // &
424	' is not allowed.'
425	CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 )
426	ENDIF
427	!
428	!-- Read string with all measured variables at this station
429	measured_variables_file = ''
430	CALL netcdf_data_input_var( measured_variables_file, &
431	char_mv // TRIM( dum ), id_vm )
432	!
433	!-- Count the number of measured variables which match with the variables
434	!-- which are allowed to be measured in PALM. Please note, for some
435	!-- NetCDF interal reasons characters end with a NULL, i.e. also empty
436	!-- characters contain a NULL. Therefore, check the strings for a Null to
437	!-- get the correct character length in order to compare them with the list
438	!-- of allowed variables.
439	vmea(l)%nvar = 0
440	DO ll = 1, SIZE( measured_variables_file )
441	IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. &
442	measured_variables_file(ll)(1:1) /= ' ') THEN
443	!
444	!-- Obtain character length of the character
445	len_char = 1
446	DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
447	.AND. measured_variables_file(ll)(len_char:len_char) /= ' ' )
448	len_char = len_char + 1
449	ENDDO
450	len_char = len_char - 1
451	!
452	!-- Now, compare the measured variable with the list of allowed
453	!-- variables.
454	DO lll= 1, SIZE( list_allowed_variables )
455	IF ( measured_variables_file(ll)(1:len_char) == &
456	TRIM( list_allowed_variables(lll) ) ) THEN
457	vmea(l)%nvar = vmea(l)%nvar + 1
458	measured_variables(vmea(l)%nvar) = &
459	measured_variables_file(ll)(1:len_char)
460	ENDIF
461	ENDDO
462	ENDIF
463	ENDDO
464	!
465	!-- Allocate array for the measured variables names for the station l.
466	ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
467
468	DO ll = 1, vmea(l)%nvar
469	vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
470	ENDDO
471	!
472	!-- In case of chemistry, check if species is considered in the modelled
473	!-- chemistry mechanism.
474	IF ( air_chemistry ) THEN
475	DO ll = 1, vmea(l)%nvar
476	chem_include = .FALSE.
477	DO n = 1, nspec
478	IF ( TRIM( vmea(l)%measured_vars_name(ll) ) == &
479	TRIM( chem_species(n)%name ) ) chem_include = .TRUE.
480	ENDDO
481	IF ( .NOT. chem_include ) THEN
482	message_string = TRIM( vmea(l)%measured_vars_name(ll) ) // &
483	' is not considered in the modelled ' // &
484	'chemistry mechanism'
485	CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
486	ENDIF
487	ENDDO
488	ENDIF
489	!
490	!-- For the actual measurement ID read the UTM coordinates. Based on these,
491	!-- define the index space on each subdomain where measurements should be
492	!-- taken. Note, the entire coordinate arrays will not be stored on data
493	!-- type as this would exceed memory requirements, particularly for
494	!-- trajectory measurements. If no variable will be virtually measured,
495	!-- skip the reading.
496	IF ( vmea(l)%nvar > 0 ) THEN
497	!
498	!-- For stationary measurements UTM coordinates are just one value and
499	!-- its dimension is "station", while for mobile measurements UTM
500	!-- coordinates are arrays. First, inquire dimension length for
501	!-- UTM coordinates.
502	IF ( vmea(l)%trajectory ) THEN
503	!
504	!-- For non-stationary measurements read the number of trajectories
505	CALL netcdf_data_input_get_dimension_length( id_vm, &
506	vmea(l)%ntraj, &
507	"traj" // &
508	TRIM( dum ) )
509	CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime, &
510	"ntime" // &
511	TRIM( dum ) )
512	!
513	!-- For stationary measurements the dimension for UTM coordinates is 1
514	ELSE
515	vmea(l)%ntraj = 1
516	dim_ntime = 1
517	ENDIF
518
519	!
520	!- Allocate array which defines individual time frame for each
521	!-- trajectory or station
522	ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
523	!
524	!-- Allocate temporary arrays for UTM and height coordinates. Note,
525	!-- on file UTM coordinates might be 1D or 2D variables
526	ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
527	ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
528	ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime) )
529	!
530	!-- Read _FillValue attributes
531	CALL netcdf_data_input_att( fill_eutm, char_fillvalue, &
532	id_vm, '', .NOT. global_attribute, '', &
533	char_eutm // TRIM( dum ) )
534	CALL netcdf_data_input_att( fill_nutm, char_fillvalue, &
535	id_vm, '', .NOT. global_attribute, '', &
536	char_nutm // TRIM( dum ) )
537	CALL netcdf_data_input_att( fill_zag, char_fillvalue, &
538	id_vm, '', .NOT. global_attribute, '', &
539	char_zag // TRIM( dum ) )
540	!
541	!-- Read UTM and height coordinates coordinates for all trajectories and
542	!-- times.
543	IF ( vmea(l)%trajectory ) THEN
544	CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm, &
545	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
546	CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm, &
547	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
548	CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), id_vm, &
549	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
550	ELSE
551	CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm )
552	CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm )
553	CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), id_vm )
554	ENDIF
555	!
556	!-- For testing:
557	e_utm = e_utm - e_utm(1,1)
558	n_utm = n_utm - n_utm(1,1)
559
560	!
561	!-- First, compute relative x- and y-coordinates with respect to the
562	!-- lower-left origin of the model domain, which is the difference
563	!-- betwen UTM coordinates.
564	! e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t)) &
565	! - init_model%origin_x
566	! n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t)) &
567	! - init_model%origin_y
568
569	!
570	!-- Based on UTM coordinates, check if the measurement station or parts
571	!-- of the trajectory is on subdomain. This case, setup grid index space
572	!-- sample these quantities.
573	meas_flag = 0
574	DO t = 1, vmea(l)%ntraj
575	!
576	!-- Determine the individual time coordinate length for each station and
577	!-- trajectory. This is required as several stations and trajectories
578	!-- are merged into one file but they do not have the same number of
579	!-- points in time, hence, missing values may occur and cannot be
580	!-- processed further.
581	vmea(l)%dim_t(t) = 0
582	DO n = 1, dim_ntime
583	IF ( e_utm(t,n) /= fill_eutm .AND. &
584	n_utm(t,n) /= fill_nutm .AND. &
585	z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n
586	ENDDO
587
588	!
589	!-- Compute grid indices relative to origin and check if these are
590	!-- on the subdomain. Note, virtual measurements will be taken also
591	!-- at grid points surrounding the station, hence, check also for
592	!-- these grid points.
593	DO n = 1, vmea(l)%dim_t(t)
594	is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
595	js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
596	!
597	!-- Is the observation point on subdomain?
598	on_pe = ( is >= nxl .AND. is <= nxr .AND. &
599	js >= nys .AND. js <= nyn )
600	!
601	!-- Check if observation coordinate is on subdomain
602	IF ( on_pe ) THEN
603	!
604	!-- Determine vertical index which correspond to the observation
605	!-- height.
606	ksurf = get_topography_top_index_ji( js, is, 's' )
607	ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
608	!
609	!-- Set mask array at the observation coordinates. Also, flag the
610	!-- surrounding coordinate points, but first check whether the
611	!-- surrounding coordinate points are on the subdomain.
612	kl = MERGE( ks-1, ks, ks-1 >= nzb .AND. ks-1 > ksurf )
613	ku = MERGE( ks+1, ks, ks+1 <= nzt+1 )
614
615	meas_flag(kl:ku,js-1:js+1,is-1:is+1) = &
616	IBSET( meas_flag(kl:ku,js-1:js+1,is-1:is+1), 0 )
617	ENDIF
618	ENDDO
619
620	ENDDO
621	!
622	!-- Based on the flag array count the number of observation points.
623	ns = 0
624	DO is = nxl-1, nxr+1
625	DO js = nys-1, nyn+1
626	DO ks = nzb, nzt+1
627	ns = ns + MERGE( 1, 0, BTEST( meas_flag(ks,js,is), 0 ) )
628	ENDDO
629	ENDDO
630	ENDDO
631	!
632	!-- Store number of observation points on subdomain and allocate index
633	!-- arrays.
634	vmea(l)%ns = ns
635	ns = 0
636
637	ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
638	ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
639	ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
640	!
641	!-- Based on the flag array store the grid indices which correspond to
642	!-- the observation coordinates.
643	DO is = nxl-1, nxr+1
644	DO js = nys-1, nyn+1
645	DO ks = nzb, nzt+1
646	IF ( BTEST( meas_flag(ks,js,is), 0 ) ) THEN
647	ns = ns + 1
648	vmea(l)%i(ns) = is
649	vmea(l)%j(ns) = js
650	vmea(l)%k(ns) = ks
651	ENDIF
652	ENDDO
653	ENDDO
654	ENDDO
655
656	! write(9,*) l, "NS: ", ns
657	! IF ( ns > 0 ) THEN
658	! write(9,*) "i ", vmea(l)%i
659	! write(9,*) "j ", vmea(l)%j
660	! write(9,*) "k ", vmea(l)%k
661	! write(9,*)
662	! ENDIF
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	!
687	!-- Dellocate flag array
688	DEALLOCATE( meas_flag )
689
690	END SUBROUTINE vm_init
691
692
693	!------------------------------------------------------------------------------!
694	! Description:
695	! ------------
696	!> Sampling of the actual quantities along the observation coordinates
697	!------------------------------------------------------------------------------!
698	SUBROUTINE vm_sampling
699
700	USE arrays_3d, &
701	ONLY: exner, pt, q, u, v, w
702
703	USE basic_constants_and_equations_mod, &
704	ONLY: pi
705
706	USE radiation_model_mod, &
707	ONLY: radiation
708
709	USE surface_mod, &
710	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
711
712	IMPLICIT NONE
713
714	CHARACTER(LEN=10) :: trimvar !< dummy for the measured variable name
715
716	INTEGER(iwp) :: i !< grid index in x-direction
717	INTEGER(iwp) :: j !< grid index in y-direction
718	INTEGER(iwp) :: k !< grid index in z-direction
719	INTEGER(iwp) :: l !< running index over the number of stations
720	INTEGER(iwp) :: m !< running index over all virtual observation coordinates
721	INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate
722	INTEGER(iwp) :: n !< running index over all measured variables at a station
723	INTEGER(iwp) :: nn !< running index over the number of chemcal species
724	!
725	!-- Loop over all stations. For each possible variable loop over all
726	!-- observation points
727	DO l = 1, nvm
728	!
729	!-- Loop over all measured variables at this station. Please note,
730	!-- velocity components are interpolated onto scalar grid.
731	DO n = 1, vmea(l)%nvar
732
733	SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
734
735	CASE ( 'theta' )
736	IF ( .NOT. neutral ) THEN
737	DO m = 1, vmea(l)%ns
738	k = vmea(l)%k(m)
739	j = vmea(l)%j(m)
740	i = vmea(l)%i(m)
741	vmea(l)%measured_vars(n,m) = pt(k,j,i)
742	ENDDO
743	ENDIF
744
745	CASE ( 'ta', 't_va' )
746	IF ( .NOT. neutral ) THEN
747	DO m = 1, vmea(l)%ns
748	k = vmea(l)%k(m)
749	j = vmea(l)%j(m)
750	i = vmea(l)%i(m)
751	vmea(l)%measured_vars(n,m) = pt(k,j,i) * exner( k )
752	ENDDO
753	ENDIF
754
755	CASE ( 'hus', 'haa' )
756	IF ( humidity ) THEN
757	DO m = 1, vmea(l)%ns
758	k = vmea(l)%k(m)
759	j = vmea(l)%j(m)
760	i = vmea(l)%i(m)
761	vmea(l)%measured_vars(n,m) = q(k,j,i)
762	ENDDO
763	ENDIF
764
765	CASE ( 'u', 'ua' )
766	DO m = 1, vmea(l)%ns
767	k = vmea(l)%k(m)
768	j = vmea(l)%j(m)
769	i = vmea(l)%i(m)
770	vmea(l)%measured_vars(n,m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
771	ENDDO
772
773	CASE ( 'v', 'va' )
774	DO m = 1, vmea(l)%ns
775	k = vmea(l)%k(m)
776	j = vmea(l)%j(m)
777	i = vmea(l)%i(m)
778	vmea(l)%measured_vars(n,m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
779	ENDDO
780
781	CASE ( 'w' )
782	DO m = 1, vmea(l)%ns
783	k = vmea(l)%k(m)
784	j = vmea(l)%j(m)
785	i = vmea(l)%i(m)
786	vmea(l)%measured_vars(n,m) = w(k,j,i) !0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
787	ENDDO
788
789	CASE ( 'wspeed' )
790	DO m = 1, vmea(l)%ns
791	k = vmea(l)%k(m)
792	j = vmea(l)%j(m)
793	i = vmea(l)%i(m)
794	vmea(l)%measured_vars(n,m) = SQRT( &
795	( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
796	( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 &
797	)
798	ENDDO
799
800	CASE ( 'wdir' )
801	DO m = 1, vmea(l)%ns
802	k = vmea(l)%k(m)
803	j = vmea(l)%j(m)
804	i = vmea(l)%i(m)
805
806	vmea(l)%measured_vars(n,m) = ATAN2( &
807	- 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ), &
808	- 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) &
809	) * 180.0_wp / pi
810	ENDDO
811	!
812	!-- MS: list of variables needs extension.
813	CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mcco', 'mcco2', 'mcbcda', &
814	'ncaa', 'mfco2', 'mfco', 'mfch4', 'mfnh3', 'mfno', &
815	'mfno2', 'mfso2', 'mfh20', 'tr03' )
816	IF ( air_chemistry ) THEN
817	DO nn = 1, nspec
818	IF ( TRIM( vmea(l)%measured_vars_name(m) ) == &
819	TRIM( chem_species(nn)%name ) ) THEN
820	DO m = 1, vmea(l)%ns
821	k = vmea(l)%k(m)
822	j = vmea(l)%j(m)
823	i = vmea(l)%i(m)
824	vmea(l)%measured_vars(n,m) = &
825	chem_species(nn)%conc(k,j,i)
826	ENDDO
827	ENDIF
828	ENDDO
829	ENDIF
830
831	CASE ( 'us' )
832	DO m = 1, vmea(l)%ns
833	!
834	!-- Surface data is only available on inner subdomains, not
835	!-- on ghost points. Hence, limit the indices.
836	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
837	j = MERGE( j , nyn, j > nyn )
838	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
839	i = MERGE( i , nxr, i > nxr )
840
841	DO mm = surf_def_h(0)%start_index(j,i), &
842	surf_def_h(0)%end_index(j,i)
843	vmea(l)%measured_vars(n,m) = surf_def_h(0)%us(mm)
844	ENDDO
845	DO mm = surf_lsm_h%start_index(j,i), &
846	surf_lsm_h%end_index(j,i)
847	vmea(l)%measured_vars(n,m) = surf_lsm_h%us(mm)
848	ENDDO
849	DO mm = surf_usm_h%start_index(j,i), &
850	surf_usm_h%end_index(j,i)
851	vmea(l)%measured_vars(n,m) = surf_usm_h%us(mm)
852	ENDDO
853	ENDDO
854
855	CASE ( 'ts' )
856	DO m = 1, vmea(l)%ns
857	!
858	!-- Surface data is only available on inner subdomains, not
859	!-- on ghost points. Hence, limit the indices.
860	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
861	j = MERGE( j , nyn, j > nyn )
862	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
863	i = MERGE( i , nxr, i > nxr )
864
865	DO mm = surf_def_h(0)%start_index(j,i), &
866	surf_def_h(0)%end_index(j,i)
867	vmea(l)%measured_vars(n,m) = surf_def_h(0)%ts(mm)
868	ENDDO
869	DO mm = surf_lsm_h%start_index(j,i), &
870	surf_lsm_h%end_index(j,i)
871	vmea(l)%measured_vars(n,m) = surf_lsm_h%ts(mm)
872	ENDDO
873	DO mm = surf_usm_h%start_index(j,i), &
874	surf_usm_h%end_index(j,i)
875	vmea(l)%measured_vars(n,m) = surf_usm_h%ts(mm)
876	ENDDO
877	ENDDO
878
879	CASE ( 'hfls' )
880	DO m = 1, vmea(l)%ns
881	!
882	!-- Surface data is only available on inner subdomains, not
883	!-- on ghost points. Hence, limit the indices.
884	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
885	j = MERGE( j , nyn, j > nyn )
886	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
887	i = MERGE( i , nxr, i > nxr )
888
889	DO mm = surf_def_h(0)%start_index(j,i), &
890	surf_def_h(0)%end_index(j,i)
891	vmea(l)%measured_vars(n,m) = surf_def_h(0)%qsws(mm)
892	ENDDO
893	DO mm = surf_lsm_h%start_index(j,i), &
894	surf_lsm_h%end_index(j,i)
895	vmea(l)%measured_vars(n,m) = surf_lsm_h%qsws(mm)
896	ENDDO
897	DO mm = surf_usm_h%start_index(j,i), &
898	surf_usm_h%end_index(j,i)
899	vmea(l)%measured_vars(n,m) = surf_usm_h%qsws(mm)
900	ENDDO
901	ENDDO
902
903	CASE ( 'hfss' )
904	DO m = 1, vmea(l)%ns
905	!
906	!-- Surface data is only available on inner subdomains, not
907	!-- on ghost points. Hence, limit the indices.
908	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
909	j = MERGE( j , nyn, j > nyn )
910	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
911	i = MERGE( i , nxr, i > nxr )
912
913	DO mm = surf_def_h(0)%start_index(j,i), &
914	surf_def_h(0)%end_index(j,i)
915	vmea(l)%measured_vars(n,m) = surf_def_h(0)%shf(mm)
916	ENDDO
917	DO mm = surf_lsm_h%start_index(j,i), &
918	surf_lsm_h%end_index(j,i)
919	vmea(l)%measured_vars(n,m) = surf_lsm_h%shf(mm)
920	ENDDO
921	DO mm = surf_usm_h%start_index(j,i), &
922	surf_usm_h%end_index(j,i)
923	vmea(l)%measured_vars(n,m) = surf_usm_h%shf(mm)
924	ENDDO
925	ENDDO
926
927	CASE ( 'rnds' )
928	IF ( radiation ) THEN
929	DO m = 1, vmea(l)%ns
930	!
931	!-- Surface data is only available on inner subdomains, not
932	!-- on ghost points. Hence, limit the indices.
933	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
934	j = MERGE( j , nyn, j > nyn )
935	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
936	i = MERGE( i , nxr, i > nxr )
937
938	DO mm = surf_lsm_h%start_index(j,i), &
939	surf_lsm_h%end_index(j,i)
940	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_net(mm)
941	ENDDO
942	DO mm = surf_usm_h%start_index(j,i), &
943	surf_usm_h%end_index(j,i)
944	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_net(mm)
945	ENDDO
946	ENDDO
947	ENDIF
948
949	CASE ( 'rsus', 'rsu' )
950	IF ( radiation ) THEN
951	DO m = 1, vmea(l)%ns
952	!
953	!-- Surface data is only available on inner subdomains, not
954	!-- on ghost points. Hence, limit the indices.
955	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
956	j = MERGE( j , nyn, j > nyn )
957	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
958	i = MERGE( i , nxr, i > nxr )
959
960	DO mm = surf_lsm_h%start_index(j,i), &
961	surf_lsm_h%end_index(j,i)
962	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_out(mm)
963	ENDDO
964	DO mm = surf_usm_h%start_index(j,i), &
965	surf_usm_h%end_index(j,i)
966	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_out(mm)
967	ENDDO
968	ENDDO
969	ENDIF
970
971	CASE ( 'rsds', 'rsd' )
972	IF ( radiation ) THEN
973	DO m = 1, vmea(l)%ns
974	!
975	!-- Surface data is only available on inner subdomains, not
976	!-- on ghost points. Hence, limit the indices.
977	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
978	j = MERGE( j , nyn, j > nyn )
979	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
980	i = MERGE( i , nxr, i > nxr )
981
982	DO mm = surf_lsm_h%start_index(j,i), &
983	surf_lsm_h%end_index(j,i)
984	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_in(mm)
985	ENDDO
986	DO mm = surf_usm_h%start_index(j,i), &
987	surf_usm_h%end_index(j,i)
988	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_in(mm)
989	ENDDO
990	ENDDO
991	ENDIF
992
993	CASE ( 'rlus', 'rlu' )
994	IF ( radiation ) THEN
995	DO m = 1, vmea(l)%ns
996	!
997	!-- Surface data is only available on inner subdomains, not
998	!-- on ghost points. Hence, limit the indices.
999	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1000	j = MERGE( j , nyn, j > nyn )
1001	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1002	i = MERGE( i , nxr, i > nxr )
1003
1004	DO mm = surf_lsm_h%start_index(j,i), &
1005	surf_lsm_h%end_index(j,i)
1006	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_out(mm)
1007	ENDDO
1008	DO mm = surf_usm_h%start_index(j,i), &
1009	surf_usm_h%end_index(j,i)
1010	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_out(mm)
1011	ENDDO
1012	ENDDO
1013	ENDIF
1014
1015	CASE ( 'rlds', 'rld' )
1016	IF ( radiation ) THEN
1017	DO m = 1, vmea(l)%ns
1018	!
1019	!-- Surface data is only available on inner subdomains, not
1020	!-- on ghost points. Hence, limit the indices.
1021	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1022	j = MERGE( j , nyn, j > nyn )
1023	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1024	i = MERGE( i , nxr, i > nxr )
1025
1026	DO mm = surf_lsm_h%start_index(j,i), &
1027	surf_lsm_h%end_index(j,i)
1028	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_in(mm)
1029	ENDDO
1030	DO mm = surf_usm_h%start_index(j,i), &
1031	surf_usm_h%end_index(j,i)
1032	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_in(mm)
1033	ENDDO
1034	ENDDO
1035	ENDIF
1036	!
1037	!-- More will follow ...
1038
1039	END SELECT
1040
1041	ENDDO
1042
1043	ENDDO
1044
1045	END SUBROUTINE vm_sampling
1046
1047
1048	END MODULE virtual_measurement_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |