Home

Context Navigation

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

Last change on this file since 236 was 198, checked in by raasch, 16 years ago
file headers updated for the next release 3.5
Property svn:keywords set to `Id`
File size: 68.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |