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

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

last commit documented

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