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

Last change on this file since 138 was 138, checked in by letzel, 14 years ago

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.

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