Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |