Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |