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

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

last commit documented

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