Home

Context Navigation

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

Last change on this file since 75 was 75, checked in by raasch, 18 years ago
preliminary update for changes concerning non-cyclic boundary conditions
Property svn:keywords set to `Id`
File size: 54.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |