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

Last change on this file since 737 was 737, checked in by suehring, 10 years ago

last commit documented

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