Home

Context Navigation

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

Last change on this file since 4566 was 4517, checked in by raasch, 5 years ago
added restart with MPI-IO for reading local arrays
Property svn:keywords set to `Id`
File size: 61.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |