source: palm/tags/release-3.3/SOURCE/prognostic_equations.f90 @ 708

Last change on this file since 708 was 98, checked in by raasch, 17 years ago

updating comments and rc-file

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