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

Last change on this file since 388 was 388, checked in by raasch, 12 years ago

in-situ AND potential density are calculated and used in the ocean version

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