Home

Context Navigation

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

Last change on this file since 18 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 53.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |