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

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