Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |