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

Last change on this file since 708 was 674, checked in by suehring, 13 years ago

last commit documented

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