source: palm/tags/release-3.4/SOURCE/prognostic_equations.f90 @ 2638

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

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:


Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:


Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

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