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

Last change on this file since 550 was 532, checked in by heinze, 14 years ago

NCL scripts allow the output of png files

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