Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |