Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |