Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |