source: palm/tags/release-3.1c/SOURCE/prognostic_equations.f90 @ 4011

Last change on this file since 4011 was 39, checked in by raasch, 17 years ago

comments prepared for 3.1c

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