Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |