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

Last change on this file since 1274 was 1258, checked in by raasch, 11 years ago

last commit documented

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