Home

Context Navigation

source: palm/trunk/SOURCE/dynamics_mod.f90 @ 4844

Last change on this file since 4844 was 4843, checked in by raasch, 4 years ago
local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates
Property svn:keywords set to `Id`
File size: 72.2 KB

Rev	Line
[4047]	1	!> @file dynamics_mod.f90
	2	!--------------------------------------------------------------------------------------------------!
	3	! This file is part of the PALM model system.
	4	!
[4626]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[4047]	8	!
[4626]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[4047]	12	!
[4626]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[4047]	15	!
[4828]	16	! Copyright 1997-2021 Leibniz Universitaet Hannover
[4047]	17	!--------------------------------------------------------------------------------------------------!
	18	!
	19	! Current revisions:
	20	! -----------------
[4626]	21	!
	22	!
[4047]	23	! Former revisions:
	24	! -----------------
	25	! $Id: dynamics_mod.f90 4843 2021-01-15 15:22:11Z raasch $
[4843]	26	! local namelist parameter added to switch off the module although the respective module namelist
	27	! appears in the namelist file
	28	!
	29	! 4842 2021-01-14 10:42:28Z raasch
[4842]	30	! reading of namelist file and actions in case of namelist errors revised so that statement labels
	31	! and goto statements are not required any more
	32	!
	33	! 4828 2021-01-05 11:21:41Z Giersch
[4768]	34	! Enable 3D data output also with 64-bit precision
	35	!
	36	! 4760 2020-10-26 13:26:47Z schwenkel
[4757]	37	! Implement relative humidity as diagnostic output quantity
	38	!
	39	! 4731 2020-10-07 13:25:11Z schwenkel
[4731]	40	! Move exchange_horiz from time_integration to modules
	41	!
	42	! 4627 2020-07-26 10:14:44Z raasch
[4627]	43	! bugfix for r4626
	44	!
	45	! 4626 2020-07-26 09:49:48Z raasch
[4626]	46	! file re-formatted to follow the PALM coding standard
	47	!
	48	! 4517 2020-05-03 14:29:30Z raasch
[4517]	49	! added restart with MPI-IO for reading local arrays
[4626]	50	!
[4517]	51	! 4505 2020-04-20 15:37:15Z schwenkel
[4505]	52	! Add flag for saturation check
[4626]	53	!
[4505]	54	! 4495 2020-04-13 20:11:20Z resler
[4495]	55	! restart data handling with MPI-IO added
[4626]	56	!
[4495]	57	! 4360 2020-01-07 11:25:50Z suehring
[4360]	58	! Bugfix for last commit.
[4626]	59	!
[4360]	60	! 4359 2019-12-30 13:36:50Z suehring
[4358]	61	! Refine post-initialization check for realistically inital values of mixing ratio. Give an error
[4626]	62	! message for faulty initial values, but only a warning in a restart run.
	63	!
[4358]	64	! 4347 2019-12-18 13:18:33Z suehring
[4347]	65	! Implement post-initialization check for realistically inital values of mixing ratio
[4626]	66	!
[4347]	67	! 4281 2019-10-29 15:15:39Z schwenkel
[4281]	68	! Moved boundary conditions in dynamics module
[4626]	69	!
[4281]	70	! 4097 2019-07-15 11:59:11Z suehring
[4097]	71	! Avoid overlong lines - limit is 132 characters per line
	72	!
	73	! 4047 2019-06-21 18:58:09Z knoop
[4047]	74	! Initial introduction of the dynamics module with only dynamics_swap_timelevel implemented
	75	!
	76	!
	77	! Description:
	78	! ------------
	79	!> This module contains the dynamics of PALM.
	80	!--------------------------------------------------------------------------------------------------!
	81	MODULE dynamics_mod
	82
	83
[4626]	84	USE arrays_3d, &
	85	ONLY: c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l, &
[4731]	86	diss, &
	87	diss_p, &
[4626]	88	dzu, &
[4731]	89	e, &
	90	e_p, &
[4626]	91	exner, &
	92	hyp, &
	93	pt, pt_1, pt_2, pt_init, pt_p, &
	94	q, q_1, q_2, q_p, &
	95	s, s_1, s_2, s_p, &
	96	u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s, &
	97	v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s, &
[4757]	98	w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s, &
	99	zu
[4047]	100
[4347]	101	USE basic_constants_and_equations_mod, &
	102	ONLY: magnus, &
	103	rd_d_rv
	104
[4626]	105	USE control_parameters, &
	106	ONLY: bc_dirichlet_l, &
	107	bc_dirichlet_s, &
	108	bc_radiation_l, &
	109	bc_radiation_n, &
	110	bc_radiation_r, &
	111	bc_radiation_s, &
	112	bc_pt_t_val, &
	113	bc_q_t_val, &
	114	bc_s_t_val, &
	115	check_realistic_q, &
	116	child_domain, &
	117	coupling_mode, &
[4731]	118	constant_diffusion, &
[4626]	119	dt_3d, &
	120	humidity, &
	121	ibc_pt_b, &
	122	ibc_pt_t, &
	123	ibc_q_b, &
	124	ibc_q_t, &
	125	ibc_s_b, &
	126	ibc_s_t, &
	127	ibc_uv_b, &
	128	ibc_uv_t, &
	129	initializing_actions, &
	130	intermediate_timestep_count, &
	131	length, &
	132	message_string, &
	133	nesting_offline, &
	134	neutral, &
	135	nudging, &
	136	passive_scalar, &
	137	restart_string, &
[4731]	138	rans_mode, &
	139	rans_tke_e, &
[4626]	140	tsc, &
[4281]	141	use_cmax
[4047]	142
[4731]	143	USE exchange_horiz_mod, &
	144	ONLY: exchange_horiz
	145
	146
[4626]	147	USE grid_variables, &
	148	ONLY: ddx, &
	149	ddy, &
	150	dx, &
[4281]	151	dy
	152
[4626]	153	USE indices, &
	154	ONLY: nbgp, &
	155	nx, &
	156	nxl, &
	157	nxlg, &
	158	nxr, &
	159	nxrg, &
	160	ny, &
	161	nys, &
	162	nysg, &
	163	nyn, &
	164	nyng, &
	165	nzb, &
[4047]	166	nzt
	167
	168	USE kinds
	169
[4281]	170	USE pegrid
	171
[4626]	172	USE pmc_interface, &
[4281]	173	ONLY : nesting_mode
	174
[4495]	175	! USE restart_data_mpi_io_mod, &
	176	! ONLY:
	177
[4626]	178	USE surface_mod, &
[4281]	179	ONLY : bc_h
	180
	181
[4047]	182	IMPLICIT NONE
	183
	184	LOGICAL :: dynamics_module_enabled = .FALSE. !<
	185
	186	SAVE
	187
	188	PRIVATE
	189
	190	!
	191	!-- Public functions
[4626]	192	PUBLIC &
	193	dynamics_parin, &
	194	dynamics_check_parameters, &
	195	dynamics_check_data_output_ts, &
	196	dynamics_check_data_output_pr, &
	197	dynamics_check_data_output, &
	198	dynamics_init_masks, &
	199	dynamics_define_netcdf_grid, &
	200	dynamics_init_arrays, &
	201	dynamics_init, &
	202	dynamics_init_checks, &
	203	dynamics_header, &
	204	dynamics_actions, &
	205	dynamics_non_advective_processes, &
	206	dynamics_exchange_horiz, &
	207	dynamics_prognostic_equations, &
	208	dynamics_boundary_conditions, &
	209	dynamics_swap_timelevel, &
	210	dynamics_3d_data_averaging, &
	211	dynamics_data_output_2d, &
	212	dynamics_data_output_3d, &
	213	dynamics_statistics, &
	214	dynamics_rrd_global, &
	215	dynamics_rrd_local, &
	216	dynamics_wrd_global, &
	217	dynamics_wrd_local, &
[4047]	218	dynamics_last_actions
	219
	220	!
	221	!-- Public parameters, constants and initial values
[4626]	222	PUBLIC &
[4047]	223	dynamics_module_enabled
	224
	225	INTERFACE dynamics_parin
	226	MODULE PROCEDURE dynamics_parin
	227	END INTERFACE dynamics_parin
	228
	229	INTERFACE dynamics_check_parameters
	230	MODULE PROCEDURE dynamics_check_parameters
	231	END INTERFACE dynamics_check_parameters
	232
	233	INTERFACE dynamics_check_data_output_ts
	234	MODULE PROCEDURE dynamics_check_data_output_ts
	235	END INTERFACE dynamics_check_data_output_ts
	236
	237	INTERFACE dynamics_check_data_output_pr
	238	MODULE PROCEDURE dynamics_check_data_output_pr
	239	END INTERFACE dynamics_check_data_output_pr
	240
	241	INTERFACE dynamics_check_data_output
	242	MODULE PROCEDURE dynamics_check_data_output
	243	END INTERFACE dynamics_check_data_output
	244
	245	INTERFACE dynamics_init_masks
	246	MODULE PROCEDURE dynamics_init_masks
	247	END INTERFACE dynamics_init_masks
	248
	249	INTERFACE dynamics_define_netcdf_grid
	250	MODULE PROCEDURE dynamics_define_netcdf_grid
	251	END INTERFACE dynamics_define_netcdf_grid
	252
	253	INTERFACE dynamics_init_arrays
	254	MODULE PROCEDURE dynamics_init_arrays
	255	END INTERFACE dynamics_init_arrays
	256
	257	INTERFACE dynamics_init
	258	MODULE PROCEDURE dynamics_init
	259	END INTERFACE dynamics_init
	260
	261	INTERFACE dynamics_init_checks
	262	MODULE PROCEDURE dynamics_init_checks
	263	END INTERFACE dynamics_init_checks
	264
	265	INTERFACE dynamics_header
	266	MODULE PROCEDURE dynamics_header
	267	END INTERFACE dynamics_header
	268
	269	INTERFACE dynamics_actions
	270	MODULE PROCEDURE dynamics_actions
	271	MODULE PROCEDURE dynamics_actions_ij
	272	END INTERFACE dynamics_actions
	273
	274	INTERFACE dynamics_non_advective_processes
	275	MODULE PROCEDURE dynamics_non_advective_processes
	276	MODULE PROCEDURE dynamics_non_advective_processes_ij
	277	END INTERFACE dynamics_non_advective_processes
	278
	279	INTERFACE dynamics_exchange_horiz
	280	MODULE PROCEDURE dynamics_exchange_horiz
	281	END INTERFACE dynamics_exchange_horiz
	282
	283	INTERFACE dynamics_prognostic_equations
	284	MODULE PROCEDURE dynamics_prognostic_equations
	285	MODULE PROCEDURE dynamics_prognostic_equations_ij
	286	END INTERFACE dynamics_prognostic_equations
	287
[4281]	288	INTERFACE dynamics_boundary_conditions
	289	MODULE PROCEDURE dynamics_boundary_conditions
	290	END INTERFACE dynamics_boundary_conditions
	291
[4047]	292	INTERFACE dynamics_swap_timelevel
	293	MODULE PROCEDURE dynamics_swap_timelevel
	294	END INTERFACE dynamics_swap_timelevel
	295
	296	INTERFACE dynamics_3d_data_averaging
	297	MODULE PROCEDURE dynamics_3d_data_averaging
	298	END INTERFACE dynamics_3d_data_averaging
	299
	300	INTERFACE dynamics_data_output_2d
	301	MODULE PROCEDURE dynamics_data_output_2d
	302	END INTERFACE dynamics_data_output_2d
	303
	304	INTERFACE dynamics_data_output_3d
	305	MODULE PROCEDURE dynamics_data_output_3d
	306	END INTERFACE dynamics_data_output_3d
	307
	308	INTERFACE dynamics_statistics
	309	MODULE PROCEDURE dynamics_statistics
	310	END INTERFACE dynamics_statistics
	311
	312	INTERFACE dynamics_rrd_global
[4495]	313	MODULE PROCEDURE dynamics_rrd_global_ftn
	314	MODULE PROCEDURE dynamics_rrd_global_mpi
[4047]	315	END INTERFACE dynamics_rrd_global
	316
	317	INTERFACE dynamics_rrd_local
[4517]	318	MODULE PROCEDURE dynamics_rrd_local_ftn
	319	MODULE PROCEDURE dynamics_rrd_local_mpi
[4047]	320	END INTERFACE dynamics_rrd_local
	321
	322	INTERFACE dynamics_wrd_global
	323	MODULE PROCEDURE dynamics_wrd_global
	324	END INTERFACE dynamics_wrd_global
	325
	326	INTERFACE dynamics_wrd_local
	327	MODULE PROCEDURE dynamics_wrd_local
	328	END INTERFACE dynamics_wrd_local
	329
	330	INTERFACE dynamics_last_actions
	331	MODULE PROCEDURE dynamics_last_actions
	332	END INTERFACE dynamics_last_actions
	333
	334
	335	CONTAINS
	336
	337
	338	!--------------------------------------------------------------------------------------------------!
	339	! Description:
	340	! ------------
	341	!> Read module-specific namelist
	342	!--------------------------------------------------------------------------------------------------!
	343	SUBROUTINE dynamics_parin
	344
[4842]	345	CHARACTER(LEN=100) :: line !< dummy string that contains the current line of the parameter
	346	!< file
	347	INTEGER(iwp) :: io_status !< status after reading the namelist file
[4047]	348
[4843]	349	LOGICAL :: switch_off_module = .FALSE. !< local namelist parameter to switch off the module
	350	!< although the respective module namelist appears in
	351	!< the namelist file
[4842]	352
[4843]	353	NAMELIST /dynamics_parameters/ switch_off_module
[4047]	354
[4843]	355
[4047]	356	!
[4843]	357	!-- For the time beeing (unless the dynamics module is further developed), set default module
	358	!-- switch to true.
[4047]	359	dynamics_module_enabled = .TRUE.
	360
[4843]	361	!
[4842]	362	!-- Move to the beginning of the namelist file and try to find and read the namelist.
	363	REWIND( 11 )
	364	READ( 11, dynamics_parameters, IOSTAT=io_status )
[4047]	365
[4842]	366	!
	367	!-- Action depending on the READ status
[4843]	368	IF ( io_status == 0 ) THEN
[4842]	369	!
[4843]	370	!-- dynamics_parameters namelist was found and read correctly.
	371	IF ( .NOT. switch_off_module ) dynamics_module_enabled = .TRUE.
	372
	373	ELSEIF ( io_status > 0 ) THEN
	374	!
[4842]	375	!-- dynamics_parameters namelist was found, but contained errors. Print an error message
	376	!-- including the line that caused the problem.
	377	BACKSPACE( 11 )
	378	READ( 11 , '(A)') line
	379	CALL parin_fail_message( 'dynamics_parameters', line )
[4047]	380
[4842]	381	ENDIF
[4047]	382
	383	END SUBROUTINE dynamics_parin
	384
	385
	386	!--------------------------------------------------------------------------------------------------!
	387	! Description:
	388	! ------------
	389	!> Check control parameters and deduce further quantities.
	390	!--------------------------------------------------------------------------------------------------!
	391	SUBROUTINE dynamics_check_parameters
	392
	393
	394	END SUBROUTINE dynamics_check_parameters
	395
	396
	397	!--------------------------------------------------------------------------------------------------!
	398	! Description:
	399	! ------------
	400	!> Set module-specific timeseries units and labels
	401	!--------------------------------------------------------------------------------------------------!
	402	SUBROUTINE dynamics_check_data_output_ts( dots_max, dots_num, dots_label, dots_unit )
	403
[4627]	404	INTEGER(iwp), INTENT(IN) :: dots_max
	405
[4626]	406	CHARACTER (LEN=*), DIMENSION(dots_max), INTENT(INOUT) :: dots_label
	407	CHARACTER (LEN=*), DIMENSION(dots_max), INTENT(INOUT) :: dots_unit
[4047]	408
	409	INTEGER(iwp), INTENT(INOUT) :: dots_num
	410
	411	!
	412	!-- Next line is to avoid compiler warning about unused variables. Please remove.
	413	IF ( dots_num == 0 .OR. dots_label(1)(1:1) == ' ' .OR. dots_unit(1)(1:1) == ' ' ) CONTINUE
	414
	415
	416	END SUBROUTINE dynamics_check_data_output_ts
	417
	418
	419	!--------------------------------------------------------------------------------------------------!
	420	! Description:
	421	! ------------
	422	!> Set the unit of module-specific profile output quantities. For those variables not recognized,
	423	!> the parameter unit is set to "illegal", which tells the calling routine that the output variable
	424	!> is not defined and leads to a program abort.
	425	!--------------------------------------------------------------------------------------------------!
	426	SUBROUTINE dynamics_check_data_output_pr( variable, var_count, unit, dopr_unit )
	427
	428
[4626]	429	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
[4047]	430	CHARACTER (LEN=*) :: unit !<
	431	CHARACTER (LEN=*) :: variable !<
	432
	433	INTEGER(iwp) :: var_count !<
	434
	435	!
	436	!-- Next line is to avoid compiler warning about unused variables. Please remove.
	437	IF ( unit(1:1) == ' ' .OR. dopr_unit(1:1) == ' ' .OR. var_count == 0 ) CONTINUE
	438
	439	SELECT CASE ( TRIM( variable ) )
	440
	441	! CASE ( 'var_name' )
	442
	443	CASE DEFAULT
	444	unit = 'illegal'
	445
	446	END SELECT
	447
	448
	449	END SUBROUTINE dynamics_check_data_output_pr
	450
	451
	452	!--------------------------------------------------------------------------------------------------!
	453	! Description:
	454	! ------------
	455	!> Set the unit of module-specific output quantities. For those variables not recognized,
	456	!> the parameter unit is set to "illegal", which tells the calling routine that the output variable
	457	!< is not defined and leads to a program abort.
	458	!--------------------------------------------------------------------------------------------------!
	459	SUBROUTINE dynamics_check_data_output( variable, unit )
	460
	461
	462	CHARACTER (LEN=*) :: unit !<
	463	CHARACTER (LEN=*) :: variable !<
	464
	465	SELECT CASE ( TRIM( variable ) )
	466
	467	! CASE ( 'u2' )
	468
	469	CASE DEFAULT
	470	unit = 'illegal'
	471
	472	END SELECT
	473
	474
	475	END SUBROUTINE dynamics_check_data_output
	476
	477
[4626]	478	!--------------------------------------------------------------------------------------------------!
[4047]	479	!
	480	! Description:
	481	! ------------
	482	!> Initialize module-specific masked output
[4626]	483	!--------------------------------------------------------------------------------------------------!
[4047]	484	SUBROUTINE dynamics_init_masks( variable, unit )
	485
	486
	487	CHARACTER (LEN=*) :: unit !<
	488	CHARACTER (LEN=*) :: variable !<
	489
	490
	491	SELECT CASE ( TRIM( variable ) )
	492
	493	! CASE ( 'u2' )
	494
	495	CASE DEFAULT
	496	unit = 'illegal'
	497
	498	END SELECT
	499
	500
	501	END SUBROUTINE dynamics_init_masks
	502
	503
	504	!--------------------------------------------------------------------------------------------------!
	505	! Description:
	506	! ------------
	507	!> Initialize module-specific arrays
	508	!--------------------------------------------------------------------------------------------------!
	509	SUBROUTINE dynamics_init_arrays
	510
	511
	512	END SUBROUTINE dynamics_init_arrays
	513
	514
	515	!--------------------------------------------------------------------------------------------------!
	516	! Description:
	517	! ------------
	518	!> Execution of module-specific initializing actions
	519	!--------------------------------------------------------------------------------------------------!
	520	SUBROUTINE dynamics_init
	521
	522
	523	END SUBROUTINE dynamics_init
	524
	525
	526	!--------------------------------------------------------------------------------------------------!
	527	! Description:
	528	! ------------
	529	!> Perform module-specific post-initialization checks
	530	!--------------------------------------------------------------------------------------------------!
	531	SUBROUTINE dynamics_init_checks
	532
[4347]	533	INTEGER(iwp) :: i !< loop index in x-direction
	534	INTEGER(iwp) :: j !< loop index in y-direction
	535	INTEGER(iwp) :: k !< loop index in z-direction
[4047]	536
[4347]	537	LOGICAL :: realistic_q = .TRUE. !< flag indicating realistic mixing ratios
	538
	539	REAL(wp) :: e_s !< saturation water vapor pressure
	540	REAL(wp) :: q_s !< saturation mixing ratio
	541	REAL(wp) :: t_l !< actual temperature
[4760]	542	REAL(wp) :: rh_check = 9999999.9_wp !< relative humidity
[4757]	543	REAL(wp) :: rh_min = 9999999.9_wp !< max relative humidity
	544	REAL(wp) :: height = 9999999.9_wp !< height of supersaturated regions
	545	REAL(wp) :: min_height = 9999999.9_wp !< height of supersaturated regions
[4347]	546
	547	!
[4626]	548	!-- Check for realistic initial mixing ratio. This must be in a realistic phyiscial range and must
	549	!-- not exceed the saturation mixing ratio by more than 2 percent. Please note, the check is
	550	!-- performed for each grid point (not just for a vertical profile), in order to cover also
[4358]	551	!-- three-dimensional initialization. Note, this check gives an error only for the initial run not
	552	!-- for a restart run. In case there are no cloud physics considered, the mixing ratio can exceed
[4626]	553	!-- the saturation moisture. This case a warning is given.
[4505]	554	IF ( humidity .AND. .NOT. neutral .AND. check_realistic_q ) THEN
[4347]	555	DO i = nxl, nxr
	556	DO j = nys, nyn
	557	DO k = nzb+1, nzt
	558	!
	559	!-- Calculate actual temperature, water vapor saturation pressure, and based on this
	560	!-- the saturation mixing ratio.
	561	t_l = exner(k) * pt(k,j,i)
	562	e_s = magnus( t_l )
	563	q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
	564
[4757]	565	IF ( q(k,j,i) > 1.02_wp * q_s ) THEN
	566	realistic_q = .FALSE.
[4760]	567	rh_check = q(k,j,i) / q_s * 100.0_wp
[4757]	568	height = zu(k)
	569	ENDIF
[4347]	570	ENDDO
	571	ENDDO
	572	ENDDO
	573	!
[4626]	574	!-- Since the check is performed locally, merge the logical flag from all mpi ranks,
[4347]	575	!-- in order to do not print the error message multiple times.
	576	#if defined( __parallel )
	577	CALL MPI_ALLREDUCE( MPI_IN_PLACE, realistic_q, 1, MPI_LOGICAL, MPI_LAND, comm2d, ierr)
[4760]	578	CALL MPI_ALLREDUCE( rh_check, rh_min, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
[4757]	579	CALL MPI_ALLREDUCE( height, min_height, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
[4347]	580	#endif
	581
[4358]	582	IF ( .NOT. realistic_q .AND. &
	583	TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
[4757]	584	WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' // &
	585	'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time'
[4347]	586	CALL message( 'dynamic_init_checks', 'PA0697', 2, 2, 0, 6, 0 )
[4358]	587	ELSEIF ( .NOT. realistic_q .AND. &
	588	TRIM( initializing_actions ) == 'read_restart_data' ) THEN
[4757]	589	WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' // &
	590	'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time'
[4358]	591	CALL message( 'dynamic_init_checks', 'PA0697', 0, 1, 0, 6, 0 )
[4347]	592	ENDIF
	593	ENDIF
	594
[4047]	595	END SUBROUTINE dynamics_init_checks
	596
	597
	598	!--------------------------------------------------------------------------------------------------!
	599	! Description:
	600	! ------------
	601	!> Set the grids on which module-specific output quantities are defined. Allowed values for
	602	!> grid_x are "x" and "xu", for grid_y "y" and "yv", and for grid_z "zu" and "zw".
	603	!--------------------------------------------------------------------------------------------------!
	604	SUBROUTINE dynamics_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
	605
	606
	607	CHARACTER (LEN=*) :: grid_x !<
	608	CHARACTER (LEN=*) :: grid_y !<
	609	CHARACTER (LEN=*) :: grid_z !<
	610	CHARACTER (LEN=*) :: variable !<
	611
	612	LOGICAL :: found !<
	613
	614
	615	SELECT CASE ( TRIM( variable ) )
	616
	617	! CASE ( 'u2' )
	618
	619	CASE DEFAULT
	620	found = .FALSE.
	621	grid_x = 'none'
	622	grid_y = 'none'
	623	grid_z = 'none'
	624
	625	END SELECT
	626
	627
	628	END SUBROUTINE dynamics_define_netcdf_grid
	629
	630
	631	!--------------------------------------------------------------------------------------------------!
	632	! Description:
	633	! ------------
	634	!> Print a header with module-specific information.
	635	!--------------------------------------------------------------------------------------------------!
	636	SUBROUTINE dynamics_header( io )
	637
	638
	639	INTEGER(iwp) :: io !<
	640
	641	!
	642	!-- If no module-specific variables are read from the namelist-file, no information will be printed.
	643	IF ( .NOT. dynamics_module_enabled ) THEN
	644	WRITE ( io, 100 )
	645	RETURN
	646	ENDIF
	647
	648	!
	649	!-- Printing the information.
	650	WRITE ( io, 110 )
	651
	652	!
	653	!-- Format-descriptors
	654	100 FORMAT (//' *** dynamic module disabled'/)
[4626]	655	110 FORMAT (//1X,78('#') &
	656	//' User-defined variables and actions:'/ &
[4047]	657	' -----------------------------------'//)
	658
	659	END SUBROUTINE dynamics_header
	660
	661
	662	!--------------------------------------------------------------------------------------------------!
	663	! Description:
	664	! ------------
	665	!> Execute module-specific actions for all grid points
	666	!--------------------------------------------------------------------------------------------------!
	667	SUBROUTINE dynamics_actions( location )
	668
	669
	670	CHARACTER (LEN=*) :: location !<
	671
	672	! INTEGER(iwp) :: i !<
	673	! INTEGER(iwp) :: j !<
	674	! INTEGER(iwp) :: k !<
	675
	676	!
	677	!-- Here the user-defined actions follow
[4626]	678	!-- No calls for single grid points are allowed at locations before and after the timestep, since
	679	!-- these calls are not within an i,j-loop
[4047]	680	SELECT CASE ( location )
	681
	682	CASE ( 'before_timestep' )
	683
	684
	685	CASE ( 'before_prognostic_equations' )
	686
	687
	688	CASE ( 'after_integration' )
	689
	690
	691	CASE ( 'after_timestep' )
	692
	693
	694	CASE ( 'u-tendency' )
	695
	696
	697	CASE ( 'v-tendency' )
	698
	699
	700	CASE ( 'w-tendency' )
	701
	702
	703	CASE ( 'pt-tendency' )
	704
	705
	706	CASE ( 'sa-tendency' )
	707
	708
	709	CASE ( 'e-tendency' )
	710
	711
	712	CASE ( 'q-tendency' )
	713
	714
	715	CASE ( 's-tendency' )
	716
	717
	718	CASE DEFAULT
	719	CONTINUE
	720
	721	END SELECT
	722
	723	END SUBROUTINE dynamics_actions
	724
	725
	726	!--------------------------------------------------------------------------------------------------!
	727	! Description:
	728	! ------------
	729	!> Execute module-specific actions for grid point i,j
	730	!--------------------------------------------------------------------------------------------------!
	731	SUBROUTINE dynamics_actions_ij( i, j, location )
	732
	733
	734	CHARACTER (LEN=*) :: location
	735
	736	INTEGER(iwp) :: i
	737	INTEGER(iwp) :: j
	738
	739	!
	740	!-- Here the user-defined actions follow
	741	SELECT CASE ( location )
	742
	743	CASE ( 'u-tendency' )
	744
[4626]	745	!
[4047]	746	!-- Next line is to avoid compiler warning about unused variables. Please remove.
	747	IF ( i + j < 0 ) CONTINUE
	748
	749	CASE ( 'v-tendency' )
	750
	751
	752	CASE ( 'w-tendency' )
	753
	754
	755	CASE ( 'pt-tendency' )
	756
	757
	758	CASE ( 'sa-tendency' )
	759
	760
	761	CASE ( 'e-tendency' )
	762
	763
	764	CASE ( 'q-tendency' )
	765
	766
	767	CASE ( 's-tendency' )
	768
	769
	770	CASE DEFAULT
	771	CONTINUE
	772
	773	END SELECT
	774
	775	END SUBROUTINE dynamics_actions_ij
	776
	777
	778	!--------------------------------------------------------------------------------------------------!
	779	! Description:
	780	! ------------
	781	!> Compute module-specific non-advective processes for all grid points
	782	!--------------------------------------------------------------------------------------------------!
	783	SUBROUTINE dynamics_non_advective_processes
	784
	785
	786
	787	END SUBROUTINE dynamics_non_advective_processes
	788
	789
	790	!--------------------------------------------------------------------------------------------------!
	791	! Description:
	792	! ------------
	793	!> Compute module-specific non-advective processes for grid points i,j
	794	!--------------------------------------------------------------------------------------------------!
	795	SUBROUTINE dynamics_non_advective_processes_ij( i, j )
	796
	797
	798	INTEGER(iwp) :: i !<
	799	INTEGER(iwp) :: j !<
	800
	801	!
	802	!-- Next line is just to avoid compiler warnings about unused variables. You may remove it.
	803	IF ( i + j < 0 ) CONTINUE
	804
	805
	806	END SUBROUTINE dynamics_non_advective_processes_ij
	807
	808
	809	!--------------------------------------------------------------------------------------------------!
	810	! Description:
	811	! ------------
	812	!> Perform module-specific horizontal boundary exchange
	813	!--------------------------------------------------------------------------------------------------!
[4731]	814	SUBROUTINE dynamics_exchange_horiz( location )
[4047]	815
[4731]	816	CHARACTER (LEN=*), INTENT(IN) :: location !< call location string
[4047]	817
[4731]	818	SELECT CASE ( location )
[4047]	819
[4731]	820	CASE ( 'before_prognostic_equation' )
	821
	822	CASE ( 'after_prognostic_equation' )
	823
	824	CALL exchange_horiz( u_p, nbgp )
	825	CALL exchange_horiz( v_p, nbgp )
	826	CALL exchange_horiz( w_p, nbgp )
	827	CALL exchange_horiz( pt_p, nbgp )
	828	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e_p, nbgp )
	829	IF ( rans_tke_e ) CALL exchange_horiz( diss_p, nbgp )
	830	IF ( humidity ) CALL exchange_horiz( q_p, nbgp )
	831	IF ( passive_scalar ) CALL exchange_horiz( s_p, nbgp )
	832
	833	CASE ( 'after_anterpolation' )
	834
	835	CALL exchange_horiz( u, nbgp )
	836	CALL exchange_horiz( v, nbgp )
	837	CALL exchange_horiz( w, nbgp )
	838	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
	839	IF ( humidity ) CALL exchange_horiz( q, nbgp )
	840	IF ( passive_scalar ) CALL exchange_horiz( s, nbgp )
	841	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e, nbgp )
	842	IF ( .NOT. constant_diffusion .AND. rans_mode .AND. rans_tke_e ) THEN
	843	CALL exchange_horiz( diss, nbgp )
	844	ENDIF
	845
	846	END SELECT
	847
[4047]	848	END SUBROUTINE dynamics_exchange_horiz
	849
	850
	851	!--------------------------------------------------------------------------------------------------!
	852	! Description:
	853	! ------------
	854	!> Compute module-specific prognostic equations for all grid points
	855	!--------------------------------------------------------------------------------------------------!
	856	SUBROUTINE dynamics_prognostic_equations
	857
	858
	859
	860	END SUBROUTINE dynamics_prognostic_equations
	861
	862
	863	!--------------------------------------------------------------------------------------------------!
	864	! Description:
	865	! ------------
	866	!> Compute module-specific prognostic equations for grid point i,j
	867	!--------------------------------------------------------------------------------------------------!
	868	SUBROUTINE dynamics_prognostic_equations_ij( i, j, i_omp_start, tn )
	869
	870
	871	INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction
[4626]	872	INTEGER(iwp), INTENT(IN) :: i_omp_start !< first loop index of i-loop in prognostic_equations
[4047]	873	INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction
	874	INTEGER(iwp), INTENT(IN) :: tn !< task number of openmp task
	875
	876	!
	877	!-- Next line is just to avoid compiler warnings about unused variables. You may remove it.
	878	IF ( i + j + i_omp_start + tn < 0 ) CONTINUE
	879
	880	END SUBROUTINE dynamics_prognostic_equations_ij
	881
	882
[4281]	883	!--------------------------------------------------------------------------------------------------!
	884	! Description:
	885	! ------------
	886	!> Compute boundary conditions of dynamics model
	887	!--------------------------------------------------------------------------------------------------!
	888	SUBROUTINE dynamics_boundary_conditions
	889
	890	IMPLICIT NONE
	891
	892	INTEGER(iwp) :: i !< grid index x direction
	893	INTEGER(iwp) :: j !< grid index y direction
	894	INTEGER(iwp) :: k !< grid index z direction
	895	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls
	896	INTEGER(iwp) :: m !< running index surface elements
	897
	898	REAL(wp) :: c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
	899	REAL(wp) :: denom !< horizontal gradient of velocity component normal to the outflow boundary
	900
	901	!
	902	!-- Bottom boundary
	903	IF ( ibc_uv_b == 1 ) THEN
	904	u_p(nzb,:,:) = u_p(nzb+1,:,:)
	905	v_p(nzb,:,:) = v_p(nzb+1,:,:)
	906	ENDIF
	907	!
	908	!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
	909	!-- of downward-facing surfaces.
	910	DO l = 0, 1
	911	!$OMP PARALLEL DO PRIVATE( i, j, k )
	912	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	913	!$ACC PRESENT(bc_h, w_p)
	914	DO m = 1, bc_h(l)%ns
	915	i = bc_h(l)%i(m)
	916	j = bc_h(l)%j(m)
	917	k = bc_h(l)%k(m)
	918	w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	919	ENDDO
	920	ENDDO
	921
	922	!
	923	!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
	924	IF ( ibc_uv_t == 0 ) THEN
	925	!$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
	926	u_p(nzt+1,:,:) = u_init(nzt+1)
	927	v_p(nzt+1,:,:) = v_init(nzt+1)
	928	!$ACC END KERNELS
	929	ELSEIF ( ibc_uv_t == 1 ) THEN
	930	u_p(nzt+1,:,:) = u_p(nzt,:,:)
	931	v_p(nzt+1,:,:) = v_p(nzt,:,:)
	932	ENDIF
	933
	934	!
	935	!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
[4626]	936	IF ( .NOT. child_domain .AND. .NOT. nesting_offline .AND. &
	937	TRIM(coupling_mode) /= 'vnested_fine' ) THEN
[4281]	938	!$ACC KERNELS PRESENT(w_p)
	939	w_p(nzt:nzt+1,:,:) = 0.0_wp !< nzt is not a prognostic level (but cf. pres)
	940	!$ACC END KERNELS
	941	ENDIF
	942
	943	!
	944	!-- Temperature at bottom and top boundary.
[4626]	945	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by the sea surface temperature
	946	!-- of the coupled ocean model.
[4281]	947	!-- Dirichlet
	948	IF ( .NOT. neutral ) THEN
	949	IF ( ibc_pt_b == 0 ) THEN
	950	DO l = 0, 1
	951	!$OMP PARALLEL DO PRIVATE( i, j, k )
	952	DO m = 1, bc_h(l)%ns
	953	i = bc_h(l)%i(m)
	954	j = bc_h(l)%j(m)
	955	k = bc_h(l)%k(m)
	956	pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
	957	ENDDO
	958	ENDDO
	959	!
	960	!-- Neumann, zero-gradient
	961	ELSEIF ( ibc_pt_b == 1 ) THEN
	962	DO l = 0, 1
	963	!$OMP PARALLEL DO PRIVATE( i, j, k )
	964	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	965	!$ACC PRESENT(bc_h, pt_p)
	966	DO m = 1, bc_h(l)%ns
	967	i = bc_h(l)%i(m)
	968	j = bc_h(l)%j(m)
	969	k = bc_h(l)%k(m)
	970	pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
	971	ENDDO
	972	ENDDO
	973	ENDIF
	974
	975	!
	976	!-- Temperature at top boundary
	977	IF ( ibc_pt_t == 0 ) THEN
	978	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
	979	!
	980	!-- In case of nudging adjust top boundary to pt which is
	981	!-- read in from NUDGING-DATA
	982	IF ( nudging ) THEN
	983	pt_p(nzt+1,:,:) = pt_init(nzt+1)
	984	ENDIF
	985	ELSEIF ( ibc_pt_t == 1 ) THEN
	986	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
	987	ELSEIF ( ibc_pt_t == 2 ) THEN
	988	!$ACC KERNELS PRESENT(pt_p, dzu)
	989	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
	990	!$ACC END KERNELS
	991	ENDIF
	992	ENDIF
	993	!
[4626]	994	!-- Boundary conditions for total water content, bottom and top boundary (see also temperature)
[4281]	995	IF ( humidity ) THEN
	996	!
	997	!-- Surface conditions for constant_humidity_flux
[4626]	998	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
	999	!-- belongs to the atmospheric grid point, therefore, set q_p at k-1
[4281]	1000	IF ( ibc_q_b == 0 ) THEN
	1001
	1002	DO l = 0, 1
	1003	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1004	DO m = 1, bc_h(l)%ns
	1005	i = bc_h(l)%i(m)
	1006	j = bc_h(l)%j(m)
	1007	k = bc_h(l)%k(m)
	1008	q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
	1009	ENDDO
	1010	ENDDO
	1011
	1012	ELSE
	1013
	1014	DO l = 0, 1
	1015	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1016	DO m = 1, bc_h(l)%ns
	1017	i = bc_h(l)%i(m)
	1018	j = bc_h(l)%j(m)
	1019	k = bc_h(l)%k(m)
	1020	q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
	1021	ENDDO
	1022	ENDDO
	1023	ENDIF
	1024	!
	1025	!-- Top boundary
	1026	IF ( ibc_q_t == 0 ) THEN
	1027	q_p(nzt+1,:,:) = q(nzt+1,:,:)
	1028	ELSEIF ( ibc_q_t == 1 ) THEN
	1029	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
	1030	ENDIF
	1031	ENDIF
	1032	!
[4626]	1033	!-- Boundary conditions for scalar, bottom and top boundary (see also temperature)
[4281]	1034	IF ( passive_scalar ) THEN
	1035	!
	1036	!-- Surface conditions for constant_humidity_flux
[4626]	1037	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
	1038	!-- belongs to the atmospheric grid point, therefore, set s_p at k-1
[4281]	1039	IF ( ibc_s_b == 0 ) THEN
	1040
	1041	DO l = 0, 1
	1042	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1043	DO m = 1, bc_h(l)%ns
	1044	i = bc_h(l)%i(m)
	1045	j = bc_h(l)%j(m)
	1046	k = bc_h(l)%k(m)
	1047	s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
	1048	ENDDO
	1049	ENDDO
	1050
	1051	ELSE
	1052
	1053	DO l = 0, 1
	1054	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1055	DO m = 1, bc_h(l)%ns
	1056	i = bc_h(l)%i(m)
	1057	j = bc_h(l)%j(m)
	1058	k = bc_h(l)%k(m)
	1059	s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
	1060	ENDDO
	1061	ENDDO
	1062	ENDIF
	1063	!
	1064	!-- Top boundary condition
	1065	IF ( ibc_s_t == 0 ) THEN
	1066	s_p(nzt+1,:,:) = s(nzt+1,:,:)
	1067	ELSEIF ( ibc_s_t == 1 ) THEN
	1068	s_p(nzt+1,:,:) = s_p(nzt,:,:)
	1069	ELSEIF ( ibc_s_t == 2 ) THEN
	1070	s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
	1071	ENDIF
	1072
	1073	ENDIF
	1074	!
[4626]	1075	!-- In case of inflow or nest boundary at the south boundary the boundary for v is at nys and in
	1076	!-- case of inflow or nest boundary at the left boundary the boundary for u is at nxl. Since in
	1077	!-- prognostic_equations (cache optimized version) these levels are handled as a prognostic level,
	1078	!-- boundary values have to be restored here.
[4281]	1079	IF ( bc_dirichlet_s ) THEN
	1080	v_p(:,nys,:) = v_p(:,nys-1,:)
	1081	ELSEIF ( bc_dirichlet_l ) THEN
	1082	u_p(:,:,nxl) = u_p(:,:,nxl-1)
	1083	ENDIF
	1084
	1085	!
[4626]	1086	!-- The same restoration for u at i=nxl and v at j=nys as above must be made in case of nest
	1087	!-- boundaries. This must not be done in case of vertical nesting mode as in that case the lateral
	1088	!-- boundaries are actually cyclic.
	1089	!-- Lateral oundary conditions for TKE and dissipation are set in tcm_boundary_conds.
[4281]	1090	IF ( nesting_mode /= 'vertical' .OR. nesting_offline ) THEN
	1091	IF ( bc_dirichlet_s ) THEN
	1092	v_p(:,nys,:) = v_p(:,nys-1,:)
	1093	ENDIF
	1094	IF ( bc_dirichlet_l ) THEN
	1095	u_p(:,:,nxl) = u_p(:,:,nxl-1)
	1096	ENDIF
	1097	ENDIF
	1098
	1099	!
	1100	!-- Lateral boundary conditions for scalar quantities at the outflow.
[4626]	1101	!-- Lateral oundary conditions for TKE and dissipation are set in tcm_boundary_conds.
[4281]	1102	IF ( bc_radiation_s ) THEN
	1103	pt_p(:,nys-1,:) = pt_p(:,nys,:)
	1104	IF ( humidity ) THEN
	1105	q_p(:,nys-1,:) = q_p(:,nys,:)
	1106	ENDIF
	1107	IF ( passive_scalar ) s_p(:,nys-1,:) = s_p(:,nys,:)
	1108	ELSEIF ( bc_radiation_n ) THEN
	1109	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
	1110	IF ( humidity ) THEN
	1111	q_p(:,nyn+1,:) = q_p(:,nyn,:)
	1112	ENDIF
	1113	IF ( passive_scalar ) s_p(:,nyn+1,:) = s_p(:,nyn,:)
	1114	ELSEIF ( bc_radiation_l ) THEN
	1115	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
	1116	IF ( humidity ) THEN
	1117	q_p(:,:,nxl-1) = q_p(:,:,nxl)
	1118	ENDIF
	1119	IF ( passive_scalar ) s_p(:,:,nxl-1) = s_p(:,:,nxl)
	1120	ELSEIF ( bc_radiation_r ) THEN
	1121	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
	1122	IF ( humidity ) THEN
	1123	q_p(:,:,nxr+1) = q_p(:,:,nxr)
	1124	ENDIF
	1125	IF ( passive_scalar ) s_p(:,:,nxr+1) = s_p(:,:,nxr)
	1126	ENDIF
	1127
	1128	!
	1129	!-- Radiation boundary conditions for the velocities at the respective outflow.
[4626]	1130	!-- The phase velocity is either assumed to the maximum phase velocity that ensures numerical
	1131	!-- stability (CFL-condition) or calculated after Orlanski(1976) and averaged along the outflow
	1132	!-- boundary.
[4281]	1133	IF ( bc_radiation_s ) THEN
	1134
	1135	IF ( use_cmax ) THEN
	1136	u_p(:,-1,:) = u(:,0,:)
	1137	v_p(:,0,:) = v(:,1,:)
	1138	w_p(:,-1,:) = w(:,0,:)
	1139	ELSEIF ( .NOT. use_cmax ) THEN
	1140
	1141	c_max = dy / dt_3d
	1142
	1143	c_u_m_l = 0.0_wp
	1144	c_v_m_l = 0.0_wp
	1145	c_w_m_l = 0.0_wp
	1146
	1147	c_u_m = 0.0_wp
	1148	c_v_m = 0.0_wp
	1149	c_w_m = 0.0_wp
	1150
	1151	!
[4626]	1152	!-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow
	1153	!-- boundary.
[4281]	1154	DO k = nzb+1, nzt+1
	1155	DO i = nxl, nxr
	1156
	1157	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
	1158
	1159	IF ( denom /= 0.0_wp ) THEN
	1160	c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
	1161	IF ( c_u(k,i) < 0.0_wp ) THEN
	1162	c_u(k,i) = 0.0_wp
	1163	ELSEIF ( c_u(k,i) > c_max ) THEN
	1164	c_u(k,i) = c_max
	1165	ENDIF
	1166	ELSE
	1167	c_u(k,i) = c_max
	1168	ENDIF
	1169
	1170	denom = v_m_s(k,1,i) - v_m_s(k,2,i)
	1171
	1172	IF ( denom /= 0.0_wp ) THEN
	1173	c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
	1174	IF ( c_v(k,i) < 0.0_wp ) THEN
	1175	c_v(k,i) = 0.0_wp
	1176	ELSEIF ( c_v(k,i) > c_max ) THEN
	1177	c_v(k,i) = c_max
	1178	ENDIF
	1179	ELSE
	1180	c_v(k,i) = c_max
	1181	ENDIF
	1182
	1183	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
	1184
	1185	IF ( denom /= 0.0_wp ) THEN
	1186	c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
	1187	IF ( c_w(k,i) < 0.0_wp ) THEN
	1188	c_w(k,i) = 0.0_wp
	1189	ELSEIF ( c_w(k,i) > c_max ) THEN
	1190	c_w(k,i) = c_max
	1191	ENDIF
	1192	ELSE
	1193	c_w(k,i) = c_max
	1194	ENDIF
	1195
	1196	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
	1197	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
	1198	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
	1199
	1200	ENDDO
	1201	ENDDO
	1202
	1203	#if defined( __parallel )
	1204	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[4626]	1205	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
	1206	ierr )
[4281]	1207	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[4626]	1208	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
	1209	ierr )
[4281]	1210	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[4626]	1211	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
	1212	ierr )
[4281]	1213	#else
	1214	c_u_m = c_u_m_l
	1215	c_v_m = c_v_m_l
	1216	c_w_m = c_w_m_l
	1217	#endif
	1218
	1219	c_u_m = c_u_m / (nx+1)
	1220	c_v_m = c_v_m / (nx+1)
	1221	c_w_m = c_w_m / (nx+1)
	1222
	1223	!
	1224	!-- Save old timelevels for the next timestep
	1225	IF ( intermediate_timestep_count == 1 ) THEN
	1226	u_m_s(:,:,:) = u(:,0:1,:)
	1227	v_m_s(:,:,:) = v(:,1:2,:)
	1228	w_m_s(:,:,:) = w(:,0:1,:)
	1229	ENDIF
	1230
	1231	!
	1232	!-- Calculate the new velocities
	1233	DO k = nzb+1, nzt+1
	1234	DO i = nxlg, nxrg
[4626]	1235	u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * &
	1236	( u(k,-1,i) - u(k,0,i) ) * ddy
[4281]	1237
[4626]	1238	v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * &
	1239	( v(k,0,i) - v(k,1,i) ) * ddy
[4281]	1240
[4626]	1241	w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * &
	1242	( w(k,-1,i) - w(k,0,i) ) * ddy
[4281]	1243	ENDDO
	1244	ENDDO
	1245
	1246	!
	1247	!-- Bottom boundary at the outflow
	1248	IF ( ibc_uv_b == 0 ) THEN
	1249	u_p(nzb,-1,:) = 0.0_wp
	1250	v_p(nzb,0,:) = 0.0_wp
	1251	ELSE
	1252	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
	1253	v_p(nzb,0,:) = v_p(nzb+1,0,:)
	1254	ENDIF
	1255	w_p(nzb,-1,:) = 0.0_wp
	1256
	1257	!
	1258	!-- Top boundary at the outflow
	1259	IF ( ibc_uv_t == 0 ) THEN
	1260	u_p(nzt+1,-1,:) = u_init(nzt+1)
	1261	v_p(nzt+1,0,:) = v_init(nzt+1)
	1262	ELSE
	1263	u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
	1264	v_p(nzt+1,0,:) = v_p(nzt,0,:)
	1265	ENDIF
	1266	w_p(nzt:nzt+1,-1,:) = 0.0_wp
	1267
	1268	ENDIF
	1269
	1270	ENDIF
	1271
	1272	IF ( bc_radiation_n ) THEN
	1273
	1274	IF ( use_cmax ) THEN
	1275	u_p(:,ny+1,:) = u(:,ny,:)
	1276	v_p(:,ny+1,:) = v(:,ny,:)
	1277	w_p(:,ny+1,:) = w(:,ny,:)
	1278	ELSEIF ( .NOT. use_cmax ) THEN
	1279
	1280	c_max = dy / dt_3d
	1281
	1282	c_u_m_l = 0.0_wp
	1283	c_v_m_l = 0.0_wp
	1284	c_w_m_l = 0.0_wp
	1285
	1286	c_u_m = 0.0_wp
	1287	c_v_m = 0.0_wp
	1288	c_w_m = 0.0_wp
	1289
	1290	!
[4626]	1291	!-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow
	1292	!-- boundary.
[4281]	1293	DO k = nzb+1, nzt+1
	1294	DO i = nxl, nxr
	1295
	1296	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
	1297
	1298	IF ( denom /= 0.0_wp ) THEN
	1299	c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
	1300	IF ( c_u(k,i) < 0.0_wp ) THEN
	1301	c_u(k,i) = 0.0_wp
	1302	ELSEIF ( c_u(k,i) > c_max ) THEN
	1303	c_u(k,i) = c_max
	1304	ENDIF
	1305	ELSE
	1306	c_u(k,i) = c_max
	1307	ENDIF
	1308
	1309	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
	1310
	1311	IF ( denom /= 0.0_wp ) THEN
	1312	c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
	1313	IF ( c_v(k,i) < 0.0_wp ) THEN
	1314	c_v(k,i) = 0.0_wp
	1315	ELSEIF ( c_v(k,i) > c_max ) THEN
	1316	c_v(k,i) = c_max
	1317	ENDIF
	1318	ELSE
	1319	c_v(k,i) = c_max
	1320	ENDIF
	1321
	1322	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
	1323
	1324	IF ( denom /= 0.0_wp ) THEN
	1325	c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
	1326	IF ( c_w(k,i) < 0.0_wp ) THEN
	1327	c_w(k,i) = 0.0_wp
	1328	ELSEIF ( c_w(k,i) > c_max ) THEN
	1329	c_w(k,i) = c_max
	1330	ENDIF
	1331	ELSE
	1332	c_w(k,i) = c_max
	1333	ENDIF
	1334
	1335	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
	1336	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
	1337	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
	1338
	1339	ENDDO
	1340	ENDDO
	1341
	1342	#if defined( __parallel )
	1343	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[4626]	1344	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
	1345	ierr )
[4281]	1346	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[4626]	1347	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
	1348	ierr )
[4281]	1349	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[4626]	1350	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
	1351	ierr )
[4281]	1352	#else
	1353	c_u_m = c_u_m_l
	1354	c_v_m = c_v_m_l
	1355	c_w_m = c_w_m_l
	1356	#endif
	1357
	1358	c_u_m = c_u_m / (nx+1)
	1359	c_v_m = c_v_m / (nx+1)
	1360	c_w_m = c_w_m / (nx+1)
	1361
	1362	!
	1363	!-- Save old timelevels for the next timestep
	1364	IF ( intermediate_timestep_count == 1 ) THEN
	1365	u_m_n(:,:,:) = u(:,ny-1:ny,:)
	1366	v_m_n(:,:,:) = v(:,ny-1:ny,:)
	1367	w_m_n(:,:,:) = w(:,ny-1:ny,:)
	1368	ENDIF
	1369
	1370	!
	1371	!-- Calculate the new velocities
	1372	DO k = nzb+1, nzt+1
	1373	DO i = nxlg, nxrg
[4626]	1374	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * &
	1375	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
[4281]	1376
[4626]	1377	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * &
	1378	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
[4281]	1379
[4626]	1380	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * &
	1381	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
[4281]	1382	ENDDO
	1383	ENDDO
	1384
	1385	!
	1386	!-- Bottom boundary at the outflow
	1387	IF ( ibc_uv_b == 0 ) THEN
	1388	u_p(nzb,ny+1,:) = 0.0_wp
	1389	v_p(nzb,ny+1,:) = 0.0_wp
	1390	ELSE
	1391	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
	1392	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
	1393	ENDIF
	1394	w_p(nzb,ny+1,:) = 0.0_wp
	1395
	1396	!
	1397	!-- Top boundary at the outflow
	1398	IF ( ibc_uv_t == 0 ) THEN
	1399	u_p(nzt+1,ny+1,:) = u_init(nzt+1)
	1400	v_p(nzt+1,ny+1,:) = v_init(nzt+1)
	1401	ELSE
	1402	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
	1403	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
	1404	ENDIF
	1405	w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
	1406
	1407	ENDIF
	1408
	1409	ENDIF
	1410
	1411	IF ( bc_radiation_l ) THEN
	1412
	1413	IF ( use_cmax ) THEN
	1414	u_p(:,:,0) = u(:,:,1)
	1415	v_p(:,:,-1) = v(:,:,0)
	1416	w_p(:,:,-1) = w(:,:,0)
	1417	ELSEIF ( .NOT. use_cmax ) THEN
	1418
	1419	c_max = dx / dt_3d
	1420
	1421	c_u_m_l = 0.0_wp
	1422	c_v_m_l = 0.0_wp
	1423	c_w_m_l = 0.0_wp
	1424
	1425	c_u_m = 0.0_wp
	1426	c_v_m = 0.0_wp
	1427	c_w_m = 0.0_wp
	1428
	1429	!
[4626]	1430	!-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow
	1431	!-- boundary.
[4281]	1432	DO k = nzb+1, nzt+1
	1433	DO j = nys, nyn
	1434
	1435	denom = u_m_l(k,j,1) - u_m_l(k,j,2)
	1436
	1437	IF ( denom /= 0.0_wp ) THEN
	1438	c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
	1439	IF ( c_u(k,j) < 0.0_wp ) THEN
	1440	c_u(k,j) = 0.0_wp
	1441	ELSEIF ( c_u(k,j) > c_max ) THEN
	1442	c_u(k,j) = c_max
	1443	ENDIF
	1444	ELSE
	1445	c_u(k,j) = c_max
	1446	ENDIF
	1447
	1448	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
	1449
	1450	IF ( denom /= 0.0_wp ) THEN
	1451	c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
	1452	IF ( c_v(k,j) < 0.0_wp ) THEN
	1453	c_v(k,j) = 0.0_wp
	1454	ELSEIF ( c_v(k,j) > c_max ) THEN
	1455	c_v(k,j) = c_max
	1456	ENDIF
	1457	ELSE
	1458	c_v(k,j) = c_max
	1459	ENDIF
	1460
	1461	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
	1462
	1463	IF ( denom /= 0.0_wp ) THEN
	1464	c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
	1465	IF ( c_w(k,j) < 0.0_wp ) THEN
	1466	c_w(k,j) = 0.0_wp
	1467	ELSEIF ( c_w(k,j) > c_max ) THEN
	1468	c_w(k,j) = c_max
	1469	ENDIF
	1470	ELSE
	1471	c_w(k,j) = c_max
	1472	ENDIF
	1473
	1474	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
	1475	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
	1476	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
	1477
	1478	ENDDO
	1479	ENDDO
	1480
	1481	#if defined( __parallel )
	1482	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[4626]	1483	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
	1484	ierr )
[4281]	1485	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[4626]	1486	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
	1487	ierr )
[4281]	1488	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[4626]	1489	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
	1490	ierr )
[4281]	1491	#else
	1492	c_u_m = c_u_m_l
	1493	c_v_m = c_v_m_l
	1494	c_w_m = c_w_m_l
	1495	#endif
	1496
	1497	c_u_m = c_u_m / (ny+1)
	1498	c_v_m = c_v_m / (ny+1)
	1499	c_w_m = c_w_m / (ny+1)
	1500
	1501	!
	1502	!-- Save old timelevels for the next timestep
	1503	IF ( intermediate_timestep_count == 1 ) THEN
	1504	u_m_l(:,:,:) = u(:,:,1:2)
	1505	v_m_l(:,:,:) = v(:,:,0:1)
	1506	w_m_l(:,:,:) = w(:,:,0:1)
	1507	ENDIF
	1508
	1509	!
	1510	!-- Calculate the new velocities
	1511	DO k = nzb+1, nzt+1
	1512	DO j = nysg, nyng
[4626]	1513	u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * &
	1514	( u(k,j,0) - u(k,j,1) ) * ddx
[4281]	1515
[4626]	1516	v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * &
	1517	( v(k,j,-1) - v(k,j,0) ) * ddx
[4281]	1518
[4626]	1519	w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * &
	1520	( w(k,j,-1) - w(k,j,0) ) * ddx
[4281]	1521	ENDDO
	1522	ENDDO
	1523
	1524	!
	1525	!-- Bottom boundary at the outflow
	1526	IF ( ibc_uv_b == 0 ) THEN
	1527	u_p(nzb,:,0) = 0.0_wp
	1528	v_p(nzb,:,-1) = 0.0_wp
	1529	ELSE
	1530	u_p(nzb,:,0) = u_p(nzb+1,:,0)
	1531	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
	1532	ENDIF
	1533	w_p(nzb,:,-1) = 0.0_wp
	1534
	1535	!
	1536	!-- Top boundary at the outflow
	1537	IF ( ibc_uv_t == 0 ) THEN
	1538	u_p(nzt+1,:,0) = u_init(nzt+1)
	1539	v_p(nzt+1,:,-1) = v_init(nzt+1)
	1540	ELSE
	1541	u_p(nzt+1,:,0) = u_p(nzt,:,0)
	1542	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
	1543	ENDIF
	1544	w_p(nzt:nzt+1,:,-1) = 0.0_wp
	1545
	1546	ENDIF
	1547
	1548	ENDIF
	1549
	1550	IF ( bc_radiation_r ) THEN
	1551
	1552	IF ( use_cmax ) THEN
	1553	u_p(:,:,nx+1) = u(:,:,nx)
	1554	v_p(:,:,nx+1) = v(:,:,nx)
	1555	w_p(:,:,nx+1) = w(:,:,nx)
	1556	ELSEIF ( .NOT. use_cmax ) THEN
	1557
	1558	c_max = dx / dt_3d
	1559
	1560	c_u_m_l = 0.0_wp
	1561	c_v_m_l = 0.0_wp
	1562	c_w_m_l = 0.0_wp
	1563
	1564	c_u_m = 0.0_wp
	1565	c_v_m = 0.0_wp
	1566	c_w_m = 0.0_wp
	1567
	1568	!
[4626]	1569	!-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow
	1570	!-- boundary.
[4281]	1571	DO k = nzb+1, nzt+1
	1572	DO j = nys, nyn
	1573
	1574	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
	1575
	1576	IF ( denom /= 0.0_wp ) THEN
	1577	c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
	1578	IF ( c_u(k,j) < 0.0_wp ) THEN
	1579	c_u(k,j) = 0.0_wp
	1580	ELSEIF ( c_u(k,j) > c_max ) THEN
	1581	c_u(k,j) = c_max
	1582	ENDIF
	1583	ELSE
	1584	c_u(k,j) = c_max
	1585	ENDIF
	1586
	1587	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
	1588
	1589	IF ( denom /= 0.0_wp ) THEN
	1590	c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
	1591	IF ( c_v(k,j) < 0.0_wp ) THEN
	1592	c_v(k,j) = 0.0_wp
	1593	ELSEIF ( c_v(k,j) > c_max ) THEN
	1594	c_v(k,j) = c_max
	1595	ENDIF
	1596	ELSE
	1597	c_v(k,j) = c_max
	1598	ENDIF
	1599
	1600	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
	1601
	1602	IF ( denom /= 0.0_wp ) THEN
	1603	c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
	1604	IF ( c_w(k,j) < 0.0_wp ) THEN
	1605	c_w(k,j) = 0.0_wp
	1606	ELSEIF ( c_w(k,j) > c_max ) THEN
	1607	c_w(k,j) = c_max
	1608	ENDIF
	1609	ELSE
	1610	c_w(k,j) = c_max
	1611	ENDIF
	1612
	1613	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
	1614	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
	1615	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
	1616
	1617	ENDDO
	1618	ENDDO
	1619
	1620	#if defined( __parallel )
	1621	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[4626]	1622	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
	1623	ierr )
[4281]	1624	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[4626]	1625	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
	1626	ierr )
[4281]	1627	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[4626]	1628	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
	1629	ierr )
[4281]	1630	#else
	1631	c_u_m = c_u_m_l
	1632	c_v_m = c_v_m_l
	1633	c_w_m = c_w_m_l
	1634	#endif
	1635
	1636	c_u_m = c_u_m / (ny+1)
	1637	c_v_m = c_v_m / (ny+1)
	1638	c_w_m = c_w_m / (ny+1)
	1639
	1640	!
	1641	!-- Save old timelevels for the next timestep
	1642	IF ( intermediate_timestep_count == 1 ) THEN
	1643	u_m_r(:,:,:) = u(:,:,nx-1:nx)
	1644	v_m_r(:,:,:) = v(:,:,nx-1:nx)
	1645	w_m_r(:,:,:) = w(:,:,nx-1:nx)
	1646	ENDIF
	1647
	1648	!
	1649	!-- Calculate the new velocities
	1650	DO k = nzb+1, nzt+1
	1651	DO j = nysg, nyng
[4626]	1652	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * &
	1653	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
[4281]	1654
[4626]	1655	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * &
	1656	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
[4281]	1657
[4626]	1658	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * &
	1659	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
[4281]	1660	ENDDO
	1661	ENDDO
	1662
	1663	!
	1664	!-- Bottom boundary at the outflow
	1665	IF ( ibc_uv_b == 0 ) THEN
	1666	u_p(nzb,:,nx+1) = 0.0_wp
	1667	v_p(nzb,:,nx+1) = 0.0_wp
	1668	ELSE
	1669	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
	1670	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
	1671	ENDIF
	1672	w_p(nzb,:,nx+1) = 0.0_wp
	1673
	1674	!
	1675	!-- Top boundary at the outflow
	1676	IF ( ibc_uv_t == 0 ) THEN
	1677	u_p(nzt+1,:,nx+1) = u_init(nzt+1)
	1678	v_p(nzt+1,:,nx+1) = v_init(nzt+1)
	1679	ELSE
	1680	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
	1681	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
	1682	ENDIF
	1683	w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
	1684
	1685	ENDIF
	1686
	1687	ENDIF
	1688
	1689	END SUBROUTINE dynamics_boundary_conditions
[4626]	1690	!--------------------------------------------------------------------------------------------------!
[4047]	1691	! Description:
	1692	! ------------
	1693	!> Swap timelevels of module-specific array pointers
[4626]	1694	!--------------------------------------------------------------------------------------------------!
[4047]	1695	SUBROUTINE dynamics_swap_timelevel ( mod_count )
	1696
	1697
	1698	INTEGER, INTENT(IN) :: mod_count
	1699
	1700
	1701	SELECT CASE ( mod_count )
	1702
	1703	CASE ( 0 )
	1704
	1705	u => u_1; u_p => u_2
	1706	v => v_1; v_p => v_2
	1707	w => w_1; w_p => w_2
	1708	IF ( .NOT. neutral ) THEN
	1709	pt => pt_1; pt_p => pt_2
	1710	ENDIF
	1711	IF ( humidity ) THEN
	1712	q => q_1; q_p => q_2
	1713	ENDIF
	1714	IF ( passive_scalar ) THEN
	1715	s => s_1; s_p => s_2
	1716	ENDIF
	1717
	1718	CASE ( 1 )
	1719
	1720	u => u_2; u_p => u_1
	1721	v => v_2; v_p => v_1
	1722	w => w_2; w_p => w_1
	1723	IF ( .NOT. neutral ) THEN
	1724	pt => pt_2; pt_p => pt_1
	1725	ENDIF
	1726	IF ( humidity ) THEN
	1727	q => q_2; q_p => q_1
	1728	ENDIF
	1729	IF ( passive_scalar ) THEN
	1730	s => s_2; s_p => s_1
	1731	ENDIF
	1732
	1733	END SELECT
	1734
	1735	END SUBROUTINE dynamics_swap_timelevel
	1736
	1737
	1738	!--------------------------------------------------------------------------------------------------!
	1739	! Description:
	1740	! ------------
[4626]	1741	!> Sum up and time-average module-specific output quantities as well as allocate the array necessary
	1742	!> for storing the average.
[4047]	1743	!--------------------------------------------------------------------------------------------------!
	1744	SUBROUTINE dynamics_3d_data_averaging( mode, variable )
	1745
	1746
	1747	CHARACTER (LEN=*) :: mode !<
	1748	CHARACTER (LEN=*) :: variable !<
	1749
	1750
	1751	IF ( mode == 'allocate' ) THEN
	1752
	1753	SELECT CASE ( TRIM( variable ) )
	1754
	1755	! CASE ( 'u2' )
	1756
	1757	CASE DEFAULT
	1758	CONTINUE
	1759
	1760	END SELECT
	1761
	1762	ELSEIF ( mode == 'sum' ) THEN
	1763
	1764	SELECT CASE ( TRIM( variable ) )
	1765
	1766	! CASE ( 'u2' )
	1767
	1768	CASE DEFAULT
	1769	CONTINUE
	1770
	1771	END SELECT
	1772
	1773	ELSEIF ( mode == 'average' ) THEN
	1774
	1775	SELECT CASE ( TRIM( variable ) )
	1776
	1777	! CASE ( 'u2' )
	1778
	1779	END SELECT
	1780
	1781	ENDIF
	1782
	1783
	1784	END SUBROUTINE dynamics_3d_data_averaging
	1785
	1786
	1787	!--------------------------------------------------------------------------------------------------!
	1788	! Description:
	1789	! ------------
[4626]	1790	!> Resorts the module-specific output quantity with indices (k,j,i) to a temporary array with
	1791	!> indices (i,j,k) and sets the grid on which it is defined.
[4047]	1792	!> Allowed values for grid are "zu" and "zw".
	1793	!--------------------------------------------------------------------------------------------------!
[4626]	1794	SUBROUTINE dynamics_data_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, &
	1795	nzt_do, fill_value )
[4047]	1796
	1797
[4626]	1798	CHARACTER (LEN=*) :: grid !<
[4047]	1799	CHARACTER (LEN=*), INTENT(IN) :: mode !< either 'xy', 'xz' or 'yz'
[4626]	1800	CHARACTER (LEN=*) :: variable !<
[4047]	1801
	1802	INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data
	1803	! INTEGER(iwp) :: i !< grid index along x-direction
	1804	! INTEGER(iwp) :: j !< grid index along y-direction
	1805	! INTEGER(iwp) :: k !< grid index along z-direction
	1806	! INTEGER(iwp) :: m !< running index surface elements
	1807	INTEGER(iwp) :: nzb_do !< lower limit of the domain (usually nzb)
	1808	INTEGER(iwp) :: nzt_do !< upper limit of the domain (usually nzt+1)
	1809
	1810	LOGICAL :: found !<
	1811	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
	1812
	1813	REAL(wp), INTENT(IN) :: fill_value
	1814
	1815	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
	1816
	1817	!
	1818	!-- Next line is just to avoid compiler warnings about unused variables. You may remove it.
	1819	IF ( two_d .AND. av + LEN( mode ) + local_pf(nxl,nys,nzb_do) + fill_value < 0.0 ) CONTINUE
	1820
	1821	found = .TRUE.
	1822
	1823	SELECT CASE ( TRIM( variable ) )
	1824
	1825	! CASE ( 'u2_xy', 'u2_xz', 'u2_yz' )
	1826
	1827	CASE DEFAULT
	1828	found = .FALSE.
	1829	grid = 'none'
	1830
	1831	END SELECT
	1832
	1833
	1834	END SUBROUTINE dynamics_data_output_2d
	1835
	1836
	1837	!--------------------------------------------------------------------------------------------------!
	1838	! Description:
	1839	! ------------
[4626]	1840	!> Resorts the module-specific output quantity with indices (k,j,i) to a temporary array with
	1841	!> indices (i,j,k).
[4047]	1842	!--------------------------------------------------------------------------------------------------!
	1843	SUBROUTINE dynamics_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
	1844
	1845
	1846	CHARACTER (LEN=*) :: variable !<
	1847
	1848	INTEGER(iwp) :: av !<
	1849	! INTEGER(iwp) :: i !<
	1850	! INTEGER(iwp) :: j !<
	1851	! INTEGER(iwp) :: k !<
	1852	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
	1853	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
	1854
	1855	LOGICAL :: found !<
	1856
	1857	REAL(wp), INTENT(IN) :: fill_value !< value for the _FillValue attribute
	1858
[4768]	1859	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
[4047]	1860
	1861	!
	1862	!-- Next line is to avoid compiler warning about unused variables. Please remove.
	1863	IF ( av + local_pf(nxl,nys,nzb_do) + fill_value < 0.0 ) CONTINUE
	1864
	1865
	1866	found = .TRUE.
	1867
	1868	SELECT CASE ( TRIM( variable ) )
	1869
	1870	! CASE ( 'u2' )
	1871
	1872	CASE DEFAULT
	1873	found = .FALSE.
	1874
	1875	END SELECT
	1876
	1877
	1878	END SUBROUTINE dynamics_data_output_3d
	1879
	1880
	1881	!--------------------------------------------------------------------------------------------------!
	1882	! Description:
	1883	! ------------
	1884	!> Calculation of module-specific statistics, i.e. horizontally averaged profiles and time series.
	1885	!> This is called for every statistic region sr, but at least for the region "total domain" (sr=0).
	1886	!--------------------------------------------------------------------------------------------------!
	1887	SUBROUTINE dynamics_statistics( mode, sr, tn )
	1888
	1889
	1890	CHARACTER (LEN=*) :: mode !<
	1891	! INTEGER(iwp) :: i !<
	1892	! INTEGER(iwp) :: j !<
	1893	! INTEGER(iwp) :: k !<
	1894	INTEGER(iwp) :: sr !<
	1895	INTEGER(iwp) :: tn !<
	1896
	1897	!
	1898	!-- Next line is to avoid compiler warning about unused variables. Please remove.
	1899	IF ( sr == 0 .OR. tn == 0 ) CONTINUE
	1900
	1901	IF ( mode == 'profiles' ) THEN
	1902
	1903	ELSEIF ( mode == 'time_series' ) THEN
	1904
	1905	ENDIF
	1906
	1907	END SUBROUTINE dynamics_statistics
	1908
	1909
	1910	!--------------------------------------------------------------------------------------------------!
	1911	! Description:
	1912	! ------------
[4495]	1913	!> Read module-specific global restart data (Fortran binary format).
[4047]	1914	!--------------------------------------------------------------------------------------------------!
[4495]	1915	SUBROUTINE dynamics_rrd_global_ftn( found )
[4047]	1916
	1917	LOGICAL, INTENT(OUT) :: found
	1918
	1919
	1920	found = .TRUE.
	1921
	1922
	1923	SELECT CASE ( restart_string(1:length) )
	1924
	1925	CASE ( 'global_paramter' )
	1926	! READ ( 13 ) global_parameter
	1927
	1928	CASE DEFAULT
	1929
	1930	found = .FALSE.
	1931
	1932	END SELECT
	1933
	1934
[4495]	1935	END SUBROUTINE dynamics_rrd_global_ftn
[4047]	1936
	1937
	1938	!--------------------------------------------------------------------------------------------------!
	1939	! Description:
	1940	! ------------
[4495]	1941	!> Read module-specific global restart data (MPI-IO).
	1942	!--------------------------------------------------------------------------------------------------!
	1943	SUBROUTINE dynamics_rrd_global_mpi
	1944
	1945	! CALL rrd_mpi_io( 'global_parameter', global_parameter )
	1946	CONTINUE
	1947
	1948	END SUBROUTINE dynamics_rrd_global_mpi
	1949
	1950
	1951	!--------------------------------------------------------------------------------------------------!
	1952	! Description:
	1953	! ------------
[4517]	1954	!> Read module-specific local restart data arrays (Fortran binary format).
[4047]	1955	!> Subdomain index limits on file are given by nxl_on_file, etc.
	1956	!> Indices nxlc, etc. indicate the range of gridpoints to be mapped from the subdomain on file (f)
	1957	!> to the subdomain of the current PE (c). They have been calculated in routine rrd_local.
	1958	!--------------------------------------------------------------------------------------------------!
[4626]	1959	SUBROUTINE dynamics_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, &
	1960	nync, nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, &
	1961	found )
[4047]	1962
	1963
	1964	INTEGER(iwp) :: k !<
	1965	INTEGER(iwp) :: nxlc !<
	1966	INTEGER(iwp) :: nxlf !<
	1967	INTEGER(iwp) :: nxl_on_file !<
	1968	INTEGER(iwp) :: nxrc !<
	1969	INTEGER(iwp) :: nxrf !<
	1970	INTEGER(iwp) :: nxr_on_file !<
	1971	INTEGER(iwp) :: nync !<
	1972	INTEGER(iwp) :: nynf !<
	1973	INTEGER(iwp) :: nyn_on_file !<
	1974	INTEGER(iwp) :: nysc !<
	1975	INTEGER(iwp) :: nysf !<
	1976	INTEGER(iwp) :: nys_on_file !<
	1977
	1978	LOGICAL, INTENT(OUT) :: found
	1979
	1980	REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !<
	1981	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<
	1982
	1983	!
	1984	!-- Next line is to avoid compiler warning about unused variables. Please remove.
[4626]	1985	IF ( k + nxlc + nxlf + nxrc + nxrf + nync + nynf + nysc + nysf + &
	1986	tmp_2d(nys_on_file,nxl_on_file) + &
[4097]	1987	tmp_3d(nzb,nys_on_file,nxl_on_file) < 0.0 ) CONTINUE
[4047]	1988	!
	1989	!-- Here the reading of user-defined restart data follows:
	1990	!-- Sample for user-defined output
	1991
	1992	found = .TRUE.
	1993
	1994	SELECT CASE ( restart_string(1:length) )
	1995
	1996	! CASE ( 'u2_av' )
	1997
	1998	CASE DEFAULT
	1999
	2000	found = .FALSE.
	2001
	2002	END SELECT
	2003
[4517]	2004	END SUBROUTINE dynamics_rrd_local_ftn
[4047]	2005
	2006
	2007	!--------------------------------------------------------------------------------------------------!
	2008	! Description:
	2009	! ------------
[4517]	2010	!> Read module-specific local restart data arrays (MPI-IO).
	2011	!--------------------------------------------------------------------------------------------------!
	2012	SUBROUTINE dynamics_rrd_local_mpi
	2013
	2014	IMPLICIT NONE
	2015
	2016	! LOGICAL :: array_found !<
	2017
	2018
	2019	! CALL rd_mpi_io_check_array( 'u2_av' , found = array_found )
	2020	! IF ( array_found ) THEN
	2021	! IF ( .NOT. ALLOCATED( u2_av ) ) ALLOCATE( u2_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2022	! CALL rrd_mpi_io( 'u2_av', u2_av )
	2023	! ENDIF
	2024
	2025	CONTINUE
	2026
	2027	END SUBROUTINE dynamics_rrd_local_mpi
	2028
	2029
	2030
	2031	!--------------------------------------------------------------------------------------------------!
	2032	! Description:
	2033	! ------------
[4047]	2034	!> Writes global module-specific restart data into binary file(s) for restart runs.
	2035	!--------------------------------------------------------------------------------------------------!
	2036	SUBROUTINE dynamics_wrd_global
	2037
	2038
	2039	END SUBROUTINE dynamics_wrd_global
	2040
	2041
	2042	!--------------------------------------------------------------------------------------------------!
	2043	! Description:
	2044	! ------------
	2045	!> Writes processor specific and module-specific restart data into binary file(s) for restart runs.
	2046	!--------------------------------------------------------------------------------------------------!
	2047	SUBROUTINE dynamics_wrd_local
	2048
	2049
	2050	END SUBROUTINE dynamics_wrd_local
	2051
	2052
	2053	!--------------------------------------------------------------------------------------------------!
	2054	! Description:
	2055	! ------------
	2056	!> Execute module-specific actions at the very end of the program.
	2057	!--------------------------------------------------------------------------------------------------!
	2058	SUBROUTINE dynamics_last_actions
	2059
	2060
	2061	END SUBROUTINE dynamics_last_actions
	2062
[4627]	2063	END MODULE dynamics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |