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

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

New:
---

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

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

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

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

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

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

Changed:


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

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

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

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

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

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

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

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

Errors:


Bugfix: summation of sums_l_l in diffusivities.

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

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

Bugfix: model_string needed a default value (combine_plot_fields)

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

Bugfix: assignment of fluxes at walls

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

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

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