Home

Context Navigation

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

Last change on this file since 2772 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 38.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |