Home

Context Navigation

prognostic_equations.f90 @ 403

Last change on this file since 403 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

Property svn:keywords set to Id

File size: 70.4 KB

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

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

Context Navigation

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

Download in other formats: