Home

Context Navigation

source: palm/trunk/SOURCE/virtual_flight_mod.f90 @ 2334

Last change on this file since 2334 was 2271, checked in by sward, 7 years ago
error messages and numbers updated
Property svn:keywords set to `Id`
File size: 37.5 KB

Line
1	!> @file virtual_flights_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
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 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: virtual_flight_mod.f90 2271 2017-06-09 12:34:55Z gronemeier $
27	! Todo added
28	!
29	! 2101 2017-01-05 16:42:31Z suehring
30	!
31	! 2000 2016-08-20 18:09:15Z knoop
32	! Forced header and separation lines into 80 columns
33	!
34	! 1960 2016-07-12 16:34:24Z suehring
35	! Separate humidity and passive scalar.
36	! Bugfix concerning labeling of timeseries.
37	!
38	! 1957 2016-07-07 10:43:48Z suehring
39	! Initial revision
40	!
41	! Description:
42	! ------------
43	!> Module for virtual flight measurements.
44	!> @todo Err msg PA0438: flight can be inside topography -> extra check?
45	!--------------------------------------------------------------------------------!
46	MODULE flight_mod
47
48	USE control_parameters, &
49	ONLY: fl_max, num_leg, num_var_fl, num_var_fl_user, virtual_flight
50
51	USE kinds
52
53	CHARACTER(LEN=6), DIMENSION(fl_max) :: leg_mode = 'cyclic' !< flight mode through the model domain, either 'cyclic' or 'return'
54
55	INTEGER(iwp) :: l !< index for flight leg
56	INTEGER(iwp) :: var_index !< index for measured variable
57
58	LOGICAL, DIMENSION(:), ALLOCATABLE :: cyclic_leg !< flag to identify fly mode
59
60	REAL(wp) :: flight_end = 9999999.9_wp !< end time of virtual flight
61	REAL(wp) :: flight_begin = 0.0_wp !< end time of virtual flight
62
63	REAL(wp), DIMENSION(fl_max) :: flight_angle = 45.0_wp !< angle determining the horizontal flight direction
64	REAL(wp), DIMENSION(fl_max) :: flight_level = 100.0_wp !< flight level
65	REAL(wp), DIMENSION(fl_max) :: max_elev_change = 0.0_wp !< maximum elevation change for the respective flight leg
66	REAL(wp), DIMENSION(fl_max) :: rate_of_climb = 0.0_wp !< rate of climb or descent
67	REAL(wp), DIMENSION(fl_max) :: speed_agl = 25.0_wp !< absolute horizontal flight speed above ground level (agl)
68	REAL(wp), DIMENSION(fl_max) :: x_start = 999999999.0_wp !< start x position
69	REAL(wp), DIMENSION(fl_max) :: x_end = 999999999.0_wp !< end x position
70	REAL(wp), DIMENSION(fl_max) :: y_start = 999999999.0_wp !< start y position
71	REAL(wp), DIMENSION(fl_max) :: y_end = 999999999.0_wp !< end y position
72
73	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_agl !< u-component of flight speed
74	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_agl !< v-component of flight speed
75	REAL(wp), DIMENSION(:), ALLOCATABLE :: w_agl !< w-component of flight speed
76	REAL(wp), DIMENSION(:), ALLOCATABLE :: x_pos !< aircraft x-position
77	REAL(wp), DIMENSION(:), ALLOCATABLE :: y_pos !< aircraft y-position
78	REAL(wp), DIMENSION(:), ALLOCATABLE :: z_pos !< aircraft z-position
79
80	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sensor_l !< measured data on local PE
81	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sensor !< measured data
82
83	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: var_u !< dummy array for possibly user-defined quantities
84
85	SAVE
86
87	PRIVATE
88
89	INTERFACE flight_header
90	MODULE PROCEDURE flight_header
91	END INTERFACE flight_header
92
93	INTERFACE flight_init
94	MODULE PROCEDURE flight_init
95	END INTERFACE flight_init
96
97	INTERFACE flight_init_output
98	MODULE PROCEDURE flight_init_output
99	END INTERFACE flight_init_output
100
101	INTERFACE flight_check_parameters
102	MODULE PROCEDURE flight_check_parameters
103	END INTERFACE flight_check_parameters
104
105	INTERFACE flight_parin
106	MODULE PROCEDURE flight_parin
107	END INTERFACE flight_parin
108
109	INTERFACE interpolate_xyz
110	MODULE PROCEDURE interpolate_xyz
111	END INTERFACE interpolate_xyz
112
113	INTERFACE flight_measurement
114	MODULE PROCEDURE flight_measurement
115	END INTERFACE flight_measurement
116
117	INTERFACE flight_skip_var_list
118	MODULE PROCEDURE flight_skip_var_list
119	END INTERFACE flight_skip_var_list
120
121	INTERFACE flight_read_restart_data
122	MODULE PROCEDURE flight_read_restart_data
123	END INTERFACE flight_read_restart_data
124
125	INTERFACE flight_write_restart_data
126	MODULE PROCEDURE flight_write_restart_data
127	END INTERFACE flight_write_restart_data
128
129	!
130	!-- Private interfaces
131	PRIVATE flight_check_parameters, flight_init_output, interpolate_xyz
132	!
133	!-- Public interfaces
134	PUBLIC flight_init, flight_header, flight_parin, flight_measurement, &
135	flight_skip_var_list, flight_read_restart_data, &
136	flight_write_restart_data
137	!
138	!-- Public variables
139	PUBLIC fl_max, sensor, x_pos, y_pos, z_pos
140
141	CONTAINS
142
143	!------------------------------------------------------------------------------!
144	! Description:
145	! ------------
146	!> Header output for flight module.
147	!------------------------------------------------------------------------------!
148	SUBROUTINE flight_header ( io )
149
150
151	IMPLICIT NONE
152
153	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
154
155	WRITE ( io, 1 )
156	WRITE ( io, 2 )
157	WRITE ( io, 3 ) num_leg
158	WRITE ( io, 4 ) flight_begin
159	WRITE ( io, 5 ) flight_end
160
161	DO l=1, num_leg
162	WRITE ( io, 6 ) l
163	WRITE ( io, 7 ) speed_agl(l)
164	WRITE ( io, 8 ) flight_level(l)
165	WRITE ( io, 9 ) max_elev_change(l)
166	WRITE ( io, 10 ) rate_of_climb(l)
167	WRITE ( io, 11 ) leg_mode(l)
168	ENDDO
169
170
171	1 FORMAT (' Virtual flights:'/ &
172	' ----------------')
173	2 FORMAT (' Output every timestep')
174	3 FORMAT (' Number of flight legs:', I3 )
175	4 FORMAT (' Begin of measurements:', F10.1 , ' s' )
176	5 FORMAT (' End of measurements:', F10.1 , ' s' )
177	6 FORMAT (' Leg', I3/, &
178	' ------' )
179	7 FORMAT (' Flight speed : ', F5.1, ' m/s' )
180	8 FORMAT (' Flight level : ', F5.1, ' m' )
181	9 FORMAT (' Maximum elevation change: ', F5.1, ' m/s' )
182	10 FORMAT (' Rate of climb / descent : ', F5.1, ' m/s' )
183	11 FORMAT (' Leg mode : ', A/ )
184
185	END SUBROUTINE flight_header
186
187	!------------------------------------------------------------------------------!
188	! Description:
189	! ------------
190	!> Reads the namelist flight_par.
191	!------------------------------------------------------------------------------!
192	SUBROUTINE flight_parin
193
194
195	IMPLICIT NONE
196
197	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
198
199	NAMELIST /flight_par/ flight_angle, flight_end, flight_begin, leg_mode,&
200	flight_level, max_elev_change, rate_of_climb, &
201	speed_agl, x_end, x_start, y_end, y_start
202
203	!
204	!-- Try to find the namelist flight_par
205	REWIND ( 11 )
206	line = ' '
207	DO WHILE ( INDEX( line, '&flight_par' ) == 0 )
208	READ ( 11, '(A)', END=10 ) line
209	ENDDO
210	BACKSPACE ( 11 )
211
212	!
213	!-- Read namelist
214	READ ( 11, flight_par )
215	!
216	!-- Set switch that virtual flights shall be carried out
217	virtual_flight = .TRUE.
218
219	10 CONTINUE
220
221	END SUBROUTINE flight_parin
222
223	!------------------------------------------------------------------------------!
224	! Description:
225	! ------------
226	!> Inititalization of required arrays, number of legs and flags. Moreover,
227	!> initialize flight speed and -direction, as well as start positions.
228	!------------------------------------------------------------------------------!
229	SUBROUTINE flight_init
230
231	USE constants, &
232	ONLY: pi
233
234	USE control_parameters, &
235	ONLY: initializing_actions
236
237	USE indices, &
238	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
239
240	IMPLICIT NONE
241
242	REAL(wp) :: distance !< distance between start and end position of a flight leg
243	!
244	!-- Determine the number of flight legs
245	l = 1
246	DO WHILE ( x_start(l) /= 999999999.0_wp .AND. l <= SIZE(x_start) )
247	l = l + 1
248	ENDDO
249	num_leg = l-1
250	!
251	!-- Check for proper parameter settings
252	CALL flight_check_parameters
253	!
254	!-- Allocate and initialize logical array for flight pattern
255	ALLOCATE( cyclic_leg(1:num_leg) )
256	!
257	!-- Initialize flags for cyclic/return legs
258	DO l = 1, num_leg
259	cyclic_leg(l) = MERGE( .TRUE., .FALSE., &
260	TRIM( leg_mode(l) ) == 'cyclic' &
261	)
262	ENDDO
263	!
264	!-- Allocate and initialize arraxs for flight position and speed. In case of
265	!-- restart runs these data are read by the routine read_flight_restart_data
266	!-- instead.
267	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
268
269	ALLOCATE( x_pos(1:num_leg), y_pos(1:num_leg ), z_pos(1:num_leg) )
270	!
271	!-- Inititalize x-, y-, and z-positions with initial start position
272	x_pos(1:num_leg) = x_start(1:num_leg)
273	y_pos(1:num_leg) = y_start(1:num_leg)
274	z_pos(1:num_leg) = flight_level(1:num_leg)
275	!
276	!-- Allocate arrays for flight-speed components
277	ALLOCATE( u_agl(1:num_leg), &
278	v_agl(1:num_leg), &
279	w_agl(1:num_leg) )
280	!
281	!-- Inititalize u-, v- and w-component.
282	DO l = 1, num_leg
283	!
284	!-- In case of return-legs, the flight direction, i.e. the horizontal
285	!-- flight-speed components, are derived from the given start/end
286	!-- positions.
287	IF ( .NOT. cyclic_leg(l) ) THEN
288	distance = SQRT( ( x_end(l) - x_start(l) )**2 &
289	+ ( y_end(l) - y_start(l) )**2 )
290	u_agl(l) = speed_agl(l) * ( x_end(l) - x_start(l) ) / distance
291	v_agl(l) = speed_agl(l) * ( y_end(l) - y_start(l) ) / distance
292	w_agl(l) = rate_of_climb(l)
293	!
294	!-- In case of cyclic-legs, flight direction is directly derived from
295	!-- the given flight angle.
296	ELSE
297	u_agl(l) = speed_agl(l) * COS( flight_angle(l) * pi / 180.0_wp )
298	v_agl(l) = speed_agl(l) * SIN( flight_angle(l) * pi / 180.0_wp )
299	w_agl(l) = rate_of_climb(l)
300	ENDIF
301
302	ENDDO
303
304	ENDIF
305	!
306	!-- Initialized data output
307	CALL flight_init_output
308	!
309	!-- Allocate array required for user-defined quantities if necessary.
310	IF ( num_var_fl_user > 0 ) &
311	ALLOCATE( var_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
312	!
313	!-- Allocate and initialize arrays containing the measured data
314	ALLOCATE( sensor_l(1:num_var_fl,1:num_leg) )
315	ALLOCATE( sensor(1:num_var_fl,1:num_leg) )
316	sensor_l = 0.0_wp
317	sensor = 0.0_wp
318
319	END SUBROUTINE flight_init
320
321
322	!------------------------------------------------------------------------------!
323	! Description:
324	! ------------
325	!> Initialization of output-variable names and units.
326	!------------------------------------------------------------------------------!
327	SUBROUTINE flight_init_output
328
329	USE control_parameters, &
330	ONLY: cloud_droplets, cloud_physics, humidity, neutral, &
331	passive_scalar
332
333	USE netcdf_interface
334
335	IMPLICIT NONE
336
337	CHARACTER(LEN=10) :: label_leg !< dummy argument to convert integer to string
338
339	INTEGER(iwp) :: i !< loop variable
340	INTEGER(iwp) :: id_pt !< identifyer for labeling
341	INTEGER(iwp) :: id_q !< identifyer for labeling
342	INTEGER(iwp) :: id_ql !< identifyer for labeling
343	INTEGER(iwp) :: id_s !< identifyer for labeling
344	INTEGER(iwp) :: id_u = 1 !< identifyer for labeling
345	INTEGER(iwp) :: id_v = 2 !< identifyer for labeling
346	INTEGER(iwp) :: id_w = 3 !< identifyer for labeling
347	INTEGER(iwp) :: k !< index variable
348
349	LOGICAL :: init = .TRUE. !< flag to distiquish calls of user_init_flight
350	!
351	!-- Define output quanities, at least three variables are measured (u,v,w)
352	num_var_fl = 3
353	IF ( .NOT. neutral ) THEN
354	num_var_fl = num_var_fl + 1
355	id_pt = num_var_fl
356	ENDIF
357	IF ( humidity ) THEN
358	num_var_fl = num_var_fl + 1
359	id_q = num_var_fl
360	ENDIF
361	IF ( cloud_physics .OR. cloud_droplets ) THEN
362	num_var_fl = num_var_fl + 1
363	id_ql = num_var_fl
364	ENDIF
365	IF ( passive_scalar ) THEN
366	num_var_fl = num_var_fl + 1
367	id_s = num_var_fl
368	ENDIF
369	!
370	!-- Write output strings for dimensions x, y, z
371	DO l=1, num_leg
372
373	IF ( l < 10 ) WRITE( label_leg, '(I1)') l
374	IF ( l >= 10 .AND. l < 100 ) WRITE( label_leg, '(I2)') l
375	IF ( l >= 100 .AND. l < 1000 ) WRITE( label_leg, '(I3)') l
376
377	dofl_dim_label_x(l) = 'x_' // TRIM( label_leg )
378	dofl_dim_label_y(l) = 'y_' // TRIM( label_leg )
379	dofl_dim_label_z(l) = 'z_' // TRIM( label_leg )
380
381	ENDDO
382
383	!
384	!-- Call user routine to initialize further variables
385	CALL user_init_flight( init )
386	!
387	!-- Write output labels and units for the quanities
388	k = 1
389	DO l=1, num_leg
390
391	IF ( l < 10 ) WRITE( label_leg, '(I1)') l
392	IF ( l >= 10 .AND. l < 100 ) WRITE( label_leg, '(I2)') l
393	IF ( l >= 100 .AND. l < 1000 ) WRITE( label_leg, '(I3)') l
394
395	label_leg = 'leg_' // TRIM(label_leg)
396	DO i=1, num_var_fl
397
398	IF ( i == id_u ) THEN
399	dofl_label(k) = TRIM( label_leg ) // '_u'
400	dofl_unit(k) = 'm/s'
401	k = k + 1
402	ELSEIF ( i == id_v ) THEN
403	dofl_label(k) = TRIM( label_leg ) // '_v'
404	dofl_unit(k) = 'm/s'
405	k = k + 1
406	ELSEIF ( i == id_w ) THEN
407	dofl_label(k) = TRIM( label_leg ) // '_w'
408	dofl_unit(k) = 'm/s'
409	k = k + 1
410	ELSEIF ( i == id_pt ) THEN
411	dofl_label(k) = TRIM( label_leg ) // '_pt'
412	dofl_unit(k) = 'K'
413	k = k + 1
414	ELSEIF ( i == id_q ) THEN
415	dofl_label(k) = TRIM( label_leg ) // '_q'
416	dofl_unit(k) = 'kg/kg'
417	k = k + 1
418	ELSEIF ( i == id_ql ) THEN
419	dofl_label(k) = TRIM( label_leg ) // '_ql'
420	dofl_unit(k) = 'kg/kg'
421	k = k + 1
422	ELSEIF ( i == id_s ) THEN
423	dofl_label(k) = TRIM( label_leg ) // '_s'
424	dofl_unit(k) = 'kg/kg'
425	k = k + 1
426	ENDIF
427	ENDDO
428
429	DO i=1, num_var_fl_user
430	CALL user_init_flight( init, k, i, label_leg )
431	ENDDO
432
433	ENDDO
434	!
435	!-- Finally, set the total number of flight-output quantities.
436	num_var_fl = num_var_fl + num_var_fl_user
437
438	END SUBROUTINE flight_init_output
439
440	!------------------------------------------------------------------------------!
441	! Description:
442	! ------------
443	!> This routine calculates the current flight positions and calls the
444	!> respective interpolation routine to measures the data at the current
445	!> flight position.
446	!------------------------------------------------------------------------------!
447	SUBROUTINE flight_measurement
448
449	USE arrays_3d, &
450	ONLY: ddzu, ddzw, pt, q, ql, s, u, v, w, zu, zw
451
452	USE control_parameters, &
453	ONLY: cloud_droplets, cloud_physics, dz, dz_stretch_level, dt_3d, &
454	humidity, neutral, passive_scalar, simulated_time
455
456	USE cpulog, &
457	ONLY: cpu_log, log_point
458
459	USE grid_variables, &
460	ONLY: ddx, ddy, dx, dy
461
462	USE indices, &
463	ONLY: nx, nxl, nxlg, nxr, nxrg, ny, nys, nysg, nyn, nyng, nzb, nzt
464
465	USE pegrid
466
467	IMPLICIT NONE
468
469	LOGICAL :: on_pe !< flag to check if current flight position is on current PE
470
471	REAL(wp) :: x !< distance between left edge of current grid box and flight position
472	REAL(wp) :: y !< distance between south edge of current grid box and flight position
473
474	INTEGER(iwp) :: i !< index of current grid box along x
475	INTEGER(iwp) :: j !< index of current grid box along y
476	INTEGER(iwp) :: n !< loop variable for number of user-defined output quantities
477
478	CALL cpu_log( log_point(65), 'flight_measurement', 'start' )
479	!
480	!-- Perform flight measurement if start time is reached.
481	IF ( simulated_time >= flight_begin .AND. &
482	simulated_time <= flight_end ) THEN
483
484	sensor_l = 0.0_wp
485	sensor = 0.0_wp
486	!
487	!-- Loop over all flight legs
488	DO l=1, num_leg
489	!
490	!-- Update location for each leg
491	x_pos(l) = x_pos(l) + u_agl(l) * dt_3d
492	y_pos(l) = y_pos(l) + v_agl(l) * dt_3d
493	z_pos(l) = z_pos(l) + w_agl(l) * dt_3d
494	!
495	!-- Check if location must be modified for return legs.
496	!-- Carry out horizontal reflection if required.
497	IF ( .NOT. cyclic_leg(l) ) THEN
498	!
499	!-- Outward flight, i.e. from start to end
500	IF ( u_agl(l) >= 0.0_wp .AND. x_pos(l) > x_end(l) ) THEN
501	x_pos(l) = 2.0_wp * x_end(l) - x_pos(l)
502	u_agl(l) = - u_agl(l)
503	!
504	!-- Return flight, i.e. from end to start
505	ELSEIF ( u_agl(l) < 0.0_wp .AND. x_pos(l) < x_start(l) ) THEN
506	x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
507	u_agl(l) = - u_agl(l)
508	ENDIF
509	!
510	!-- Outward flight, i.e. from start to end
511	IF ( v_agl(l) >= 0.0_wp .AND. y_pos(l) > y_end(l) ) THEN
512	y_pos(l) = 2.0_wp * y_end(l) - y_pos(l)
513	v_agl(l) = - v_agl(l)
514	!
515	!-- Return flight, i.e. from end to start
516	ELSEIF ( v_agl(l) < 0.0_wp .AND. y_pos(l) < y_start(l) ) THEN
517	y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
518	v_agl(l) = - v_agl(l)
519	ENDIF
520	!
521	!-- Check if flight position is out of the model domain and apply
522	!-- cyclic conditions if required
523	ELSEIF ( cyclic_leg(l) ) THEN
524	!
525	!-- Check if aircraft leaves the model domain at the right boundary
526	IF ( ( flight_angle(l) >= 0.0_wp .AND. &
527	flight_angle(l) <= 90.0_wp ) .OR. &
528	( flight_angle(l) >= 270.0_wp .AND. &
529	flight_angle(l) <= 360.0_wp ) ) THEN
530	IF ( x_pos(l) >= ( nx + 0.5_wp ) * dx ) &
531	x_pos(l) = x_pos(l) - ( nx + 1 ) * dx
532	!
533	!-- Check if aircraft leaves the model domain at the left boundary
534	ELSEIF ( flight_angle(l) > 90.0_wp .AND. &
535	flight_angle(l) < 270.0_wp ) THEN
536	IF ( x_pos(l) < -0.5_wp * dx ) &
537	x_pos(l) = ( nx + 1 ) * dx + x_pos(l)
538	ENDIF
539	!
540	!-- Check if aircraft leaves the model domain at the north boundary
541	IF ( flight_angle(l) >= 0.0_wp .AND. &
542	flight_angle(l) <= 180.0_wp ) THEN
543	IF ( y_pos(l) >= ( ny + 0.5_wp ) * dy ) &
544	y_pos(l) = y_pos(l) - ( ny + 1 ) * dy
545	!
546	!-- Check if aircraft leaves the model domain at the south boundary
547	ELSEIF ( flight_angle(l) > 180.0_wp .AND. &
548	flight_angle(l) < 360.0_wp ) THEN
549	IF ( y_pos(l) < -0.5_wp * dy ) &
550	y_pos(l) = ( ny + 1 ) * dy + y_pos(l)
551	ENDIF
552
553	ENDIF
554	!
555	!-- Check if maximum elevation change is already reached. If required
556	!-- reflect vertically.
557	IF ( rate_of_climb(l) /= 0.0_wp ) THEN
558	!
559	!-- First check if aircraft is too high
560	IF ( w_agl(l) > 0.0_wp .AND. &
561	z_pos(l) - flight_level(l) > max_elev_change(l) ) THEN
562	z_pos(l) = 2.0_wp * ( flight_level(l) + max_elev_change(l) )&
563	- z_pos(l)
564	w_agl(l) = - w_agl(l)
565	!
566	!-- Check if aircraft is too low
567	ELSEIF ( w_agl(l) < 0.0_wp .AND. z_pos(l) < flight_level(l) ) THEN
568	z_pos(l) = 2.0_wp * flight_level(l) - z_pos(l)
569	w_agl(l) = - w_agl(l)
570	ENDIF
571
572	ENDIF
573	!
574	!-- Determine grid indices for flight position along x- and y-direction.
575	!-- Please note, there is a special treatment for the index
576	!-- along z-direction, which is due to vertical grid stretching.
577	i = ( x_pos(l) + 0.5_wp * dx ) * ddx
578	j = ( y_pos(l) + 0.5_wp * dy ) * ddy
579	!
580	!-- Check if indices are on current PE
581	on_pe = ( i >= nxl .AND. i <= nxr .AND. &
582	j >= nys .AND. j <= nyn )
583
584	IF ( on_pe ) THEN
585
586	var_index = 1
587	!
588	!-- Recalculate indices, as u is shifted by -0.5*dx.
589	i = x_pos(l) * ddx
590	j = ( y_pos(l) + 0.5_wp * dy ) * ddy
591	!
592	!-- Calculate distance from left and south grid-box edges.
593	x = x_pos(l) - ( 0.5_wp - i ) * dx
594	y = y_pos(l) - j * dy
595	!
596	!-- Interpolate u-component onto current flight position.
597	CALL interpolate_xyz( u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
598	var_index = var_index + 1
599	!
600	!-- Recalculate indices, as v is shifted by -0.5*dy.
601	i = ( x_pos(l) + 0.5_wp * dx ) * ddx
602	j = y_pos(l) * ddy
603
604	x = x_pos(l) - i * dx
605	y = y_pos(l) - ( 0.5_wp - j ) * dy
606	CALL interpolate_xyz( v, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
607	var_index = var_index + 1
608	!
609	!-- Interpolate w and scalar quantities. Recalculate indices.
610	i = ( x_pos(l) + 0.5_wp * dx ) * ddx
611	j = ( y_pos(l) + 0.5_wp * dy ) * ddy
612	x = x_pos(l) - i * dx
613	y = y_pos(l) - j * dy
614	!
615	!-- Interpolate w-velocity component.
616	CALL interpolate_xyz( w, zw, ddzw, 0.0_wp, x, y, var_index, j, i )
617	var_index = var_index + 1
618	!
619	!-- Potential temerature
620	IF ( .NOT. neutral ) THEN
621	CALL interpolate_xyz( pt, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
622	var_index = var_index + 1
623	ENDIF
624	!
625	!-- Humidity
626	IF ( humidity ) THEN
627	CALL interpolate_xyz( q, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
628	var_index = var_index + 1
629	ENDIF
630	!
631	!-- Liquid water content
632	IF ( cloud_physics .OR. cloud_droplets ) THEN
633	CALL interpolate_xyz( ql, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
634	var_index = var_index + 1
635	ENDIF
636	!
637	!-- Passive scalar
638	IF ( passive_scalar ) THEN
639	CALL interpolate_xyz( s, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
640	var_index = var_index + 1
641	ENDIF
642	!
643	!-- Treat user-defined variables if required
644	DO n = 1, num_var_fl_user
645	CALL user_flight( var_u, n )
646	CALL interpolate_xyz( var_u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
647	var_index = var_index + 1
648	ENDDO
649	ENDIF
650
651	ENDDO
652	!
653	!-- Write local data on global array.
654	#if defined( __parallel )
655	CALL MPI_ALLREDUCE(sensor_l(1,1), sensor(1,1), &
656	num_var_fl*num_leg, MPI_REAL, MPI_SUM, &
657	comm2d, ierr )
658	#else
659	sensor = sensor_l
660	#endif
661	ENDIF
662
663	CALL cpu_log( log_point(65), 'flight_measurement', 'stop' )
664
665	END SUBROUTINE flight_measurement
666
667	!------------------------------------------------------------------------------!
668	! Description:
669	! ------------
670	!> This subroutine bi-linearly interpolates the respective data onto the current
671	!> flight position.
672	!------------------------------------------------------------------------------!
673	SUBROUTINE interpolate_xyz( var, z_uw, ddz_uw, fac, x, y, var_ind, j, i )
674
675	USE control_parameters, &
676	ONLY: dz, dz_stretch_level
677
678	USE grid_variables, &
679	ONLY: dx, dy
680
681	USE indices, &
682	ONLY: nzb, nzt, nxlg, nxrg, nysg, nyng
683
684	IMPLICIT NONE
685
686	INTEGER(iwp) :: i !< index along x
687	INTEGER(iwp) :: j !< index along y
688	INTEGER(iwp) :: k !< index along z
689	INTEGER(iwp) :: k1 !< dummy variable
690	INTEGER(iwp) :: var_ind !< index variable for output quantity
691
692	REAL(wp) :: aa !< dummy argument for horizontal interpolation
693	REAL(wp) :: bb !< dummy argument for horizontal interpolation
694	REAL(wp) :: cc !< dummy argument for horizontal interpolation
695	REAL(wp) :: dd !< dummy argument for horizontal interpolation
696	REAL(wp) :: gg !< dummy argument for horizontal interpolation
697	REAL(wp) :: fac !< flag to indentify if quantity is on zu or zw level
698	REAL(wp) :: var_int !< horizontal interpolated variable at current position
699	REAL(wp) :: var_int_l !< horizontal interpolated variable at k-level
700	REAL(wp) :: var_int_u !< horizontal interpolated variable at (k+1)-level
701	REAL(wp) :: x !< distance between left edge of current grid box and flight position
702	REAL(wp) :: y !< distance between south edge of current grid box and flight position
703
704	REAL(wp), DIMENSION(1:nzt+1) :: ddz_uw !< inverse vertical grid spacing
705	REAL(wp), DIMENSION(nzb:nzt+1) :: z_uw !< height level
706
707	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< treted quantity
708	!
709	!-- Calculate interpolation coefficients
710	aa = x2 + y2
711	bb = ( dx - x )2 + y2
712	cc = x2 + ( dy - y )2
713	dd = ( dx - x )2 + ( dy - y )2
714	gg = aa + bb + cc + dd
715	!
716	!-- Obtain vertical index. Special treatment for grid index along z-direction
717	!-- if flight position is above the vertical grid-stretching level.
718	!-- fac=1 if variable is on scalar grid level, fac=0 for w-component.
719	IF ( z_pos(l) < dz_stretch_level ) THEN
720	k = ( z_pos(l) + fac * 0.5_wp * dz ) / dz
721	ELSE
722	!
723	!-- Search for k-index
724	DO k1=nzb, nzt
725	IF ( z_pos(l) >= z_uw(k1) .AND. z_pos(l) < z_uw(k1+1) ) THEN
726	k = k1
727	EXIT
728	ENDIF
729	ENDDO
730	ENDIF
731	!
732	!-- (x,y)-interpolation at k-level
733	var_int_l = ( ( gg - aa ) * var(k,j,i) + &
734	( gg - bb ) * var(k,j,i+1) + &
735	( gg - cc ) * var(k,j+1,i) + &
736	( gg - dd ) * var(k,j+1,i+1) &
737	) / ( 3.0_wp * gg )
738	!
739	!-- (x,y)-interpolation on (k+1)-level
740	var_int_u = ( ( gg - aa ) * var(k+1,j,i) + &
741	( gg - bb ) * var(k+1,j,i+1) + &
742	( gg - cc ) * var(k+1,j+1,i) + &
743	( gg - dd ) * var(k+1,j+1,i+1) &
744	) / ( 3.0_wp * gg )
745	!
746	!-- z-interpolation onto current flight postion
747	var_int = var_int_l &
748	+ ( z_pos(l) - z_uw(k) ) * ddz_uw(k+1) &
749	* (var_int_u - var_int_l )
750	!
751	!-- Store on local data array
752	sensor_l(var_ind,l) = var_int
753
754	END SUBROUTINE interpolate_xyz
755
756	!------------------------------------------------------------------------------!
757	! Description:
758	! ------------
759	!> Perform parameter checks.
760	!------------------------------------------------------------------------------!
761	SUBROUTINE flight_check_parameters
762
763	USE arrays_3d, &
764	ONLY: zu
765
766	USE control_parameters, &
767	ONLY: bc_lr_cyc, bc_ns_cyc, dz, message_string
768
769	USE grid_variables, &
770	ONLY: dx, dy
771
772	USE indices, &
773	ONLY: nx, ny, nz
774
775	USE netcdf_interface, &
776	ONLY: netcdf_data_format
777
778	IMPLICIT NONE
779
780	!
781	!-- Check if start positions are properly set.
782	DO l=1, num_leg
783	IF ( x_start(l) < 0.0_wp .OR. x_start(l) > ( nx + 1 ) * dx ) THEN
784	message_string = 'Start x position is outside the model domain'
785	CALL message( 'flight_check_parameters', 'PA0431', 1, 2, 0, 6, 0 )
786	ENDIF
787	IF ( y_start(l) < 0.0_wp .OR. y_start(l) > ( ny + 1 ) * dy ) THEN
788	message_string = 'Start y position is outside the model domain'
789	CALL message( 'flight_check_parameters', 'PA0432', 1, 2, 0, 6, 0 )
790	ENDIF
791
792	ENDDO
793	!
794	!-- Check for leg mode
795	DO l=1, num_leg
796	!
797	!-- Check if leg mode matches the overall lateral model boundary
798	!-- conditions.
799	IF ( TRIM( leg_mode(l) ) == 'cyclic' ) THEN
800	IF ( .NOT. bc_lr_cyc .OR. .NOT. bc_ns_cyc ) THEN
801	message_string = 'Cyclic flight leg does not match ' // &
802	'lateral boundary condition'
803	CALL message( 'flight_check_parameters', 'PA0433', 1, 2, 0, 6, 0 )
804	ENDIF
805	!
806	!-- Check if end-positions are inside the model domain in case of
807	!.. return-legs.
808	ELSEIF ( TRIM( leg_mode(l) ) == 'return' ) THEN
809	IF ( x_end(l) > ( nx + 1 ) * dx .OR. &
810	y_end(l) > ( ny + 1 ) * dx ) THEN
811	message_string = 'Flight leg or parts of it are outside ' // &
812	'the model domain'
813	CALL message( 'flight_check_parameters', 'PA0434', 1, 2, 0, 6, 0 )
814	ENDIF
815	ELSE
816	message_string = 'Unknown flight mode'
817	CALL message( 'flight_check_parameters', 'PA0435', 1, 2, 0, 6, 0 )
818	ENDIF
819
820	ENDDO
821	!
822	!-- Check if start and end positions are properly set in case of return legs.
823	DO l=1, num_leg
824
825	IF ( x_start(l) > x_end(l) .AND. leg_mode(l) == 'return' ) THEN
826	message_string = 'x_start position must be <= x_end ' // &
827	'position for return legs'
828	CALL message( 'flight_check_parameters', 'PA0436', 1, 2, 0, 6, 0 )
829	ENDIF
830	IF ( y_start(l) > y_end(l) .AND. leg_mode(l) == 'return' ) THEN
831	message_string = 'y_start position must be <= y_end ' // &
832	'position for return legs'
833	CALL message( 'flight_check_parameters', 'PA0437', 1, 2, 0, 6, 0 )
834	ENDIF
835	ENDDO
836	!
837	!-- Check if given flight object remains inside model domain if a rate of
838	!-- climb / descent is prescribed.
839	DO l=1, num_leg
840	IF ( flight_level(l) + max_elev_change(l) > zu(nz) .OR. &
841	flight_level(l) + max_elev_change(l) <= 0.0_wp ) THEN
842	message_string = 'Flight level is outside the model domain '
843	CALL message( 'flight_check_parameters', 'PA0438', 1, 2, 0, 6, 0 )
844	ENDIF
845	ENDDO
846	!
847	!-- Check for appropriate NetCDF format. Definition of more than one
848	!-- unlimited dimension is unfortunately only possible with NetCDF4/HDF5.
849	IF ( netcdf_data_format <= 2 ) THEN
850	message_string = 'netcdf_data_format must be > 2'
851	CALL message( 'flight_check_parameters', 'PA0439', 1, 2, 0, 6, 0 )
852	ENDIF
853
854
855	END SUBROUTINE flight_check_parameters
856
857	!------------------------------------------------------------------------------!
858	! Description:
859	! ------------
860	!> Skipping the flight-module variables from restart-file (binary format).
861	!------------------------------------------------------------------------------!
862	SUBROUTINE flight_skip_var_list
863
864	IMPLICIT NONE
865
866	CHARACTER (LEN=1) :: cdum
867	CHARACTER (LEN=30) :: variable_chr
868
869	READ ( 13 ) variable_chr
870	DO WHILE ( TRIM( variable_chr ) /= '* end flight *' )
871	READ ( 13 ) cdum
872	READ ( 13 ) variable_chr
873	ENDDO
874
875	END SUBROUTINE flight_skip_var_list
876
877	!------------------------------------------------------------------------------!
878	! Description:
879	! ------------
880	!> This routine reads the respective restart data.
881	!------------------------------------------------------------------------------!
882	SUBROUTINE flight_read_restart_data
883
884
885	IMPLICIT NONE
886
887	CHARACTER (LEN=30) :: variable_chr !< dummy variable to read string
888
889
890	READ ( 13 ) variable_chr
891	DO WHILE ( TRIM( variable_chr ) /= '* end flight *' )
892
893	SELECT CASE ( TRIM( variable_chr ) )
894
895	CASE ( 'u_agl' )
896	IF ( .NOT. ALLOCATED( u_agl ) ) ALLOCATE( u_agl(1:num_leg) )
897	READ ( 13 ) u_agl
898	CASE ( 'v_agl' )
899	IF ( .NOT. ALLOCATED( v_agl ) ) ALLOCATE( v_agl(1:num_leg) )
900	READ ( 13 ) v_agl
901	CASE ( 'w_agl' )
902	IF ( .NOT. ALLOCATED( w_agl ) ) ALLOCATE( w_agl(1:num_leg) )
903	READ ( 13 ) w_agl
904	CASE ( 'x_pos' )
905	IF ( .NOT. ALLOCATED( x_pos ) ) ALLOCATE( x_pos(1:num_leg) )
906	READ ( 13 ) x_pos
907	CASE ( 'y_pos' )
908	IF ( .NOT. ALLOCATED( y_pos ) ) ALLOCATE( y_pos(1:num_leg) )
909	READ ( 13 ) y_pos
910	CASE ( 'z_pos' )
911	IF ( .NOT. ALLOCATED( z_pos ) ) ALLOCATE( z_pos(1:num_leg) )
912	READ ( 13 ) z_pos
913
914	END SELECT
915
916	READ ( 13 ) variable_chr
917
918	ENDDO
919
920	END SUBROUTINE flight_read_restart_data
921
922	!------------------------------------------------------------------------------!
923	! Description:
924	! ------------
925	!> This routine writes the respective restart data.
926	!------------------------------------------------------------------------------!
927	SUBROUTINE flight_write_restart_data
928
929	IMPLICIT NONE
930
931	WRITE ( 14 ) 'u_agl '
932	WRITE ( 14 ) u_agl
933	WRITE ( 14 ) 'v_agl '
934	WRITE ( 14 ) v_agl
935	WRITE ( 14 ) 'w_agl '
936	WRITE ( 14 ) w_agl
937	WRITE ( 14 ) 'x_pos '
938	WRITE ( 14 ) x_pos
939	WRITE ( 14 ) 'y_pos '
940	WRITE ( 14 ) y_pos
941	WRITE ( 14 ) 'z_pos '
942	WRITE ( 14 ) z_pos
943
944	WRITE ( 14 ) '* end flight * '
945
946	END SUBROUTINE flight_write_restart_data
947
948
949	END MODULE flight_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |