source: palm/tags/release-3.1b/SOURCE/prognostic_equations.f90 @ 4355

Last change on this file since 4355 was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

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