Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |