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

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

New:
---
compare_palm_logs is additionally compiled with mbuild -u (Makefile in trunk/UTIL)

make options (mopts) to be set by configuration file implemented (mrun, mbuild)

humidity=.T. is now usable for runs with topography. wall_humidityflux and
wall_scalarflux are the corresponding new parin arrays.
(check_parameters, init_3d_model, parin)

Large scale vertical motion (subsidence/ascent) can be applied to the
prognostic equation for the potential temperature. (check_parameters, header,
Makefile, modules, parin, prognostic_equations, read_var_list, subsidence,
write_var_list)

A simple method for installing and running palm (with limited features)
has been added. (Makefile, palm_simple_install, palm_simple_run)

Changed:


2d-decomposition is default for Cray-XT machines (init_pegrid)

var_ts is replaced by dots_max (modules,init_3d_model)

Every cloud droplet has now an own weighting factor and can be deleted due to
collisions. Condensation and collision of cloud droplets are adjusted
accordingly. (advec_particles)

Collision efficiency for large cloud droplets has changed according to table of
Rogers and Yau. (collision_efficiency)

Errors:


Bugfix for generating serial jobs (subjob)

Bugfix: index problem concerning gradient_level indices removed (header)

Dimension of array stat in cascade change to prevent type problems with
mpi2 libraries (poisfft_hybrid)

Loop was split to make runs reproducible when using ifort compiler.
(disturb_field)

Bugfix: exchange of ghost points for prho included (time_integration)

Bugfix: calculation of time-averaged surface heatfluxes (sum_up_3d_data)

Bugfix: calculation of precipitation_rate (calc_precipitation)

Bugfix: initial data assignments to some dvrp arrays changed due to error
messages from gfortran compiler (modules)

Bugfix: calculation of cloud droplet velocity (advec_particles)

Bugfix: transfer of particles at south/left edge (advec_particles)

Bugfix: calculation of collision_efficiency (collision_efficiency)

Bugfix: initialisation of var_mod (subsidence)

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