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

Last change on this file since 987 was 979, checked in by fricke, 12 years ago

last commit documented

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