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

Last change on this file since 673 was 673, checked in by suehring, 11 years ago

Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.

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