Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |