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

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

preliminary update for changes concerning non-cyclic boundary conditions

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