Home

Context Navigation

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

Last change on this file since 668 was 668, checked in by suehring, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 76.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |