Home

Context Navigation

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

Last change on this file since 673 was 673, checked in by suehring, 13 years ago
Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.
Property svn:keywords set to `Id`
File size: 75.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |