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

Last change on this file since 847 was 786, checked in by raasch, 13 years ago

last commit documented

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