Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |