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

Last change on this file since 234 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

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