source: palm/tags/release-3.8/SOURCE/prognostic_equations.f90 @ 4343

Last change on this file since 4343 was 710, checked in by raasch, 13 years ago

last commit documented

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