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

Last change on this file since 978 was 978, checked in by fricke, 9 years ago

merge fricke branch back into trunk

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