Home

Context Navigation

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

Last change on this file since 4001 was 3885, checked in by kanani, 6 years ago
restructure/add location/debug messages
Property svn:keywords set to `Id`
File size: 40.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |