Home

Context Navigation

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

Last change on this file since 95 was 95, checked in by raasch, 18 years ago
further preliminary uncomplete changes for ocean version
Property svn:keywords set to `Id`
File size: 62.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |