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

Last change on this file since 95 was 95, checked in by raasch, 14 years ago

further preliminary uncomplete changes for ocean version

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