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

Last change on this file since 2718 was 2718, checked in by maronga, 6 years ago

deleting of deprecated files; headers updated where needed

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