Home

Context Navigation

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

Last change on this file since 193 was 153, checked in by steinfeld, 17 years ago
Update for the plant canopy model
Property svn:keywords set to `Id`
File size: 68.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |