Home

Context Navigation

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

Last change on this file since 3372 was 3355, checked in by knoop, 6 years ago
removed calc_radiation.f90 and related cloud_top_radiation namelist parameter (functionality replaced by radiation model)
Property svn:keywords set to `Id`
File size: 91.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |