Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |