Home

Context Navigation

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

Last change on this file since 2287 was 2271, checked in by sward, 7 years ago
error messages and numbers updated
Property svn:keywords set to `Id`
File size: 37.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |