Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |