Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |