Home

Context Navigation

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

Last change on this file since 513 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 71.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |