Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |