Home

Context Navigation

prognostic_equations.f90 @ 86

Last change on this file since 86 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 54.6 KB

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

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

Context Navigation

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

Download in other formats: