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

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

last commit documented

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