Home

Context Navigation

source: palm/trunk/SOURCE/prognostic_equations.f90 @ 991

Last change on this file since 991 was 979, checked in by fricke, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 80.3 KB

Rev	Line
[736]	1	MODULE prognostic_equations_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Current revisions:
	5	! -----------------
[979]	6	!
	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: prognostic_equations.f90 979 2012-08-09 08:50:11Z heinze $
	11	!
	12	! 978 2012-08-09 08:28:32Z fricke
[978]	13	! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
	14	! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
	15	! boundary in case of non-cyclic lateral boundaries
	16	! Bugfix: first thread index changes for WS-scheme at the inflow
[736]	17	!
[941]	18	! 940 2012-07-09 14:31:00Z raasch
	19	! temperature equation can be switched off
	20	!
[786]	21	! 785 2011-11-28 09:47:19Z raasch
	22	! new factor rdf_sc allows separate Rayleigh damping of scalars
	23	!
[737]	24	! 736 2011-08-17 14:13:26Z suehring
	25	! Bugfix: determination of first thread index i for WS-scheme
	26	!
[736]	27	! 709 2011-03-30 09:31:40Z raasch
	28	! formatting adjustments
	29	!
	30	! 673 2011-01-18 16:19:48Z suehring
	31	! Consideration of the pressure gradient (steered by tsc(4)) during the time
	32	! integration removed.
	33	!
	34	! 667 2010-12-23 12:06:00Z suehring/gryschka
	35	! Calls of the advection routines with WS5 added.
	36	! Calls of ws_statistics added to set the statistical arrays to zero after each
	37	! time step.
	38	!
	39	! 531 2010-04-21 06:47:21Z heinze
	40	! add call of subsidence in the equation for humidity / passive scalar
	41	!
	42	! 411 2009-12-11 14:15:58Z heinze
	43	! add call of subsidence in the equation for potential temperature
	44	!
	45	! 388 2009-09-23 09:40:33Z raasch
	46	! prho is used instead of rho in diffusion_e,
	47	! external pressure gradient
	48	!
	49	! 153 2008-03-19 09:41:30Z steinfeld
	50	! add call of plant_canopy_model in the prognostic equation for
	51	! the potential temperature and for the passive scalar
	52	!
	53	! 138 2007-11-28 10:03:58Z letzel
	54	! add call of subroutines that evaluate the canopy drag terms,
	55	! add wall_*flux to parameter list of calls of diffusion_s
	56	!
	57	! 106 2007-08-16 14:30:26Z raasch
	58	! +uswst, vswst as arguments in calls of diffusion_u\|v,
	59	! loops for u and v are starting from index nxlu, nysv, respectively (needed
	60	! for non-cyclic boundary conditions)
	61	!
	62	! 97 2007-06-21 08:23:15Z raasch
	63	! prognostic equation for salinity, density is calculated from equation of
	64	! state for seawater and is used for calculation of buoyancy,
	65	! +eqn_state_seawater_mod
	66	! diffusion_e is called with argument rho in case of ocean runs,
	67	! new argument zw in calls of diffusion_e, new argument pt_/prho_reference
	68	! in calls of buoyancy and diffusion_e, calc_mean_pt_profile renamed
	69	! calc_mean_profile
	70	!
	71	! 75 2007-03-22 09:54:05Z raasch
	72	! checking for negative q and limiting for positive values,
	73	! z0 removed from arguments in calls of diffusion_u/v/w, uxrp, vynp eliminated,
	74	! subroutine names changed to .._noopt, .._cache, and .._vector,
	75	! moisture renamed humidity, Bott-Chlond-scheme can be used in the
	76	! _vector-version
	77	!
	78	! 19 2007-02-23 04:53:48Z raasch
	79	! Calculation of e, q, and pt extended for gridpoint nzt,
	80	! handling of given temperature/humidity/scalar fluxes at top surface
	81	!
	82	! RCS Log replace by Id keyword, revision history cleaned up
	83	!
	84	! Revision 1.21 2006/08/04 15:01:07 raasch
	85	! upstream scheme can be forced to be used for tke (use_upstream_for_tke)
	86	! regardless of the timestep scheme used for the other quantities,
	87	! new argument diss in call of diffusion_e
	88	!
	89	! Revision 1.1 2000/04/13 14:56:27 schroeter
	90	! Initial revision
	91	!
	92	!
	93	! Description:
	94	! ------------
	95	! Solving the prognostic equations.
	96	!------------------------------------------------------------------------------!
	97
	98	USE arrays_3d
	99	USE control_parameters
	100	USE cpulog
	101	USE eqn_state_seawater_mod
	102	USE grid_variables
	103	USE indices
	104	USE interfaces
	105	USE pegrid
	106	USE pointer_interfaces
	107	USE statistics
	108	USE advec_ws
	109	USE advec_s_pw_mod
	110	USE advec_s_up_mod
	111	USE advec_u_pw_mod
	112	USE advec_u_up_mod
	113	USE advec_v_pw_mod
	114	USE advec_v_up_mod
	115	USE advec_w_pw_mod
	116	USE advec_w_up_mod
	117	USE buoyancy_mod
	118	USE calc_precipitation_mod
	119	USE calc_radiation_mod
	120	USE coriolis_mod
	121	USE diffusion_e_mod
	122	USE diffusion_s_mod
	123	USE diffusion_u_mod
	124	USE diffusion_v_mod
	125	USE diffusion_w_mod
	126	USE impact_of_latent_heat_mod
	127	USE plant_canopy_model_mod
	128	USE production_e_mod
	129	USE subsidence_mod
	130	USE user_actions_mod
	131
	132
	133	PRIVATE
	134	PUBLIC prognostic_equations_noopt, prognostic_equations_cache, &
	135	prognostic_equations_vector
	136
	137	INTERFACE prognostic_equations_noopt
	138	MODULE PROCEDURE prognostic_equations_noopt
	139	END INTERFACE prognostic_equations_noopt
	140
	141	INTERFACE prognostic_equations_cache
	142	MODULE PROCEDURE prognostic_equations_cache
	143	END INTERFACE prognostic_equations_cache
	144
	145	INTERFACE prognostic_equations_vector
	146	MODULE PROCEDURE prognostic_equations_vector
	147	END INTERFACE prognostic_equations_vector
	148
	149
	150	CONTAINS
	151
	152
	153	SUBROUTINE prognostic_equations_noopt
	154
	155	!------------------------------------------------------------------------------!
	156	! Version with single loop optimization
	157	!
	158	! (Optimized over each single prognostic equation.)
	159	!------------------------------------------------------------------------------!
	160
	161	IMPLICIT NONE
	162
	163	CHARACTER (LEN=9) :: time_to_string
	164	INTEGER :: i, i_omp_start, j, k, tn = 0
	165	REAL :: sat, sbt
	166
	167	!
	168	!-- Calculate those variables needed in the tendency terms which need
	169	!-- global communication
[940]	170	IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 )
	171	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	172	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
[736]	173	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
	174	intermediate_timestep_count == 1 ) CALL ws_statistics
	175
	176	!
	177	!-- u-velocity component
	178	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	179
	180	!
	181	!-- u-tendency terms with communication
	182	IF ( momentum_advec == 'ups-scheme' ) THEN
	183	tend = 0.0
	184	CALL advec_u_ups
	185	ENDIF
	186
	187	!
	188	!-- u-tendency terms with no communication
	189	i_omp_start = nxlu
	190	DO i = nxlu, nxr
	191	DO j = nys, nyn
	192	!
	193	!-- Tendency terms
	194	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	195	tend(:,j,i) = 0.0
	196	IF ( ws_scheme_mom ) THEN
	197	CALL advec_u_ws( i, j, i_omp_start, tn )
	198	ELSE
	199	CALL advec_u_pw( i, j )
	200	ENDIF
	201
	202	ELSE
	203	IF ( momentum_advec /= 'ups-scheme' ) THEN
	204	tend(:,j,i) = 0.0
	205	CALL advec_u_up( i, j )
	206	ENDIF
	207	ENDIF
	208	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[978]	209	CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend, u_m, usws_m, &
	210	uswst_m, v_m, w_m )
[736]	211	ELSE
[978]	212	CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, uswst, &
	213	v, w )
[736]	214	ENDIF
	215	CALL coriolis( i, j, 1 )
[940]	216	IF ( sloping_surface .AND. .NOT. neutral ) THEN
	217	CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
	218	ENDIF
[736]	219
	220	!
	221	!-- Drag by plant canopy
	222	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 1 )
	223
	224	!
	225	!-- External pressure gradient
	226	IF ( dp_external ) THEN
	227	DO k = dp_level_ind_b+1, nzt
	228	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	229	ENDDO
	230	ENDIF
	231
	232	CALL user_actions( i, j, 'u-tendency' )
	233
	234	!
	235	!-- Prognostic equation for u-velocity component
	236	DO k = nzb_u_inner(j,i)+1, nzt
	237	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
	238	dt_3d * ( &
	239	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
	240	) - &
	241	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	242	ENDDO
	243
	244	!
	245	!-- Calculate tendencies for the next Runge-Kutta step
	246	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	247	IF ( intermediate_timestep_count == 1 ) THEN
	248	DO k = nzb_u_inner(j,i)+1, nzt
	249	tu_m(k,j,i) = tend(k,j,i)
	250	ENDDO
	251	ELSEIF ( intermediate_timestep_count < &
	252	intermediate_timestep_count_max ) THEN
	253	DO k = nzb_u_inner(j,i)+1, nzt
	254	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	255	ENDDO
	256	ENDIF
	257	ENDIF
	258
	259	ENDDO
	260	ENDDO
	261
	262	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	263
	264	!
	265	!-- v-velocity component
	266	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	267
	268	!
	269	!-- v-tendency terms with communication
	270	IF ( momentum_advec == 'ups-scheme' ) THEN
	271	tend = 0.0
	272	CALL advec_v_ups
	273	ENDIF
	274
	275	!
	276	!-- v-tendency terms with no communication
	277	i_omp_start = nxl
	278	DO i = nxl, nxr
	279	DO j = nysv, nyn
	280	!
	281	!-- Tendency terms
	282	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	283	tend(:,j,i) = 0.0
	284	IF ( ws_scheme_mom ) THEN
	285	CALL advec_v_ws( i, j, i_omp_start, tn )
	286	ELSE
	287	CALL advec_v_pw( i, j )
	288	ENDIF
	289
	290	ELSE
	291	IF ( momentum_advec /= 'ups-scheme' ) THEN
	292	tend(:,j,i) = 0.0
	293	CALL advec_v_up( i, j )
	294	ENDIF
	295	ENDIF
	296	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[978]	297	CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
	298	vsws_m, vswst_m, w_m )
[736]	299	ELSE
[978]	300	CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, &
	301	vswst, w )
[736]	302	ENDIF
	303	CALL coriolis( i, j, 2 )
	304
	305	!
	306	!-- Drag by plant canopy
	307	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 2 )
	308
	309	!
	310	!-- External pressure gradient
	311	IF ( dp_external ) THEN
	312	DO k = dp_level_ind_b+1, nzt
	313	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	314	ENDDO
	315	ENDIF
	316
	317	CALL user_actions( i, j, 'v-tendency' )
	318
	319	!
	320	!-- Prognostic equation for v-velocity component
	321	DO k = nzb_v_inner(j,i)+1, nzt
	322	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
	323	dt_3d * ( &
	324	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
	325	) - &
	326	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	327	ENDDO
	328
	329	!
	330	!-- Calculate tendencies for the next Runge-Kutta step
	331	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	332	IF ( intermediate_timestep_count == 1 ) THEN
	333	DO k = nzb_v_inner(j,i)+1, nzt
	334	tv_m(k,j,i) = tend(k,j,i)
	335	ENDDO
	336	ELSEIF ( intermediate_timestep_count < &
	337	intermediate_timestep_count_max ) THEN
	338	DO k = nzb_v_inner(j,i)+1, nzt
	339	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	340	ENDDO
	341	ENDIF
	342	ENDIF
	343
	344	ENDDO
	345	ENDDO
	346
	347	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	348
	349	!
	350	!-- w-velocity component
	351	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	352
	353	!
	354	!-- w-tendency terms with communication
	355	IF ( momentum_advec == 'ups-scheme' ) THEN
	356	tend = 0.0
	357	CALL advec_w_ups
	358	ENDIF
	359
	360	!
	361	!-- w-tendency terms with no communication
	362	DO i = nxl, nxr
	363	DO j = nys, nyn
	364	!
	365	!-- Tendency terms
	366	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	367	tend(:,j,i) = 0.0
	368	IF ( ws_scheme_mom ) THEN
	369	CALL advec_w_ws( i, j, i_omp_start, tn )
	370	ELSE
	371	CALL advec_w_pw( i, j )
	372	ENDIF
	373
	374	ELSE
	375	IF ( momentum_advec /= 'ups-scheme' ) THEN
	376	tend(:,j,i) = 0.0
	377	CALL advec_w_up( i, j )
	378	ENDIF
	379	ENDIF
	380	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[978]	381	CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, w_m )
[736]	382	ELSE
[978]	383	CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w )
[736]	384	ENDIF
	385	CALL coriolis( i, j, 3 )
[940]	386
	387	IF ( .NOT. neutral ) THEN
	388	IF ( ocean ) THEN
	389	CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
[736]	390	ELSE
[940]	391	IF ( .NOT. humidity ) THEN
	392	CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
	393	ELSE
	394	CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
	395	ENDIF
[736]	396	ENDIF
	397	ENDIF
	398
	399	!
	400	!-- Drag by plant canopy
	401	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 3 )
	402
	403	CALL user_actions( i, j, 'w-tendency' )
	404
	405	!
	406	!-- Prognostic equation for w-velocity component
	407	DO k = nzb_w_inner(j,i)+1, nzt-1
	408	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
	409	dt_3d * ( &
	410	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
	411	) - &
	412	tsc(5) * rdf(k) * w(k,j,i)
	413	ENDDO
	414
	415	!
	416	!-- Calculate tendencies for the next Runge-Kutta step
	417	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	418	IF ( intermediate_timestep_count == 1 ) THEN
	419	DO k = nzb_w_inner(j,i)+1, nzt-1
	420	tw_m(k,j,i) = tend(k,j,i)
	421	ENDDO
	422	ELSEIF ( intermediate_timestep_count < &
	423	intermediate_timestep_count_max ) THEN
	424	DO k = nzb_w_inner(j,i)+1, nzt-1
	425	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	426	ENDDO
	427	ENDIF
	428	ENDIF
	429
	430	ENDDO
	431	ENDDO
	432
	433	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	434
	435	!
[940]	436	!-- If required, compute prognostic equation for potential temperature
	437	IF ( .NOT. neutral ) THEN
[736]	438
[940]	439	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
[736]	440
[940]	441	!
	442	!-- pt-tendency terms with communication
	443	sat = tsc(1)
	444	sbt = tsc(2)
	445	IF ( scalar_advec == 'bc-scheme' ) THEN
	446
	447	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	448	!
	449	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	450	!-- switched on. Thus:
	451	sat = 1.0
	452	sbt = 1.0
	453	ENDIF
[736]	454	tend = 0.0
[940]	455	CALL advec_s_bc( pt, 'pt' )
	456	ELSE
	457	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
	458	tend = 0.0
	459	CALL advec_s_ups( pt, 'pt' )
	460	ENDIF
[736]	461	ENDIF
	462
[940]	463	!
	464	!-- pt-tendency terms with no communication
	465	DO i = nxl, nxr
	466	DO j = nys, nyn
	467	!
	468	!-- Tendency terms
	469	IF ( scalar_advec == 'bc-scheme' ) THEN
	470	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
	471	wall_heatflux, tend )
[736]	472	ELSE
[940]	473	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	474	THEN
[736]	475	tend(:,j,i) = 0.0
[940]	476	IF ( ws_scheme_sca ) THEN
	477	CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
	478	diss_s_pt, flux_l_pt, diss_l_pt, &
	479	i_omp_start, tn )
	480	ELSE
	481	CALL advec_s_pw( i, j, pt )
	482	ENDIF
	483	ELSE
	484	IF ( scalar_advec /= 'ups-scheme' ) THEN
	485	tend(:,j,i) = 0.0
	486	CALL advec_s_up( i, j, pt )
	487	ENDIF
[736]	488	ENDIF
[940]	489	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	490	THEN
	491	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
	492	tswst_m, wall_heatflux, tend )
	493	ELSE
	494	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
	495	wall_heatflux, tend )
	496	ENDIF
[736]	497	ENDIF
[940]	498
	499	!
	500	!-- If required compute heating/cooling due to long wave radiation
	501	!-- processes
	502	IF ( radiation ) THEN
	503	CALL calc_radiation( i, j )
[736]	504	ENDIF
	505
[940]	506	!
	507	!-- If required compute impact of latent heat due to precipitation
	508	IF ( precipitation ) THEN
	509	CALL impact_of_latent_heat( i, j )
	510	ENDIF
[736]	511
[940]	512	!
	513	!-- Consideration of heat sources within the plant canopy
	514	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	515	CALL plant_canopy_model( i, j, 4 )
	516	ENDIF
[736]	517
[940]	518	!
	519	!-- If required compute influence of large-scale subsidence/ascent
	520	IF ( large_scale_subsidence ) THEN
	521	CALL subsidence( i, j, tend, pt, pt_init )
	522	ENDIF
[736]	523
[940]	524	CALL user_actions( i, j, 'pt-tendency' )
[736]	525
[940]	526	!
	527	!-- Prognostic equation for potential temperature
	528	DO k = nzb_s_inner(j,i)+1, nzt
	529	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
[978]	530	dt_3d * ( &
[940]	531	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
[978]	532	) - &
	533	tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
	534	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
[940]	535	ENDDO
[736]	536
[940]	537	!
	538	!-- Calculate tendencies for the next Runge-Kutta step
	539	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	540	IF ( intermediate_timestep_count == 1 ) THEN
	541	DO k = nzb_s_inner(j,i)+1, nzt
	542	tpt_m(k,j,i) = tend(k,j,i)
	543	ENDDO
	544	ELSEIF ( intermediate_timestep_count < &
	545	intermediate_timestep_count_max ) THEN
	546	DO k = nzb_s_inner(j,i)+1, nzt
	547	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	548	5.3125 * tpt_m(k,j,i)
	549	ENDDO
	550	ENDIF
[736]	551	ENDIF
	552
[940]	553	ENDDO
[736]	554	ENDDO
	555
[940]	556	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
[736]	557
[940]	558	ENDIF
	559
[736]	560	!
	561	!-- If required, compute prognostic equation for salinity
	562	IF ( ocean ) THEN
	563
	564	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	565
	566	!
	567	!-- sa-tendency terms with communication
	568	sat = tsc(1)
	569	sbt = tsc(2)
	570	IF ( scalar_advec == 'bc-scheme' ) THEN
	571
	572	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	573	!
	574	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	575	!-- switched on. Thus:
	576	sat = 1.0
	577	sbt = 1.0
	578	ENDIF
	579	tend = 0.0
	580	CALL advec_s_bc( sa, 'sa' )
	581	ELSE
	582	IF ( tsc(2) /= 2.0 ) THEN
	583	IF ( scalar_advec == 'ups-scheme' ) THEN
	584	tend = 0.0
	585	CALL advec_s_ups( sa, 'sa' )
	586	ENDIF
	587	ENDIF
	588	ENDIF
	589
	590	!
	591	!-- sa terms with no communication
	592	DO i = nxl, nxr
	593	DO j = nys, nyn
	594	!
	595	!-- Tendency-terms
	596	IF ( scalar_advec == 'bc-scheme' ) THEN
	597	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	598	wall_salinityflux, tend )
	599	ELSE
	600	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	601	tend(:,j,i) = 0.0
	602	IF ( ws_scheme_sca ) THEN
	603	CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, &
	604	diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
	605	ELSE
	606	CALL advec_s_pw( i, j, sa )
	607	ENDIF
	608
	609	ELSE
	610	IF ( scalar_advec /= 'ups-scheme' ) THEN
	611	tend(:,j,i) = 0.0
	612	CALL advec_s_up( i, j, sa )
	613	ENDIF
	614	ENDIF
	615	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	616	wall_salinityflux, tend )
	617	ENDIF
	618
	619	CALL user_actions( i, j, 'sa-tendency' )
	620
	621	!
	622	!-- Prognostic equation for salinity
	623	DO k = nzb_s_inner(j,i)+1, nzt
	624	sa_p(k,j,i) = sat * sa(k,j,i) + &
	625	dt_3d * ( &
	626	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	627	) - &
[785]	628	tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) )
[736]	629	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	630	ENDDO
	631
	632	!
	633	!-- Calculate tendencies for the next Runge-Kutta step
	634	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	635	IF ( intermediate_timestep_count == 1 ) THEN
	636	DO k = nzb_s_inner(j,i)+1, nzt
	637	tsa_m(k,j,i) = tend(k,j,i)
	638	ENDDO
	639	ELSEIF ( intermediate_timestep_count < &
	640	intermediate_timestep_count_max ) THEN
	641	DO k = nzb_s_inner(j,i)+1, nzt
	642	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	643	5.3125 * tsa_m(k,j,i)
	644	ENDDO
	645	ENDIF
	646	ENDIF
	647
	648	!
	649	!-- Calculate density by the equation of state for seawater
	650	CALL eqn_state_seawater( i, j )
	651
	652	ENDDO
	653	ENDDO
	654
	655	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	656
	657	ENDIF
	658
	659	!
	660	!-- If required, compute prognostic equation for total water content / scalar
	661	IF ( humidity .OR. passive_scalar ) THEN
	662
	663	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	664
	665	!
	666	!-- Scalar/q-tendency terms with communication
	667	sat = tsc(1)
	668	sbt = tsc(2)
	669	IF ( scalar_advec == 'bc-scheme' ) THEN
	670
	671	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	672	!
	673	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	674	!-- switched on. Thus:
	675	sat = 1.0
	676	sbt = 1.0
	677	ENDIF
	678	tend = 0.0
	679	CALL advec_s_bc( q, 'q' )
	680	ELSE
	681	IF ( tsc(2) /= 2.0 ) THEN
	682	IF ( scalar_advec == 'ups-scheme' ) THEN
	683	tend = 0.0
	684	CALL advec_s_ups( q, 'q' )
	685	ENDIF
	686	ENDIF
	687	ENDIF
	688
	689	!
	690	!-- Scalar/q-tendency terms with no communication
	691	DO i = nxl, nxr
	692	DO j = nys, nyn
	693	!
	694	!-- Tendency-terms
	695	IF ( scalar_advec == 'bc-scheme' ) THEN
	696	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
	697	wall_qflux, tend )
	698	ELSE
	699	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	700	tend(:,j,i) = 0.0
	701	IF ( ws_scheme_sca ) THEN
	702	CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
	703	diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
	704	ELSE
	705	CALL advec_s_pw( i, j, q )
	706	ENDIF
	707	ELSE
	708	IF ( scalar_advec /= 'ups-scheme' ) THEN
	709	tend(:,j,i) = 0.0
	710	CALL advec_s_up( i, j, q )
	711	ENDIF
	712	ENDIF
	713	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	714	THEN
	715	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
	716	qswst_m, wall_qflux, tend )
	717	ELSE
	718	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
	719	wall_qflux, tend )
	720	ENDIF
	721	ENDIF
	722
	723	!
	724	!-- If required compute decrease of total water content due to
	725	!-- precipitation
	726	IF ( precipitation ) THEN
	727	CALL calc_precipitation( i, j )
	728	ENDIF
	729
	730	!
	731	!-- Sink or source of scalar concentration due to canopy elements
	732	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
	733
	734	!
	735	!-- If required compute influence of large-scale subsidence/ascent
[940]	736	IF ( large_scale_subsidence ) THEN
	737	CALL subsidence( i, j, tend, q, q_init )
[736]	738	ENDIF
	739
	740	CALL user_actions( i, j, 'q-tendency' )
	741
	742	!
	743	!-- Prognostic equation for total water content / scalar
	744	DO k = nzb_s_inner(j,i)+1, nzt
	745	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
	746	dt_3d * ( &
	747	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	748	) - &
[785]	749	tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) )
[736]	750	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	751	ENDDO
	752
	753	!
	754	!-- Calculate tendencies for the next Runge-Kutta step
	755	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	756	IF ( intermediate_timestep_count == 1 ) THEN
	757	DO k = nzb_s_inner(j,i)+1, nzt
	758	tq_m(k,j,i) = tend(k,j,i)
	759	ENDDO
	760	ELSEIF ( intermediate_timestep_count < &
	761	intermediate_timestep_count_max ) THEN
	762	DO k = nzb_s_inner(j,i)+1, nzt
	763	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	764	ENDDO
	765	ENDIF
	766	ENDIF
	767
	768	ENDDO
	769	ENDDO
	770
	771	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	772
	773	ENDIF
	774
	775	!
	776	!-- If required, compute prognostic equation for turbulent kinetic
	777	!-- energy (TKE)
	778	IF ( .NOT. constant_diffusion ) THEN
	779
	780	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	781
	782	!
	783	!-- TKE-tendency terms with communication
	784	CALL production_e_init
	785
	786	sat = tsc(1)
	787	sbt = tsc(2)
	788	IF ( .NOT. use_upstream_for_tke ) THEN
	789	IF ( scalar_advec == 'bc-scheme' ) THEN
	790
	791	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	792	!
	793	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	794	!-- switched on. Thus:
	795	sat = 1.0
	796	sbt = 1.0
	797	ENDIF
	798	tend = 0.0
	799	CALL advec_s_bc( e, 'e' )
	800	ELSE
	801	IF ( tsc(2) /= 2.0 ) THEN
	802	IF ( scalar_advec == 'ups-scheme' ) THEN
	803	tend = 0.0
	804	CALL advec_s_ups( e, 'e' )
	805	ENDIF
	806	ENDIF
	807	ENDIF
	808	ENDIF
	809
	810	!
	811	!-- TKE-tendency terms with no communication
	812	DO i = nxl, nxr
	813	DO j = nys, nyn
	814	!
	815	!-- Tendency-terms
	816	IF ( scalar_advec == 'bc-scheme' .AND. &
	817	.NOT. use_upstream_for_tke ) THEN
	818	IF ( .NOT. humidity ) THEN
	819	IF ( ocean ) THEN
	820	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	821	l_grid, prho, prho_reference, rif, &
	822	tend, zu, zw )
	823	ELSE
	824	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	825	l_grid, pt, pt_reference, rif, tend, &
	826	zu, zw )
	827	ENDIF
	828	ELSE
	829	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	830	l_grid, vpt, pt_reference, rif, tend, zu, &
	831	zw )
	832	ENDIF
	833	ELSE
	834	IF ( use_upstream_for_tke ) THEN
	835	tend(:,j,i) = 0.0
	836	CALL advec_s_up( i, j, e )
	837	ELSE
	838	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	839	THEN
	840	tend(:,j,i) = 0.0
	841	IF ( ws_scheme_sca ) THEN
	842	CALL advec_s_ws( i, j, e, 'e', flux_s_e, &
	843	diss_s_e, flux_l_e, diss_l_e, i_omp_start, tn )
	844	ELSE
	845	CALL advec_s_pw( i, j, e )
	846	ENDIF
	847	ELSE
	848	IF ( scalar_advec /= 'ups-scheme' ) THEN
	849	tend(:,j,i) = 0.0
	850	CALL advec_s_up( i, j, e )
	851	ENDIF
	852	ENDIF
	853	ENDIF
	854	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	855	THEN
	856	IF ( .NOT. humidity ) THEN
	857	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	858	km_m, l_grid, pt_m, pt_reference, &
	859	rif_m, tend, zu, zw )
	860	ELSE
	861	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	862	km_m, l_grid, vpt_m, pt_reference, &
	863	rif_m, tend, zu, zw )
	864	ENDIF
	865	ELSE
	866	IF ( .NOT. humidity ) THEN
	867	IF ( ocean ) THEN
	868	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	869	km, l_grid, prho, prho_reference, &
	870	rif, tend, zu, zw )
	871	ELSE
	872	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	873	km, l_grid, pt, pt_reference, rif, &
	874	tend, zu, zw )
	875	ENDIF
	876	ELSE
	877	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	878	l_grid, vpt, pt_reference, rif, tend, &
	879	zu, zw )
	880	ENDIF
	881	ENDIF
	882	ENDIF
	883	CALL production_e( i, j )
	884
	885	!
	886	!-- Additional sink term for flows through plant canopies
	887	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
	888
	889	CALL user_actions( i, j, 'e-tendency' )
	890
	891	!
	892	!-- Prognostic equation for TKE.
	893	!-- Eliminate negative TKE values, which can occur due to numerical
	894	!-- reasons in the course of the integration. In such cases the old TKE
	895	!-- value is reduced by 90%.
	896	DO k = nzb_s_inner(j,i)+1, nzt
	897	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
	898	dt_3d * ( &
	899	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	900	)
	901	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	902	ENDDO
	903
	904	!
	905	!-- Calculate tendencies for the next Runge-Kutta step
	906	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	907	IF ( intermediate_timestep_count == 1 ) THEN
	908	DO k = nzb_s_inner(j,i)+1, nzt
	909	te_m(k,j,i) = tend(k,j,i)
	910	ENDDO
	911	ELSEIF ( intermediate_timestep_count < &
	912	intermediate_timestep_count_max ) THEN
	913	DO k = nzb_s_inner(j,i)+1, nzt
	914	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	915	ENDDO
	916	ENDIF
	917	ENDIF
	918
	919	ENDDO
	920	ENDDO
	921
	922	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	923
	924	ENDIF
	925
	926
	927	END SUBROUTINE prognostic_equations_noopt
	928
	929
	930	SUBROUTINE prognostic_equations_cache
	931
	932	!------------------------------------------------------------------------------!
	933	! Version with one optimized loop over all equations. It is only allowed to
	934	! be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
	935	!
	936	! Here the calls of most subroutines are embedded in two DO loops over i and j,
	937	! so communication between CPUs is not allowed (does not make sense) within
	938	! these loops.
	939	!
	940	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
	941	!------------------------------------------------------------------------------!
	942
	943	IMPLICIT NONE
	944
	945	CHARACTER (LEN=9) :: time_to_string
	946	INTEGER :: i, i_omp_start, j, k, omp_get_thread_num, tn = 0
	947	LOGICAL :: loop_start
	948
	949
	950	!
	951	!-- Time measurement can only be performed for the whole set of equations
	952	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
	953
	954
	955	!
	956	!-- Calculate those variables needed in the tendency terms which need
	957	!-- global communication
[940]	958	IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 )
	959	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	960	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
[736]	961	IF ( .NOT. constant_diffusion ) CALL production_e_init
	962	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
	963	intermediate_timestep_count == 1 ) CALL ws_statistics
	964
	965	!
	966	!-- Loop over all prognostic equations
	967	!$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
	968
	969	!$ tn = omp_get_thread_num()
	970	loop_start = .TRUE.
	971	!$OMP DO
	972	DO i = nxl, nxr
	973
	974	!
	975	!-- Store the first loop index. It differs for each thread and is required
	976	!-- later in advec_ws
	977	IF ( loop_start ) THEN
	978	loop_start = .FALSE.
	979	i_omp_start = i
	980	ENDIF
	981
	982	DO j = nys, nyn
	983	!
	984	!-- Tendency terms for u-velocity component
	985	IF ( .NOT. outflow_l .OR. i > nxl ) THEN
	986
	987	tend(:,j,i) = 0.0
	988	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	989	IF ( ws_scheme_mom ) THEN
[978]	990	IF ( ( inflow_l .OR. outflow_l ) .AND. i_omp_start == nxl ) THEN
[736]	991	! CALL local_diss( i, j, u) ! dissipation control
	992	CALL advec_u_ws( i, j, i_omp_start + 1, tn )
	993	ELSE
	994	CALL advec_u_ws( i, j, i_omp_start, tn )
	995	ENDIF
	996	ELSE
	997	CALL advec_u_pw( i, j )
	998	ENDIF
	999	ELSE
	1000	CALL advec_u_up( i, j )
	1001	ENDIF
	1002	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	1003	THEN
[978]	1004	CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend, u_m, &
	1005	usws_m, uswst_m, v_m, w_m )
[736]	1006	ELSE
[978]	1007	CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, &
	1008	uswst, v, w )
[736]	1009	ENDIF
	1010	CALL coriolis( i, j, 1 )
[940]	1011	IF ( sloping_surface .AND. .NOT. neutral ) THEN
	1012	CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
	1013	ENDIF
[736]	1014
	1015	!
	1016	!-- Drag by plant canopy
	1017	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 1 )
	1018
	1019	!
	1020	!-- External pressure gradient
	1021	IF ( dp_external ) THEN
	1022	DO k = dp_level_ind_b+1, nzt
	1023	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	1024	ENDDO
	1025	ENDIF
	1026
	1027	CALL user_actions( i, j, 'u-tendency' )
	1028
	1029	!
	1030	!-- Prognostic equation for u-velocity component
	1031	DO k = nzb_u_inner(j,i)+1, nzt
	1032	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
	1033	dt_3d * ( &
	1034	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
	1035	) - &
	1036	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	1037	ENDDO
	1038
	1039	!
	1040	!-- Calculate tendencies for the next Runge-Kutta step
	1041	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1042	IF ( intermediate_timestep_count == 1 ) THEN
	1043	DO k = nzb_u_inner(j,i)+1, nzt
	1044	tu_m(k,j,i) = tend(k,j,i)
	1045	ENDDO
	1046	ELSEIF ( intermediate_timestep_count < &
	1047	intermediate_timestep_count_max ) THEN
	1048	DO k = nzb_u_inner(j,i)+1, nzt
	1049	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	1050	ENDDO
	1051	ENDIF
	1052	ENDIF
	1053
	1054	ENDIF
	1055
	1056	!
	1057	!-- Tendency terms for v-velocity component
	1058	IF ( .NOT. outflow_s .OR. j > nys ) THEN
	1059
	1060	tend(:,j,i) = 0.0
	1061	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1062	IF ( ws_scheme_mom ) THEN
	1063	! CALL local_diss( i, j, v)
	1064	CALL advec_v_ws( i, j, i_omp_start, tn )
	1065	ELSE
	1066	CALL advec_v_pw( i, j )
	1067	ENDIF
	1068	ELSE
	1069	CALL advec_v_up( i, j )
	1070	ENDIF
	1071	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	1072	THEN
[978]	1073	CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
	1074	vsws_m, vswst_m, w_m )
[736]	1075	ELSE
[978]	1076	CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, &
	1077	vswst, w )
[736]	1078	ENDIF
	1079	CALL coriolis( i, j, 2 )
	1080
	1081	!
	1082	!-- Drag by plant canopy
	1083	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 2 )
	1084
	1085	!
	1086	!-- External pressure gradient
	1087	IF ( dp_external ) THEN
	1088	DO k = dp_level_ind_b+1, nzt
	1089	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1090	ENDDO
	1091	ENDIF
	1092
	1093	CALL user_actions( i, j, 'v-tendency' )
	1094
	1095	!
	1096	!-- Prognostic equation for v-velocity component
	1097	DO k = nzb_v_inner(j,i)+1, nzt
	1098	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
	1099	dt_3d * ( &
	1100	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
	1101	) - &
	1102	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	1103	ENDDO
	1104
	1105	!
	1106	!-- Calculate tendencies for the next Runge-Kutta step
	1107	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1108	IF ( intermediate_timestep_count == 1 ) THEN
	1109	DO k = nzb_v_inner(j,i)+1, nzt
	1110	tv_m(k,j,i) = tend(k,j,i)
	1111	ENDDO
	1112	ELSEIF ( intermediate_timestep_count < &
	1113	intermediate_timestep_count_max ) THEN
	1114	DO k = nzb_v_inner(j,i)+1, nzt
	1115	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	1116	ENDDO
	1117	ENDIF
	1118	ENDIF
	1119
	1120	ENDIF
	1121
	1122	!
	1123	!-- Tendency terms for w-velocity component
	1124	tend(:,j,i) = 0.0
	1125	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1126	IF ( ws_scheme_mom ) THEN
	1127	! CALL local_diss( i, j, w)
	1128	CALL advec_w_ws( i, j, i_omp_start, tn )
	1129	ELSE
	1130	CALL advec_w_pw( i, j )
	1131	END IF
	1132	ELSE
	1133	CALL advec_w_up( i, j )
	1134	ENDIF
	1135	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	1136	THEN
[978]	1137	CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
	1138	w_m )
[736]	1139	ELSE
[978]	1140	CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w )
[736]	1141	ENDIF
	1142	CALL coriolis( i, j, 3 )
[940]	1143
	1144	IF ( .NOT. neutral ) THEN
	1145	IF ( ocean ) THEN
	1146	CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
[736]	1147	ELSE
[940]	1148	IF ( .NOT. humidity ) THEN
	1149	CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
	1150	ELSE
	1151	CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
	1152	ENDIF
[736]	1153	ENDIF
	1154	ENDIF
	1155
	1156	!
	1157	!-- Drag by plant canopy
	1158	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 3 )
	1159
	1160	CALL user_actions( i, j, 'w-tendency' )
	1161
	1162	!
	1163	!-- Prognostic equation for w-velocity component
	1164	DO k = nzb_w_inner(j,i)+1, nzt-1
	1165	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
	1166	dt_3d * ( &
	1167	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
	1168	) - &
	1169	tsc(5) * rdf(k) * w(k,j,i)
	1170	ENDDO
	1171
	1172	!
	1173	!-- Calculate tendencies for the next Runge-Kutta step
	1174	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1175	IF ( intermediate_timestep_count == 1 ) THEN
	1176	DO k = nzb_w_inner(j,i)+1, nzt-1
	1177	tw_m(k,j,i) = tend(k,j,i)
	1178	ENDDO
	1179	ELSEIF ( intermediate_timestep_count < &
	1180	intermediate_timestep_count_max ) THEN
	1181	DO k = nzb_w_inner(j,i)+1, nzt-1
	1182	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1183	ENDDO
	1184	ENDIF
	1185	ENDIF
	1186
	1187	!
[940]	1188	!-- If required, compute prognostic equation for potential temperature
	1189	IF ( .NOT. neutral ) THEN
	1190	!
	1191	!-- Tendency terms for potential temperature
	1192	tend(:,j,i) = 0.0
	1193	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1194	IF ( ws_scheme_sca ) THEN
	1195	CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
	1196	flux_l_pt, diss_l_pt, i_omp_start, tn )
	1197	ELSE
	1198	CALL advec_s_pw( i, j, pt )
	1199	ENDIF
	1200	ELSE
	1201	CALL advec_s_up( i, j, pt )
	1202	ENDIF
	1203	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	1204	THEN
	1205	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
	1206	tswst_m, wall_heatflux, tend )
	1207	ELSE
	1208	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
	1209	wall_heatflux, tend )
	1210	ENDIF
[736]	1211
	1212	!
[940]	1213	!-- If required compute heating/cooling due to long wave radiation
	1214	!-- processes
	1215	IF ( radiation ) THEN
	1216	CALL calc_radiation( i, j )
	1217	ENDIF
[736]	1218
	1219	!
[940]	1220	!-- If required compute impact of latent heat due to precipitation
	1221	IF ( precipitation ) THEN
	1222	CALL impact_of_latent_heat( i, j )
	1223	ENDIF
[736]	1224
	1225	!
[940]	1226	!-- Consideration of heat sources within the plant canopy
	1227	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	1228	CALL plant_canopy_model( i, j, 4 )
	1229	ENDIF
[736]	1230
[940]	1231	!
	1232	!-- If required, compute influence of large-scale subsidence/ascent
	1233	IF ( large_scale_subsidence ) THEN
	1234	CALL subsidence( i, j, tend, pt, pt_init )
	1235	ENDIF
[736]	1236
	1237
[940]	1238	CALL user_actions( i, j, 'pt-tendency' )
[736]	1239
	1240	!
[940]	1241	!-- Prognostic equation for potential temperature
	1242	DO k = nzb_s_inner(j,i)+1, nzt
	1243	pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
	1244	dt_3d * ( &
	1245	tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
[978]	1246	) - &
	1247	tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
	1248	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
[940]	1249	ENDDO
[736]	1250
	1251	!
[940]	1252	!-- Calculate tendencies for the next Runge-Kutta step
	1253	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1254	IF ( intermediate_timestep_count == 1 ) THEN
	1255	DO k = nzb_s_inner(j,i)+1, nzt
	1256	tpt_m(k,j,i) = tend(k,j,i)
	1257	ENDDO
	1258	ELSEIF ( intermediate_timestep_count < &
	1259	intermediate_timestep_count_max ) THEN
	1260	DO k = nzb_s_inner(j,i)+1, nzt
	1261	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1262	5.3125 * tpt_m(k,j,i)
	1263	ENDDO
	1264	ENDIF
[736]	1265	ENDIF
[940]	1266
[736]	1267	ENDIF
	1268
	1269	!
	1270	!-- If required, compute prognostic equation for salinity
	1271	IF ( ocean ) THEN
	1272
	1273	!
	1274	!-- Tendency-terms for salinity
	1275	tend(:,j,i) = 0.0
	1276	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	1277	THEN
	1278	IF ( ws_scheme_sca ) THEN
	1279	! CALL local_diss( i, j, sa )
	1280	CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, &
	1281	diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
	1282	ELSE
	1283	CALL advec_s_pw( i, j, sa )
	1284	ENDIF
	1285	ELSE
	1286	CALL advec_s_up( i, j, sa )
	1287	ENDIF
	1288	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	1289	wall_salinityflux, tend )
	1290
	1291	CALL user_actions( i, j, 'sa-tendency' )
	1292
	1293	!
	1294	!-- Prognostic equation for salinity
	1295	DO k = nzb_s_inner(j,i)+1, nzt
	1296	sa_p(k,j,i) = tsc(1) * sa(k,j,i) + &
	1297	dt_3d * ( &
	1298	tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	1299	) - &
[785]	1300	tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) )
[736]	1301	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1302	ENDDO
	1303
	1304	!
	1305	!-- Calculate tendencies for the next Runge-Kutta step
	1306	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1307	IF ( intermediate_timestep_count == 1 ) THEN
	1308	DO k = nzb_s_inner(j,i)+1, nzt
	1309	tsa_m(k,j,i) = tend(k,j,i)
	1310	ENDDO
	1311	ELSEIF ( intermediate_timestep_count < &
	1312	intermediate_timestep_count_max ) THEN
	1313	DO k = nzb_s_inner(j,i)+1, nzt
	1314	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1315	5.3125 * tsa_m(k,j,i)
	1316	ENDDO
	1317	ENDIF
	1318	ENDIF
	1319
	1320	!
	1321	!-- Calculate density by the equation of state for seawater
	1322	CALL eqn_state_seawater( i, j )
	1323
	1324	ENDIF
	1325
	1326	!
	1327	!-- If required, compute prognostic equation for total water content /
	1328	!-- scalar
	1329	IF ( humidity .OR. passive_scalar ) THEN
	1330
	1331	!
	1332	!-- Tendency-terms for total water content / scalar
	1333	tend(:,j,i) = 0.0
	1334	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	1335	THEN
	1336	IF ( ws_scheme_sca ) THEN
	1337	! CALL local_diss( i, j, q )
	1338	CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
	1339	diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
	1340	ELSE
	1341	CALL advec_s_pw( i, j, q )
	1342	ENDIF
	1343	ELSE
	1344	CALL advec_s_up( i, j, q )
	1345	ENDIF
	1346	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	1347	THEN
	1348	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
	1349	qswst_m, wall_qflux, tend )
	1350	ELSE
	1351	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
	1352	wall_qflux, tend )
	1353	ENDIF
	1354
	1355	!
	1356	!-- If required compute decrease of total water content due to
	1357	!-- precipitation
	1358	IF ( precipitation ) THEN
	1359	CALL calc_precipitation( i, j )
	1360	ENDIF
	1361
	1362	!
	1363	!-- Sink or source of scalar concentration due to canopy elements
	1364	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
	1365
	1366	!-- If required compute influence of large-scale subsidence/ascent
[940]	1367	IF ( large_scale_subsidence ) THEN
	1368	CALL subsidence( i, j, tend, q, q_init )
[736]	1369	ENDIF
	1370
	1371	CALL user_actions( i, j, 'q-tendency' )
	1372
	1373	!
	1374	!-- Prognostic equation for total water content / scalar
	1375	DO k = nzb_s_inner(j,i)+1, nzt
	1376	q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
	1377	dt_3d * ( &
	1378	tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	1379	) - &
[785]	1380	tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) )
[736]	1381	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	1382	ENDDO
	1383
	1384	!
	1385	!-- Calculate tendencies for the next Runge-Kutta step
	1386	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1387	IF ( intermediate_timestep_count == 1 ) THEN
	1388	DO k = nzb_s_inner(j,i)+1, nzt
	1389	tq_m(k,j,i) = tend(k,j,i)
	1390	ENDDO
	1391	ELSEIF ( intermediate_timestep_count < &
	1392	intermediate_timestep_count_max ) THEN
	1393	DO k = nzb_s_inner(j,i)+1, nzt
	1394	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1395	5.3125 * tq_m(k,j,i)
	1396	ENDDO
	1397	ENDIF
	1398	ENDIF
	1399
	1400	ENDIF
	1401
	1402	!
	1403	!-- If required, compute prognostic equation for turbulent kinetic
	1404	!-- energy (TKE)
	1405	IF ( .NOT. constant_diffusion ) THEN
	1406
	1407	!
	1408	!-- Tendency-terms for TKE
	1409	tend(:,j,i) = 0.0
	1410	IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	1411	.AND. .NOT. use_upstream_for_tke ) THEN
	1412	IF ( ws_scheme_sca ) THEN
	1413	! CALL local_diss( i, j, e )
	1414	CALL advec_s_ws( i, j, e, 'e', flux_s_e, &
	1415	diss_s_e, flux_l_e, diss_l_e , i_omp_start, tn )
	1416	ELSE
	1417	CALL advec_s_pw( i, j, e )
	1418	ENDIF
	1419	ELSE
	1420	CALL advec_s_up( i, j, e )
	1421	ENDIF
	1422	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	1423	THEN
	1424	IF ( .NOT. humidity ) THEN
	1425	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	1426	km_m, l_grid, pt_m, pt_reference, &
	1427	rif_m, tend, zu, zw )
	1428	ELSE
	1429	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	1430	km_m, l_grid, vpt_m, pt_reference, &
	1431	rif_m, tend, zu, zw )
	1432	ENDIF
	1433	ELSE
	1434	IF ( .NOT. humidity ) THEN
	1435	IF ( ocean ) THEN
	1436	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	1437	km, l_grid, prho, prho_reference, &
	1438	rif, tend, zu, zw )
	1439	ELSE
	1440	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	1441	km, l_grid, pt, pt_reference, rif, &
	1442	tend, zu, zw )
	1443	ENDIF
	1444	ELSE
	1445	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	1446	l_grid, vpt, pt_reference, rif, tend, &
	1447	zu, zw )
	1448	ENDIF
	1449	ENDIF
	1450	CALL production_e( i, j )
	1451
	1452	!
	1453	!-- Additional sink term for flows through plant canopies
	1454	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
	1455
	1456	CALL user_actions( i, j, 'e-tendency' )
	1457
	1458	!
	1459	!-- Prognostic equation for TKE.
	1460	!-- Eliminate negative TKE values, which can occur due to numerical
	1461	!-- reasons in the course of the integration. In such cases the old
	1462	!-- TKE value is reduced by 90%.
	1463	DO k = nzb_s_inner(j,i)+1, nzt
	1464	e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
	1465	dt_3d * ( &
	1466	tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	1467	)
	1468	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	1469	ENDDO
	1470
	1471	!
	1472	!-- Calculate tendencies for the next Runge-Kutta step
	1473	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1474	IF ( intermediate_timestep_count == 1 ) THEN
	1475	DO k = nzb_s_inner(j,i)+1, nzt
	1476	te_m(k,j,i) = tend(k,j,i)
	1477	ENDDO
	1478	ELSEIF ( intermediate_timestep_count < &
	1479	intermediate_timestep_count_max ) THEN
	1480	DO k = nzb_s_inner(j,i)+1, nzt
	1481	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1482	5.3125 * te_m(k,j,i)
	1483	ENDDO
	1484	ENDIF
	1485	ENDIF
	1486
	1487	ENDIF ! TKE equation
	1488
	1489	ENDDO
	1490	ENDDO
	1491	!$OMP END PARALLEL
	1492
	1493	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
	1494
	1495
	1496	END SUBROUTINE prognostic_equations_cache
	1497
	1498
	1499	SUBROUTINE prognostic_equations_vector
	1500
	1501	!------------------------------------------------------------------------------!
	1502	! Version for vector machines
	1503	!------------------------------------------------------------------------------!
	1504
	1505	IMPLICIT NONE
	1506
	1507	CHARACTER (LEN=9) :: time_to_string
	1508	INTEGER :: i, j, k
	1509	REAL :: sat, sbt
	1510
	1511	!
	1512	!-- Calculate those variables needed in the tendency terms which need
	1513	!-- global communication
[940]	1514	IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 )
	1515	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	1516	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
[736]	1517	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
	1518	intermediate_timestep_count == 1 ) CALL ws_statistics
	1519
	1520	!
	1521	!-- u-velocity component
	1522	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	1523
	1524	!
	1525	!-- u-tendency terms with communication
	1526	IF ( momentum_advec == 'ups-scheme' ) THEN
	1527	tend = 0.0
	1528	CALL advec_u_ups
	1529	ENDIF
	1530
	1531	!
	1532	!-- u-tendency terms with no communication
	1533	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1534	tend = 0.0
	1535	IF ( ws_scheme_mom ) THEN
	1536	CALL advec_u_ws
	1537	ELSE
	1538	CALL advec_u_pw
	1539	ENDIF
	1540	ELSE
	1541	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1542	tend = 0.0
	1543	CALL advec_u_up
	1544	ENDIF
	1545	ENDIF
	1546	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[978]	1547	CALL diffusion_u( ddzu, ddzw, km_m, tend, u_m, usws_m, uswst_m, &
	1548	v_m, w_m )
[736]	1549	ELSE
[978]	1550	CALL diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w )
[736]	1551	ENDIF
	1552	CALL coriolis( 1 )
[940]	1553	IF ( sloping_surface .AND. .NOT. neutral ) THEN
	1554	CALL buoyancy( pt, pt_reference, 1, 4 )
	1555	ENDIF
[736]	1556
	1557	!
	1558	!-- Drag by plant canopy
	1559	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
	1560
	1561	!
	1562	!-- External pressure gradient
	1563	IF ( dp_external ) THEN
	1564	DO i = nxlu, nxr
	1565	DO j = nys, nyn
	1566	DO k = dp_level_ind_b+1, nzt
	1567	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	1568	ENDDO
	1569	ENDDO
	1570	ENDDO
	1571	ENDIF
	1572
	1573	CALL user_actions( 'u-tendency' )
	1574
	1575	!
	1576	!-- Prognostic equation for u-velocity component
	1577	DO i = nxlu, nxr
	1578	DO j = nys, nyn
	1579	DO k = nzb_u_inner(j,i)+1, nzt
	1580	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
	1581	dt_3d * ( &
	1582	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
	1583	) - &
	1584	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	1585	ENDDO
	1586	ENDDO
	1587	ENDDO
	1588
	1589	!
	1590	!-- Calculate tendencies for the next Runge-Kutta step
	1591	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1592	IF ( intermediate_timestep_count == 1 ) THEN
	1593	DO i = nxlu, nxr
	1594	DO j = nys, nyn
	1595	DO k = nzb_u_inner(j,i)+1, nzt
	1596	tu_m(k,j,i) = tend(k,j,i)
	1597	ENDDO
	1598	ENDDO
	1599	ENDDO
	1600	ELSEIF ( intermediate_timestep_count < &
	1601	intermediate_timestep_count_max ) THEN
	1602	DO i = nxlu, nxr
	1603	DO j = nys, nyn
	1604	DO k = nzb_u_inner(j,i)+1, nzt
	1605	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	1606	ENDDO
	1607	ENDDO
	1608	ENDDO
	1609	ENDIF
	1610	ENDIF
	1611
	1612	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	1613
	1614	!
	1615	!-- v-velocity component
	1616	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	1617
	1618	!
	1619	!-- v-tendency terms with communication
	1620	IF ( momentum_advec == 'ups-scheme' ) THEN
	1621	tend = 0.0
	1622	CALL advec_v_ups
	1623	ENDIF
	1624
	1625	!
	1626	!-- v-tendency terms with no communication
	1627	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1628	tend = 0.0
	1629	IF ( ws_scheme_mom ) THEN
	1630	CALL advec_v_ws
	1631	ELSE
	1632	CALL advec_v_pw
	1633	END IF
	1634	ELSE
	1635	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1636	tend = 0.0
	1637	CALL advec_v_up
	1638	ENDIF
	1639	ENDIF
	1640	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[978]	1641	CALL diffusion_v( ddzu, ddzw, km_m, tend, u_m, v_m, vsws_m, vswst_m, &
	1642	w_m )
[736]	1643	ELSE
[978]	1644	CALL diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w )
[736]	1645	ENDIF
	1646	CALL coriolis( 2 )
	1647
	1648	!
	1649	!-- Drag by plant canopy
	1650	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
	1651
	1652	!
	1653	!-- External pressure gradient
	1654	IF ( dp_external ) THEN
	1655	DO i = nxl, nxr
	1656	DO j = nysv, nyn
	1657	DO k = dp_level_ind_b+1, nzt
	1658	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1659	ENDDO
	1660	ENDDO
	1661	ENDDO
	1662	ENDIF
	1663
	1664	CALL user_actions( 'v-tendency' )
	1665
	1666	!
	1667	!-- Prognostic equation for v-velocity component
	1668	DO i = nxl, nxr
	1669	DO j = nysv, nyn
	1670	DO k = nzb_v_inner(j,i)+1, nzt
	1671	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
	1672	dt_3d * ( &
	1673	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
	1674	) - &
	1675	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	1676	ENDDO
	1677	ENDDO
	1678	ENDDO
	1679
	1680	!
	1681	!-- Calculate tendencies for the next Runge-Kutta step
	1682	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1683	IF ( intermediate_timestep_count == 1 ) THEN
	1684	DO i = nxl, nxr
	1685	DO j = nysv, nyn
	1686	DO k = nzb_v_inner(j,i)+1, nzt
	1687	tv_m(k,j,i) = tend(k,j,i)
	1688	ENDDO
	1689	ENDDO
	1690	ENDDO
	1691	ELSEIF ( intermediate_timestep_count < &
	1692	intermediate_timestep_count_max ) THEN
	1693	DO i = nxl, nxr
	1694	DO j = nysv, nyn
	1695	DO k = nzb_v_inner(j,i)+1, nzt
	1696	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	1697	ENDDO
	1698	ENDDO
	1699	ENDDO
	1700	ENDIF
	1701	ENDIF
	1702
	1703	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1704
	1705	!
	1706	!-- w-velocity component
	1707	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1708
	1709	!
	1710	!-- w-tendency terms with communication
	1711	IF ( momentum_advec == 'ups-scheme' ) THEN
	1712	tend = 0.0
	1713	CALL advec_w_ups
	1714	ENDIF
	1715
	1716	!
	1717	!-- w-tendency terms with no communication
	1718	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1719	tend = 0.0
	1720	IF ( ws_scheme_mom ) THEN
	1721	CALL advec_w_ws
	1722	ELSE
	1723	CALL advec_w_pw
	1724	ENDIF
	1725	ELSE
	1726	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1727	tend = 0.0
	1728	CALL advec_w_up
	1729	ENDIF
	1730	ENDIF
	1731	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[978]	1732	CALL diffusion_w( ddzu, ddzw, km_m, tend, u_m, v_m, w_m )
[736]	1733	ELSE
[978]	1734	CALL diffusion_w( ddzu, ddzw, km, tend, u, v, w )
[736]	1735	ENDIF
	1736	CALL coriolis( 3 )
[940]	1737
	1738	IF ( .NOT. neutral ) THEN
	1739	IF ( ocean ) THEN
	1740	CALL buoyancy( rho, rho_reference, 3, 64 )
[736]	1741	ELSE
[940]	1742	IF ( .NOT. humidity ) THEN
	1743	CALL buoyancy( pt, pt_reference, 3, 4 )
	1744	ELSE
	1745	CALL buoyancy( vpt, pt_reference, 3, 44 )
	1746	ENDIF
[736]	1747	ENDIF
	1748	ENDIF
	1749
	1750	!
	1751	!-- Drag by plant canopy
	1752	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
	1753
	1754	CALL user_actions( 'w-tendency' )
	1755
	1756	!
	1757	!-- Prognostic equation for w-velocity component
	1758	DO i = nxl, nxr
	1759	DO j = nys, nyn
	1760	DO k = nzb_w_inner(j,i)+1, nzt-1
	1761	w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
	1762	dt_3d * ( &
	1763	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
	1764	) - &
	1765	tsc(5) * rdf(k) * w(k,j,i)
	1766	ENDDO
	1767	ENDDO
	1768	ENDDO
	1769
	1770	!
	1771	!-- Calculate tendencies for the next Runge-Kutta step
	1772	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1773	IF ( intermediate_timestep_count == 1 ) THEN
	1774	DO i = nxl, nxr
	1775	DO j = nys, nyn
	1776	DO k = nzb_w_inner(j,i)+1, nzt-1
	1777	tw_m(k,j,i) = tend(k,j,i)
	1778	ENDDO
	1779	ENDDO
	1780	ENDDO
	1781	ELSEIF ( intermediate_timestep_count < &
	1782	intermediate_timestep_count_max ) THEN
	1783	DO i = nxl, nxr
	1784	DO j = nys, nyn
	1785	DO k = nzb_w_inner(j,i)+1, nzt-1
	1786	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1787	ENDDO
	1788	ENDDO
	1789	ENDDO
	1790	ENDIF
	1791	ENDIF
	1792
	1793	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1794
[940]	1795
[736]	1796	!
[940]	1797	!-- If required, compute prognostic equation for potential temperature
	1798	IF ( .NOT. neutral ) THEN
[736]	1799
[940]	1800	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1801
[736]	1802	!
[940]	1803	!-- pt-tendency terms with communication
	1804	sat = tsc(1)
	1805	sbt = tsc(2)
	1806	IF ( scalar_advec == 'bc-scheme' ) THEN
[736]	1807
[940]	1808	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[736]	1809	!
[940]	1810	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1811	!-- switched on. Thus:
	1812	sat = 1.0
	1813	sbt = 1.0
	1814	ENDIF
[736]	1815	tend = 0.0
[940]	1816	CALL advec_s_bc( pt, 'pt' )
[736]	1817	ELSE
[940]	1818	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
[736]	1819	tend = 0.0
[940]	1820	CALL advec_s_ups( pt, 'pt' )
[736]	1821	ENDIF
	1822	ENDIF
[940]	1823
	1824	!
	1825	!-- pt-tendency terms with no communication
	1826	IF ( scalar_advec == 'bc-scheme' ) THEN
[736]	1827	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
	1828	tend )
[940]	1829	ELSE
	1830	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1831	tend = 0.0
	1832	IF ( ws_scheme_sca ) THEN
	1833	CALL advec_s_ws( pt, 'pt' )
	1834	ELSE
	1835	CALL advec_s_pw( pt )
	1836	ENDIF
	1837	ELSE
	1838	IF ( scalar_advec /= 'ups-scheme' ) THEN
	1839	tend = 0.0
	1840	CALL advec_s_up( pt )
	1841	ENDIF
	1842	ENDIF
	1843	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	1844	CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, &
	1845	wall_heatflux, tend )
	1846	ELSE
	1847	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
	1848	tend )
	1849	ENDIF
[736]	1850	ENDIF
	1851
	1852	!
[940]	1853	!-- If required compute heating/cooling due to long wave radiation processes
	1854	IF ( radiation ) THEN
	1855	CALL calc_radiation
	1856	ENDIF
[736]	1857
	1858	!
[940]	1859	!-- If required compute impact of latent heat due to precipitation
	1860	IF ( precipitation ) THEN
	1861	CALL impact_of_latent_heat
	1862	ENDIF
[736]	1863
	1864	!
[940]	1865	!-- Consideration of heat sources within the plant canopy
	1866	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	1867	CALL plant_canopy_model( 4 )
	1868	ENDIF
[736]	1869
[940]	1870	!
	1871	!-- If required compute influence of large-scale subsidence/ascent
	1872	IF ( large_scale_subsidence ) THEN
	1873	CALL subsidence( tend, pt, pt_init )
	1874	ENDIF
[736]	1875
[940]	1876	CALL user_actions( 'pt-tendency' )
[736]	1877
	1878	!
[940]	1879	!-- Prognostic equation for potential temperature
	1880	DO i = nxl, nxr
	1881	DO j = nys, nyn
	1882	DO k = nzb_s_inner(j,i)+1, nzt
[978]	1883	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
	1884	dt_3d * ( &
	1885	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i)&
	1886	) - &
	1887	tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
	1888	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
[940]	1889	ENDDO
[736]	1890	ENDDO
	1891	ENDDO
	1892
	1893	!
[940]	1894	!-- Calculate tendencies for the next Runge-Kutta step
	1895	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1896	IF ( intermediate_timestep_count == 1 ) THEN
	1897	DO i = nxl, nxr
	1898	DO j = nys, nyn
	1899	DO k = nzb_s_inner(j,i)+1, nzt
	1900	tpt_m(k,j,i) = tend(k,j,i)
	1901	ENDDO
[736]	1902	ENDDO
	1903	ENDDO
[940]	1904	ELSEIF ( intermediate_timestep_count < &
	1905	intermediate_timestep_count_max ) THEN
	1906	DO i = nxl, nxr
	1907	DO j = nys, nyn
	1908	DO k = nzb_s_inner(j,i)+1, nzt
	1909	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1910	5.3125 * tpt_m(k,j,i)
	1911	ENDDO
[736]	1912	ENDDO
	1913	ENDDO
[940]	1914	ENDIF
[736]	1915	ENDIF
[940]	1916
	1917	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1918
[736]	1919	ENDIF
	1920
	1921	!
	1922	!-- If required, compute prognostic equation for salinity
	1923	IF ( ocean ) THEN
	1924
	1925	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	1926
	1927	!
	1928	!-- sa-tendency terms with communication
	1929	sat = tsc(1)
	1930	sbt = tsc(2)
	1931	IF ( scalar_advec == 'bc-scheme' ) THEN
	1932
	1933	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1934	!
	1935	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1936	!-- switched on. Thus:
	1937	sat = 1.0
	1938	sbt = 1.0
	1939	ENDIF
	1940	tend = 0.0
	1941	CALL advec_s_bc( sa, 'sa' )
	1942	ELSE
	1943	IF ( tsc(2) /= 2.0 ) THEN
	1944	IF ( scalar_advec == 'ups-scheme' ) THEN
	1945	tend = 0.0
	1946	CALL advec_s_ups( sa, 'sa' )
	1947	ENDIF
	1948	ENDIF
	1949	ENDIF
	1950
	1951	!
	1952	!-- sa-tendency terms with no communication
	1953	IF ( scalar_advec == 'bc-scheme' ) THEN
	1954	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
	1955	wall_salinityflux, tend )
	1956	ELSE
	1957	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1958	tend = 0.0
	1959	IF ( ws_scheme_sca ) THEN
	1960	CALL advec_s_ws( sa, 'sa' )
	1961	ELSE
	1962	CALL advec_s_pw( sa )
	1963	ENDIF
	1964	ELSE
	1965	IF ( scalar_advec /= 'ups-scheme' ) THEN
	1966	tend = 0.0
	1967	CALL advec_s_up( sa )
	1968	ENDIF
	1969	ENDIF
	1970	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
	1971	wall_salinityflux, tend )
	1972	ENDIF
	1973
	1974	CALL user_actions( 'sa-tendency' )
	1975
	1976	!
	1977	!-- Prognostic equation for salinity
	1978	DO i = nxl, nxr
	1979	DO j = nys, nyn
	1980	DO k = nzb_s_inner(j,i)+1, nzt
	1981	sa_p(k,j,i) = sat * sa(k,j,i) + &
	1982	dt_3d * ( &
	1983	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	1984	) - &
[785]	1985	tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) )
[736]	1986	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1987	ENDDO
	1988	ENDDO
	1989	ENDDO
	1990
	1991	!
	1992	!-- Calculate tendencies for the next Runge-Kutta step
	1993	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1994	IF ( intermediate_timestep_count == 1 ) THEN
	1995	DO i = nxl, nxr
	1996	DO j = nys, nyn
	1997	DO k = nzb_s_inner(j,i)+1, nzt
	1998	tsa_m(k,j,i) = tend(k,j,i)
	1999	ENDDO
	2000	ENDDO
	2001	ENDDO
	2002	ELSEIF ( intermediate_timestep_count < &
	2003	intermediate_timestep_count_max ) THEN
	2004	DO i = nxl, nxr
	2005	DO j = nys, nyn
	2006	DO k = nzb_s_inner(j,i)+1, nzt
	2007	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	2008	5.3125 * tsa_m(k,j,i)
	2009	ENDDO
	2010	ENDDO
	2011	ENDDO
	2012	ENDIF
	2013	ENDIF
	2014
	2015	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	2016
	2017	!
	2018	!-- Calculate density by the equation of state for seawater
	2019	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	2020	CALL eqn_state_seawater
	2021	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	2022
	2023	ENDIF
	2024
	2025	!
	2026	!-- If required, compute prognostic equation for total water content / scalar
	2027	IF ( humidity .OR. passive_scalar ) THEN
	2028
	2029	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	2030
	2031	!
	2032	!-- Scalar/q-tendency terms with communication
	2033	sat = tsc(1)
	2034	sbt = tsc(2)
	2035	IF ( scalar_advec == 'bc-scheme' ) THEN
	2036
	2037	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2038	!
	2039	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	2040	!-- switched on. Thus:
	2041	sat = 1.0
	2042	sbt = 1.0
	2043	ENDIF
	2044	tend = 0.0
	2045	CALL advec_s_bc( q, 'q' )
	2046	ELSE
	2047	IF ( tsc(2) /= 2.0 ) THEN
	2048	IF ( scalar_advec == 'ups-scheme' ) THEN
	2049	tend = 0.0
	2050	CALL advec_s_ups( q, 'q' )
	2051	ENDIF
	2052	ENDIF
	2053	ENDIF
	2054
	2055	!
	2056	!-- Scalar/q-tendency terms with no communication
	2057	IF ( scalar_advec == 'bc-scheme' ) THEN
	2058	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, wall_qflux, tend )
	2059	ELSE
	2060	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	2061	tend = 0.0
	2062	IF ( ws_scheme_sca ) THEN
	2063	CALL advec_s_ws( q, 'q' )
	2064	ELSE
	2065	CALL advec_s_pw( q )
	2066	ENDIF
	2067	ELSE
	2068	IF ( scalar_advec /= 'ups-scheme' ) THEN
	2069	tend = 0.0
	2070	CALL advec_s_up( q )
	2071	ENDIF
	2072	ENDIF
	2073	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	2074	CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, &
	2075	wall_qflux, tend )
	2076	ELSE
	2077	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, &
	2078	wall_qflux, tend )
	2079	ENDIF
	2080	ENDIF
	2081
	2082	!
	2083	!-- If required compute decrease of total water content due to
	2084	!-- precipitation
	2085	IF ( precipitation ) THEN
	2086	CALL calc_precipitation
	2087	ENDIF
	2088
	2089	!
	2090	!-- Sink or source of scalar concentration due to canopy elements
	2091	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
	2092
	2093	!
	2094	!-- If required compute influence of large-scale subsidence/ascent
[940]	2095	IF ( large_scale_subsidence ) THEN
	2096	CALL subsidence( tend, q, q_init )
[736]	2097	ENDIF
	2098
	2099	CALL user_actions( 'q-tendency' )
	2100
	2101	!
	2102	!-- Prognostic equation for total water content / scalar
	2103	DO i = nxl, nxr
	2104	DO j = nys, nyn
	2105	DO k = nzb_s_inner(j,i)+1, nzt
	2106	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
	2107	dt_3d * ( &
	2108	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	2109	) - &
[785]	2110	tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) )
[736]	2111	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	2112	ENDDO
	2113	ENDDO
	2114	ENDDO
	2115
	2116	!
	2117	!-- Calculate tendencies for the next Runge-Kutta step
	2118	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2119	IF ( intermediate_timestep_count == 1 ) THEN
	2120	DO i = nxl, nxr
	2121	DO j = nys, nyn
	2122	DO k = nzb_s_inner(j,i)+1, nzt
	2123	tq_m(k,j,i) = tend(k,j,i)
	2124	ENDDO
	2125	ENDDO
	2126	ENDDO
	2127	ELSEIF ( intermediate_timestep_count < &
	2128	intermediate_timestep_count_max ) THEN
	2129	DO i = nxl, nxr
	2130	DO j = nys, nyn
	2131	DO k = nzb_s_inner(j,i)+1, nzt
	2132	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	2133	ENDDO
	2134	ENDDO
	2135	ENDDO
	2136	ENDIF
	2137	ENDIF
	2138
	2139	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	2140
	2141	ENDIF
	2142
	2143	!
	2144	!-- If required, compute prognostic equation for turbulent kinetic
	2145	!-- energy (TKE)
	2146	IF ( .NOT. constant_diffusion ) THEN
	2147
	2148	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	2149
	2150	!
	2151	!-- TKE-tendency terms with communication
	2152	CALL production_e_init
	2153
	2154	sat = tsc(1)
	2155	sbt = tsc(2)
	2156	IF ( .NOT. use_upstream_for_tke ) THEN
	2157	IF ( scalar_advec == 'bc-scheme' ) THEN
	2158
	2159	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2160	!
	2161	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	2162	!-- switched on. Thus:
	2163	sat = 1.0
	2164	sbt = 1.0
	2165	ENDIF
	2166	tend = 0.0
	2167	CALL advec_s_bc( e, 'e' )
	2168	ELSE
	2169	IF ( tsc(2) /= 2.0 ) THEN
	2170	IF ( scalar_advec == 'ups-scheme' ) THEN
	2171	tend = 0.0
	2172	CALL advec_s_ups( e, 'e' )
	2173	ENDIF
	2174	ENDIF
	2175	ENDIF
	2176	ENDIF
	2177
	2178	!
	2179	!-- TKE-tendency terms with no communication
	2180	IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) &
	2181	THEN
	2182	IF ( .NOT. humidity ) THEN
	2183	IF ( ocean ) THEN
	2184	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	2185	prho, prho_reference, rif, tend, zu, zw )
	2186	ELSE
	2187	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
	2188	pt_reference, rif, tend, zu, zw )
	2189	ENDIF
	2190	ELSE
	2191	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
	2192	pt_reference, rif, tend, zu, zw )
	2193	ENDIF
	2194	ELSE
	2195	IF ( use_upstream_for_tke ) THEN
	2196	tend = 0.0
	2197	CALL advec_s_up( e )
	2198	ELSE
	2199	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	2200	tend = 0.0
	2201	IF ( ws_scheme_sca ) THEN
	2202	CALL advec_s_ws( e, 'e' )
	2203	ELSE
	2204	CALL advec_s_pw( e )
	2205	ENDIF
	2206	ELSE
	2207	IF ( scalar_advec /= 'ups-scheme' ) THEN
	2208	tend = 0.0
	2209	CALL advec_s_up( e )
	2210	ENDIF
	2211	ENDIF
	2212	ENDIF
	2213	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	2214	IF ( .NOT. humidity ) THEN
	2215	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
	2216	pt_m, pt_reference, rif_m, tend, zu, zw )
	2217	ELSE
	2218	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
	2219	vpt_m, pt_reference, rif_m, tend, zu, zw )
	2220	ENDIF
	2221	ELSE
	2222	IF ( .NOT. humidity ) THEN
	2223	IF ( ocean ) THEN
	2224	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	2225	prho, prho_reference, rif, tend, zu, zw )
	2226	ELSE
	2227	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	2228	pt, pt_reference, rif, tend, zu, zw )
	2229	ENDIF
	2230	ELSE
	2231	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
	2232	pt_reference, rif, tend, zu, zw )
	2233	ENDIF
	2234	ENDIF
	2235	ENDIF
	2236	CALL production_e
	2237
	2238	!
	2239	!-- Additional sink term for flows through plant canopies
	2240	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
	2241	CALL user_actions( 'e-tendency' )
	2242
	2243	!
	2244	!-- Prognostic equation for TKE.
	2245	!-- Eliminate negative TKE values, which can occur due to numerical
	2246	!-- reasons in the course of the integration. In such cases the old TKE
	2247	!-- value is reduced by 90%.
	2248	DO i = nxl, nxr
	2249	DO j = nys, nyn
	2250	DO k = nzb_s_inner(j,i)+1, nzt
	2251	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
	2252	dt_3d * ( &
	2253	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	2254	)
	2255	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	2256	ENDDO
	2257	ENDDO
	2258	ENDDO
	2259
	2260	!
	2261	!-- Calculate tendencies for the next Runge-Kutta step
	2262	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2263	IF ( intermediate_timestep_count == 1 ) THEN
	2264	DO i = nxl, nxr
	2265	DO j = nys, nyn
	2266	DO k = nzb_s_inner(j,i)+1, nzt
	2267	te_m(k,j,i) = tend(k,j,i)
	2268	ENDDO
	2269	ENDDO
	2270	ENDDO
	2271	ELSEIF ( intermediate_timestep_count < &
	2272	intermediate_timestep_count_max ) THEN
	2273	DO i = nxl, nxr
	2274	DO j = nys, nyn
	2275	DO k = nzb_s_inner(j,i)+1, nzt
	2276	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	2277	ENDDO
	2278	ENDDO
	2279	ENDDO
	2280	ENDIF
	2281	ENDIF
	2282
	2283	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	2284
	2285	ENDIF
	2286
	2287
	2288	END SUBROUTINE prognostic_equations_vector
	2289
	2290
	2291	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |