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

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