Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |