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

Last change on this file since 136 was 129, checked in by letzel, 16 years ago

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

  • it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

Bugfix: assignment of fluxes at walls

Updates to documentation:
chapter_4.2.html#mode_dvrp
chapter_3.5.4.html#time_series

Default for mrun options -q and -n is "sla3" for lctit. Queues bes1 and bes2
removed. DOC/misc/Tsubame.html updated.

Modified default paths (/work/...) for lctit in .mrun.config.default

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