Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |