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

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

preliminary changes for radiation conditions

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