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

Last change on this file since 411 was 411, checked in by heinze, 12 years ago

Large scale vertical motion (subsidence/ascent) can be applied to the prognostic equation for the potential temperature

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