source: palm/tags/release-3.9/SOURCE/prognostic_equations.f90 @ 1038

Last change on this file since 1038 was 1037, checked in by raasch, 11 years ago

last commit documented

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