Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |