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

Last change on this file since 57 was 57, checked in by raasch, 15 years ago

preliminary update of further changes, advec_particles is not running!

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