Home

Context Navigation

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

Last change on this file since 4351 was 4346, checked in by motisi, 5 years ago
Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0
Property svn:keywords set to `Id`
File size: 59.3 KB

Rev	Line
[1873]	1	!> @file prognostic_equations.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1875]	4	!
[2000]	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.
[1875]	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	!
[3655]	17	! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1875]	19	!
	20	! Current revisions:
	21	! ------------------
[4110]	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: prognostic_equations.f90 4346 2019-12-18 11:55:56Z oliver.maas $
[4346]	27	! Introduction of wall_flags_total_0, which currently sets bits based on static
	28	! topography information used in wall_flags_static_0
	29	!
	30	! 4329 2019-12-10 15:46:36Z motisi
[4329]	31	! Renamed wall_flags_0 to wall_flags_static_0
	32	!
	33	! 4182 2019-08-22 15:20:23Z scharf
[4182]	34	! Corrected "Former revisions" section
	35	!
	36	! 4110 2019-07-22 17:05:21Z suehring
[4109]	37	! pass integer flag array to WS scalar advection routine which is now necessary
	38	! as the flags may differ for scalars, e.g. pt can be cyclic while chemical
	39	! species may be non-cyclic. Further, pass boundary flags.
[2156]	40	!
[4110]	41	! 4109 2019-07-22 17:00:34Z suehring
[4079]	42	! Application of monotonic flux limiter for the vertical scalar advection
	43	! up to the topography top (only for the cache-optimized version at the
	44	! moment). Please note, at the moment the limiter is only applied for passive
	45	! scalars.
	46	!
	47	! 4048 2019-06-21 21:00:21Z knoop
[4048]	48	! Moved tcm_prognostic_equations to module_interface
	49	!
	50	! 3987 2019-05-22 09:52:13Z kanani
[3987]	51	! Introduce alternative switch for debug output during timestepping
	52	!
	53	! 3956 2019-05-07 12:32:52Z monakurppa
[3956]	54	! Removed salsa calls.
	55	!
	56	! 3931 2019-04-24 16:34:28Z schwenkel
[3929]	57	! Correct/complete module_interface introduction for chemistry model
	58	!
	59	! 3899 2019-04-16 14:05:27Z monakurppa
[3899]	60	! Corrections in the OpenMP version of salsa
[3929]	61	!
	62	! 3887 2019 -04-12 08:47:41Z schwenkel
[3887]	63	! Implicit Bugfix for chemistry model, loop for non_transport_physics over
	64	! ghost points is avoided. Instead introducing module_interface_exchange_horiz.
	65	!
	66	! 3885 2019-04-11 11:29:34Z kanani
[3885]	67	! Changes related to global restructuring of location messages and introduction
	68	! of additional debug messages
	69	!
	70	! 3881 2019-04-10 09:31:22Z suehring
[3881]	71	! Bugfix in OpenMP directive
	72	!
	73	! 3880 2019-04-08 21:43:02Z knoop
[3875]	74	! Moved wtm_tendencies to module_interface_actions
	75	!
	76	! 3874 2019-04-08 16:53:48Z knoop
[3874]	77	! Added non_transport_physics module interfaces and moved bcm code into it
	78	!
	79	! 3872 2019-04-08 15:03:06Z knoop
[3870]	80	! Moving prognostic equations of bcm into bulk_cloud_model_mod
	81	!
	82	! 3864 2019-04-05 09:01:56Z monakurppa
[3864]	83	! Modifications made for salsa:
	84	! - salsa_prognostic_equations moved to salsa_mod (and the call to
	85	! module_interface_mod)
	86	! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
	87	! ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
	88	!
	89	! 3840 2019-03-29 10:35:52Z knoop
[3879]	90	! added USE chem_gasphase_mod for nspec, nspec and spc_names
[3833]	91	!
	92	! 3820 2019-03-27 11:53:41Z forkel
[3820]	93	! renamed do_depo to deposition_dry (ecc)
	94	!
	95	! 3797 2019-03-15 11:15:38Z forkel
[3797]	96	! Call chem_integegrate in OpenMP loop (ketelsen)
	97	!
	98	!
	99	! 3771 2019-02-28 12:19:33Z raasch
[3771]	100	! preprocessor directivs fro rrtmg added
	101	!
	102	! 3761 2019-02-25 15:31:42Z raasch
[3761]	103	! unused variable removed
	104	!
	105	! 3719 2019-02-06 13:10:18Z kanani
[3719]	106	! Cleaned up chemistry cpu measurements
	107	!
	108	! 3684 2019-01-20 20:20:58Z knoop
[3634]	109	! OpenACC port for SPEC
[3458]	110	!
[4182]	111	! Revision 1.1 2000/04/13 14:56:27 schroeter
	112	! Initial revision
	113	!
	114	!
[1875]	115	! Description:
	116	! ------------
	117	!> Solving the prognostic equations.
	118	!------------------------------------------------------------------------------!
	119	MODULE prognostic_equations_mod
	120
[3294]	121	USE advec_s_bc_mod, &
	122	ONLY: advec_s_bc
[2155]	123
[3294]	124	USE advec_s_pw_mod, &
	125	ONLY: advec_s_pw
	126
	127	USE advec_s_up_mod, &
	128	ONLY: advec_s_up
	129
	130	USE advec_u_pw_mod, &
	131	ONLY: advec_u_pw
	132
	133	USE advec_u_up_mod, &
	134	ONLY: advec_u_up
	135
	136	USE advec_v_pw_mod, &
	137	ONLY: advec_v_pw
	138
	139	USE advec_v_up_mod, &
	140	ONLY: advec_v_up
	141
	142	USE advec_w_pw_mod, &
	143	ONLY: advec_w_pw
	144
	145	USE advec_w_up_mod, &
	146	ONLY: advec_w_up
	147
	148	USE advec_ws, &
	149	ONLY: advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
	150
[1875]	151	USE arrays_3d, &
[3870]	152	ONLY: diss_l_e, diss_l_pt, diss_l_q, &
	153	diss_l_s, diss_l_sa, diss_s_e, &
	154	diss_s_pt, diss_s_q, diss_s_s, diss_s_sa, &
	155	e, e_p, flux_s_e, flux_s_pt, flux_s_q, &
	156	flux_s_s, flux_s_sa, flux_l_e, &
	157	flux_l_pt, flux_l_q, flux_l_s, &
	158	flux_l_sa, pt, ptdf_x, ptdf_y, pt_init, &
	159	pt_p, prho, q, q_init, q_p, rdf, rdf_sc, &
	160	ref_state, rho_ocean, s, s_init, s_p, tend, te_m, &
	161	tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u, &
[3294]	162	ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
[2155]	163
[3294]	164	USE buoyancy_mod, &
	165	ONLY: buoyancy
[3864]	166
[1875]	167	USE control_parameters, &
[4109]	168	ONLY: bc_dirichlet_l, &
	169	bc_dirichlet_n, &
	170	bc_dirichlet_r, &
	171	bc_dirichlet_s, &
	172	bc_radiation_l, &
	173	bc_radiation_n, &
	174	bc_radiation_r, &
	175	bc_radiation_s, &
	176	constant_diffusion, &
[3987]	177	debug_output_timestep, &
[2696]	178	dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, &
[3182]	179	humidity, intermediate_timestep_count, &
[1875]	180	intermediate_timestep_count_max, large_scale_forcing, &
[4079]	181	large_scale_subsidence, &
	182	monotonic_limiter_z, &
	183	neutral, nudging, &
[3294]	184	ocean_mode, passive_scalar, plant_canopy, pt_reference, &
[1875]	185	scalar_advec, scalar_advec, simulated_time, sloping_surface, &
[2232]	186	timestep_scheme, tsc, use_subsidence_tendencies, &
[2563]	187	use_upstream_for_tke, wind_turbine, ws_scheme_mom, &
[3467]	188	ws_scheme_sca, urban_surface, land_surface, &
[3582]	189	time_since_reference_point, salsa
[1875]	190
[3294]	191	USE coriolis_mod, &
	192	ONLY: coriolis
	193
[1875]	194	USE cpulog, &
[2696]	195	ONLY: cpu_log, log_point, log_point_s
[1875]	196
	197	USE diffusion_s_mod, &
[2118]	198	ONLY: diffusion_s
[1875]	199
	200	USE diffusion_u_mod, &
[2118]	201	ONLY: diffusion_u
[1875]	202
	203	USE diffusion_v_mod, &
[2118]	204	ONLY: diffusion_v
[1875]	205
	206	USE diffusion_w_mod, &
[2118]	207	ONLY: diffusion_w
[1875]	208
[3294]	209	USE indices, &
[4109]	210	ONLY: advc_flags_s, &
	211	nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
[4346]	212	nzb, nzt, wall_flags_total_0
[3294]	213
[1875]	214	USE kinds
	215
[2320]	216	USE lsf_nudging_mod, &
	217	ONLY: ls_advec, nudge
[1875]	218
[3684]	219	USE module_interface, &
[3837]	220	ONLY: module_interface_actions, &
[3931]	221	module_interface_non_advective_processes, &
[3887]	222	module_interface_exchange_horiz, &
[3837]	223	module_interface_prognostic_equations
[3684]	224
[3294]	225	USE ocean_mod, &
[3840]	226	ONLY: stokes_drift_terms, stokes_force, &
[3302]	227	wave_breaking, wave_breaking_term
[3294]	228
[1875]	229	USE plant_canopy_model_mod, &
[2746]	230	ONLY: cthf, pcm_tendency
[1875]	231
[3771]	232	#if defined( __rrtmg )
[1875]	233	USE radiation_model_mod, &
[1976]	234	ONLY: radiation, radiation_tendency, &
[1875]	235	skip_time_do_radiation
[3771]	236	#endif
[3864]	237
[1875]	238	USE statistics, &
	239	ONLY: hom
	240
	241	USE subsidence_mod, &
	242	ONLY: subsidence
	243
[3294]	244	USE surface_mod, &
	245	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	246	surf_usm_v
	247
[3874]	248	IMPLICIT NONE
[1914]	249
[1875]	250	PRIVATE
[2118]	251	PUBLIC prognostic_equations_cache, prognostic_equations_vector
[1875]	252
	253	INTERFACE prognostic_equations_cache
	254	MODULE PROCEDURE prognostic_equations_cache
	255	END INTERFACE prognostic_equations_cache
	256
	257	INTERFACE prognostic_equations_vector
	258	MODULE PROCEDURE prognostic_equations_vector
	259	END INTERFACE prognostic_equations_vector
	260
	261
	262	CONTAINS
	263
	264
	265	!------------------------------------------------------------------------------!
	266	! Description:
	267	! ------------
	268	!> Version with one optimized loop over all equations. It is only allowed to
	269	!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
	270	!>
	271	!> Here the calls of most subroutines are embedded in two DO loops over i and j,
	272	!> so communication between CPUs is not allowed (does not make sense) within
	273	!> these loops.
	274	!>
	275	!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
	276	!------------------------------------------------------------------------------!
[2155]	277
[1875]	278	SUBROUTINE prognostic_equations_cache
	279
	280
	281	INTEGER(iwp) :: i !<
	282	INTEGER(iwp) :: i_omp_start !<
	283	INTEGER(iwp) :: j !<
	284	INTEGER(iwp) :: k !<
[3241]	285	!$ INTEGER(iwp) :: omp_get_thread_num !<
[1875]	286	INTEGER(iwp) :: tn = 0 !<
[2155]	287
[1875]	288	LOGICAL :: loop_start !<
	289
	290
[3987]	291	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_cache', 'start' )
[1875]	292	!
	293	!-- Time measurement can only be performed for the whole set of equations
	294	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
	295
[3878]	296	!$OMP PARALLEL PRIVATE (i,j)
	297	!$OMP DO
[3887]	298	DO i = nxl, nxr
	299	DO j = nys, nyn
[1875]	300	!
[3931]	301	!-- Calculate non advective processes for all other modules
	302	CALL module_interface_non_advective_processes( i, j )
[3878]	303	ENDDO
	304	ENDDO
[3887]	305	!
[3931]	306	!-- Module Inferface for exchange horiz after non_advective_processes but before
[3956]	307	!-- advection. Therefore, non_advective_processes must not run for ghost points.
	308	!$OMP END PARALLEL
[3887]	309	CALL module_interface_exchange_horiz()
[2696]	310	!
[2192]	311	!-- Loop over all prognostic equations
[3881]	312	!$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn)
[2192]	313
	314	!$ tn = omp_get_thread_num()
	315	loop_start = .TRUE.
	316
	317	!$OMP DO
[1875]	318	DO i = nxl, nxr
	319
	320	!
	321	!-- Store the first loop index. It differs for each thread and is required
	322	!-- later in advec_ws
	323	IF ( loop_start ) THEN
	324	loop_start = .FALSE.
[2155]	325	i_omp_start = i
[1875]	326	ENDIF
	327
	328	DO j = nys, nyn
	329	!
[3022]	330	!-- Tendency terms for u-velocity component. Please note, in case of
	331	!-- non-cyclic boundary conditions the grid point i=0 is excluded from
[3899]	332	!-- the prognostic equations for the u-component.
[3022]	333	IF ( i >= nxlu ) THEN
[1875]	334
	335	tend(:,j,i) = 0.0_wp
	336	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	337	IF ( ws_scheme_mom ) THEN
	338	CALL advec_u_ws( i, j, i_omp_start, tn )
[2155]	339	ELSE
[1875]	340	CALL advec_u_pw( i, j )
[2155]	341	ENDIF
[1875]	342	ELSE
	343	CALL advec_u_up( i, j )
	344	ENDIF
	345	CALL diffusion_u( i, j )
	346	CALL coriolis( i, j, 1 )
	347	IF ( sloping_surface .AND. .NOT. neutral ) THEN
	348	CALL buoyancy( i, j, pt, 1 )
	349	ENDIF
	350
	351	!
	352	!-- Drag by plant canopy
	353	IF ( plant_canopy ) CALL pcm_tendency( i, j, 1 )
	354
	355	!
	356	!-- External pressure gradient
	357	IF ( dp_external ) THEN
	358	DO k = dp_level_ind_b+1, nzt
	359	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	360	ENDDO
	361	ENDIF
	362
	363	!
	364	!-- Nudging
	365	IF ( nudging ) CALL nudge( i, j, simulated_time, 'u' )
	366
[1914]	367	!
[3302]	368	!-- Effect of Stokes drift (in ocean mode only)
	369	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 1 )
	370
[3684]	371	CALL module_interface_actions( i, j, 'u-tendency' )
[1875]	372	!
	373	!-- Prognostic equation for u-velocity component
[2232]	374	DO k = nzb+1, nzt
	375
	376	u_p(k,j,i) = u(k,j,i) + ( dt_3d * &
	377	( tsc(2) * tend(k,j,i) + &
	378	tsc(3) * tu_m(k,j,i) ) &
	379	- tsc(5) * rdf(k) &
	380	* ( u(k,j,i) - u_init(k) ) &
	381	) * MERGE( 1.0_wp, 0.0_wp, &
[4346]	382	BTEST( wall_flags_total_0(k,j,i), 1 )&
[2232]	383	)
[1875]	384	ENDDO
	385
	386	!
[3302]	387	!-- Add turbulence generated by wave breaking (in ocean mode only)
	388	IF ( wave_breaking .AND. &
	389	intermediate_timestep_count == intermediate_timestep_count_max )&
	390	THEN
	391	CALL wave_breaking_term( i, j, 1 )
	392	ENDIF
	393
	394	!
[1875]	395	!-- Calculate tendencies for the next Runge-Kutta step
	396	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	397	IF ( intermediate_timestep_count == 1 ) THEN
[2232]	398	DO k = nzb+1, nzt
[1875]	399	tu_m(k,j,i) = tend(k,j,i)
	400	ENDDO
	401	ELSEIF ( intermediate_timestep_count < &
	402	intermediate_timestep_count_max ) THEN
[2232]	403	DO k = nzb+1, nzt
	404	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
	405	+ 5.3125_wp * tu_m(k,j,i)
[1875]	406	ENDDO
	407	ENDIF
	408	ENDIF
	409
	410	ENDIF
	411	!
[3022]	412	!-- Tendency terms for v-velocity component. Please note, in case of
	413	!-- non-cyclic boundary conditions the grid point j=0 is excluded from
	414	!-- the prognostic equations for the v-component. !--
	415	IF ( j >= nysv ) THEN
[1875]	416
	417	tend(:,j,i) = 0.0_wp
	418	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	419	IF ( ws_scheme_mom ) THEN
	420	CALL advec_v_ws( i, j, i_omp_start, tn )
[2155]	421	ELSE
[1875]	422	CALL advec_v_pw( i, j )
	423	ENDIF
	424	ELSE
	425	CALL advec_v_up( i, j )
	426	ENDIF
	427	CALL diffusion_v( i, j )
	428	CALL coriolis( i, j, 2 )
	429
	430	!
	431	!-- Drag by plant canopy
[2155]	432	IF ( plant_canopy ) CALL pcm_tendency( i, j, 2 )
[1875]	433
	434	!
	435	!-- External pressure gradient
	436	IF ( dp_external ) THEN
	437	DO k = dp_level_ind_b+1, nzt
	438	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	439	ENDDO
	440	ENDIF
	441
	442	!
	443	!-- Nudging
	444	IF ( nudging ) CALL nudge( i, j, simulated_time, 'v' )
	445
[1914]	446	!
[3302]	447	!-- Effect of Stokes drift (in ocean mode only)
	448	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 2 )
	449
[3684]	450	CALL module_interface_actions( i, j, 'v-tendency' )
[1875]	451	!
	452	!-- Prognostic equation for v-velocity component
[2232]	453	DO k = nzb+1, nzt
	454	v_p(k,j,i) = v(k,j,i) + ( dt_3d * &
	455	( tsc(2) * tend(k,j,i) + &
	456	tsc(3) * tv_m(k,j,i) ) &
	457	- tsc(5) * rdf(k) &
	458	* ( v(k,j,i) - v_init(k) )&
	459	) * MERGE( 1.0_wp, 0.0_wp, &
[4346]	460	BTEST( wall_flags_total_0(k,j,i), 2 )&
[2232]	461	)
[1875]	462	ENDDO
	463
	464	!
[3302]	465	!-- Add turbulence generated by wave breaking (in ocean mode only)
	466	IF ( wave_breaking .AND. &
	467	intermediate_timestep_count == intermediate_timestep_count_max )&
	468	THEN
	469	CALL wave_breaking_term( i, j, 2 )
	470	ENDIF
	471
	472	!
[1875]	473	!-- Calculate tendencies for the next Runge-Kutta step
	474	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	475	IF ( intermediate_timestep_count == 1 ) THEN
[2232]	476	DO k = nzb+1, nzt
[1875]	477	tv_m(k,j,i) = tend(k,j,i)
	478	ENDDO
	479	ELSEIF ( intermediate_timestep_count < &
	480	intermediate_timestep_count_max ) THEN
[2232]	481	DO k = nzb+1, nzt
	482	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
	483	+ 5.3125_wp * tv_m(k,j,i)
[1875]	484	ENDDO
	485	ENDIF
	486	ENDIF
	487
	488	ENDIF
	489
	490	!
	491	!-- Tendency terms for w-velocity component
	492	tend(:,j,i) = 0.0_wp
	493	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	494	IF ( ws_scheme_mom ) THEN
	495	CALL advec_w_ws( i, j, i_omp_start, tn )
[2155]	496	ELSE
[1875]	497	CALL advec_w_pw( i, j )
	498	END IF
	499	ELSE
	500	CALL advec_w_up( i, j )
	501	ENDIF
	502	CALL diffusion_w( i, j )
	503	CALL coriolis( i, j, 3 )
	504
	505	IF ( .NOT. neutral ) THEN
[3294]	506	IF ( ocean_mode ) THEN
[2031]	507	CALL buoyancy( i, j, rho_ocean, 3 )
[1875]	508	ELSE
	509	IF ( .NOT. humidity ) THEN
	510	CALL buoyancy( i, j, pt, 3 )
	511	ELSE
	512	CALL buoyancy( i, j, vpt, 3 )
	513	ENDIF
	514	ENDIF
	515	ENDIF
	516
	517	!
	518	!-- Drag by plant canopy
	519	IF ( plant_canopy ) CALL pcm_tendency( i, j, 3 )
	520
[1914]	521	!
[3302]	522	!-- Effect of Stokes drift (in ocean mode only)
	523	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 3 )
	524
[3684]	525	CALL module_interface_actions( i, j, 'w-tendency' )
[1875]	526	!
	527	!-- Prognostic equation for w-velocity component
[2232]	528	DO k = nzb+1, nzt-1
	529	w_p(k,j,i) = w(k,j,i) + ( dt_3d * &
	530	( tsc(2) * tend(k,j,i) + &
[1875]	531	tsc(3) * tw_m(k,j,i) ) &
[2232]	532	- tsc(5) * rdf(k) * w(k,j,i) &
	533	) * MERGE( 1.0_wp, 0.0_wp, &
[4346]	534	BTEST( wall_flags_total_0(k,j,i), 3 )&
[2232]	535	)
[1875]	536	ENDDO
	537
	538	!
	539	!-- Calculate tendencies for the next Runge-Kutta step
	540	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	541	IF ( intermediate_timestep_count == 1 ) THEN
[2232]	542	DO k = nzb+1, nzt-1
[1875]	543	tw_m(k,j,i) = tend(k,j,i)
	544	ENDDO
	545	ELSEIF ( intermediate_timestep_count < &
	546	intermediate_timestep_count_max ) THEN
[2232]	547	DO k = nzb+1, nzt-1
	548	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
	549	+ 5.3125_wp * tw_m(k,j,i)
[1875]	550	ENDDO
	551	ENDIF
	552	ENDIF
	553
	554	!
	555	!-- If required, compute prognostic equation for potential temperature
	556	IF ( .NOT. neutral ) THEN
	557	!
	558	!-- Tendency terms for potential temperature
	559	tend(:,j,i) = 0.0_wp
	560	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	561	IF ( ws_scheme_sca ) THEN
[4109]	562	CALL advec_s_ws( advc_flags_s, &
	563	i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
	564	flux_l_pt, diss_l_pt, i_omp_start, tn, &
	565	bc_dirichlet_l .OR. bc_radiation_l, &
	566	bc_dirichlet_n .OR. bc_radiation_n, &
	567	bc_dirichlet_r .OR. bc_radiation_r, &
	568	bc_dirichlet_s .OR. bc_radiation_s )
[1875]	569	ELSE
	570	CALL advec_s_pw( i, j, pt )
	571	ENDIF
	572	ELSE
	573	CALL advec_s_up( i, j, pt )
	574	ENDIF
[2232]	575	CALL diffusion_s( i, j, pt, &
	576	surf_def_h(0)%shf, surf_def_h(1)%shf, &
	577	surf_def_h(2)%shf, &
	578	surf_lsm_h%shf, surf_usm_h%shf, &
	579	surf_def_v(0)%shf, surf_def_v(1)%shf, &
	580	surf_def_v(2)%shf, surf_def_v(3)%shf, &
	581	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
	582	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
	583	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
	584	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
[1875]	585
	586	!
	587	!-- Consideration of heat sources within the plant canopy
[3014]	588	IF ( plant_canopy .AND. &
	589	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
[1875]	590	CALL pcm_tendency( i, j, 4 )
	591	ENDIF
	592
	593	!
	594	!-- Large scale advection
	595	IF ( large_scale_forcing ) THEN
	596	CALL ls_advec( i, j, simulated_time, 'pt' )
[2155]	597	ENDIF
[1875]	598
	599	!
	600	!-- Nudging
[2155]	601	IF ( nudging ) CALL nudge( i, j, simulated_time, 'pt' )
[1875]	602
	603	!
	604	!-- If required, compute effect of large-scale subsidence/ascent
	605	IF ( large_scale_subsidence .AND. &
	606	.NOT. use_subsidence_tendencies ) THEN
	607	CALL subsidence( i, j, tend, pt, pt_init, 2 )
	608	ENDIF
	609
[3771]	610	#if defined( __rrtmg )
[1875]	611	!
	612	!-- If required, add tendency due to radiative heating/cooling
[1976]	613	IF ( radiation .AND. &
[1875]	614	simulated_time > skip_time_do_radiation ) THEN
	615	CALL radiation_tendency ( i, j, tend )
	616	ENDIF
[3771]	617	#endif
[1875]	618
[3684]	619	CALL module_interface_actions( i, j, 'pt-tendency' )
[1875]	620	!
	621	!-- Prognostic equation for potential temperature
[2232]	622	DO k = nzb+1, nzt
	623	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * &
	624	( tsc(2) * tend(k,j,i) + &
[1875]	625	tsc(3) * tpt_m(k,j,i) ) &
[2232]	626	- tsc(5) &
	627	* ( pt(k,j,i) - pt_init(k) ) &
	628	* ( rdf_sc(k) + ptdf_x(i) &
	629	+ ptdf_y(j) ) &
	630	) &
	631	* MERGE( 1.0_wp, 0.0_wp, &
[4346]	632	BTEST( wall_flags_total_0(k,j,i), 0 )&
[2232]	633	)
[1875]	634	ENDDO
	635
	636	!
	637	!-- Calculate tendencies for the next Runge-Kutta step
	638	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	639	IF ( intermediate_timestep_count == 1 ) THEN
[2232]	640	DO k = nzb+1, nzt
[1875]	641	tpt_m(k,j,i) = tend(k,j,i)
	642	ENDDO
	643	ELSEIF ( intermediate_timestep_count < &
	644	intermediate_timestep_count_max ) THEN
[2232]	645	DO k = nzb+1, nzt
	646	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	647	5.3125_wp * tpt_m(k,j,i)
[1875]	648	ENDDO
	649	ENDIF
	650	ENDIF
	651
	652	ENDIF
	653
	654	!
[1960]	655	!-- If required, compute prognostic equation for total water content
	656	IF ( humidity ) THEN
[1875]	657
	658	!
	659	!-- Tendency-terms for total water content / scalar
	660	tend(:,j,i) = 0.0_wp
[4109]	661	IF ( timestep_scheme(1:5) == 'runge' ) &
[1875]	662	THEN
	663	IF ( ws_scheme_sca ) THEN
[4109]	664	CALL advec_s_ws( advc_flags_s, &
	665	i, j, q, 'q', flux_s_q, &
	666	diss_s_q, flux_l_q, diss_l_q, &
	667	i_omp_start, tn, &
	668	bc_dirichlet_l .OR. bc_radiation_l, &
	669	bc_dirichlet_n .OR. bc_radiation_n, &
	670	bc_dirichlet_r .OR. bc_radiation_r, &
	671	bc_dirichlet_s .OR. bc_radiation_s )
[2155]	672	ELSE
[1875]	673	CALL advec_s_pw( i, j, q )
	674	ENDIF
	675	ELSE
	676	CALL advec_s_up( i, j, q )
	677	ENDIF
[2232]	678	CALL diffusion_s( i, j, q, &
	679	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
	680	surf_def_h(2)%qsws, &
	681	surf_lsm_h%qsws, surf_usm_h%qsws, &
	682	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
	683	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
	684	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
	685	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
	686	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
	687	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
[1875]	688
	689	!
[1960]	690	!-- Sink or source of humidity due to canopy elements
[1875]	691	IF ( plant_canopy ) CALL pcm_tendency( i, j, 5 )
	692
	693	!
	694	!-- Large scale advection
	695	IF ( large_scale_forcing ) THEN
	696	CALL ls_advec( i, j, simulated_time, 'q' )
	697	ENDIF
	698
	699	!
	700	!-- Nudging
[2155]	701	IF ( nudging ) CALL nudge( i, j, simulated_time, 'q' )
[1875]	702
	703	!
	704	!-- If required compute influence of large-scale subsidence/ascent
	705	IF ( large_scale_subsidence .AND. &
	706	.NOT. use_subsidence_tendencies ) THEN
	707	CALL subsidence( i, j, tend, q, q_init, 3 )
	708	ENDIF
	709
[3684]	710	CALL module_interface_actions( i, j, 'q-tendency' )
[1875]	711
	712	!
	713	!-- Prognostic equation for total water content / scalar
[2232]	714	DO k = nzb+1, nzt
	715	q_p(k,j,i) = q(k,j,i) + ( dt_3d * &
	716	( tsc(2) * tend(k,j,i) + &
[1875]	717	tsc(3) * tq_m(k,j,i) ) &
[2232]	718	- tsc(5) * rdf_sc(k) * &
	719	( q(k,j,i) - q_init(k) ) &
	720	) &
	721	* MERGE( 1.0_wp, 0.0_wp, &
[4346]	722	BTEST( wall_flags_total_0(k,j,i), 0 )&
[2232]	723	)
[1875]	724	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
	725	ENDDO
	726
	727	!
	728	!-- Calculate tendencies for the next Runge-Kutta step
	729	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	730	IF ( intermediate_timestep_count == 1 ) THEN
[2232]	731	DO k = nzb+1, nzt
[1875]	732	tq_m(k,j,i) = tend(k,j,i)
	733	ENDDO
	734	ELSEIF ( intermediate_timestep_count < &
	735	intermediate_timestep_count_max ) THEN
[2232]	736	DO k = nzb+1, nzt
	737	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	738	5.3125_wp * tq_m(k,j,i)
[1875]	739	ENDDO
	740	ENDIF
	741	ENDIF
	742
	743	ENDIF
[2155]	744
[1960]	745	!
	746	!-- If required, compute prognostic equation for scalar
	747	IF ( passive_scalar ) THEN
	748	!
	749	!-- Tendency-terms for total water content / scalar
	750	tend(:,j,i) = 0.0_wp
[4109]	751	IF ( timestep_scheme(1:5) == 'runge' ) &
[1960]	752	THEN
	753	IF ( ws_scheme_sca ) THEN
[4079]	754	!
	755	!-- For scalar advection apply monotonic flux limiter near
	756	!-- topography.
[4109]	757	CALL advec_s_ws( advc_flags_s, &
	758	i, j, s, 's', flux_s_s, &
[4079]	759	diss_s_s, flux_l_s, diss_l_s, i_omp_start, &
[4109]	760	tn, &
	761	bc_dirichlet_l .OR. bc_radiation_l, &
	762	bc_dirichlet_n .OR. bc_radiation_n, &
	763	bc_dirichlet_r .OR. bc_radiation_r, &
	764	bc_dirichlet_s .OR. bc_radiation_s, &
	765	monotonic_limiter_z )
[2155]	766	ELSE
[1960]	767	CALL advec_s_pw( i, j, s )
	768	ENDIF
	769	ELSE
	770	CALL advec_s_up( i, j, s )
	771	ENDIF
[2232]	772	CALL diffusion_s( i, j, s, &
	773	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
	774	surf_def_h(2)%ssws, &
	775	surf_lsm_h%ssws, surf_usm_h%ssws, &
	776	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
	777	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
	778	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
	779	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
	780	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
	781	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
[1875]	782
	783	!
[1960]	784	!-- Sink or source of scalar concentration due to canopy elements
	785	IF ( plant_canopy ) CALL pcm_tendency( i, j, 7 )
	786
	787	!
	788	!-- Large scale advection, still need to be extended for scalars
	789	! IF ( large_scale_forcing ) THEN
	790	! CALL ls_advec( i, j, simulated_time, 's' )
	791	! ENDIF
	792
	793	!
	794	!-- Nudging, still need to be extended for scalars
[2155]	795	! IF ( nudging ) CALL nudge( i, j, simulated_time, 's' )
[1960]	796
	797	!
	798	!-- If required compute influence of large-scale subsidence/ascent.
[2155]	799	!-- Note, the last argument is of no meaning in this case, as it is
	800	!-- only used in conjunction with large_scale_forcing, which is to
[1960]	801	!-- date not implemented for scalars.
	802	IF ( large_scale_subsidence .AND. &
	803	.NOT. use_subsidence_tendencies .AND. &
	804	.NOT. large_scale_forcing ) THEN
	805	CALL subsidence( i, j, tend, s, s_init, 3 )
	806	ENDIF
	807
[3684]	808	CALL module_interface_actions( i, j, 's-tendency' )
[1960]	809
	810	!
	811	!-- Prognostic equation for scalar
[2232]	812	DO k = nzb+1, nzt
	813	s_p(k,j,i) = s(k,j,i) + ( dt_3d * &
	814	( tsc(2) * tend(k,j,i) + &
	815	tsc(3) * ts_m(k,j,i) ) &
	816	- tsc(5) * rdf_sc(k) &
	817	* ( s(k,j,i) - s_init(k) )&
	818	) &
	819	* MERGE( 1.0_wp, 0.0_wp, &
[4346]	820	BTEST( wall_flags_total_0(k,j,i), 0 )&
[2232]	821	)
[1960]	822	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
	823	ENDDO
	824
	825	!
	826	!-- Calculate tendencies for the next Runge-Kutta step
	827	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	828	IF ( intermediate_timestep_count == 1 ) THEN
[2232]	829	DO k = nzb+1, nzt
[1960]	830	ts_m(k,j,i) = tend(k,j,i)
	831	ENDDO
	832	ELSEIF ( intermediate_timestep_count < &
	833	intermediate_timestep_count_max ) THEN
[2232]	834	DO k = nzb+1, nzt
	835	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	836	5.3125_wp * ts_m(k,j,i)
[1960]	837	ENDDO
	838	ENDIF
	839	ENDIF
	840
[2155]	841	ENDIF
[1960]	842	!
[3837]	843	!-- Calculate prognostic equations for all other modules
	844	CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
[1875]	845
[3294]	846	ENDDO ! loop over j
	847	ENDDO ! loop over i
[2192]	848	!$OMP END PARALLEL
[1875]	849
[2232]	850
[1875]	851	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
	852
[3987]	853	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_cache', 'end' )
[1875]	854
	855	END SUBROUTINE prognostic_equations_cache
	856
	857
	858	!------------------------------------------------------------------------------!
	859	! Description:
	860	! ------------
	861	!> Version for vector machines
	862	!------------------------------------------------------------------------------!
[2155]	863
[1875]	864	SUBROUTINE prognostic_equations_vector
	865
	866
[2815]	867	INTEGER(iwp) :: i !<
	868	INTEGER(iwp) :: j !<
	869	INTEGER(iwp) :: k !<
[1875]	870
	871	REAL(wp) :: sbt !<
	872
[3885]	873
[3987]	874	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_vector', 'start' )
[3467]	875	!
[3931]	876	!-- Calculate non advective processes for all other modules
	877	CALL module_interface_non_advective_processes
[1875]	878	!
[3931]	879	!-- Module Inferface for exchange horiz after non_advective_processes but before
[3956]	880	!-- advection. Therefore, non_advective_processes must not run for ghost points.
[3887]	881	CALL module_interface_exchange_horiz()
	882	!
[1875]	883	!-- u-velocity component
	884	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	885
[3634]	886	!$ACC KERNELS PRESENT(tend)
[1875]	887	tend = 0.0_wp
[3634]	888	!$ACC END KERNELS
[1875]	889	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	890	IF ( ws_scheme_mom ) THEN
	891	CALL advec_u_ws
	892	ELSE
	893	CALL advec_u_pw
	894	ENDIF
	895	ELSE
	896	CALL advec_u_up
	897	ENDIF
	898	CALL diffusion_u
	899	CALL coriolis( 1 )
	900	IF ( sloping_surface .AND. .NOT. neutral ) THEN
	901	CALL buoyancy( pt, 1 )
	902	ENDIF
	903
	904	!
	905	!-- Drag by plant canopy
	906	IF ( plant_canopy ) CALL pcm_tendency( 1 )
	907
	908	!
	909	!-- External pressure gradient
	910	IF ( dp_external ) THEN
	911	DO i = nxlu, nxr
	912	DO j = nys, nyn
	913	DO k = dp_level_ind_b+1, nzt
	914	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	915	ENDDO
	916	ENDDO
	917	ENDDO
	918	ENDIF
	919
	920	!
	921	!-- Nudging
	922	IF ( nudging ) CALL nudge( simulated_time, 'u' )
	923
[1914]	924	!
[3302]	925	!-- Effect of Stokes drift (in ocean mode only)
	926	IF ( stokes_force ) CALL stokes_drift_terms( 1 )
	927
[3684]	928	CALL module_interface_actions( 'u-tendency' )
[1875]	929
	930	!
	931	!-- Prognostic equation for u-velocity component
[3634]	932	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
[4346]	933	!$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_total_0) &
[3634]	934	!$ACC PRESENT(tsc(2:5)) &
	935	!$ACC PRESENT(u_p)
[1875]	936	DO i = nxlu, nxr
	937	DO j = nys, nyn
[2232]	938	DO k = nzb+1, nzt
	939	u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
	940	tsc(3) * tu_m(k,j,i) ) &
	941	- tsc(5) * rdf(k) * &
	942	( u(k,j,i) - u_init(k) ) &
	943	) * MERGE( 1.0_wp, 0.0_wp, &
[4346]	944	BTEST( wall_flags_total_0(k,j,i), 1 ) &
[2232]	945	)
[1875]	946	ENDDO
	947	ENDDO
	948	ENDDO
	949
	950	!
[3302]	951	!-- Add turbulence generated by wave breaking (in ocean mode only)
	952	IF ( wave_breaking .AND. &
	953	intermediate_timestep_count == intermediate_timestep_count_max ) &
	954	THEN
	955	CALL wave_breaking_term( 1 )
	956	ENDIF
	957
	958	!
[1875]	959	!-- Calculate tendencies for the next Runge-Kutta step
	960	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	961	IF ( intermediate_timestep_count == 1 ) THEN
[3634]	962	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	963	!$ACC PRESENT(tend, tu_m)
[1875]	964	DO i = nxlu, nxr
	965	DO j = nys, nyn
[2232]	966	DO k = nzb+1, nzt
[1875]	967	tu_m(k,j,i) = tend(k,j,i)
	968	ENDDO
	969	ENDDO
	970	ENDDO
	971	ELSEIF ( intermediate_timestep_count < &
	972	intermediate_timestep_count_max ) THEN
[3634]	973	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	974	!$ACC PRESENT(tend, tu_m)
[1875]	975	DO i = nxlu, nxr
	976	DO j = nys, nyn
[2232]	977	DO k = nzb+1, nzt
	978	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
	979	+ 5.3125_wp * tu_m(k,j,i)
[1875]	980	ENDDO
	981	ENDDO
	982	ENDDO
	983	ENDIF
	984	ENDIF
	985
	986	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	987
	988	!
	989	!-- v-velocity component
	990	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	991
[3634]	992	!$ACC KERNELS PRESENT(tend)
[1875]	993	tend = 0.0_wp
[3634]	994	!$ACC END KERNELS
[1875]	995	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	996	IF ( ws_scheme_mom ) THEN
	997	CALL advec_v_ws
[2155]	998	ELSE
[1875]	999	CALL advec_v_pw
	1000	END IF
	1001	ELSE
	1002	CALL advec_v_up
	1003	ENDIF
	1004	CALL diffusion_v
	1005	CALL coriolis( 2 )
	1006
	1007	!
	1008	!-- Drag by plant canopy
	1009	IF ( plant_canopy ) CALL pcm_tendency( 2 )
	1010
	1011	!
	1012	!-- External pressure gradient
	1013	IF ( dp_external ) THEN
	1014	DO i = nxl, nxr
	1015	DO j = nysv, nyn
	1016	DO k = dp_level_ind_b+1, nzt
	1017	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1018	ENDDO
	1019	ENDDO
	1020	ENDDO
	1021	ENDIF
	1022
	1023	!
	1024	!-- Nudging
	1025	IF ( nudging ) CALL nudge( simulated_time, 'v' )
	1026
[1914]	1027	!
[3302]	1028	!-- Effect of Stokes drift (in ocean mode only)
	1029	IF ( stokes_force ) CALL stokes_drift_terms( 2 )
	1030
[3684]	1031	CALL module_interface_actions( 'v-tendency' )
[1875]	1032
	1033	!
	1034	!-- Prognostic equation for v-velocity component
[3634]	1035	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
[4346]	1036	!$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_total_0) &
[3634]	1037	!$ACC PRESENT(tsc(2:5)) &
	1038	!$ACC PRESENT(v_p)
[1875]	1039	DO i = nxl, nxr
	1040	DO j = nysv, nyn
[2232]	1041	DO k = nzb+1, nzt
	1042	v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
	1043	tsc(3) * tv_m(k,j,i) ) &
	1044	- tsc(5) * rdf(k) * &
	1045	( v(k,j,i) - v_init(k) ) &
	1046	) * MERGE( 1.0_wp, 0.0_wp, &
[4346]	1047	BTEST( wall_flags_total_0(k,j,i), 2 )&
[2232]	1048	)
[1875]	1049	ENDDO
	1050	ENDDO
	1051	ENDDO
	1052
	1053	!
[3302]	1054	!-- Add turbulence generated by wave breaking (in ocean mode only)
	1055	IF ( wave_breaking .AND. &
	1056	intermediate_timestep_count == intermediate_timestep_count_max ) &
	1057	THEN
	1058	CALL wave_breaking_term( 2 )
	1059	ENDIF
	1060
	1061	!
[1875]	1062	!-- Calculate tendencies for the next Runge-Kutta step
	1063	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1064	IF ( intermediate_timestep_count == 1 ) THEN
[3634]	1065	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	1066	!$ACC PRESENT(tend, tv_m)
[1875]	1067	DO i = nxl, nxr
	1068	DO j = nysv, nyn
[2232]	1069	DO k = nzb+1, nzt
[1875]	1070	tv_m(k,j,i) = tend(k,j,i)
	1071	ENDDO
	1072	ENDDO
	1073	ENDDO
	1074	ELSEIF ( intermediate_timestep_count < &
	1075	intermediate_timestep_count_max ) THEN
[3634]	1076	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	1077	!$ACC PRESENT(tend, tv_m)
[1875]	1078	DO i = nxl, nxr
	1079	DO j = nysv, nyn
[2232]	1080	DO k = nzb+1, nzt
	1081	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
	1082	+ 5.3125_wp * tv_m(k,j,i)
[1875]	1083	ENDDO
	1084	ENDDO
	1085	ENDDO
	1086	ENDIF
	1087	ENDIF
	1088
	1089	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1090
	1091	!
	1092	!-- w-velocity component
	1093	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1094
[3634]	1095	!$ACC KERNELS PRESENT(tend)
[1875]	1096	tend = 0.0_wp
[3634]	1097	!$ACC END KERNELS
[1875]	1098	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1099	IF ( ws_scheme_mom ) THEN
	1100	CALL advec_w_ws
	1101	ELSE
	1102	CALL advec_w_pw
	1103	ENDIF
	1104	ELSE
	1105	CALL advec_w_up
	1106	ENDIF
	1107	CALL diffusion_w
	1108	CALL coriolis( 3 )
	1109
	1110	IF ( .NOT. neutral ) THEN
[3294]	1111	IF ( ocean_mode ) THEN
[2031]	1112	CALL buoyancy( rho_ocean, 3 )
[1875]	1113	ELSE
	1114	IF ( .NOT. humidity ) THEN
	1115	CALL buoyancy( pt, 3 )
	1116	ELSE
	1117	CALL buoyancy( vpt, 3 )
	1118	ENDIF
	1119	ENDIF
	1120	ENDIF
	1121
	1122	!
	1123	!-- Drag by plant canopy
	1124	IF ( plant_canopy ) CALL pcm_tendency( 3 )
	1125
[1914]	1126	!
[3302]	1127	!-- Effect of Stokes drift (in ocean mode only)
	1128	IF ( stokes_force ) CALL stokes_drift_terms( 3 )
	1129
[3684]	1130	CALL module_interface_actions( 'w-tendency' )
[1875]	1131
	1132	!
	1133	!-- Prognostic equation for w-velocity component
[3634]	1134	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
[4346]	1135	!$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_total_0) &
[3634]	1136	!$ACC PRESENT(tsc(2:5)) &
	1137	!$ACC PRESENT(w_p)
[1875]	1138	DO i = nxl, nxr
	1139	DO j = nys, nyn
[2232]	1140	DO k = nzb+1, nzt-1
	1141	w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
	1142	tsc(3) * tw_m(k,j,i) ) &
	1143	- tsc(5) * rdf(k) * w(k,j,i) &
	1144	) * MERGE( 1.0_wp, 0.0_wp, &
[4346]	1145	BTEST( wall_flags_total_0(k,j,i), 3 )&
[2232]	1146	)
[1875]	1147	ENDDO
	1148	ENDDO
	1149	ENDDO
	1150
	1151	!
	1152	!-- Calculate tendencies for the next Runge-Kutta step
	1153	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1154	IF ( intermediate_timestep_count == 1 ) THEN
[3634]	1155	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	1156	!$ACC PRESENT(tend, tw_m)
[1875]	1157	DO i = nxl, nxr
	1158	DO j = nys, nyn
[2232]	1159	DO k = nzb+1, nzt-1
[1875]	1160	tw_m(k,j,i) = tend(k,j,i)
	1161	ENDDO
	1162	ENDDO
	1163	ENDDO
	1164	ELSEIF ( intermediate_timestep_count < &
	1165	intermediate_timestep_count_max ) THEN
[3634]	1166	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	1167	!$ACC PRESENT(tend, tw_m)
[1875]	1168	DO i = nxl, nxr
	1169	DO j = nys, nyn
[2232]	1170	DO k = nzb+1, nzt-1
	1171	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
	1172	+ 5.3125_wp * tw_m(k,j,i)
[1875]	1173	ENDDO
	1174	ENDDO
	1175	ENDDO
	1176	ENDIF
	1177	ENDIF
	1178
	1179	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1180
	1181
	1182	!
	1183	!-- If required, compute prognostic equation for potential temperature
	1184	IF ( .NOT. neutral ) THEN
	1185
	1186	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1187
	1188	!
	1189	!-- pt-tendency terms with communication
	1190	sbt = tsc(2)
	1191	IF ( scalar_advec == 'bc-scheme' ) THEN
	1192
	1193	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1194	!
	1195	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1196	sbt = 1.0_wp
	1197	ENDIF
	1198	tend = 0.0_wp
	1199	CALL advec_s_bc( pt, 'pt' )
	1200
	1201	ENDIF
	1202
	1203	!
	1204	!-- pt-tendency terms with no communication
	1205	IF ( scalar_advec /= 'bc-scheme' ) THEN
[3634]	1206	!$ACC KERNELS PRESENT(tend)
[1875]	1207	tend = 0.0_wp
[3634]	1208	!$ACC END KERNELS
[1875]	1209	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1210	IF ( ws_scheme_sca ) THEN
[4109]	1211	CALL advec_s_ws( advc_flags_s, pt, 'pt', &
	1212	bc_dirichlet_l .OR. bc_radiation_l, &
	1213	bc_dirichlet_n .OR. bc_radiation_n, &
	1214	bc_dirichlet_r .OR. bc_radiation_r, &
	1215	bc_dirichlet_s .OR. bc_radiation_s )
[1875]	1216	ELSE
	1217	CALL advec_s_pw( pt )
	1218	ENDIF
	1219	ELSE
	1220	CALL advec_s_up( pt )
	1221	ENDIF
	1222	ENDIF
	1223
[2232]	1224	CALL diffusion_s( pt, &
	1225	surf_def_h(0)%shf, surf_def_h(1)%shf, &
	1226	surf_def_h(2)%shf, &
	1227	surf_lsm_h%shf, surf_usm_h%shf, &
	1228	surf_def_v(0)%shf, surf_def_v(1)%shf, &
	1229	surf_def_v(2)%shf, surf_def_v(3)%shf, &
	1230	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
	1231	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
	1232	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
	1233	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
[1875]	1234
	1235	!
	1236	!-- Consideration of heat sources within the plant canopy
[3014]	1237	IF ( plant_canopy .AND. &
	1238	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
[1875]	1239	CALL pcm_tendency( 4 )
	1240	ENDIF
	1241
	1242	!
	1243	!-- Large scale advection
	1244	IF ( large_scale_forcing ) THEN
	1245	CALL ls_advec( simulated_time, 'pt' )
	1246	ENDIF
	1247
	1248	!
	1249	!-- Nudging
[2155]	1250	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
[1875]	1251
	1252	!
	1253	!-- If required compute influence of large-scale subsidence/ascent
	1254	IF ( large_scale_subsidence .AND. &
	1255	.NOT. use_subsidence_tendencies ) THEN
	1256	CALL subsidence( tend, pt, pt_init, 2 )
	1257	ENDIF
	1258
[3771]	1259	#if defined( __rrtmg )
[1875]	1260	!
	1261	!-- If required, add tendency due to radiative heating/cooling
[1976]	1262	IF ( radiation .AND. &
[1875]	1263	simulated_time > skip_time_do_radiation ) THEN
	1264	CALL radiation_tendency ( tend )
	1265	ENDIF
[3771]	1266	#endif
[1875]	1267
[3684]	1268	CALL module_interface_actions( 'pt-tendency' )
[1875]	1269
	1270	!
	1271	!-- Prognostic equation for potential temperature
[3634]	1272	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
[4346]	1273	!$ACC PRESENT(pt, tend, tpt_m, wall_flags_total_0) &
[3634]	1274	!$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
	1275	!$ACC PRESENT(tsc(3:5)) &
	1276	!$ACC PRESENT(pt_p)
[1875]	1277	DO i = nxl, nxr
	1278	DO j = nys, nyn
[2232]	1279	DO k = nzb+1, nzt
	1280	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
	1281	tsc(3) * tpt_m(k,j,i) ) &
	1282	- tsc(5) * &
	1283	( pt(k,j,i) - pt_init(k) ) *&
	1284	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
	1285	) &
	1286	* MERGE( 1.0_wp, 0.0_wp, &
[4346]	1287	BTEST( wall_flags_total_0(k,j,i), 0 ) &
[2232]	1288	)
[1875]	1289	ENDDO
	1290	ENDDO
	1291	ENDDO
	1292	!
	1293	!-- Calculate tendencies for the next Runge-Kutta step
	1294	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1295	IF ( intermediate_timestep_count == 1 ) THEN
[3634]	1296	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	1297	!$ACC PRESENT(tend, tpt_m)
[1875]	1298	DO i = nxl, nxr
	1299	DO j = nys, nyn
[2232]	1300	DO k = nzb+1, nzt
[1875]	1301	tpt_m(k,j,i) = tend(k,j,i)
	1302	ENDDO
	1303	ENDDO
	1304	ENDDO
	1305	ELSEIF ( intermediate_timestep_count < &
	1306	intermediate_timestep_count_max ) THEN
[3634]	1307	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	1308	!$ACC PRESENT(tend, tpt_m)
[1875]	1309	DO i = nxl, nxr
	1310	DO j = nys, nyn
[2232]	1311	DO k = nzb+1, nzt
	1312	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	1313	5.3125_wp * tpt_m(k,j,i)
[1875]	1314	ENDDO
	1315	ENDDO
	1316	ENDDO
	1317	ENDIF
	1318	ENDIF
	1319
	1320	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1321
	1322	ENDIF
	1323
	1324	!
[1960]	1325	!-- If required, compute prognostic equation for total water content
	1326	IF ( humidity ) THEN
[1875]	1327
[1960]	1328	CALL cpu_log( log_point(29), 'q-equation', 'start' )
[1875]	1329
	1330	!
	1331	!-- Scalar/q-tendency terms with communication
	1332	sbt = tsc(2)
	1333	IF ( scalar_advec == 'bc-scheme' ) THEN
	1334
	1335	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1336	!
	1337	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1338	sbt = 1.0_wp
	1339	ENDIF
	1340	tend = 0.0_wp
	1341	CALL advec_s_bc( q, 'q' )
	1342
	1343	ENDIF
	1344
	1345	!
	1346	!-- Scalar/q-tendency terms with no communication
	1347	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1348	tend = 0.0_wp
	1349	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1350	IF ( ws_scheme_sca ) THEN
[4109]	1351	CALL advec_s_ws( advc_flags_s, q, 'q', &
	1352	bc_dirichlet_l .OR. bc_radiation_l, &
	1353	bc_dirichlet_n .OR. bc_radiation_n, &
	1354	bc_dirichlet_r .OR. bc_radiation_r, &
	1355	bc_dirichlet_s .OR. bc_radiation_s )
[1875]	1356	ELSE
	1357	CALL advec_s_pw( q )
	1358	ENDIF
	1359	ELSE
	1360	CALL advec_s_up( q )
	1361	ENDIF
	1362	ENDIF
	1363
[2232]	1364	CALL diffusion_s( q, &
	1365	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
	1366	surf_def_h(2)%qsws, &
	1367	surf_lsm_h%qsws, surf_usm_h%qsws, &
	1368	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
	1369	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
	1370	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
	1371	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
	1372	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
	1373	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
[2155]	1374
[1875]	1375	!
[1960]	1376	!-- Sink or source of humidity due to canopy elements
[1875]	1377	IF ( plant_canopy ) CALL pcm_tendency( 5 )
	1378
	1379	!
	1380	!-- Large scale advection
	1381	IF ( large_scale_forcing ) THEN
	1382	CALL ls_advec( simulated_time, 'q' )
	1383	ENDIF
	1384
	1385	!
	1386	!-- Nudging
[2155]	1387	IF ( nudging ) CALL nudge( simulated_time, 'q' )
[1875]	1388
	1389	!
	1390	!-- If required compute influence of large-scale subsidence/ascent
	1391	IF ( large_scale_subsidence .AND. &
	1392	.NOT. use_subsidence_tendencies ) THEN
	1393	CALL subsidence( tend, q, q_init, 3 )
	1394	ENDIF
	1395
[3684]	1396	CALL module_interface_actions( 'q-tendency' )
[1875]	1397
	1398	!
[1960]	1399	!-- Prognostic equation for total water content
[1875]	1400	DO i = nxl, nxr
	1401	DO j = nys, nyn
[2232]	1402	DO k = nzb+1, nzt
	1403	q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
	1404	tsc(3) * tq_m(k,j,i) ) &
	1405	- tsc(5) * rdf_sc(k) * &
	1406	( q(k,j,i) - q_init(k) ) &
	1407	) * MERGE( 1.0_wp, 0.0_wp, &
[4346]	1408	BTEST( wall_flags_total_0(k,j,i), 0 ) &
[2232]	1409	)
[1875]	1410	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
	1411	ENDDO
	1412	ENDDO
	1413	ENDDO
	1414
	1415	!
	1416	!-- Calculate tendencies for the next Runge-Kutta step
	1417	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1418	IF ( intermediate_timestep_count == 1 ) THEN
	1419	DO i = nxl, nxr
	1420	DO j = nys, nyn
[2232]	1421	DO k = nzb+1, nzt
[1875]	1422	tq_m(k,j,i) = tend(k,j,i)
	1423	ENDDO
	1424	ENDDO
	1425	ENDDO
	1426	ELSEIF ( intermediate_timestep_count < &
	1427	intermediate_timestep_count_max ) THEN
	1428	DO i = nxl, nxr
	1429	DO j = nys, nyn
[2232]	1430	DO k = nzb+1, nzt
	1431	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
	1432	+ 5.3125_wp * tq_m(k,j,i)
[1875]	1433	ENDDO
	1434	ENDDO
	1435	ENDDO
	1436	ENDIF
	1437	ENDIF
	1438
[1960]	1439	CALL cpu_log( log_point(29), 'q-equation', 'stop' )
[1875]	1440
	1441	ENDIF
[1960]	1442	!
	1443	!-- If required, compute prognostic equation for scalar
	1444	IF ( passive_scalar ) THEN
[1875]	1445
[1960]	1446	CALL cpu_log( log_point(66), 's-equation', 'start' )
	1447
[1875]	1448	!
[1960]	1449	!-- Scalar/q-tendency terms with communication
	1450	sbt = tsc(2)
	1451	IF ( scalar_advec == 'bc-scheme' ) THEN
	1452
	1453	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1454	!
	1455	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1456	sbt = 1.0_wp
	1457	ENDIF
	1458	tend = 0.0_wp
	1459	CALL advec_s_bc( s, 's' )
	1460
	1461	ENDIF
	1462
	1463	!
	1464	!-- Scalar/q-tendency terms with no communication
	1465	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1466	tend = 0.0_wp
	1467	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1468	IF ( ws_scheme_sca ) THEN
[4109]	1469	CALL advec_s_ws( advc_flags_s, s, 's', &
	1470	bc_dirichlet_l .OR. bc_radiation_l, &
	1471	bc_dirichlet_n .OR. bc_radiation_n, &
	1472	bc_dirichlet_r .OR. bc_radiation_r, &
	1473	bc_dirichlet_s .OR. bc_radiation_s )
[1960]	1474	ELSE
	1475	CALL advec_s_pw( s )
	1476	ENDIF
	1477	ELSE
	1478	CALL advec_s_up( s )
	1479	ENDIF
	1480	ENDIF
	1481
[2232]	1482	CALL diffusion_s( s, &
	1483	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
	1484	surf_def_h(2)%ssws, &
	1485	surf_lsm_h%ssws, surf_usm_h%ssws, &
	1486	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
	1487	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
	1488	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
	1489	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
	1490	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
	1491	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
[2155]	1492
[1960]	1493	!
	1494	!-- Sink or source of humidity due to canopy elements
	1495	IF ( plant_canopy ) CALL pcm_tendency( 7 )
	1496
	1497	!
	1498	!-- Large scale advection. Not implemented for scalars so far.
	1499	! IF ( large_scale_forcing ) THEN
	1500	! CALL ls_advec( simulated_time, 'q' )
	1501	! ENDIF
	1502
	1503	!
	1504	!-- Nudging. Not implemented for scalars so far.
[2155]	1505	! IF ( nudging ) CALL nudge( simulated_time, 'q' )
[1960]	1506
	1507	!
	1508	!-- If required compute influence of large-scale subsidence/ascent.
	1509	!-- Not implemented for scalars so far.
	1510	IF ( large_scale_subsidence .AND. &
	1511	.NOT. use_subsidence_tendencies .AND. &
	1512	.NOT. large_scale_forcing ) THEN
	1513	CALL subsidence( tend, s, s_init, 3 )
	1514	ENDIF
	1515
[3684]	1516	CALL module_interface_actions( 's-tendency' )
[1960]	1517
	1518	!
	1519	!-- Prognostic equation for total water content
	1520	DO i = nxl, nxr
	1521	DO j = nys, nyn
[2232]	1522	DO k = nzb+1, nzt
	1523	s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
	1524	tsc(3) * ts_m(k,j,i) ) &
	1525	- tsc(5) * rdf_sc(k) * &
	1526	( s(k,j,i) - s_init(k) ) &
	1527	) &
	1528	* MERGE( 1.0_wp, 0.0_wp, &
[4346]	1529	BTEST( wall_flags_total_0(k,j,i), 0 ) &
[2232]	1530	)
[1960]	1531	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
	1532	ENDDO
	1533	ENDDO
	1534	ENDDO
	1535
	1536	!
	1537	!-- Calculate tendencies for the next Runge-Kutta step
	1538	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1539	IF ( intermediate_timestep_count == 1 ) THEN
	1540	DO i = nxl, nxr
	1541	DO j = nys, nyn
[2232]	1542	DO k = nzb+1, nzt
[1960]	1543	ts_m(k,j,i) = tend(k,j,i)
	1544	ENDDO
	1545	ENDDO
	1546	ENDDO
	1547	ELSEIF ( intermediate_timestep_count < &
	1548	intermediate_timestep_count_max ) THEN
	1549	DO i = nxl, nxr
	1550	DO j = nys, nyn
[2232]	1551	DO k = nzb+1, nzt
	1552	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
	1553	+ 5.3125_wp * ts_m(k,j,i)
[1960]	1554	ENDDO
	1555	ENDDO
	1556	ENDDO
	1557	ENDIF
	1558	ENDIF
	1559
	1560	CALL cpu_log( log_point(66), 's-equation', 'stop' )
	1561
	1562	ENDIF
[3294]	1563	!
[3837]	1564	!-- Calculate prognostic equations for all other modules
	1565	CALL module_interface_prognostic_equations()
[1875]	1566
[3987]	1567	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_vector', 'end' )
[3885]	1568
[1875]	1569	END SUBROUTINE prognostic_equations_vector
	1570
	1571
	1572	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |