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

Last change on this file since 153 was 153, checked in by steinfeld, 14 years ago

Update for the plant canopy model

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