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

Last change on this file since 407 was 407, checked in by maronga, 12 years ago

humidity for non-flat topography implemented, re-adjustments for lcxt4

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