Home

Context Navigation

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

Last change on this file since 4656 was 4649, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 65.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |