Home

Context Navigation

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

Last change on this file since 2157 was 2156, checked in by hoffmann, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 66.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |