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

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

preliminary version of modified boundary conditions at top

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