Home

Context Navigation

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

Last change on this file since 2865 was 2776, checked in by Giersch, 7 years ago
Skipping of module related restart data changed + adapting synthetic turbulence generator to current restart procedure
Property svn:keywords set to `Id`
File size: 38.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |