Home

Context Navigation

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

Last change on this file since 2239 was 2233, checked in by suehring, 7 years ago
last commit documented
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 2233 2017-05-30 18:08:54Z suehring $
27	!
28	! 2232 2017-05-30 17:47:52Z suehring
29	! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
30	! which is realized directly in diffusion_s now.
31	!
32	! 2192 2017-03-22 04:14:10Z raasch
33	! Bugfix for misplaced and missing openMP directives from r2155
34	!
35	! 2155 2017-02-21 09:57:40Z hoffmann
36	! Bugfix in the calculation of microphysical quantities on ghost points.
37	!
38	! 2118 2017-01-17 16:38:49Z raasch
39	! OpenACC version of subroutine removed
40	!
41	! 2031 2016-10-21 15:11:58Z knoop
42	! renamed variable rho to rho_ocean
43	!
44	! 2011 2016-09-19 17:29:57Z kanani
45	! Flag urban_surface is now defined in module control_parameters.
46	!
47	! 2007 2016-08-24 15:47:17Z kanani
48	! Added pt tendency calculation based on energy balance at urban surfaces
49	! (new urban surface model)
50	!
51	! 2000 2016-08-20 18:09:15Z knoop
52	! Forced header and separation lines into 80 columns
53	!
54	! 1976 2016-07-27 13:28:04Z maronga
55	! Simplied calls to radiation model
56	!
57	! 1960 2016-07-12 16:34:24Z suehring
58	! Separate humidity and passive scalar
59	!
60	! 1914 2016-05-26 14:44:07Z witha
61	! Added calls for wind turbine model
62	!
63	! 1873 2016-04-18 14:50:06Z maronga
64	! Module renamed (removed _mod)
65	!
66	! 1850 2016-04-08 13:29:27Z maronga
67	! Module renamed
68	!
69	! 1826 2016-04-07 12:01:39Z maronga
70	! Renamed canopy model calls.
71	!
72	! 1822 2016-04-07 07:49:42Z hoffmann
73	! Kessler microphysics scheme moved to microphysics.
74	!
75	! 1757 2016-02-22 15:49:32Z maronga
76	!
77	! 1691 2015-10-26 16:17:44Z maronga
78	! Added optional model spin-up without radiation / land surface model calls.
79	! Formatting corrections.
80	!
81	! 1682 2015-10-07 23:56:08Z knoop
82	! Code annotations made doxygen readable
83	!
84	! 1585 2015-04-30 07:05:52Z maronga
85	! Added call for temperature tendency calculation due to radiative flux divergence
86	!
87	! 1517 2015-01-07 19:12:25Z hoffmann
88	! advec_s_bc_mod addded, since advec_s_bc is now a module
89	!
90	! 1496 2014-12-02 17:25:50Z maronga
91	! Renamed "radiation" -> "cloud_top_radiation"
92	!
93	! 1484 2014-10-21 10:53:05Z kanani
94	! Changes due to new module structure of the plant canopy model:
95	! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
96	! Removed double-listing of use_upstream_for_tke in ONLY-list of module
97	! control_parameters
98	!
99	! 1409 2014-05-23 12:11:32Z suehring
100	! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
101	! This ensures that left-hand side fluxes are also calculated for nxl in that
102	! case, even though the solution at nxl is overwritten in boundary_conds()
103	!
104	! 1398 2014-05-07 11:15:00Z heinze
105	! Rayleigh-damping for horizontal velocity components changed: instead of damping
106	! against ug and vg, damping against u_init and v_init is used to allow for a
107	! homogenized treatment in case of nudging
108	!
109	! 1380 2014-04-28 12:40:45Z heinze
110	! Change order of calls for scalar prognostic quantities:
111	! ls_advec -> nudging -> subsidence since initial profiles
112	!
113	! 1374 2014-04-25 12:55:07Z raasch
114	! missing variables added to ONLY lists
115	!
116	! 1365 2014-04-22 15:03:56Z boeske
117	! Calls of ls_advec for large scale advection added,
118	! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
119	! new argument ls_index added to the calls of subsidence
120	! +ls_index
121	!
122	! 1361 2014-04-16 15:17:48Z hoffmann
123	! Two-moment microphysics moved to the start of prognostic equations. This makes
124	! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
125	! Additionally, it is allowed to call the microphysics just once during the time
126	! step (not at each sub-time step).
127	!
128	! Two-moment cloud physics added for vector and accelerator optimization.
129	!
130	! 1353 2014-04-08 15:21:23Z heinze
131	! REAL constants provided with KIND-attribute
132	!
133	! 1337 2014-03-25 15:11:48Z heinze
134	! Bugfix: REAL constants provided with KIND-attribute
135	!
136	! 1332 2014-03-25 11:59:43Z suehring
137	! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
138	!
139	! 1330 2014-03-24 17:29:32Z suehring
140	! In case of SGS-particle velocity advection of TKE is also allowed with
141	! dissipative 5th-order scheme.
142	!
143	! 1320 2014-03-20 08:40:49Z raasch
144	! ONLY-attribute added to USE-statements,
145	! kind-parameters added to all INTEGER and REAL declaration statements,
146	! kinds are defined in new module kinds,
147	! old module precision_kind is removed,
148	! revision history before 2012 removed,
149	! comment fields (!:) to be used for variable explanations added to
150	! all variable declaration statements
151	!
152	! 1318 2014-03-17 13:35:16Z raasch
153	! module interfaces removed
154	!
155	! 1257 2013-11-08 15:18:40Z raasch
156	! openacc loop vector clauses removed, independent clauses added
157	!
158	! 1246 2013-11-01 08:59:45Z heinze
159	! enable nudging also for accelerator version
160	!
161	! 1241 2013-10-30 11:36:58Z heinze
162	! usage of nudging enabled (so far not implemented for accelerator version)
163	!
164	! 1179 2013-06-14 05:57:58Z raasch
165	! two arguments removed from routine buoyancy, ref_state updated on device
166	!
167	! 1128 2013-04-12 06:19:32Z raasch
168	! those parts requiring global communication moved to time_integration,
169	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
170	! j_north
171	!
172	! 1115 2013-03-26 18:16:16Z hoffmann
173	! optimized cloud physics: calculation of microphysical tendencies transfered
174	! to microphysics.f90; qr and nr are only calculated if precipitation is required
175	!
176	! 1111 2013-03-08 23:54:10Z raasch
177	! update directives for prognostic quantities removed
178	!
179	! 1106 2013-03-04 05:31:38Z raasch
180	! small changes in code formatting
181	!
182	! 1092 2013-02-02 11:24:22Z raasch
183	! unused variables removed
184	!
185	! 1053 2012-11-13 17:11:03Z hoffmann
186	! implementation of two new prognostic equations for rain drop concentration (nr)
187	! and rain water content (qr)
188	!
189	! currently, only available for cache loop optimization
190	!
191	! 1036 2012-10-22 13:43:42Z raasch
192	! code put under GPL (PALM 3.9)
193	!
194	! 1019 2012-09-28 06:46:45Z raasch
195	! non-optimized version of prognostic_equations removed
196	!
197	! 1015 2012-09-27 09:23:24Z raasch
198	! new branch prognostic_equations_acc
199	! OpenACC statements added + code changes required for GPU optimization
200	!
201	! 1001 2012-09-13 14:08:46Z raasch
202	! all actions concerning leapfrog- and upstream-spline-scheme removed
203	!
204	! 978 2012-08-09 08:28:32Z fricke
205	! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
206	! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
207	! boundary in case of non-cyclic lateral boundaries
208	! Bugfix: first thread index changes for WS-scheme at the inflow
209	!
210	! 940 2012-07-09 14:31:00Z raasch
211	! temperature equation can be switched off
212	!
213	! Revision 1.1 2000/04/13 14:56:27 schroeter
214	! Initial revision
215	!
216	!
217	! Description:
218	! ------------
219	!> Solving the prognostic equations.
220	!------------------------------------------------------------------------------!
221	MODULE prognostic_equations_mod
222
223
224
225	USE arrays_3d, &
226	ONLY: diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, &
227	diss_l_s, diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q, &
228	diss_s_qr, diss_s_s, diss_s_sa, e, e_p, flux_s_e, flux_s_nr, &
229	flux_s_pt, flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e, &
230	flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, &
231	nr, nr_p, pt, ptdf_x, ptdf_y, pt_init, pt_p, &
232	prho, q, q_init, q_p, qr, qr_p, rdf, &
233	rdf_sc, ref_state, rho_ocean, s, s_init, s_p, sa, sa_init, sa_p,&
234	tend, te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tu_m, tv_m, &
235	tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
236
237	USE control_parameters, &
238	ONLY: call_microphysics_at_all_substeps, cloud_physics, &
239	cloud_top_radiation, constant_diffusion, dp_external, &
240	dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity, &
241	inflow_l, intermediate_timestep_count, &
242	intermediate_timestep_count_max, large_scale_forcing, &
243	large_scale_subsidence, microphysics_seifert, &
244	microphysics_sat_adjust, neutral, nudging, ocean, outflow_l, &
245	outflow_s, passive_scalar, prho_reference, prho_reference, &
246	prho_reference, pt_reference, pt_reference, pt_reference, &
247	scalar_advec, scalar_advec, simulated_time, sloping_surface, &
248	timestep_scheme, tsc, use_subsidence_tendencies, &
249	use_upstream_for_tke, ws_scheme_mom, ws_scheme_sca
250
251	USE cpulog, &
252	ONLY: cpu_log, log_point
253
254	USE eqn_state_seawater_mod, &
255	ONLY: eqn_state_seawater
256
257	USE indices, &
258	ONLY: nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, nzb, &
259	nzt, wall_flags_0
260
261	USE advec_ws, &
262	ONLY: advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
263
264	USE advec_s_bc_mod, &
265	ONLY: advec_s_bc
266
267	USE advec_s_pw_mod, &
268	ONLY: advec_s_pw
269
270	USE advec_s_up_mod, &
271	ONLY: advec_s_up
272
273	USE advec_u_pw_mod, &
274	ONLY: advec_u_pw
275
276	USE advec_u_up_mod, &
277	ONLY: advec_u_up
278
279	USE advec_v_pw_mod, &
280	ONLY: advec_v_pw
281
282	USE advec_v_up_mod, &
283	ONLY: advec_v_up
284
285	USE advec_w_pw_mod, &
286	ONLY: advec_w_pw
287
288	USE advec_w_up_mod, &
289	ONLY: advec_w_up
290
291	USE buoyancy_mod, &
292	ONLY: buoyancy
293
294	USE calc_radiation_mod, &
295	ONLY: calc_radiation
296
297	USE coriolis_mod, &
298	ONLY: coriolis
299
300	USE diffusion_e_mod, &
301	ONLY: diffusion_e
302
303	USE diffusion_s_mod, &
304	ONLY: diffusion_s
305
306	USE diffusion_u_mod, &
307	ONLY: diffusion_u
308
309	USE diffusion_v_mod, &
310	ONLY: diffusion_v
311
312	USE diffusion_w_mod, &
313	ONLY: diffusion_w
314
315	USE kinds
316
317	USE ls_forcing_mod, &
318	ONLY: ls_advec
319
320	USE microphysics_mod, &
321	ONLY: microphysics_control
322
323	USE nudge_mod, &
324	ONLY: nudge
325
326	USE plant_canopy_model_mod, &
327	ONLY: cthf, plant_canopy, pcm_tendency
328
329	USE production_e_mod, &
330	ONLY: production_e
331
332	USE radiation_model_mod, &
333	ONLY: radiation, radiation_tendency, &
334	skip_time_do_radiation
335
336	USE statistics, &
337	ONLY: hom
338
339	USE subsidence_mod, &
340	ONLY: subsidence
341
342	USE user_actions_mod, &
343	ONLY: user_actions
344
345	USE surface_mod, &
346	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
347	surf_usm_v
348
349	USE wind_turbine_model_mod, &
350	ONLY: wind_turbine, wtm_tendencies
351
352
353	PRIVATE
354	PUBLIC prognostic_equations_cache, prognostic_equations_vector
355
356	INTERFACE prognostic_equations_cache
357	MODULE PROCEDURE prognostic_equations_cache
358	END INTERFACE prognostic_equations_cache
359
360	INTERFACE prognostic_equations_vector
361	MODULE PROCEDURE prognostic_equations_vector
362	END INTERFACE prognostic_equations_vector
363
364
365	CONTAINS
366
367
368	!------------------------------------------------------------------------------!
369	! Description:
370	! ------------
371	!> Version with one optimized loop over all equations. It is only allowed to
372	!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
373	!>
374	!> Here the calls of most subroutines are embedded in two DO loops over i and j,
375	!> so communication between CPUs is not allowed (does not make sense) within
376	!> these loops.
377	!>
378	!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
379	!------------------------------------------------------------------------------!
380
381	SUBROUTINE prognostic_equations_cache
382
383
384	IMPLICIT NONE
385
386	INTEGER(iwp) :: i !<
387	INTEGER(iwp) :: i_omp_start !<
388	INTEGER(iwp) :: j !<
389	INTEGER(iwp) :: k !<
390	INTEGER(iwp) :: omp_get_thread_num !<
391	INTEGER(iwp) :: tn = 0 !<
392
393	LOGICAL :: loop_start !<
394
395
396	!
397	!-- Time measurement can only be performed for the whole set of equations
398	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
399
400	!
401	!-- If required, calculate cloud microphysics
402	IF ( cloud_physics .AND. .NOT. microphysics_sat_adjust .AND. &
403	( intermediate_timestep_count == 1 .OR. &
404	call_microphysics_at_all_substeps ) ) &
405	THEN
406	!$OMP PARALLEL private (i,j)
407	!$OMP DO
408	!$OMP PARALLEL PRIVATE (i, j)
409	!$OMP DO
410	DO i = nxlg, nxrg
411	DO j = nysg, nyng
412	CALL microphysics_control( i, j )
413	ENDDO
414	ENDDO
415	!$OMP END PARALLEL
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 |