Home

Context Navigation

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

Last change on this file since 589 was 532, checked in by heinze, 15 years ago
NCL scripts allow the output of png files
Property svn:keywords set to `Id`
File size: 71.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |