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

Last change on this file since 404 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
Line 
1 MODULE prognostic_equations_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: prognostic_equations.f90 392 2009-09-24 10:39:14Z raasch $
11!
12! 388 2009-09-23 09:40:33Z raasch
13! prho is used instead of rho in diffusion_e,
14! external pressure gradient
15!
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!
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!
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!
29! 97 2007-06-21 08:23:15Z raasch
30! prognostic equation for salinity, density is calculated from equation of
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
36! calc_mean_profile
37!
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!
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!
49! RCS Log replace by Id keyword, revision history cleaned up
50!
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! ------------
62! Solving the prognostic equations.
63!------------------------------------------------------------------------------!
64
65    USE arrays_3d
66    USE control_parameters
67    USE cpulog
68    USE eqn_state_seawater_mod
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
94    USE plant_canopy_model_mod
95    USE production_e_mod
96    USE user_actions_mod
97
98
99    PRIVATE
100    PUBLIC prognostic_equations_noopt, prognostic_equations_cache, &
101           prognostic_equations_vector
102
103    INTERFACE prognostic_equations_noopt
104       MODULE PROCEDURE prognostic_equations_noopt
105    END INTERFACE prognostic_equations_noopt
106
107    INTERFACE prognostic_equations_cache
108       MODULE PROCEDURE prognostic_equations_cache
109    END INTERFACE prognostic_equations_cache
110
111    INTERFACE prognostic_equations_vector
112       MODULE PROCEDURE prognostic_equations_vector
113    END INTERFACE prognostic_equations_vector
114
115
116 CONTAINS
117
118
119 SUBROUTINE prognostic_equations_noopt
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
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 )
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
153    DO  i = nxlu, nxr
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,   &
168                               usws_m, uswst_m, v_m, w_m )
169          ELSE
170             CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
171                               uswst, v, w )
172          ENDIF
173          CALL coriolis( i, j, 1 )
174          IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
175
176!
177!--       Drag by plant canopy
178          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 1 )
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
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
235       DO  j = nysv, nyn
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, &
249                               v_m, vsws_m, vswst_m, w_m )
250          ELSE
251             CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v,  &
252                               vsws, vswst, w )
253          ENDIF
254          CALL coriolis( i, j, 2 )
255
256!
257!--       Drag by plant canopy
258          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 2 )     
259
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
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,  &
329                               tend, u_m, v_m, w_m )
330          ELSE
331             CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y,    &
332                               tend, u, v, w )
333          ENDIF
334          CALL coriolis( i, j, 3 )
335          IF ( ocean )  THEN
336             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
337          ELSE
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
343          ENDIF
344
345!
346!--       Drag by plant canopy
347          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 3 )
348
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
388    sat = tsc(1)
389    sbt = tsc(2)
390    IF ( scalar_advec == 'bc-scheme' )  THEN
391
392       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
393!
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
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
415             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
416                               wall_heatflux, tend )
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
429                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
430                                  tswst_m, wall_heatflux, tend )
431             ELSE
432                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
433                                  wall_heatflux, tend )
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
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
456          CALL user_actions( i, j, 'pt-tendency' )
457
458!
459!--       Prognostic equation for potential temperature
460          DO  k = nzb_s_inner(j,i)+1, nzt
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
472                DO  k = nzb_s_inner(j,i)+1, nzt
473                   tpt_m(k,j,i) = tend(k,j,i)
474                ENDDO
475             ELSEIF ( intermediate_timestep_count < &
476                      intermediate_timestep_count_max )  THEN
477                DO  k = nzb_s_inner(j,i)+1, nzt
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!
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, &
526                                  wall_salinityflux, tend )
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, &
538                                  wall_salinityflux, tend )
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
570!
571!--          Calculate density by the equation of state for seawater
572             CALL eqn_state_seawater( i, j )
573
574          ENDDO
575       ENDDO
576
577       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
578
579    ENDIF
580
581!
582!-- If required, compute prognostic equation for total water content / scalar
583    IF ( humidity  .OR.  passive_scalar )  THEN
584
585       CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
586
587!
588!--    Scalar/q-tendency terms with communication
589       sat = tsc(1)
590       sbt = tsc(2)
591       IF ( scalar_advec == 'bc-scheme' )  THEN
592
593          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
594!
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
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
618                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
619                                  wall_qflux, tend )
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, &
633                                     qswst_m, wall_qflux, tend )
634                ELSE
635                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
636                                     wall_qflux, tend )
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
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
651             CALL user_actions( i, j, 'q-tendency' )
652
653!
654!--          Prognostic equation for total water content / scalar
655             DO  k = nzb_s_inner(j,i)+1, nzt
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) )
661                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
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
668                   DO  k = nzb_s_inner(j,i)+1, nzt
669                      tq_m(k,j,i) = tend(k,j,i)
670                   ENDDO
671                ELSEIF ( intermediate_timestep_count < &
672                         intermediate_timestep_count_max )  THEN
673                   DO  k = nzb_s_inner(j,i)+1, nzt
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
696
697       sat = tsc(1)
698       sbt = tsc(2)
699       IF ( .NOT. use_upstream_for_tke )  THEN
700          IF ( scalar_advec == 'bc-scheme' )  THEN
701
702             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
703!
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
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
729                IF ( .NOT. humidity )  THEN
730                   IF ( ocean )  THEN
731                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
732                                        l_grid, prho, prho_reference, rif,     &
733                                        tend, zu, zw )
734                   ELSE
735                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
736                                        l_grid, pt, pt_reference, rif, tend,   &
737                                        zu, zw )
738                   ENDIF
739                ELSE
740                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,     &
741                                     l_grid, vpt, pt_reference, rif, tend, zu, &
742                                     zw )
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
762                   IF ( .NOT. humidity )  THEN
763                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
764                                        km_m, l_grid, pt_m, pt_reference,   &
765                                        rif_m, tend, zu, zw )
766                   ELSE
767                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
768                                        km_m, l_grid, vpt_m, pt_reference,  &
769                                        rif_m, tend, zu, zw )
770                   ENDIF
771                ELSE
772                   IF ( .NOT. humidity )  THEN
773                      IF ( ocean )  THEN
774                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
775                                           km, l_grid, prho, prho_reference,  &
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
782                   ELSE
783                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
784                                        l_grid, vpt, pt_reference, rif, tend, &
785                                        zu, zw )
786                   ENDIF
787                ENDIF
788             ENDIF
789             CALL production_e( i, j )
790
791!
792!--          Additional sink term for flows through plant canopies
793             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 6 )         
794
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%.
802             DO  k = nzb_s_inner(j,i)+1, nzt
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
814                   DO  k = nzb_s_inner(j,i)+1, nzt
815                      te_m(k,j,i) = tend(k,j,i)
816                   ENDDO
817                ELSEIF ( intermediate_timestep_count < &
818                         intermediate_timestep_count_max )  THEN
819                   DO  k = nzb_s_inner(j,i)+1, nzt
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
833 END SUBROUTINE prognostic_equations_noopt
834
835
836 SUBROUTINE prognostic_equations_cache
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!
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.
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
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 )
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
873    DO  i = nxl, nxr
874       DO  j = nys, nyn
875!
876!--       Tendency terms for u-velocity component
877          IF ( .NOT. outflow_l  .OR.  i > nxl )  THEN
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,  &
888                                  u_m, usws_m, uswst_m, v_m, w_m )
889             ELSE
890                CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
891                                  usws, uswst, v, w )
892             ENDIF
893             CALL coriolis( i, j, 1 )
894             IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, &
895                                                    4 )
896
897!
898!--          Drag by plant canopy
899             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 1 )
900
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
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
941          IF ( .NOT. outflow_s  .OR.  j > nys )  THEN
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,     &
952                                  u_m, v_m, vsws_m, vswst_m, w_m )
953             ELSE
954                CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
955                                  vsws, vswst, w )
956             ENDIF
957             CALL coriolis( i, j, 2 )
958
959!
960!--          Drag by plant canopy
961             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 2 )       
962
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
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
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
1019             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
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
1027
1028!
1029!--       Drag by plant canopy
1030          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 3 )
1031
1032          CALL user_actions( i, j, 'w-tendency' )
1033
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, &
1071                               tswst_m, wall_heatflux, tend )
1072          ELSE
1073             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
1074                               wall_heatflux, tend )
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
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
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
1130             tend(:,j,i) = 0.0
1131             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
1132             THEN
1133                CALL advec_s_pw( i, j, sa )
1134             ELSE
1135                CALL advec_s_up( i, j, sa )
1136             ENDIF
1137             CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
1138                               wall_salinityflux, tend )
1139
1140             CALL user_actions( i, j, 'sa-tendency' )
1141
1142!
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)
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
1157                   DO  k = nzb_s_inner(j,i)+1, nzt
1158                      tsa_m(k,j,i) = tend(k,j,i)
1159                   ENDDO
1160                ELSEIF ( intermediate_timestep_count < &
1161                         intermediate_timestep_count_max )  THEN
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)
1165                   ENDDO
1166                ENDIF
1167             ENDIF
1168
1169!
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
1182             tend(:,j,i) = 0.0
1183             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
1184             THEN
1185                CALL advec_s_pw( i, j, q )
1186             ELSE
1187                CALL advec_s_up( i, j, q )
1188             ENDIF
1189             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
1190             THEN
1191                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
1192                                  qswst_m, wall_qflux, tend )
1193             ELSE
1194                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
1195                                  wall_qflux, tend )
1196             ENDIF
1197       
1198!
1199!--          If required compute decrease of total water content due to
1200!--          precipitation
1201             IF ( precipitation )  THEN
1202                CALL calc_precipitation( i, j )
1203             ENDIF
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
1210             CALL user_actions( i, j, 'q-tendency' )
1211
1212!
1213!--          Prognostic equation for total water content / scalar
1214             DO  k = nzb_s_inner(j,i)+1, nzt
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)
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
1227                   DO  k = nzb_s_inner(j,i)+1, nzt
1228                      tq_m(k,j,i) = tend(k,j,i)
1229                   ENDDO
1230                ELSEIF ( intermediate_timestep_count < &
1231                         intermediate_timestep_count_max )  THEN
1232                   DO  k = nzb_s_inner(j,i)+1, nzt
1233                      tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1234                                     5.3125 * tq_m(k,j,i)
1235                   ENDDO
1236                ENDIF
1237             ENDIF
1238
1239          ENDIF
1240
1241!
1242!--       If required, compute prognostic equation for turbulent kinetic
1243!--       energy (TKE)
1244          IF ( .NOT. constant_diffusion )  THEN
1245
1246!
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 )
1254             ENDIF
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 )
1261                ELSE
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 )
1265                ENDIF
1266             ELSE
1267                IF ( .NOT. humidity )  THEN
1268                   IF ( ocean )  THEN
1269                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
1270                                        km, l_grid, prho, prho_reference,  &
1271                                        rif, tend, zu, zw )
1272                   ELSE
1273                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
1274                                        km, l_grid, pt, pt_reference, rif, &
1275                                        tend, zu, zw )
1276                   ENDIF
1277                ELSE
1278                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
1279                                     l_grid, vpt, pt_reference, rif, tend, &
1280                                     zu, zw )
1281                ENDIF
1282             ENDIF
1283             CALL production_e( i, j )
1284
1285!
1286!--          Additional sink term for flows through plant canopies
1287             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 6 ) 
1288
1289             CALL user_actions( i, j, 'e-tendency' )
1290
1291!
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
1303
1304!
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
1317                ENDIF
1318             ENDIF
1319
1320          ENDIF   ! TKE equation
1321
1322       ENDDO
1323    ENDDO
1324!$OMP END PARALLEL
1325
1326    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1327
1328
1329 END SUBROUTINE prognostic_equations_cache
1330
1331
1332 SUBROUTINE prognostic_equations_vector
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
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 )
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
1374       CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
1375                         uswst_m, v_m, w_m )
1376    ELSE
1377       CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
1378    ENDIF
1379    CALL coriolis( 1 )
1380    IF ( sloping_surface )  CALL buoyancy( pt, pt_reference, 1, 4 )
1381
1382!
1383!-- Drag by plant canopy
1384    IF ( plant_canopy )  CALL plant_canopy_model( 1 )
1385
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
1398    CALL user_actions( 'u-tendency' )
1399
1400!
1401!-- Prognostic equation for u-velocity component
1402    DO  i = nxlu, nxr
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
1419          DO  i = nxlu, nxr
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
1428          DO  i = nxlu, nxr
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, &
1464                         vswst_m, w_m )
1465    ELSE
1466       CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
1467    ENDIF
1468    CALL coriolis( 2 )
1469
1470!
1471!-- Drag by plant canopy
1472    IF ( plant_canopy )  CALL plant_canopy_model( 2 )
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
1486    CALL user_actions( 'v-tendency' )
1487
1488!
1489!-- Prognostic equation for v-velocity component
1490    DO  i = nxl, nxr
1491       DO  j = nysv, nyn
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
1508             DO  j = nysv, nyn
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
1517             DO  j = nysv, nyn
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, &
1552                         v_m, w_m )
1553    ELSE
1554       CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
1555    ENDIF
1556    CALL coriolis( 3 )
1557    IF ( ocean )  THEN
1558       CALL buoyancy( rho, rho_reference, 3, 64 )
1559    ELSE
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
1565    ENDIF
1566
1567!
1568!-- Drag by plant canopy
1569    IF ( plant_canopy )  CALL plant_canopy_model( 3 )
1570
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
1619    sat = tsc(1)
1620    sbt = tsc(2)
1621    IF ( scalar_advec == 'bc-scheme' )  THEN
1622
1623       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1624!
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
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
1642       CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
1643                         tend )
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
1655          CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, &
1656                            wall_heatflux, tend )
1657       ELSE
1658          CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
1659                            tend )
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
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
1683    CALL user_actions( 'pt-tendency' )
1684
1685!
1686!-- Prognostic equation for potential temperature
1687    DO  i = nxl, nxr
1688       DO  j = nys, nyn
1689          DO  k = nzb_s_inner(j,i)+1, nzt
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
1705                DO  k = nzb_s_inner(j,i)+1, nzt
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
1714                DO  k = nzb_s_inner(j,i)+1, nzt
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!
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!
1755!--    sa-tendency terms with no communication
1756       IF ( scalar_advec == 'bc-scheme' )  THEN
1757          CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
1758                            wall_salinityflux, tend )
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
1769          CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
1770                            wall_salinityflux, tend )
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
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
1822    ENDIF
1823
1824!
1825!-- If required, compute prognostic equation for total water content / scalar
1826    IF ( humidity  .OR.  passive_scalar )  THEN
1827
1828       CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1829
1830!
1831!--    Scalar/q-tendency terms with communication
1832       sat = tsc(1)
1833       sbt = tsc(2)
1834       IF ( scalar_advec == 'bc-scheme' )  THEN
1835
1836          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1837!
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
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
1857          CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, wall_qflux, tend )
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
1869             CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, &
1870                               wall_qflux, tend )
1871          ELSE
1872             CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, &
1873                               wall_qflux, tend )
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
1883
1884!
1885!--    Sink or source of scalar concentration due to canopy elements
1886       IF ( plant_canopy ) CALL plant_canopy_model( 5 )
1887
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
1894             DO  k = nzb_s_inner(j,i)+1, nzt
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) )
1900                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
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
1911                   DO  k = nzb_s_inner(j,i)+1, nzt
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
1920                   DO  k = nzb_s_inner(j,i)+1, nzt
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
1942
1943       sat = tsc(1)
1944       sbt = tsc(2)
1945       IF ( .NOT. use_upstream_for_tke )  THEN
1946          IF ( scalar_advec == 'bc-scheme' )  THEN
1947
1948             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1949!
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
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
1971          IF ( .NOT. humidity )  THEN
1972             IF ( ocean )  THEN
1973                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
1974                                  prho, prho_reference, rif, tend, zu, zw )
1975             ELSE
1976                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
1977                                  pt_reference, rif, tend, zu, zw )
1978             ENDIF
1979          ELSE
1980             CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1981                               pt_reference, rif, tend, zu, zw )
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
1999             IF ( .NOT. humidity )  THEN
2000                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
2001                                  pt_m, pt_reference, rif_m, tend, zu, zw )
2002             ELSE
2003                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
2004                                  vpt_m, pt_reference, rif_m, tend, zu, zw )
2005             ENDIF
2006          ELSE
2007             IF ( .NOT. humidity )  THEN
2008                IF ( ocean )  THEN
2009                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
2010                                     prho, prho_reference, rif, tend, zu, zw )
2011                ELSE
2012                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
2013                                     pt, pt_reference, rif, tend, zu, zw )
2014                ENDIF
2015             ELSE
2016                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
2017                                  pt_reference, rif, tend, zu, zw )
2018             ENDIF
2019          ENDIF
2020       ENDIF
2021       CALL production_e
2022
2023!
2024!--    Additional sink term for flows through plant canopies
2025       IF ( plant_canopy )  CALL plant_canopy_model( 6 )
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
2035             DO  k = nzb_s_inner(j,i)+1, nzt
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
2051                   DO  k = nzb_s_inner(j,i)+1, nzt
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
2060                   DO  k = nzb_s_inner(j,i)+1, nzt
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
2073 END SUBROUTINE prognostic_equations_vector
2074
2075
2076 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.