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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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