Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |