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

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

preliminary update of bugfixes and extensions for non-cyclic BCs

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