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

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

Initial repository layout and content

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