Home

Context Navigation

prognostic_equations.f90 @ 139

Last change on this file since 139 was 139, checked in by raasch, 16 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:

Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:

Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

Property svn:keywords set to Id

File size: 67.7 KB

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

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

Context Navigation

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

Download in other formats: