Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |