Home

Context Navigation

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

Last change on this file since 2563 was 2563, checked in by Giersch, 7 years ago
Restart runs with the usage of the wind turbine model are possible now. Further small at reading/writing restart data
Property svn:keywords set to `Id`
File size: 37.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |