Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |