Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |