Home

Context Navigation

virtual_flight_mod.f90 @ 4879

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

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

Context Navigation

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

Download in other formats: