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

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

preliminary uncomplete changes for ocean version

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