source: palm/tags/release-3.2b/SOURCE/prognostic_equations.f90 @ 1320

Last change on this file since 1320 was 77, checked in by raasch, 17 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

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