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

Last change on this file since 1242 was 1242, checked in by heinze, 10 years ago

Last commmit documented

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