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

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

bugs fixed for particle code and bc-scheme

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