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

Last change on this file since 97 was 97, checked in by raasch, 14 years ago

New:
---
ocean version including prognostic equation for salinity and equation of state for seawater. Routine buoyancy can be used with both temperature and density.
+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux

advec_s_bc, average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, netcdf, parin, production_e, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary

New:
eqn_state_seawater, init_ocean

Changed:


inipar-parameter use_pt_reference renamed use_reference

hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile

format adjustments for the ocean version (run_control)

advec_particles, buoyancy, calc_liquid_water_content, check_parameters, diffusion_e, diffusivities, header, init_cloud_physics, modules, production_e, prognostic_equations, run_control

Errors:


Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).

Bugfix: error in boundary condition for TKE removed (advec_s_bc)

advec_s_bc, diffusion_e, prognostic_equations

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