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

Last change on this file since 669 was 669, checked in by gryschka, 11 years ago

last commit documented

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