Home

Context Navigation

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

Last change on this file since 4760 was 4760, checked in by schwenkel, 4 years ago
rename variable
Property svn:keywords set to `Id`
File size: 71.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |