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

Last change on this file since 785 was 785, checked in by raasch, 10 years ago

scalar quantities can be excluded from Rayleigh damping; bugfix for long lines in configuration file with more than 300 characters

  • Property svn:keywords set to Id
File size: 78.6 KB
Line 
1 MODULE prognostic_equations_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! new factor rdf_sc allows separate Rayleigh damping of scalars
7!
8! Former revisions:
9! -----------------
10! $Id: prognostic_equations.f90 785 2011-11-28 09:47:19Z raasch $
11!
12! 736 2011-08-17 14:13:26Z suehring
13! Bugfix: determination of first thread index i for WS-scheme
14!
15! 709 2011-03-30 09:31:40Z raasch
16! formatting adjustments
17!
18! 673 2011-01-18 16:19:48Z suehring
19! Consideration of the pressure gradient (steered by tsc(4)) during the time
20! integration removed.
21!
22! 667 2010-12-23 12:06:00Z suehring/gryschka
23! Calls of the advection routines with WS5 added.
24! Calls of ws_statistics added to set the statistical arrays to zero after each
25! time step.
26!
27! 531 2010-04-21 06:47:21Z heinze
28! add call of subsidence in the equation for humidity / passive scalar
29!
30! 411 2009-12-11 14:15:58Z heinze
31! add call of subsidence in the equation for potential temperature
32!
33! 388 2009-09-23 09:40:33Z raasch
34! prho is used instead of rho in diffusion_e,
35! external pressure gradient
36!
37! 153 2008-03-19 09:41:30Z steinfeld
38! add call of plant_canopy_model in the prognostic equation for
39! the potential temperature and for the passive scalar
40!
41! 138 2007-11-28 10:03:58Z letzel
42! add call of subroutines that evaluate the canopy drag terms,
43! add wall_*flux to parameter list of calls of diffusion_s
44!
45! 106 2007-08-16 14:30:26Z raasch
46! +uswst, vswst as arguments in calls of diffusion_u|v,
47! loops for u and v are starting from index nxlu, nysv, respectively (needed
48! for non-cyclic boundary conditions)
49!
50! 97 2007-06-21 08:23:15Z raasch
51! prognostic equation for salinity, density is calculated from equation of
52! state for seawater and is used for calculation of buoyancy,
53! +eqn_state_seawater_mod
54! diffusion_e is called with argument rho in case of ocean runs,
55! new argument zw in calls of diffusion_e, new argument pt_/prho_reference
56! in calls of buoyancy and diffusion_e, calc_mean_pt_profile renamed
57! calc_mean_profile
58!
59! 75 2007-03-22 09:54:05Z raasch
60! checking for negative q and limiting for positive values,
61! z0 removed from arguments in calls of diffusion_u/v/w, uxrp, vynp eliminated,
62! subroutine names changed to .._noopt, .._cache, and .._vector,
63! moisture renamed humidity, Bott-Chlond-scheme can be used in the
64! _vector-version
65!
66! 19 2007-02-23 04:53:48Z raasch
67! Calculation of e, q, and pt extended for gridpoint nzt,
68! handling of given temperature/humidity/scalar fluxes at top surface
69!
70! RCS Log replace by Id keyword, revision history cleaned up
71!
72! Revision 1.21  2006/08/04 15:01:07  raasch
73! upstream scheme can be forced to be used for tke (use_upstream_for_tke)
74! regardless of the timestep scheme used for the other quantities,
75! new argument diss in call of diffusion_e
76!
77! Revision 1.1  2000/04/13 14:56:27  schroeter
78! Initial revision
79!
80!
81! Description:
82! ------------
83! Solving the prognostic equations.
84!------------------------------------------------------------------------------!
85
86    USE arrays_3d
87    USE control_parameters
88    USE cpulog
89    USE eqn_state_seawater_mod
90    USE grid_variables
91    USE indices
92    USE interfaces
93    USE pegrid
94    USE pointer_interfaces
95    USE statistics
96    USE advec_ws
97    USE advec_s_pw_mod
98    USE advec_s_up_mod
99    USE advec_u_pw_mod
100    USE advec_u_up_mod
101    USE advec_v_pw_mod
102    USE advec_v_up_mod
103    USE advec_w_pw_mod
104    USE advec_w_up_mod
105    USE buoyancy_mod
106    USE calc_precipitation_mod
107    USE calc_radiation_mod
108    USE coriolis_mod
109    USE diffusion_e_mod
110    USE diffusion_s_mod
111    USE diffusion_u_mod
112    USE diffusion_v_mod
113    USE diffusion_w_mod
114    USE impact_of_latent_heat_mod
115    USE plant_canopy_model_mod
116    USE production_e_mod
117    USE subsidence_mod
118    USE user_actions_mod
119
120
121    PRIVATE
122    PUBLIC prognostic_equations_noopt, prognostic_equations_cache, &
123           prognostic_equations_vector
124
125    INTERFACE prognostic_equations_noopt
126       MODULE PROCEDURE prognostic_equations_noopt
127    END INTERFACE prognostic_equations_noopt
128
129    INTERFACE prognostic_equations_cache
130       MODULE PROCEDURE prognostic_equations_cache
131    END INTERFACE prognostic_equations_cache
132
133    INTERFACE prognostic_equations_vector
134       MODULE PROCEDURE prognostic_equations_vector
135    END INTERFACE prognostic_equations_vector
136
137
138 CONTAINS
139
140
141 SUBROUTINE prognostic_equations_noopt
142
143!------------------------------------------------------------------------------!
144! Version with single loop optimization
145!
146! (Optimized over each single prognostic equation.)
147!------------------------------------------------------------------------------!
148
149    IMPLICIT NONE
150
151    CHARACTER (LEN=9) ::  time_to_string
152    INTEGER ::  i, i_omp_start, j, k, tn = 0
153    REAL    ::  sat, sbt
154
155!
156!-- Calculate those variables needed in the tendency terms which need
157!-- global communication
158    CALL calc_mean_profile( pt, 4 )
159    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
160    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
161    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
162         intermediate_timestep_count == 1 )  CALL ws_statistics
163
164!
165!-- u-velocity component
166    CALL cpu_log( log_point(5), 'u-equation', 'start' )
167
168!
169!-- u-tendency terms with communication
170    IF ( momentum_advec == 'ups-scheme' )  THEN
171       tend = 0.0
172       CALL advec_u_ups
173    ENDIF
174
175!
176!-- u-tendency terms with no communication
177    i_omp_start = nxlu
178    DO  i = nxlu, nxr
179       DO  j = nys, nyn
180!
181!--       Tendency terms
182          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
183             tend(:,j,i) = 0.0
184             IF ( ws_scheme_mom )  THEN
185                 CALL advec_u_ws( i, j, i_omp_start, tn )
186             ELSE
187                 CALL advec_u_pw( i, j )
188             ENDIF
189
190          ELSE
191             IF ( momentum_advec /= 'ups-scheme' )  THEN
192                tend(:,j,i) = 0.0
193                CALL advec_u_up( i, j )
194             ENDIF
195          ENDIF
196          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
197             CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m,   &
198                               usws_m, uswst_m, v_m, w_m )
199          ELSE
200             CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
201                               uswst, v, w )
202          ENDIF
203          CALL coriolis( i, j, 1 )
204          IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
205
206!
207!--       Drag by plant canopy
208          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 1 )
209
210!
211!--       External pressure gradient
212          IF ( dp_external )  THEN
213             DO  k = dp_level_ind_b+1, nzt
214                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
215             ENDDO
216          ENDIF
217
218          CALL user_actions( i, j, 'u-tendency' )
219
220!
221!--       Prognostic equation for u-velocity component
222          DO  k = nzb_u_inner(j,i)+1, nzt
223             u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) +    &
224                          dt_3d * (                                            &
225                                   tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
226                                  ) -                                          &
227                           tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
228          ENDDO
229
230!
231!--       Calculate tendencies for the next Runge-Kutta step
232          IF ( timestep_scheme(1:5) == 'runge' )  THEN
233             IF ( intermediate_timestep_count == 1 )  THEN
234                DO  k = nzb_u_inner(j,i)+1, nzt
235                   tu_m(k,j,i) = tend(k,j,i)
236                ENDDO
237             ELSEIF ( intermediate_timestep_count < &
238                      intermediate_timestep_count_max )  THEN
239                DO  k = nzb_u_inner(j,i)+1, nzt
240                   tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
241                ENDDO
242             ENDIF
243          ENDIF
244
245       ENDDO
246    ENDDO
247
248    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
249
250!
251!-- v-velocity component
252    CALL cpu_log( log_point(6), 'v-equation', 'start' )
253
254!
255!-- v-tendency terms with communication
256    IF ( momentum_advec == 'ups-scheme' )  THEN
257       tend = 0.0
258       CALL advec_v_ups
259    ENDIF
260
261!
262!-- v-tendency terms with no communication
263    i_omp_start = nxl
264    DO  i = nxl, nxr
265       DO  j = nysv, nyn
266!
267!--       Tendency terms
268          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
269             tend(:,j,i) = 0.0
270             IF ( ws_scheme_mom )  THEN
271                 CALL advec_v_ws( i, j, i_omp_start, tn )
272             ELSE
273                 CALL advec_v_pw( i, j )
274             ENDIF
275
276          ELSE
277             IF ( momentum_advec /= 'ups-scheme' )  THEN
278                tend(:,j,i) = 0.0
279                CALL advec_v_up( i, j )
280             ENDIF
281          ENDIF
282          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
283             CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
284                               v_m, vsws_m, vswst_m, w_m )
285          ELSE
286             CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v,  &
287                               vsws, vswst, w )
288          ENDIF
289          CALL coriolis( i, j, 2 )
290
291!
292!--       Drag by plant canopy
293          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 2 )     
294
295!
296!--       External pressure gradient
297          IF ( dp_external )  THEN
298             DO  k = dp_level_ind_b+1, nzt
299                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
300             ENDDO
301          ENDIF
302
303          CALL user_actions( i, j, 'v-tendency' )
304
305!
306!--       Prognostic equation for v-velocity component
307          DO  k = nzb_v_inner(j,i)+1, nzt
308             v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) +    &
309                          dt_3d * (                                            &
310                                   tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
311                                  ) -                                          &
312                          tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
313          ENDDO
314
315!
316!--       Calculate tendencies for the next Runge-Kutta step
317          IF ( timestep_scheme(1:5) == 'runge' )  THEN
318             IF ( intermediate_timestep_count == 1 )  THEN
319                DO  k = nzb_v_inner(j,i)+1, nzt
320                   tv_m(k,j,i) = tend(k,j,i)
321                ENDDO
322             ELSEIF ( intermediate_timestep_count < &
323                      intermediate_timestep_count_max )  THEN
324                DO  k = nzb_v_inner(j,i)+1, nzt
325                   tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
326                ENDDO
327             ENDIF
328          ENDIF
329
330       ENDDO
331    ENDDO
332
333    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
334
335!
336!-- w-velocity component
337    CALL cpu_log( log_point(7), 'w-equation', 'start' )
338
339!
340!-- w-tendency terms with communication
341    IF ( momentum_advec == 'ups-scheme' )  THEN
342       tend = 0.0
343       CALL advec_w_ups
344    ENDIF
345
346!
347!-- w-tendency terms with no communication
348    DO  i = nxl, nxr
349       DO  j = nys, nyn
350!
351!--       Tendency terms
352          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
353             tend(:,j,i) = 0.0
354             IF ( ws_scheme_mom )  THEN
355                 CALL advec_w_ws( i, j, i_omp_start, tn )
356             ELSE 
357                 CALL advec_w_pw( i, j )
358             ENDIF
359
360          ELSE
361             IF ( momentum_advec /= 'ups-scheme' )  THEN
362                tend(:,j,i) = 0.0
363                CALL advec_w_up( i, j )
364             ENDIF
365          ENDIF
366          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
367             CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y,  &
368                               tend, u_m, v_m, w_m )
369          ELSE
370             CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y,    &
371                               tend, u, v, w )
372          ENDIF
373          CALL coriolis( i, j, 3 )
374          IF ( ocean )  THEN
375             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
376          ELSE
377             IF ( .NOT. humidity )  THEN
378                CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
379             ELSE
380                CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
381             ENDIF
382          ENDIF
383
384!
385!--       Drag by plant canopy
386          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 3 )
387
388          CALL user_actions( i, j, 'w-tendency' )
389
390!
391!--       Prognostic equation for w-velocity component
392          DO  k = nzb_w_inner(j,i)+1, nzt-1
393             w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) +    &
394                          dt_3d * (                                            &
395                                   tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
396                                  ) -                                          &
397                          tsc(5) * rdf(k) * w(k,j,i)
398          ENDDO
399
400!
401!--       Calculate tendencies for the next Runge-Kutta step
402          IF ( timestep_scheme(1:5) == 'runge' )  THEN
403             IF ( intermediate_timestep_count == 1 )  THEN
404                DO  k = nzb_w_inner(j,i)+1, nzt-1
405                   tw_m(k,j,i) = tend(k,j,i)
406                ENDDO
407             ELSEIF ( intermediate_timestep_count < &
408                      intermediate_timestep_count_max )  THEN
409                DO  k = nzb_w_inner(j,i)+1, nzt-1
410                   tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
411                ENDDO
412             ENDIF
413          ENDIF
414
415       ENDDO
416    ENDDO
417
418    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
419
420!
421!-- potential temperature
422    CALL cpu_log( log_point(13), 'pt-equation', 'start' )
423
424!
425!-- pt-tendency terms with communication
426    sat = tsc(1)
427    sbt = tsc(2)
428    IF ( scalar_advec == 'bc-scheme' )  THEN
429
430       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
431!
432!--       Bott-Chlond scheme always uses Euler time step when leapfrog is
433!--       switched on. Thus:
434          sat = 1.0
435          sbt = 1.0
436       ENDIF
437       tend = 0.0
438       CALL advec_s_bc( pt, 'pt' )
439    ELSE
440       IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
441          tend = 0.0
442          CALL advec_s_ups( pt, 'pt' )
443       ENDIF
444    ENDIF
445
446!
447!-- pt-tendency terms with no communication
448    DO  i = nxl, nxr
449       DO  j = nys, nyn
450!
451!--       Tendency terms
452          IF ( scalar_advec == 'bc-scheme' )  THEN
453             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
454                               wall_heatflux, tend )
455          ELSE
456             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
457                tend(:,j,i) = 0.0
458                IF ( ws_scheme_sca )  THEN
459                   CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
460                                 diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
461                ELSE
462                    CALL advec_s_pw( i, j, pt )
463                ENDIF
464             ELSE
465                IF ( scalar_advec /= 'ups-scheme' )  THEN
466                   tend(:,j,i) = 0.0
467                   CALL advec_s_up( i, j, pt )
468                ENDIF
469             ENDIF
470             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
471             THEN
472                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
473                                  tswst_m, wall_heatflux, tend )
474             ELSE
475                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
476                                  wall_heatflux, tend )
477             ENDIF
478          ENDIF
479
480!
481!--       If required compute heating/cooling due to long wave radiation
482!--       processes
483          IF ( radiation )  THEN
484             CALL calc_radiation( i, j )
485          ENDIF
486
487!
488!--       If required compute impact of latent heat due to precipitation
489          IF ( precipitation )  THEN
490             CALL impact_of_latent_heat( i, j )
491          ENDIF
492
493!
494!--       Consideration of heat sources within the plant canopy
495          IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
496             CALL plant_canopy_model( i, j, 4 )
497          ENDIF
498
499!
500!--       If required compute influence of large-scale subsidence/ascent
501          IF ( large_scale_subsidence ) THEN
502             CALL subsidence ( i, j, tend, pt, pt_init )
503          ENDIF
504
505          CALL user_actions( i, j, 'pt-tendency' )
506
507!
508!--       Prognostic equation for potential temperature
509          DO  k = nzb_s_inner(j,i)+1, nzt
510             pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
511                           dt_3d * (                                           &
512                                     sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
513                                   ) -                                         &
514                           tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
515          ENDDO
516
517!
518!--       Calculate tendencies for the next Runge-Kutta step
519          IF ( timestep_scheme(1:5) == 'runge' )  THEN
520             IF ( intermediate_timestep_count == 1 )  THEN
521                DO  k = nzb_s_inner(j,i)+1, nzt
522                   tpt_m(k,j,i) = tend(k,j,i)
523                ENDDO
524             ELSEIF ( intermediate_timestep_count < &
525                      intermediate_timestep_count_max )  THEN
526                DO  k = nzb_s_inner(j,i)+1, nzt
527                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
528                ENDDO
529             ENDIF
530          ENDIF
531
532       ENDDO
533    ENDDO
534
535    CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
536
537!
538!-- If required, compute prognostic equation for salinity
539    IF ( ocean )  THEN
540
541       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
542
543!
544!--    sa-tendency terms with communication
545       sat = tsc(1)
546       sbt = tsc(2)
547       IF ( scalar_advec == 'bc-scheme' )  THEN
548
549          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
550!
551!--          Bott-Chlond scheme always uses Euler time step when leapfrog is
552!--          switched on. Thus:
553             sat = 1.0
554             sbt = 1.0
555          ENDIF
556          tend = 0.0
557          CALL advec_s_bc( sa, 'sa' )
558       ELSE
559          IF ( tsc(2) /= 2.0 )  THEN
560             IF ( scalar_advec == 'ups-scheme' )  THEN
561                tend = 0.0
562                CALL advec_s_ups( sa, 'sa' )
563             ENDIF
564          ENDIF
565       ENDIF
566
567!
568!--    sa terms with no communication
569       DO  i = nxl, nxr
570          DO  j = nys, nyn
571!
572!--          Tendency-terms
573             IF ( scalar_advec == 'bc-scheme' )  THEN
574                CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
575                                  wall_salinityflux, tend )
576             ELSE
577                IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) THEN
578                   tend(:,j,i) = 0.0
579                   IF ( ws_scheme_sca )  THEN
580                       CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
581                                    diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
582                   ELSE
583                       CALL advec_s_pw( i, j, sa )
584                   ENDIF
585
586                ELSE
587                   IF ( scalar_advec /= 'ups-scheme' )  THEN
588                      tend(:,j,i) = 0.0
589                      CALL advec_s_up( i, j, sa )
590                   ENDIF
591                ENDIF
592                CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
593                                  wall_salinityflux, tend )
594             ENDIF
595       
596             CALL user_actions( i, j, 'sa-tendency' )
597
598!
599!--          Prognostic equation for salinity
600             DO  k = nzb_s_inner(j,i)+1, nzt
601                sa_p(k,j,i) = sat * sa(k,j,i) +                                &
602                              dt_3d * (                                        &
603                                     sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
604                                      ) -                                      &
605                              tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) )
606                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
607             ENDDO
608
609!
610!--          Calculate tendencies for the next Runge-Kutta step
611             IF ( timestep_scheme(1:5) == 'runge' )  THEN
612                IF ( intermediate_timestep_count == 1 )  THEN
613                   DO  k = nzb_s_inner(j,i)+1, nzt
614                      tsa_m(k,j,i) = tend(k,j,i)
615                   ENDDO
616                ELSEIF ( intermediate_timestep_count < &
617                         intermediate_timestep_count_max )  THEN
618                   DO  k = nzb_s_inner(j,i)+1, nzt
619                      tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
620                                      5.3125 * tsa_m(k,j,i)
621                   ENDDO
622                ENDIF
623             ENDIF
624
625!
626!--          Calculate density by the equation of state for seawater
627             CALL eqn_state_seawater( i, j )
628
629          ENDDO
630       ENDDO
631
632       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
633
634    ENDIF
635
636!
637!-- If required, compute prognostic equation for total water content / scalar
638    IF ( humidity  .OR.  passive_scalar )  THEN
639
640       CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
641
642!
643!--    Scalar/q-tendency terms with communication
644       sat = tsc(1)
645       sbt = tsc(2)
646       IF ( scalar_advec == 'bc-scheme' )  THEN
647
648          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
649!
650!--          Bott-Chlond scheme always uses Euler time step when leapfrog is
651!--          switched on. Thus:
652             sat = 1.0
653             sbt = 1.0
654          ENDIF
655          tend = 0.0
656          CALL advec_s_bc( q, 'q' )
657       ELSE
658          IF ( tsc(2) /= 2.0 )  THEN
659             IF ( scalar_advec == 'ups-scheme' )  THEN
660                tend = 0.0
661                CALL advec_s_ups( q, 'q' )
662             ENDIF
663          ENDIF
664       ENDIF
665
666!
667!--    Scalar/q-tendency terms with no communication
668       DO  i = nxl, nxr
669          DO  j = nys, nyn
670!
671!--          Tendency-terms
672             IF ( scalar_advec == 'bc-scheme' )  THEN
673                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
674                                  wall_qflux, tend )
675             ELSE
676                IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) THEN
677                   tend(:,j,i) = 0.0
678                   IF ( ws_scheme_sca )  THEN
679                       CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
680                                   diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
681                   ELSE
682                       CALL advec_s_pw( i, j, q )
683                   ENDIF
684                ELSE
685                   IF ( scalar_advec /= 'ups-scheme' )  THEN
686                      tend(:,j,i) = 0.0
687                      CALL advec_s_up( i, j, q )
688                   ENDIF
689                ENDIF
690                IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
691                THEN
692                   CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
693                                     qswst_m, wall_qflux, tend )
694                ELSE
695                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
696                                     wall_qflux, tend )
697                ENDIF
698             ENDIF
699       
700!
701!--          If required compute decrease of total water content due to
702!--          precipitation
703             IF ( precipitation )  THEN
704                CALL calc_precipitation( i, j )
705             ENDIF
706
707!
708!--          Sink or source of scalar concentration due to canopy elements
709             IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
710             
711!
712!--          If required compute influence of large-scale subsidence/ascent
713             IF ( large_scale_subsidence ) THEN
714                CALL subsidence ( i, j, tend, q, q_init )
715             ENDIF
716
717             CALL user_actions( i, j, 'q-tendency' )
718
719!
720!--          Prognostic equation for total water content / scalar
721             DO  k = nzb_s_inner(j,i)+1, nzt
722                q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) +       &
723                             dt_3d * (                                         &
724                                      sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
725                                     ) -                                       &
726                             tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) )
727                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
728             ENDDO
729
730!
731!--          Calculate tendencies for the next Runge-Kutta step
732             IF ( timestep_scheme(1:5) == 'runge' )  THEN
733                IF ( intermediate_timestep_count == 1 )  THEN
734                   DO  k = nzb_s_inner(j,i)+1, nzt
735                      tq_m(k,j,i) = tend(k,j,i)
736                   ENDDO
737                ELSEIF ( intermediate_timestep_count < &
738                         intermediate_timestep_count_max )  THEN
739                   DO  k = nzb_s_inner(j,i)+1, nzt
740                      tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
741                   ENDDO
742                ENDIF
743             ENDIF
744
745          ENDDO
746       ENDDO
747
748       CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
749
750    ENDIF
751
752!
753!-- If required, compute prognostic equation for turbulent kinetic
754!-- energy (TKE)
755    IF ( .NOT. constant_diffusion )  THEN
756
757       CALL cpu_log( log_point(16), 'tke-equation', 'start' )
758
759!
760!--    TKE-tendency terms with communication
761       CALL production_e_init
762
763       sat = tsc(1)
764       sbt = tsc(2)
765       IF ( .NOT. use_upstream_for_tke )  THEN
766          IF ( scalar_advec == 'bc-scheme' )  THEN
767
768             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
769!
770!--             Bott-Chlond scheme always uses Euler time step when leapfrog is
771!--             switched on. Thus:
772                sat = 1.0
773                sbt = 1.0
774             ENDIF
775             tend = 0.0
776             CALL advec_s_bc( e, 'e' )
777          ELSE
778             IF ( tsc(2) /= 2.0 )  THEN
779                IF ( scalar_advec == 'ups-scheme' )  THEN
780                   tend = 0.0
781                   CALL advec_s_ups( e, 'e' )
782                ENDIF
783             ENDIF
784          ENDIF
785       ENDIF
786
787!
788!--    TKE-tendency terms with no communication
789       DO  i = nxl, nxr
790          DO  j = nys, nyn
791!
792!--          Tendency-terms
793             IF ( scalar_advec == 'bc-scheme'  .AND.  &
794                  .NOT. use_upstream_for_tke )  THEN
795                IF ( .NOT. humidity )  THEN
796                   IF ( ocean )  THEN
797                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
798                                        l_grid, prho, prho_reference, rif,     &
799                                        tend, zu, zw )
800                   ELSE
801                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
802                                        l_grid, pt, pt_reference, rif, tend,   &
803                                        zu, zw )
804                   ENDIF
805                ELSE
806                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,     &
807                                     l_grid, vpt, pt_reference, rif, tend, zu, &
808                                     zw )
809                ENDIF
810             ELSE
811                IF ( use_upstream_for_tke )  THEN
812                   tend(:,j,i) = 0.0
813                   CALL advec_s_up( i, j, e )
814                ELSE
815                   IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
816                   THEN
817                      tend(:,j,i) = 0.0
818                      IF ( ws_scheme_sca )  THEN
819                          CALL advec_s_ws( i, j, e, 'e', flux_s_e, & 
820                                      diss_s_e, flux_l_e, diss_l_e, i_omp_start, tn )
821                      ELSE
822                          CALL advec_s_pw( i, j, e )
823                      ENDIF
824                   ELSE
825                      IF ( scalar_advec /= 'ups-scheme' )  THEN
826                         tend(:,j,i) = 0.0
827                         CALL advec_s_up( i, j, e )
828                      ENDIF
829                   ENDIF
830                ENDIF
831                IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
832                THEN
833                   IF ( .NOT. humidity )  THEN
834                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
835                                        km_m, l_grid, pt_m, pt_reference,   &
836                                        rif_m, tend, zu, zw )
837                   ELSE
838                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
839                                        km_m, l_grid, vpt_m, pt_reference,  &
840                                        rif_m, tend, zu, zw )
841                   ENDIF
842                ELSE
843                   IF ( .NOT. humidity )  THEN
844                      IF ( ocean )  THEN
845                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
846                                           km, l_grid, prho, prho_reference,  &
847                                           rif, tend, zu, zw )
848                      ELSE
849                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
850                                           km, l_grid, pt, pt_reference, rif, &
851                                           tend, zu, zw )
852                      ENDIF
853                   ELSE
854                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
855                                        l_grid, vpt, pt_reference, rif, tend, &
856                                        zu, zw )
857                   ENDIF
858                ENDIF
859             ENDIF
860             CALL production_e( i, j )
861
862!
863!--          Additional sink term for flows through plant canopies
864             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 6 )         
865
866             CALL user_actions( i, j, 'e-tendency' )
867
868!
869!--          Prognostic equation for TKE.
870!--          Eliminate negative TKE values, which can occur due to numerical
871!--          reasons in the course of the integration. In such cases the old TKE
872!--          value is reduced by 90%.
873             DO  k = nzb_s_inner(j,i)+1, nzt
874                e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) +       &
875                             dt_3d * (                                         &
876                                      sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
877                                     )
878                IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
879             ENDDO
880
881!
882!--          Calculate tendencies for the next Runge-Kutta step
883             IF ( timestep_scheme(1:5) == 'runge' )  THEN
884                IF ( intermediate_timestep_count == 1 )  THEN
885                   DO  k = nzb_s_inner(j,i)+1, nzt
886                      te_m(k,j,i) = tend(k,j,i)
887                   ENDDO
888                ELSEIF ( intermediate_timestep_count < &
889                         intermediate_timestep_count_max )  THEN
890                   DO  k = nzb_s_inner(j,i)+1, nzt
891                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
892                   ENDDO
893                ENDIF
894             ENDIF
895
896          ENDDO
897       ENDDO
898
899       CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
900
901    ENDIF
902
903
904 END SUBROUTINE prognostic_equations_noopt
905
906
907 SUBROUTINE prognostic_equations_cache
908
909!------------------------------------------------------------------------------!
910! Version with one optimized loop over all equations. It is only allowed to
911! be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
912!
913! Here the calls of most subroutines are embedded in two DO loops over i and j,
914! so communication between CPUs is not allowed (does not make sense) within
915! these loops.
916!
917! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
918!------------------------------------------------------------------------------!
919
920    IMPLICIT NONE
921
922    CHARACTER (LEN=9) ::  time_to_string
923    INTEGER ::  i, i_omp_start, j, k, omp_get_thread_num, tn = 0
924    LOGICAL ::  loop_start
925
926
927!
928!-- Time measurement can only be performed for the whole set of equations
929    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
930
931
932!
933!-- Calculate those variables needed in the tendency terms which need
934!-- global communication
935    CALL calc_mean_profile( pt, 4 )
936    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
937    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
938    IF ( .NOT. constant_diffusion )  CALL production_e_init
939    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
940         intermediate_timestep_count == 1 )  CALL ws_statistics
941
942!
943!-- Loop over all prognostic equations
944!$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
945
946!$  tn = omp_get_thread_num()
947    loop_start = .TRUE.
948!$OMP DO
949    DO  i = nxl, nxr
950
951!
952!--    Store the first loop index. It differs for each thread and is required
953!--    later in advec_ws
954       IF ( loop_start )  THEN
955          loop_start  = .FALSE.
956          i_omp_start = i
957       ENDIF
958 
959       DO  j = nys, nyn
960!
961!--       Tendency terms for u-velocity component
962          IF ( .NOT. outflow_l  .OR.  i > nxl )  THEN
963
964             tend(:,j,i) = 0.0
965             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
966                IF ( ws_scheme_mom )  THEN
967                   IF ( outflow_l .AND. i_omp_start == nxl )  THEN 
968!                     CALL local_diss( i, j, u)    ! dissipation control
969                      CALL advec_u_ws( i, j, i_omp_start + 1, tn )
970                   ELSE
971                      CALL advec_u_ws( i, j, i_omp_start, tn )
972                   ENDIF
973                ELSE
974                   CALL advec_u_pw( i, j )
975                ENDIF
976             ELSE
977                CALL advec_u_up( i, j )
978             ENDIF
979             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
980             THEN
981                CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend,  &
982                                  u_m, usws_m, uswst_m, v_m, w_m )
983             ELSE
984                CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
985                                  usws, uswst, v, w )
986             ENDIF
987             CALL coriolis( i, j, 1 )
988             IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, &
989                                                    4 )
990
991!
992!--          Drag by plant canopy
993             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 1 )
994
995!
996!--          External pressure gradient
997             IF ( dp_external )  THEN
998                DO  k = dp_level_ind_b+1, nzt
999                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1000                ENDDO
1001             ENDIF
1002
1003             CALL user_actions( i, j, 'u-tendency' )
1004
1005!
1006!--          Prognostic equation for u-velocity component
1007             DO  k = nzb_u_inner(j,i)+1, nzt
1008                u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
1009                             dt_3d * (                                         &
1010                                   tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
1011                                     ) -                                       &
1012                             tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
1013             ENDDO
1014
1015!
1016!--          Calculate tendencies for the next Runge-Kutta step
1017             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1018                IF ( intermediate_timestep_count == 1 )  THEN
1019                   DO  k = nzb_u_inner(j,i)+1, nzt
1020                      tu_m(k,j,i) = tend(k,j,i)
1021                   ENDDO
1022                ELSEIF ( intermediate_timestep_count < &
1023                         intermediate_timestep_count_max )  THEN
1024                   DO  k = nzb_u_inner(j,i)+1, nzt
1025                      tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
1026                   ENDDO
1027                ENDIF
1028             ENDIF
1029
1030          ENDIF
1031
1032!
1033!--       Tendency terms for v-velocity component
1034          IF ( .NOT. outflow_s  .OR.  j > nys )  THEN
1035
1036             tend(:,j,i) = 0.0
1037             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
1038                IF ( ws_scheme_mom )  THEN
1039                 !   CALL local_diss( i, j, v)
1040                    CALL advec_v_ws( i, j, i_omp_start, tn )
1041                ELSE
1042                    CALL advec_v_pw( i, j )
1043                ENDIF
1044             ELSE
1045                CALL advec_v_up( i, j )
1046             ENDIF
1047             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
1048             THEN
1049                CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend,     &
1050                                  u_m, v_m, vsws_m, vswst_m, w_m )
1051             ELSE
1052                CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
1053                                  vsws, vswst, w )
1054             ENDIF
1055             CALL coriolis( i, j, 2 )
1056
1057!
1058!--          Drag by plant canopy
1059             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 2 )       
1060
1061!
1062!--          External pressure gradient
1063             IF ( dp_external )  THEN
1064                DO  k = dp_level_ind_b+1, nzt
1065                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1066                ENDDO
1067             ENDIF
1068
1069             CALL user_actions( i, j, 'v-tendency' )
1070
1071!
1072!--          Prognostic equation for v-velocity component
1073             DO  k = nzb_v_inner(j,i)+1, nzt
1074                v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
1075                             dt_3d * (                                         &
1076                                   tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
1077                                     ) -                                       &
1078                             tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1079             ENDDO
1080
1081!
1082!--          Calculate tendencies for the next Runge-Kutta step
1083             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1084                IF ( intermediate_timestep_count == 1 )  THEN
1085                   DO  k = nzb_v_inner(j,i)+1, nzt
1086                      tv_m(k,j,i) = tend(k,j,i)
1087                   ENDDO
1088                ELSEIF ( intermediate_timestep_count < &
1089                         intermediate_timestep_count_max )  THEN
1090                   DO  k = nzb_v_inner(j,i)+1, nzt
1091                      tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1092                   ENDDO
1093                ENDIF
1094             ENDIF
1095
1096          ENDIF
1097
1098!
1099!--       Tendency terms for w-velocity component
1100          tend(:,j,i) = 0.0
1101          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
1102             IF ( ws_scheme_mom )  THEN
1103             !   CALL local_diss( i, j, w)
1104                CALL advec_w_ws( i, j, i_omp_start, tn )
1105             ELSE
1106                CALL advec_w_pw( i, j )
1107             END IF
1108          ELSE
1109             CALL advec_w_up( i, j )
1110          ENDIF
1111          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
1112          THEN
1113             CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,          &
1114                               km_damp_y, tend, u_m, v_m, w_m )
1115          ELSE
1116             CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
1117                               tend, u, v, w )
1118          ENDIF
1119          CALL coriolis( i, j, 3 )
1120          IF ( ocean )  THEN
1121             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
1122          ELSE
1123             IF ( .NOT. humidity )  THEN
1124                CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
1125             ELSE
1126                CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
1127             ENDIF
1128          ENDIF
1129
1130!
1131!--       Drag by plant canopy
1132          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 3 )
1133
1134          CALL user_actions( i, j, 'w-tendency' )
1135
1136!
1137!--       Prognostic equation for w-velocity component
1138          DO  k = nzb_w_inner(j,i)+1, nzt-1
1139             w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
1140                          dt_3d * (                                         &
1141                                tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
1142                                  ) -                                       &
1143                          tsc(5) * rdf(k) * w(k,j,i)
1144          ENDDO
1145
1146!
1147!--       Calculate tendencies for the next Runge-Kutta step
1148          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1149             IF ( intermediate_timestep_count == 1 )  THEN
1150                DO  k = nzb_w_inner(j,i)+1, nzt-1
1151                   tw_m(k,j,i) = tend(k,j,i)
1152                ENDDO
1153             ELSEIF ( intermediate_timestep_count < &
1154                      intermediate_timestep_count_max )  THEN
1155                DO  k = nzb_w_inner(j,i)+1, nzt-1
1156                   tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1157                ENDDO
1158             ENDIF
1159          ENDIF
1160
1161!
1162!--       Tendency terms for potential temperature
1163          tend(:,j,i) = 0.0
1164          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
1165                IF ( ws_scheme_sca )  THEN
1166       !            CALL local_diss( i, j, pt )
1167                   CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
1168                             diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
1169                ELSE
1170                   CALL advec_s_pw( i, j, pt )
1171                ENDIF
1172          ELSE
1173             CALL advec_s_up( i, j, pt )
1174          ENDIF
1175          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
1176          THEN
1177             CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
1178                               tswst_m, wall_heatflux, tend )
1179          ELSE
1180             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
1181                               wall_heatflux, tend )
1182          ENDIF
1183
1184!
1185!--       If required compute heating/cooling due to long wave radiation
1186!--       processes
1187          IF ( radiation )  THEN
1188             CALL calc_radiation( i, j )
1189          ENDIF
1190
1191!
1192!--       If required compute impact of latent heat due to precipitation
1193          IF ( precipitation )  THEN
1194             CALL impact_of_latent_heat( i, j )
1195          ENDIF
1196
1197!
1198!--       Consideration of heat sources within the plant canopy
1199          IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
1200             CALL plant_canopy_model( i, j, 4 )
1201          ENDIF
1202
1203
1204!--       If required compute influence of large-scale subsidence/ascent
1205          IF ( large_scale_subsidence ) THEN
1206             CALL subsidence ( i, j, tend, pt, pt_init )
1207          ENDIF
1208
1209
1210          CALL user_actions( i, j, 'pt-tendency' )
1211
1212!
1213!--       Prognostic equation for potential temperature
1214          DO  k = nzb_s_inner(j,i)+1, nzt
1215             pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
1216                           dt_3d * (                                        &
1217                               tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1218                                   ) -                                      &
1219                           tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
1220          ENDDO
1221
1222!
1223!--       Calculate tendencies for the next Runge-Kutta step
1224          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1225             IF ( intermediate_timestep_count == 1 )  THEN
1226                DO  k = nzb_s_inner(j,i)+1, nzt
1227                   tpt_m(k,j,i) = tend(k,j,i)
1228                ENDDO
1229             ELSEIF ( intermediate_timestep_count < &
1230                      intermediate_timestep_count_max )  THEN
1231                DO  k = nzb_s_inner(j,i)+1, nzt
1232                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1233                                   5.3125 * tpt_m(k,j,i)
1234                ENDDO
1235             ENDIF
1236          ENDIF
1237
1238!
1239!--       If required, compute prognostic equation for salinity
1240          IF ( ocean )  THEN
1241
1242!
1243!--          Tendency-terms for salinity
1244             tend(:,j,i) = 0.0
1245             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
1246             THEN
1247                IF ( ws_scheme_sca )  THEN
1248            !        CALL local_diss( i, j, sa )
1249                    CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
1250                                diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn  )
1251                ELSE
1252                    CALL advec_s_pw( i, j, sa )
1253                ENDIF
1254             ELSE
1255                CALL advec_s_up( i, j, sa )
1256             ENDIF
1257             CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
1258                               wall_salinityflux, tend )
1259
1260             CALL user_actions( i, j, 'sa-tendency' )
1261
1262!
1263!--          Prognostic equation for salinity
1264             DO  k = nzb_s_inner(j,i)+1, nzt
1265                sa_p(k,j,i) = tsc(1) * sa(k,j,i) +                          &
1266                              dt_3d * (                                     &
1267                               tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
1268                                      ) -                                   &
1269                             tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) )
1270                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
1271             ENDDO
1272
1273!
1274!--          Calculate tendencies for the next Runge-Kutta step
1275             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1276                IF ( intermediate_timestep_count == 1 )  THEN
1277                   DO  k = nzb_s_inner(j,i)+1, nzt
1278                      tsa_m(k,j,i) = tend(k,j,i)
1279                   ENDDO
1280                ELSEIF ( intermediate_timestep_count < &
1281                         intermediate_timestep_count_max )  THEN
1282                   DO  k = nzb_s_inner(j,i)+1, nzt
1283                      tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1284                                      5.3125 * tsa_m(k,j,i)
1285                   ENDDO
1286                ENDIF
1287             ENDIF
1288
1289!
1290!--          Calculate density by the equation of state for seawater
1291             CALL eqn_state_seawater( i, j )
1292
1293          ENDIF
1294
1295!
1296!--       If required, compute prognostic equation for total water content /
1297!--       scalar
1298          IF ( humidity  .OR.  passive_scalar )  THEN
1299
1300!
1301!--          Tendency-terms for total water content / scalar
1302             tend(:,j,i) = 0.0
1303             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
1304             THEN
1305                IF ( ws_scheme_sca )  THEN
1306          !         CALL local_diss( i, j, q )
1307                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 
1308                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
1309                ELSE
1310                   CALL advec_s_pw( i, j, q )
1311                ENDIF
1312             ELSE
1313                CALL advec_s_up( i, j, q )
1314             ENDIF
1315             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
1316             THEN
1317                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
1318                                  qswst_m, wall_qflux, tend )
1319             ELSE
1320                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
1321                                  wall_qflux, tend )
1322             ENDIF
1323       
1324!
1325!--          If required compute decrease of total water content due to
1326!--          precipitation
1327             IF ( precipitation )  THEN
1328                CALL calc_precipitation( i, j )
1329             ENDIF
1330
1331!
1332!--          Sink or source of scalar concentration due to canopy elements
1333             IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
1334
1335!--          If required compute influence of large-scale subsidence/ascent
1336             IF ( large_scale_subsidence ) THEN
1337                CALL subsidence ( i, j, tend, q, q_init )
1338             ENDIF
1339
1340             CALL user_actions( i, j, 'q-tendency' )
1341
1342!
1343!--          Prognostic equation for total water content / scalar
1344             DO  k = nzb_s_inner(j,i)+1, nzt
1345                q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
1346                             dt_3d * (                                      &
1347                                tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
1348                                     ) -                                    &
1349                             tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) )
1350                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
1351             ENDDO
1352
1353!
1354!--          Calculate tendencies for the next Runge-Kutta step
1355             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1356                IF ( intermediate_timestep_count == 1 )  THEN
1357                   DO  k = nzb_s_inner(j,i)+1, nzt
1358                      tq_m(k,j,i) = tend(k,j,i)
1359                   ENDDO
1360                ELSEIF ( intermediate_timestep_count < &
1361                         intermediate_timestep_count_max )  THEN
1362                   DO  k = nzb_s_inner(j,i)+1, nzt
1363                      tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1364                                     5.3125 * tq_m(k,j,i)
1365                   ENDDO
1366                ENDIF
1367             ENDIF
1368
1369          ENDIF
1370
1371!
1372!--       If required, compute prognostic equation for turbulent kinetic
1373!--       energy (TKE)
1374          IF ( .NOT. constant_diffusion )  THEN
1375
1376!
1377!--          Tendency-terms for TKE
1378             tend(:,j,i) = 0.0
1379             IF ( ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  &
1380                 .AND.  .NOT. use_upstream_for_tke )  THEN
1381                 IF ( ws_scheme_sca )  THEN
1382                 !    CALL local_diss( i, j, e )
1383                     CALL advec_s_ws( i, j, e, 'e', flux_s_e, &
1384                                diss_s_e, flux_l_e, diss_l_e , i_omp_start, tn )
1385                 ELSE
1386                     CALL advec_s_pw( i, j, e )
1387                 ENDIF
1388             ELSE
1389                CALL advec_s_up( i, j, e )
1390             ENDIF
1391             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
1392             THEN
1393                IF ( .NOT. humidity )  THEN
1394                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
1395                                     km_m, l_grid, pt_m, pt_reference,   &
1396                                     rif_m, tend, zu, zw )
1397                ELSE
1398                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
1399                                     km_m, l_grid, vpt_m, pt_reference,  &
1400                                     rif_m, tend, zu, zw )
1401                ENDIF
1402             ELSE
1403                IF ( .NOT. humidity )  THEN
1404                   IF ( ocean )  THEN
1405                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
1406                                        km, l_grid, prho, prho_reference,  &
1407                                        rif, tend, zu, zw )
1408                   ELSE
1409                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
1410                                        km, l_grid, pt, pt_reference, rif, &
1411                                        tend, zu, zw )
1412                   ENDIF
1413                ELSE
1414                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
1415                                     l_grid, vpt, pt_reference, rif, tend, &
1416                                     zu, zw )
1417                ENDIF
1418             ENDIF
1419             CALL production_e( i, j )
1420
1421!
1422!--          Additional sink term for flows through plant canopies
1423             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 6 ) 
1424
1425             CALL user_actions( i, j, 'e-tendency' )
1426
1427!
1428!--          Prognostic equation for TKE.
1429!--          Eliminate negative TKE values, which can occur due to numerical
1430!--          reasons in the course of the integration. In such cases the old
1431!--          TKE value is reduced by 90%.
1432             DO  k = nzb_s_inner(j,i)+1, nzt
1433                e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
1434                             dt_3d * (                                      &
1435                                tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
1436                                     )
1437                IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
1438             ENDDO
1439
1440!
1441!--          Calculate tendencies for the next Runge-Kutta step
1442             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1443                IF ( intermediate_timestep_count == 1 )  THEN
1444                   DO  k = nzb_s_inner(j,i)+1, nzt
1445                      te_m(k,j,i) = tend(k,j,i)
1446                   ENDDO
1447                ELSEIF ( intermediate_timestep_count < &
1448                         intermediate_timestep_count_max )  THEN
1449                   DO  k = nzb_s_inner(j,i)+1, nzt
1450                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1451                                     5.3125 * te_m(k,j,i)
1452                   ENDDO
1453                ENDIF
1454             ENDIF
1455
1456          ENDIF   ! TKE equation
1457
1458       ENDDO
1459    ENDDO
1460!$OMP END PARALLEL
1461
1462    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1463
1464
1465 END SUBROUTINE prognostic_equations_cache
1466
1467
1468 SUBROUTINE prognostic_equations_vector
1469
1470!------------------------------------------------------------------------------!
1471! Version for vector machines
1472!------------------------------------------------------------------------------!
1473
1474    IMPLICIT NONE
1475
1476    CHARACTER (LEN=9) ::  time_to_string
1477    INTEGER ::  i, j, k
1478    REAL    ::  sat, sbt
1479
1480!
1481!-- Calculate those variables needed in the tendency terms which need
1482!-- global communication
1483    CALL calc_mean_profile( pt, 4 )
1484    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
1485    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
1486    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
1487         intermediate_timestep_count == 1 )  CALL ws_statistics
1488
1489!
1490!-- u-velocity component
1491    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1492
1493!
1494!-- u-tendency terms with communication
1495    IF ( momentum_advec == 'ups-scheme' )  THEN
1496       tend = 0.0
1497       CALL advec_u_ups
1498    ENDIF
1499
1500!
1501!-- u-tendency terms with no communication
1502    IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
1503       tend = 0.0
1504       IF ( ws_scheme_mom )  THEN
1505          CALL advec_u_ws
1506       ELSE
1507          CALL advec_u_pw
1508       ENDIF
1509    ELSE
1510       IF ( momentum_advec /= 'ups-scheme' )  THEN
1511          tend = 0.0
1512          CALL advec_u_up
1513       ENDIF
1514    ENDIF
1515    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
1516       CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
1517                         uswst_m, v_m, w_m )
1518    ELSE
1519       CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
1520    ENDIF
1521    CALL coriolis( 1 )
1522    IF ( sloping_surface )  CALL buoyancy( pt, pt_reference, 1, 4 )
1523
1524!
1525!-- Drag by plant canopy
1526    IF ( plant_canopy )  CALL plant_canopy_model( 1 )
1527
1528!
1529!-- External pressure gradient
1530    IF ( dp_external )  THEN
1531       DO  i = nxlu, nxr
1532          DO  j = nys, nyn
1533             DO  k = dp_level_ind_b+1, nzt
1534                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1535             ENDDO
1536          ENDDO
1537       ENDDO
1538    ENDIF
1539
1540    CALL user_actions( 'u-tendency' )
1541
1542!
1543!-- Prognostic equation for u-velocity component
1544    DO  i = nxlu, nxr
1545       DO  j = nys, nyn
1546          DO  k = nzb_u_inner(j,i)+1, nzt
1547             u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) +    &
1548                          dt_3d * (                                            &
1549                                   tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
1550                                  ) -                                          &
1551                          tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
1552          ENDDO
1553       ENDDO
1554    ENDDO
1555
1556!
1557!-- Calculate tendencies for the next Runge-Kutta step
1558    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1559       IF ( intermediate_timestep_count == 1 )  THEN
1560          DO  i = nxlu, nxr
1561             DO  j = nys, nyn
1562                DO  k = nzb_u_inner(j,i)+1, nzt
1563                   tu_m(k,j,i) = tend(k,j,i)
1564                ENDDO
1565             ENDDO
1566          ENDDO
1567       ELSEIF ( intermediate_timestep_count < &
1568                intermediate_timestep_count_max )  THEN
1569          DO  i = nxlu, nxr
1570             DO  j = nys, nyn
1571                DO  k = nzb_u_inner(j,i)+1, nzt
1572                   tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
1573                ENDDO
1574             ENDDO
1575          ENDDO
1576       ENDIF
1577    ENDIF
1578
1579    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1580
1581!
1582!-- v-velocity component
1583    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1584
1585!
1586!-- v-tendency terms with communication
1587    IF ( momentum_advec == 'ups-scheme' )  THEN
1588       tend = 0.0
1589       CALL advec_v_ups
1590    ENDIF
1591
1592!
1593!-- v-tendency terms with no communication
1594    IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
1595       tend = 0.0
1596       IF ( ws_scheme_mom )  THEN
1597          CALL advec_v_ws
1598       ELSE
1599          CALL advec_v_pw
1600       END IF
1601    ELSE
1602       IF ( momentum_advec /= 'ups-scheme' )  THEN
1603          tend = 0.0
1604          CALL advec_v_up
1605       ENDIF
1606    ENDIF
1607    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
1608       CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
1609                         vswst_m, w_m )
1610    ELSE
1611       CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
1612    ENDIF
1613    CALL coriolis( 2 )
1614
1615!
1616!-- Drag by plant canopy
1617    IF ( plant_canopy )  CALL plant_canopy_model( 2 )
1618
1619!
1620!-- External pressure gradient
1621    IF ( dp_external )  THEN
1622       DO  i = nxl, nxr
1623          DO  j = nysv, nyn
1624             DO  k = dp_level_ind_b+1, nzt
1625                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1626             ENDDO
1627          ENDDO
1628       ENDDO
1629    ENDIF
1630
1631    CALL user_actions( 'v-tendency' )
1632
1633!
1634!-- Prognostic equation for v-velocity component
1635    DO  i = nxl, nxr
1636       DO  j = nysv, nyn
1637          DO  k = nzb_v_inner(j,i)+1, nzt
1638             v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) +    &
1639                          dt_3d * (                                            &
1640                                   tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
1641                                  ) -                                          &
1642                          tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1643          ENDDO
1644       ENDDO
1645    ENDDO
1646
1647!
1648!-- Calculate tendencies for the next Runge-Kutta step
1649    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1650       IF ( intermediate_timestep_count == 1 )  THEN
1651          DO  i = nxl, nxr
1652             DO  j = nysv, nyn
1653                DO  k = nzb_v_inner(j,i)+1, nzt
1654                   tv_m(k,j,i) = tend(k,j,i)
1655                ENDDO
1656             ENDDO
1657          ENDDO
1658       ELSEIF ( intermediate_timestep_count < &
1659                intermediate_timestep_count_max )  THEN
1660          DO  i = nxl, nxr
1661             DO  j = nysv, nyn
1662                DO  k = nzb_v_inner(j,i)+1, nzt
1663                   tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1664                ENDDO
1665             ENDDO
1666          ENDDO
1667       ENDIF
1668    ENDIF
1669
1670    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1671
1672!
1673!-- w-velocity component
1674    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1675
1676!
1677!-- w-tendency terms with communication
1678    IF ( momentum_advec == 'ups-scheme' )  THEN
1679       tend = 0.0
1680       CALL advec_w_ups
1681    ENDIF
1682
1683!
1684!-- w-tendency terms with no communication
1685    IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
1686       tend = 0.0
1687       IF ( ws_scheme_mom )  THEN
1688          CALL advec_w_ws
1689       ELSE
1690          CALL advec_w_pw
1691       ENDIF
1692    ELSE
1693       IF ( momentum_advec /= 'ups-scheme' )  THEN
1694          tend = 0.0
1695          CALL advec_w_up
1696       ENDIF
1697    ENDIF
1698    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
1699       CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
1700                         v_m, w_m )
1701    ELSE
1702       CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
1703    ENDIF
1704    CALL coriolis( 3 )
1705    IF ( ocean )  THEN
1706       CALL buoyancy( rho, rho_reference, 3, 64 )
1707    ELSE
1708       IF ( .NOT. humidity )  THEN
1709          CALL buoyancy( pt, pt_reference, 3, 4 )
1710       ELSE
1711          CALL buoyancy( vpt, pt_reference, 3, 44 )
1712       ENDIF
1713    ENDIF
1714
1715!
1716!-- Drag by plant canopy
1717    IF ( plant_canopy )  CALL plant_canopy_model( 3 )
1718
1719    CALL user_actions( 'w-tendency' )
1720
1721!
1722!-- Prognostic equation for w-velocity component
1723    DO  i = nxl, nxr
1724       DO  j = nys, nyn
1725          DO  k = nzb_w_inner(j,i)+1, nzt-1
1726             w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) +      &
1727                          dt_3d * (                                            &
1728                                   tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
1729                                  ) -                                          &
1730                          tsc(5) * rdf(k) * w(k,j,i)
1731          ENDDO
1732       ENDDO
1733    ENDDO
1734
1735!
1736!-- Calculate tendencies for the next Runge-Kutta step
1737    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1738       IF ( intermediate_timestep_count == 1 )  THEN
1739          DO  i = nxl, nxr
1740             DO  j = nys, nyn
1741                DO  k = nzb_w_inner(j,i)+1, nzt-1
1742                   tw_m(k,j,i) = tend(k,j,i)
1743                ENDDO
1744             ENDDO
1745          ENDDO
1746       ELSEIF ( intermediate_timestep_count < &
1747                intermediate_timestep_count_max )  THEN
1748          DO  i = nxl, nxr
1749             DO  j = nys, nyn
1750                DO  k = nzb_w_inner(j,i)+1, nzt-1
1751                   tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1752                ENDDO
1753             ENDDO
1754          ENDDO
1755       ENDIF
1756    ENDIF
1757
1758    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1759
1760!
1761!-- potential temperature
1762    CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1763
1764!
1765!-- pt-tendency terms with communication
1766    sat = tsc(1)
1767    sbt = tsc(2)
1768    IF ( scalar_advec == 'bc-scheme' )  THEN
1769
1770       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1771!
1772!--       Bott-Chlond scheme always uses Euler time step when leapfrog is
1773!--       switched on. Thus:
1774          sat = 1.0
1775          sbt = 1.0
1776       ENDIF
1777       tend = 0.0
1778       CALL advec_s_bc( pt, 'pt' )
1779    ELSE
1780       IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
1781          tend = 0.0
1782          CALL advec_s_ups( pt, 'pt' )
1783       ENDIF
1784    ENDIF
1785         
1786!
1787!-- pt-tendency terms with no communication
1788    IF ( scalar_advec == 'bc-scheme' )  THEN
1789       CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
1790                         tend )
1791    ELSE
1792       IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
1793          tend = 0.0
1794          IF ( ws_scheme_sca )  THEN
1795             CALL advec_s_ws( pt, 'pt' )
1796          ELSE
1797             CALL advec_s_pw( pt )
1798          ENDIF
1799       ELSE
1800          IF ( scalar_advec /= 'ups-scheme' )  THEN
1801             tend = 0.0
1802             CALL advec_s_up( pt )
1803          ENDIF
1804       ENDIF
1805       IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
1806          CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, &
1807                            wall_heatflux, tend )
1808       ELSE
1809          CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
1810                            tend )
1811       ENDIF
1812    ENDIF
1813
1814!
1815!-- If required compute heating/cooling due to long wave radiation
1816!-- processes
1817    IF ( radiation )  THEN
1818       CALL calc_radiation
1819    ENDIF
1820
1821!
1822!-- If required compute impact of latent heat due to precipitation
1823    IF ( precipitation )  THEN
1824       CALL impact_of_latent_heat
1825    ENDIF
1826
1827!
1828!-- Consideration of heat sources within the plant canopy
1829    IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
1830       CALL plant_canopy_model( 4 )
1831    ENDIF
1832
1833!--If required compute influence of large-scale subsidence/ascent
1834   IF ( large_scale_subsidence ) THEN
1835      CALL subsidence ( tend, pt, pt_init )
1836   ENDIF
1837
1838    CALL user_actions( 'pt-tendency' )
1839
1840!
1841!-- Prognostic equation for potential temperature
1842    DO  i = nxl, nxr
1843       DO  j = nys, nyn
1844          DO  k = nzb_s_inner(j,i)+1, nzt
1845             pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +       &
1846                           dt_3d * (                                           &
1847                                     sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1848                                   ) -                                         &
1849                           tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
1850          ENDDO
1851       ENDDO
1852    ENDDO
1853
1854!
1855!-- Calculate tendencies for the next Runge-Kutta step
1856    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1857       IF ( intermediate_timestep_count == 1 )  THEN
1858          DO  i = nxl, nxr
1859             DO  j = nys, nyn
1860                DO  k = nzb_s_inner(j,i)+1, nzt
1861                   tpt_m(k,j,i) = tend(k,j,i)
1862                ENDDO
1863             ENDDO
1864          ENDDO
1865       ELSEIF ( intermediate_timestep_count < &
1866                intermediate_timestep_count_max )  THEN
1867          DO  i = nxl, nxr
1868             DO  j = nys, nyn
1869                DO  k = nzb_s_inner(j,i)+1, nzt
1870                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
1871                ENDDO
1872             ENDDO
1873          ENDDO
1874       ENDIF
1875    ENDIF
1876
1877    CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1878
1879!
1880!-- If required, compute prognostic equation for salinity
1881    IF ( ocean )  THEN
1882
1883       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1884
1885!
1886!--    sa-tendency terms with communication
1887       sat = tsc(1)
1888       sbt = tsc(2)
1889       IF ( scalar_advec == 'bc-scheme' )  THEN
1890
1891          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1892!
1893!--          Bott-Chlond scheme always uses Euler time step when leapfrog is
1894!--          switched on. Thus:
1895             sat = 1.0
1896             sbt = 1.0
1897          ENDIF
1898          tend = 0.0
1899          CALL advec_s_bc( sa, 'sa' )
1900       ELSE
1901          IF ( tsc(2) /= 2.0 )  THEN
1902             IF ( scalar_advec == 'ups-scheme' )  THEN
1903                tend = 0.0
1904                CALL advec_s_ups( sa, 'sa' )
1905             ENDIF
1906          ENDIF
1907       ENDIF
1908
1909!
1910!--    sa-tendency terms with no communication
1911       IF ( scalar_advec == 'bc-scheme' )  THEN
1912          CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
1913                            wall_salinityflux, tend )
1914       ELSE
1915          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
1916             tend = 0.0
1917             IF ( ws_scheme_sca )  THEN
1918                 CALL advec_s_ws( sa, 'sa' )
1919             ELSE
1920                 CALL advec_s_pw( sa )
1921             ENDIF
1922          ELSE
1923             IF ( scalar_advec /= 'ups-scheme' )  THEN
1924                tend = 0.0
1925                CALL advec_s_up( sa )
1926             ENDIF
1927          ENDIF
1928          CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
1929                            wall_salinityflux, tend )
1930       ENDIF
1931       
1932       CALL user_actions( 'sa-tendency' )
1933
1934!
1935!--    Prognostic equation for salinity
1936       DO  i = nxl, nxr
1937          DO  j = nys, nyn
1938             DO  k = nzb_s_inner(j,i)+1, nzt
1939                sa_p(k,j,i) = sat * sa(k,j,i) +                                &
1940                              dt_3d * (                                        &
1941                                     sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
1942                                      ) -                                      &
1943                              tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) )
1944                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
1945             ENDDO
1946          ENDDO
1947       ENDDO
1948
1949!
1950!--    Calculate tendencies for the next Runge-Kutta step
1951       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1952          IF ( intermediate_timestep_count == 1 )  THEN
1953             DO  i = nxl, nxr
1954                DO  j = nys, nyn
1955                   DO  k = nzb_s_inner(j,i)+1, nzt
1956                      tsa_m(k,j,i) = tend(k,j,i)
1957                   ENDDO
1958                ENDDO
1959             ENDDO
1960          ELSEIF ( intermediate_timestep_count < &
1961                   intermediate_timestep_count_max )  THEN
1962             DO  i = nxl, nxr
1963                DO  j = nys, nyn
1964                   DO  k = nzb_s_inner(j,i)+1, nzt
1965                      tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1966                                      5.3125 * tsa_m(k,j,i)
1967                   ENDDO
1968                ENDDO
1969             ENDDO
1970          ENDIF
1971       ENDIF
1972
1973       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1974
1975!
1976!--    Calculate density by the equation of state for seawater
1977       CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1978       CALL eqn_state_seawater
1979       CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1980
1981    ENDIF
1982
1983!
1984!-- If required, compute prognostic equation for total water content / scalar
1985    IF ( humidity  .OR.  passive_scalar )  THEN
1986
1987       CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1988
1989!
1990!--    Scalar/q-tendency terms with communication
1991       sat = tsc(1)
1992       sbt = tsc(2)
1993       IF ( scalar_advec == 'bc-scheme' )  THEN
1994
1995          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1996!
1997!--          Bott-Chlond scheme always uses Euler time step when leapfrog is
1998!--          switched on. Thus:
1999             sat = 1.0
2000             sbt = 1.0
2001          ENDIF
2002          tend = 0.0
2003          CALL advec_s_bc( q, 'q' )
2004       ELSE
2005          IF ( tsc(2) /= 2.0 )  THEN
2006             IF ( scalar_advec == 'ups-scheme' )  THEN
2007                tend = 0.0
2008                CALL advec_s_ups( q, 'q' )
2009             ENDIF
2010          ENDIF
2011       ENDIF
2012
2013!
2014!--    Scalar/q-tendency terms with no communication
2015       IF ( scalar_advec == 'bc-scheme' )  THEN
2016          CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, wall_qflux, tend )
2017       ELSE
2018          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
2019             tend = 0.0
2020             IF ( ws_scheme_sca )  THEN
2021                CALL advec_s_ws( q, 'q' )
2022             ELSE
2023                CALL advec_s_pw( q )
2024             ENDIF
2025          ELSE
2026             IF ( scalar_advec /= 'ups-scheme' )  THEN
2027                tend = 0.0
2028                CALL advec_s_up( q )
2029             ENDIF
2030          ENDIF
2031          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
2032             CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, &
2033                               wall_qflux, tend )
2034          ELSE
2035             CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, &
2036                               wall_qflux, tend )
2037          ENDIF
2038       ENDIF
2039       
2040!
2041!--    If required compute decrease of total water content due to
2042!--    precipitation
2043       IF ( precipitation )  THEN
2044          CALL calc_precipitation
2045       ENDIF
2046
2047!
2048!--    Sink or source of scalar concentration due to canopy elements
2049       IF ( plant_canopy ) CALL plant_canopy_model( 5 )
2050       
2051!
2052!--    If required compute influence of large-scale subsidence/ascent
2053       IF ( large_scale_subsidence ) THEN
2054         CALL subsidence ( tend, q, q_init )
2055       ENDIF
2056
2057       CALL user_actions( 'q-tendency' )
2058
2059!
2060!--    Prognostic equation for total water content / scalar
2061       DO  i = nxl, nxr
2062          DO  j = nys, nyn
2063             DO  k = nzb_s_inner(j,i)+1, nzt
2064                q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) +       &
2065                             dt_3d * (                                         &
2066                                      sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
2067                                     ) -                                       &
2068                             tsc(5) * rdf_sc(k) * ( q(k,j,i) - q_init(k) )
2069                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
2070             ENDDO
2071          ENDDO
2072       ENDDO
2073
2074!
2075!--    Calculate tendencies for the next Runge-Kutta step
2076       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2077          IF ( intermediate_timestep_count == 1 )  THEN
2078             DO  i = nxl, nxr
2079                DO  j = nys, nyn
2080                   DO  k = nzb_s_inner(j,i)+1, nzt
2081                      tq_m(k,j,i) = tend(k,j,i)
2082                   ENDDO
2083                ENDDO
2084             ENDDO
2085          ELSEIF ( intermediate_timestep_count < &
2086                   intermediate_timestep_count_max )  THEN
2087             DO  i = nxl, nxr
2088                DO  j = nys, nyn
2089                   DO  k = nzb_s_inner(j,i)+1, nzt
2090                      tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
2091                   ENDDO
2092                ENDDO
2093             ENDDO
2094          ENDIF
2095       ENDIF
2096
2097       CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
2098
2099    ENDIF
2100
2101!
2102!-- If required, compute prognostic equation for turbulent kinetic
2103!-- energy (TKE)
2104    IF ( .NOT. constant_diffusion )  THEN
2105
2106       CALL cpu_log( log_point(16), 'tke-equation', 'start' )
2107
2108!
2109!--    TKE-tendency terms with communication
2110       CALL production_e_init
2111
2112       sat = tsc(1)
2113       sbt = tsc(2)
2114       IF ( .NOT. use_upstream_for_tke )  THEN
2115          IF ( scalar_advec == 'bc-scheme' )  THEN
2116
2117             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2118!
2119!--             Bott-Chlond scheme always uses Euler time step when leapfrog is
2120!--             switched on. Thus:
2121                sat = 1.0
2122                sbt = 1.0
2123             ENDIF
2124             tend = 0.0
2125             CALL advec_s_bc( e, 'e' )
2126          ELSE
2127             IF ( tsc(2) /= 2.0 )  THEN
2128                IF ( scalar_advec == 'ups-scheme' )  THEN
2129                   tend = 0.0
2130                   CALL advec_s_ups( e, 'e' )
2131                ENDIF
2132             ENDIF
2133          ENDIF
2134       ENDIF
2135
2136!
2137!--    TKE-tendency terms with no communication
2138       IF ( scalar_advec == 'bc-scheme'  .AND.  .NOT. use_upstream_for_tke )  &
2139       THEN
2140          IF ( .NOT. humidity )  THEN
2141             IF ( ocean )  THEN
2142                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
2143                                  prho, prho_reference, rif, tend, zu, zw )
2144             ELSE
2145                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
2146                                  pt_reference, rif, tend, zu, zw )
2147             ENDIF
2148          ELSE
2149             CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
2150                               pt_reference, rif, tend, zu, zw )
2151          ENDIF
2152       ELSE
2153          IF ( use_upstream_for_tke )  THEN
2154             tend = 0.0
2155             CALL advec_s_up( e )
2156          ELSE
2157             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
2158                tend = 0.0
2159                IF ( ws_scheme_sca )  THEN
2160                   CALL advec_s_ws( e, 'e' )
2161                ELSE
2162                   CALL advec_s_pw( e )
2163                ENDIF
2164             ELSE
2165                IF ( scalar_advec /= 'ups-scheme' )  THEN
2166                   tend = 0.0
2167                   CALL advec_s_up( e )
2168                ENDIF
2169             ENDIF
2170          ENDIF
2171          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
2172             IF ( .NOT. humidity )  THEN
2173                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
2174                                  pt_m, pt_reference, rif_m, tend, zu, zw )
2175             ELSE
2176                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
2177                                  vpt_m, pt_reference, rif_m, tend, zu, zw )
2178             ENDIF
2179          ELSE
2180             IF ( .NOT. humidity )  THEN
2181                IF ( ocean )  THEN
2182                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
2183                                     prho, prho_reference, rif, tend, zu, zw )
2184                ELSE
2185                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
2186                                     pt, pt_reference, rif, tend, zu, zw )
2187                ENDIF
2188             ELSE
2189                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
2190                                  pt_reference, rif, tend, zu, zw )
2191             ENDIF
2192          ENDIF
2193       ENDIF
2194       CALL production_e
2195
2196!
2197!--    Additional sink term for flows through plant canopies
2198       IF ( plant_canopy )  CALL plant_canopy_model( 6 )
2199       CALL user_actions( 'e-tendency' )
2200
2201!
2202!--    Prognostic equation for TKE.
2203!--    Eliminate negative TKE values, which can occur due to numerical
2204!--    reasons in the course of the integration. In such cases the old TKE
2205!--    value is reduced by 90%.
2206       DO  i = nxl, nxr
2207          DO  j = nys, nyn
2208             DO  k = nzb_s_inner(j,i)+1, nzt
2209                e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) +       &
2210                             dt_3d * (                                         &
2211                                      sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
2212                                     )
2213                IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
2214             ENDDO
2215          ENDDO
2216       ENDDO
2217
2218!
2219!--    Calculate tendencies for the next Runge-Kutta step
2220       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2221          IF ( intermediate_timestep_count == 1 )  THEN
2222             DO  i = nxl, nxr
2223                DO  j = nys, nyn
2224                   DO  k = nzb_s_inner(j,i)+1, nzt
2225                      te_m(k,j,i) = tend(k,j,i)
2226                   ENDDO
2227                ENDDO
2228             ENDDO
2229          ELSEIF ( intermediate_timestep_count < &
2230                   intermediate_timestep_count_max )  THEN
2231             DO  i = nxl, nxr
2232                DO  j = nys, nyn
2233                   DO  k = nzb_s_inner(j,i)+1, nzt
2234                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
2235                   ENDDO
2236                ENDDO
2237             ENDDO
2238          ENDIF
2239       ENDIF
2240
2241       CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
2242
2243    ENDIF
2244
2245
2246 END SUBROUTINE prognostic_equations_vector
2247
2248
2249 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.