Home

Context Navigation

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

Last change on this file since 29 was 19, checked in by raasch, 18 years ago
preliminary version of modified boundary conditions at top
Property svn:keywords set to `Id`
File size: 53.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |