Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |