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

Last change on this file since 96 was 96, checked in by raasch, 14 years ago

more preliminary uncomplete changes for ocean version

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