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

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

preliminary changes concerning update of BC-scheme

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