Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |