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

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

last commit documented

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