Home

Context Navigation

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

Last change on this file since 765 was 737, checked in by suehring, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 78.5 KB

Rev	Line
[736]	1	MODULE prognostic_equations_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Current revisions:
	5	! -----------------
	6	!
	7	! Former revisions:
	8	! -----------------
	9	! $Id: prognostic_equations.f90 737 2011-08-17 14:17:56Z raasch $
	10	!
[737]	11	! 736 2011-08-17 14:13:26Z suehring
	12	! Bugfix: determination of first thread index i for WS-scheme
	13	!
[736]	14	! 709 2011-03-30 09:31:40Z raasch
	15	! formatting adjustments
	16	!
	17	! 673 2011-01-18 16:19:48Z suehring
	18	! Consideration of the pressure gradient (steered by tsc(4)) during the time
	19	! integration removed.
	20	!
	21	! 667 2010-12-23 12:06:00Z suehring/gryschka
	22	! Calls of the advection routines with WS5 added.
	23	! Calls of ws_statistics added to set the statistical arrays to zero after each
	24	! time step.
	25	!
	26	! 531 2010-04-21 06:47:21Z heinze
	27	! add call of subsidence in the equation for humidity / passive scalar
	28	!
	29	! 411 2009-12-11 14:15:58Z heinze
	30	! add call of subsidence in the equation for potential temperature
	31	!
	32	! 388 2009-09-23 09:40:33Z raasch
	33	! prho is used instead of rho in diffusion_e,
	34	! external pressure gradient
	35	!
	36	! 153 2008-03-19 09:41:30Z steinfeld
	37	! add call of plant_canopy_model in the prognostic equation for
	38	! the potential temperature and for the passive scalar
	39	!
	40	! 138 2007-11-28 10:03:58Z letzel
	41	! add call of subroutines that evaluate the canopy drag terms,
	42	! add wall_*flux to parameter list of calls of diffusion_s
	43	!
	44	! 106 2007-08-16 14:30:26Z raasch
	45	! +uswst, vswst as arguments in calls of diffusion_u\|v,
	46	! loops for u and v are starting from index nxlu, nysv, respectively (needed
	47	! for non-cyclic boundary conditions)
	48	!
	49	! 97 2007-06-21 08:23:15Z raasch
	50	! prognostic equation for salinity, density is calculated from equation of
	51	! state for seawater and is used for calculation of buoyancy,
	52	! +eqn_state_seawater_mod
	53	! diffusion_e is called with argument rho in case of ocean runs,
	54	! new argument zw in calls of diffusion_e, new argument pt_/prho_reference
	55	! in calls of buoyancy and diffusion_e, calc_mean_pt_profile renamed
	56	! calc_mean_profile
	57	!
	58	! 75 2007-03-22 09:54:05Z raasch
	59	! checking for negative q and limiting for positive values,
	60	! z0 removed from arguments in calls of diffusion_u/v/w, uxrp, vynp eliminated,
	61	! subroutine names changed to .._noopt, .._cache, and .._vector,
	62	! moisture renamed humidity, Bott-Chlond-scheme can be used in the
	63	! _vector-version
	64	!
	65	! 19 2007-02-23 04:53:48Z raasch
	66	! Calculation of e, q, and pt extended for gridpoint nzt,
	67	! handling of given temperature/humidity/scalar fluxes at top surface
	68	!
	69	! RCS Log replace by Id keyword, revision history cleaned up
	70	!
	71	! Revision 1.21 2006/08/04 15:01:07 raasch
	72	! upstream scheme can be forced to be used for tke (use_upstream_for_tke)
	73	! regardless of the timestep scheme used for the other quantities,
	74	! new argument diss in call of diffusion_e
	75	!
	76	! Revision 1.1 2000/04/13 14:56:27 schroeter
	77	! Initial revision
	78	!
	79	!
	80	! Description:
	81	! ------------
	82	! Solving the prognostic equations.
	83	!------------------------------------------------------------------------------!
	84
	85	USE arrays_3d
	86	USE control_parameters
	87	USE cpulog
	88	USE eqn_state_seawater_mod
	89	USE grid_variables
	90	USE indices
	91	USE interfaces
	92	USE pegrid
	93	USE pointer_interfaces
	94	USE statistics
	95	USE advec_ws
	96	USE advec_s_pw_mod
	97	USE advec_s_up_mod
	98	USE advec_u_pw_mod
	99	USE advec_u_up_mod
	100	USE advec_v_pw_mod
	101	USE advec_v_up_mod
	102	USE advec_w_pw_mod
	103	USE advec_w_up_mod
	104	USE buoyancy_mod
	105	USE calc_precipitation_mod
	106	USE calc_radiation_mod
	107	USE coriolis_mod
	108	USE diffusion_e_mod
	109	USE diffusion_s_mod
	110	USE diffusion_u_mod
	111	USE diffusion_v_mod
	112	USE diffusion_w_mod
	113	USE impact_of_latent_heat_mod
	114	USE plant_canopy_model_mod
	115	USE production_e_mod
	116	USE subsidence_mod
	117	USE user_actions_mod
	118
	119
	120	PRIVATE
	121	PUBLIC prognostic_equations_noopt, prognostic_equations_cache, &
	122	prognostic_equations_vector
	123
	124	INTERFACE prognostic_equations_noopt
	125	MODULE PROCEDURE prognostic_equations_noopt
	126	END INTERFACE prognostic_equations_noopt
	127
	128	INTERFACE prognostic_equations_cache
	129	MODULE PROCEDURE prognostic_equations_cache
	130	END INTERFACE prognostic_equations_cache
	131
	132	INTERFACE prognostic_equations_vector
	133	MODULE PROCEDURE prognostic_equations_vector
	134	END INTERFACE prognostic_equations_vector
	135
	136
	137	CONTAINS
	138
	139
	140	SUBROUTINE prognostic_equations_noopt
	141
	142	!------------------------------------------------------------------------------!
	143	! Version with single loop optimization
	144	!
	145	! (Optimized over each single prognostic equation.)
	146	!------------------------------------------------------------------------------!
	147
	148	IMPLICIT NONE
	149
	150	CHARACTER (LEN=9) :: time_to_string
	151	INTEGER :: i, i_omp_start, j, k, tn = 0
	152	REAL :: sat, sbt
	153
	154	!
	155	!-- Calculate those variables needed in the tendency terms which need
	156	!-- global communication
	157	CALL calc_mean_profile( pt, 4 )
	158	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	159	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
	160	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
	161	intermediate_timestep_count == 1 ) CALL ws_statistics
	162
	163	!
	164	!-- u-velocity component
	165	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	166
	167	!
	168	!-- u-tendency terms with communication
	169	IF ( momentum_advec == 'ups-scheme' ) THEN
	170	tend = 0.0
	171	CALL advec_u_ups
	172	ENDIF
	173
	174	!
	175	!-- u-tendency terms with no communication
	176	i_omp_start = nxlu
	177	DO i = nxlu, nxr
	178	DO j = nys, nyn
	179	!
	180	!-- Tendency terms
	181	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	182	tend(:,j,i) = 0.0
	183	IF ( ws_scheme_mom ) THEN
	184	CALL advec_u_ws( i, j, i_omp_start, tn )
	185	ELSE
	186	CALL advec_u_pw( i, j )
	187	ENDIF
	188
	189	ELSE
	190	IF ( momentum_advec /= 'ups-scheme' ) THEN
	191	tend(:,j,i) = 0.0
	192	CALL advec_u_up( i, j )
	193	ENDIF
	194	ENDIF
	195	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	196	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m, &
	197	usws_m, uswst_m, v_m, w_m )
	198	ELSE
	199	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
	200	uswst, v, w )
	201	ENDIF
	202	CALL coriolis( i, j, 1 )
	203	IF ( sloping_surface ) CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
	204
	205	!
	206	!-- Drag by plant canopy
	207	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 1 )
	208
	209	!
	210	!-- External pressure gradient
	211	IF ( dp_external ) THEN
	212	DO k = dp_level_ind_b+1, nzt
	213	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	214	ENDDO
	215	ENDIF
	216
	217	CALL user_actions( i, j, 'u-tendency' )
	218
	219	!
	220	!-- Prognostic equation for u-velocity component
	221	DO k = nzb_u_inner(j,i)+1, nzt
	222	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
	223	dt_3d * ( &
	224	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
	225	) - &
	226	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	227	ENDDO
	228
	229	!
	230	!-- Calculate tendencies for the next Runge-Kutta step
	231	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	232	IF ( intermediate_timestep_count == 1 ) THEN
	233	DO k = nzb_u_inner(j,i)+1, nzt
	234	tu_m(k,j,i) = tend(k,j,i)
	235	ENDDO
	236	ELSEIF ( intermediate_timestep_count < &
	237	intermediate_timestep_count_max ) THEN
	238	DO k = nzb_u_inner(j,i)+1, nzt
	239	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	240	ENDDO
	241	ENDIF
	242	ENDIF
	243
	244	ENDDO
	245	ENDDO
	246
	247	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	248
	249	!
	250	!-- v-velocity component
	251	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	252
	253	!
	254	!-- v-tendency terms with communication
	255	IF ( momentum_advec == 'ups-scheme' ) THEN
	256	tend = 0.0
	257	CALL advec_v_ups
	258	ENDIF
	259
	260	!
	261	!-- v-tendency terms with no communication
	262	i_omp_start = nxl
	263	DO i = nxl, nxr
	264	DO j = nysv, nyn
	265	!
	266	!-- Tendency terms
	267	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	268	tend(:,j,i) = 0.0
	269	IF ( ws_scheme_mom ) THEN
	270	CALL advec_v_ws( i, j, i_omp_start, tn )
	271	ELSE
	272	CALL advec_v_pw( i, j )
	273	ENDIF
	274
	275	ELSE
	276	IF ( momentum_advec /= 'ups-scheme' ) THEN
	277	tend(:,j,i) = 0.0
	278	CALL advec_v_up( i, j )
	279	ENDIF
	280	ENDIF
	281	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	282	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
	283	v_m, vsws_m, vswst_m, w_m )
	284	ELSE
	285	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
	286	vsws, vswst, w )
	287	ENDIF
	288	CALL coriolis( i, j, 2 )
	289
	290	!
	291	!-- Drag by plant canopy
	292	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 2 )
	293
	294	!
	295	!-- External pressure gradient
	296	IF ( dp_external ) THEN
	297	DO k = dp_level_ind_b+1, nzt
	298	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	299	ENDDO
	300	ENDIF
	301
	302	CALL user_actions( i, j, 'v-tendency' )
	303
	304	!
	305	!-- Prognostic equation for v-velocity component
	306	DO k = nzb_v_inner(j,i)+1, nzt
	307	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
	308	dt_3d * ( &
	309	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
	310	) - &
	311	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	312	ENDDO
	313
	314	!
	315	!-- Calculate tendencies for the next Runge-Kutta step
	316	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	317	IF ( intermediate_timestep_count == 1 ) THEN
	318	DO k = nzb_v_inner(j,i)+1, nzt
	319	tv_m(k,j,i) = tend(k,j,i)
	320	ENDDO
	321	ELSEIF ( intermediate_timestep_count < &
	322	intermediate_timestep_count_max ) THEN
	323	DO k = nzb_v_inner(j,i)+1, nzt
	324	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	325	ENDDO
	326	ENDIF
	327	ENDIF
	328
	329	ENDDO
	330	ENDDO
	331
	332	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	333
	334	!
	335	!-- w-velocity component
	336	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	337
	338	!
	339	!-- w-tendency terms with communication
	340	IF ( momentum_advec == 'ups-scheme' ) THEN
	341	tend = 0.0
	342	CALL advec_w_ups
	343	ENDIF
	344
	345	!
	346	!-- w-tendency terms with no communication
	347	DO i = nxl, nxr
	348	DO j = nys, nyn
	349	!
	350	!-- Tendency terms
	351	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	352	tend(:,j,i) = 0.0
	353	IF ( ws_scheme_mom ) THEN
	354	CALL advec_w_ws( i, j, i_omp_start, tn )
	355	ELSE
	356	CALL advec_w_pw( i, j )
	357	ENDIF
	358
	359	ELSE
	360	IF ( momentum_advec /= 'ups-scheme' ) THEN
	361	tend(:,j,i) = 0.0
	362	CALL advec_w_up( i, j )
	363	ENDIF
	364	ENDIF
	365	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	366	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y, &
	367	tend, u_m, v_m, w_m )
	368	ELSE
	369	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
	370	tend, u, v, w )
	371	ENDIF
	372	CALL coriolis( i, j, 3 )
	373	IF ( ocean ) THEN
	374	CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
	375	ELSE
	376	IF ( .NOT. humidity ) THEN
	377	CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
	378	ELSE
	379	CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
	380	ENDIF
	381	ENDIF
	382
	383	!
	384	!-- Drag by plant canopy
	385	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 3 )
	386
	387	CALL user_actions( i, j, 'w-tendency' )
	388
	389	!
	390	!-- Prognostic equation for w-velocity component
	391	DO k = nzb_w_inner(j,i)+1, nzt-1
	392	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
	393	dt_3d * ( &
	394	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
	395	) - &
	396	tsc(5) * rdf(k) * w(k,j,i)
	397	ENDDO
	398
	399	!
	400	!-- Calculate tendencies for the next Runge-Kutta step
	401	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	402	IF ( intermediate_timestep_count == 1 ) THEN
	403	DO k = nzb_w_inner(j,i)+1, nzt-1
	404	tw_m(k,j,i) = tend(k,j,i)
	405	ENDDO
	406	ELSEIF ( intermediate_timestep_count < &
	407	intermediate_timestep_count_max ) THEN
	408	DO k = nzb_w_inner(j,i)+1, nzt-1
	409	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	410	ENDDO
	411	ENDIF
	412	ENDIF
	413
	414	ENDDO
	415	ENDDO
	416
	417	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	418
	419	!
	420	!-- potential temperature
	421	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	422
	423	!
	424	!-- pt-tendency terms with communication
	425	sat = tsc(1)
	426	sbt = tsc(2)
	427	IF ( scalar_advec == 'bc-scheme' ) THEN
	428
	429	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	430	!
	431	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	432	!-- switched on. Thus:
	433	sat = 1.0
	434	sbt = 1.0
	435	ENDIF
	436	tend = 0.0
	437	CALL advec_s_bc( pt, 'pt' )
	438	ELSE
	439	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
	440	tend = 0.0
	441	CALL advec_s_ups( pt, 'pt' )
	442	ENDIF
	443	ENDIF
	444
	445	!
	446	!-- pt-tendency terms with no communication
	447	DO i = nxl, nxr
	448	DO j = nys, nyn
	449	!
	450	!-- Tendency terms
	451	IF ( scalar_advec == 'bc-scheme' ) THEN
	452	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
	453	wall_heatflux, tend )
	454	ELSE
	455	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	456	tend(:,j,i) = 0.0
	457	IF ( ws_scheme_sca ) THEN
	458	CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
	459	diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
	460	ELSE
	461	CALL advec_s_pw( i, j, pt )
	462	ENDIF
	463	ELSE
	464	IF ( scalar_advec /= 'ups-scheme' ) THEN
	465	tend(:,j,i) = 0.0
	466	CALL advec_s_up( i, j, pt )
	467	ENDIF
	468	ENDIF
	469	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	470	THEN
	471	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
	472	tswst_m, wall_heatflux, tend )
	473	ELSE
	474	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
	475	wall_heatflux, tend )
	476	ENDIF
	477	ENDIF
	478
	479	!
	480	!-- If required compute heating/cooling due to long wave radiation
	481	!-- processes
	482	IF ( radiation ) THEN
	483	CALL calc_radiation( i, j )
	484	ENDIF
	485
	486	!
	487	!-- If required compute impact of latent heat due to precipitation
	488	IF ( precipitation ) THEN
	489	CALL impact_of_latent_heat( i, j )
	490	ENDIF
	491
	492	!
	493	!-- Consideration of heat sources within the plant canopy
	494	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	495	CALL plant_canopy_model( i, j, 4 )
	496	ENDIF
	497
	498	!
	499	!-- If required compute influence of large-scale subsidence/ascent
	500	IF ( large_scale_subsidence ) THEN
	501	CALL subsidence ( i, j, tend, pt, pt_init )
	502	ENDIF
	503
	504	CALL user_actions( i, j, 'pt-tendency' )
	505
	506	!
	507	!-- Prognostic equation for potential temperature
	508	DO k = nzb_s_inner(j,i)+1, nzt
	509	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
	510	dt_3d * ( &
	511	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
	512	) - &
	513	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
	514	ENDDO
	515
	516	!
	517	!-- Calculate tendencies for the next Runge-Kutta step
	518	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	519	IF ( intermediate_timestep_count == 1 ) THEN
	520	DO k = nzb_s_inner(j,i)+1, nzt
	521	tpt_m(k,j,i) = tend(k,j,i)
	522	ENDDO
	523	ELSEIF ( intermediate_timestep_count < &
	524	intermediate_timestep_count_max ) THEN
	525	DO k = nzb_s_inner(j,i)+1, nzt
	526	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
	527	ENDDO
	528	ENDIF
	529	ENDIF
	530
	531	ENDDO
	532	ENDDO
	533
	534	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	535
	536	!
	537	!-- If required, compute prognostic equation for salinity
	538	IF ( ocean ) THEN
	539
	540	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	541
	542	!
	543	!-- sa-tendency terms with communication
	544	sat = tsc(1)
	545	sbt = tsc(2)
	546	IF ( scalar_advec == 'bc-scheme' ) THEN
	547
	548	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	549	!
	550	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	551	!-- switched on. Thus:
	552	sat = 1.0
	553	sbt = 1.0
	554	ENDIF
	555	tend = 0.0
	556	CALL advec_s_bc( sa, 'sa' )
	557	ELSE
	558	IF ( tsc(2) /= 2.0 ) THEN
	559	IF ( scalar_advec == 'ups-scheme' ) THEN
	560	tend = 0.0
	561	CALL advec_s_ups( sa, 'sa' )
	562	ENDIF
	563	ENDIF
	564	ENDIF
	565
	566	!
	567	!-- sa terms with no communication
	568	DO i = nxl, nxr
	569	DO j = nys, nyn
	570	!
	571	!-- Tendency-terms
	572	IF ( scalar_advec == 'bc-scheme' ) THEN
	573	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	574	wall_salinityflux, tend )
	575	ELSE
	576	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	577	tend(:,j,i) = 0.0
	578	IF ( ws_scheme_sca ) THEN
	579	CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, &
	580	diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
	581	ELSE
	582	CALL advec_s_pw( i, j, sa )
	583	ENDIF
	584
	585	ELSE
	586	IF ( scalar_advec /= 'ups-scheme' ) THEN
	587	tend(:,j,i) = 0.0
	588	CALL advec_s_up( i, j, sa )
	589	ENDIF
	590	ENDIF
	591	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	592	wall_salinityflux, tend )
	593	ENDIF
	594
	595	CALL user_actions( i, j, 'sa-tendency' )
	596
	597	!
	598	!-- Prognostic equation for salinity
	599	DO k = nzb_s_inner(j,i)+1, nzt
	600	sa_p(k,j,i) = sat * sa(k,j,i) + &
	601	dt_3d * ( &
	602	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	603	) - &
	604	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
	605	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	606	ENDDO
	607
	608	!
	609	!-- Calculate tendencies for the next Runge-Kutta step
	610	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	611	IF ( intermediate_timestep_count == 1 ) THEN
	612	DO k = nzb_s_inner(j,i)+1, nzt
	613	tsa_m(k,j,i) = tend(k,j,i)
	614	ENDDO
	615	ELSEIF ( intermediate_timestep_count < &
	616	intermediate_timestep_count_max ) THEN
	617	DO k = nzb_s_inner(j,i)+1, nzt
	618	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	619	5.3125 * tsa_m(k,j,i)
	620	ENDDO
	621	ENDIF
	622	ENDIF
	623
	624	!
	625	!-- Calculate density by the equation of state for seawater
	626	CALL eqn_state_seawater( i, j )
	627
	628	ENDDO
	629	ENDDO
	630
	631	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	632
	633	ENDIF
	634
	635	!
	636	!-- If required, compute prognostic equation for total water content / scalar
	637	IF ( humidity .OR. passive_scalar ) THEN
	638
	639	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	640
	641	!
	642	!-- Scalar/q-tendency terms with communication
	643	sat = tsc(1)
	644	sbt = tsc(2)
	645	IF ( scalar_advec == 'bc-scheme' ) THEN
	646
	647	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	648	!
	649	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	650	!-- switched on. Thus:
	651	sat = 1.0
	652	sbt = 1.0
	653	ENDIF
	654	tend = 0.0
	655	CALL advec_s_bc( q, 'q' )
	656	ELSE
	657	IF ( tsc(2) /= 2.0 ) THEN
	658	IF ( scalar_advec == 'ups-scheme' ) THEN
	659	tend = 0.0
	660	CALL advec_s_ups( q, 'q' )
	661	ENDIF
	662	ENDIF
	663	ENDIF
	664
	665	!
	666	!-- Scalar/q-tendency terms with no communication
	667	DO i = nxl, nxr
	668	DO j = nys, nyn
	669	!
	670	!-- Tendency-terms
	671	IF ( scalar_advec == 'bc-scheme' ) THEN
	672	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
	673	wall_qflux, tend )
	674	ELSE
	675	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	676	tend(:,j,i) = 0.0
	677	IF ( ws_scheme_sca ) THEN
	678	CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
	679	diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
	680	ELSE
	681	CALL advec_s_pw( i, j, q )
	682	ENDIF
	683	ELSE
	684	IF ( scalar_advec /= 'ups-scheme' ) THEN
	685	tend(:,j,i) = 0.0
	686	CALL advec_s_up( i, j, q )
	687	ENDIF
	688	ENDIF
	689	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	690	THEN
	691	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
	692	qswst_m, wall_qflux, tend )
	693	ELSE
	694	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
	695	wall_qflux, tend )
	696	ENDIF
	697	ENDIF
	698
	699	!
	700	!-- If required compute decrease of total water content due to
	701	!-- precipitation
	702	IF ( precipitation ) THEN
	703	CALL calc_precipitation( i, j )
	704	ENDIF
	705
	706	!
	707	!-- Sink or source of scalar concentration due to canopy elements
	708	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
	709
	710	!
	711	!-- If required compute influence of large-scale subsidence/ascent
	712	IF ( large_scale_subsidence ) THEN
	713	CALL subsidence ( i, j, tend, q, q_init )
	714	ENDIF
	715
	716	CALL user_actions( i, j, 'q-tendency' )
	717
	718	!
	719	!-- Prognostic equation for total water content / scalar
	720	DO k = nzb_s_inner(j,i)+1, nzt
	721	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
	722	dt_3d * ( &
	723	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	724	) - &
	725	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
	726	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	727	ENDDO
	728
	729	!
	730	!-- Calculate tendencies for the next Runge-Kutta step
	731	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	732	IF ( intermediate_timestep_count == 1 ) THEN
	733	DO k = nzb_s_inner(j,i)+1, nzt
	734	tq_m(k,j,i) = tend(k,j,i)
	735	ENDDO
	736	ELSEIF ( intermediate_timestep_count < &
	737	intermediate_timestep_count_max ) THEN
	738	DO k = nzb_s_inner(j,i)+1, nzt
	739	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	740	ENDDO
	741	ENDIF
	742	ENDIF
	743
	744	ENDDO
	745	ENDDO
	746
	747	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	748
	749	ENDIF
	750
	751	!
	752	!-- If required, compute prognostic equation for turbulent kinetic
	753	!-- energy (TKE)
	754	IF ( .NOT. constant_diffusion ) THEN
	755
	756	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	757
	758	!
	759	!-- TKE-tendency terms with communication
	760	CALL production_e_init
	761
	762	sat = tsc(1)
	763	sbt = tsc(2)
	764	IF ( .NOT. use_upstream_for_tke ) THEN
	765	IF ( scalar_advec == 'bc-scheme' ) THEN
	766
	767	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	768	!
	769	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	770	!-- switched on. Thus:
	771	sat = 1.0
	772	sbt = 1.0
	773	ENDIF
	774	tend = 0.0
	775	CALL advec_s_bc( e, 'e' )
	776	ELSE
	777	IF ( tsc(2) /= 2.0 ) THEN
	778	IF ( scalar_advec == 'ups-scheme' ) THEN
	779	tend = 0.0
	780	CALL advec_s_ups( e, 'e' )
	781	ENDIF
	782	ENDIF
	783	ENDIF
	784	ENDIF
	785
	786	!
	787	!-- TKE-tendency terms with no communication
	788	DO i = nxl, nxr
	789	DO j = nys, nyn
	790	!
	791	!-- Tendency-terms
	792	IF ( scalar_advec == 'bc-scheme' .AND. &
	793	.NOT. use_upstream_for_tke ) THEN
	794	IF ( .NOT. humidity ) THEN
	795	IF ( ocean ) THEN
	796	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	797	l_grid, prho, prho_reference, rif, &
	798	tend, zu, zw )
	799	ELSE
	800	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	801	l_grid, pt, pt_reference, rif, tend, &
	802	zu, zw )
	803	ENDIF
	804	ELSE
	805	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	806	l_grid, vpt, pt_reference, rif, tend, zu, &
	807	zw )
	808	ENDIF
	809	ELSE
	810	IF ( use_upstream_for_tke ) THEN
	811	tend(:,j,i) = 0.0
	812	CALL advec_s_up( i, j, e )
	813	ELSE
	814	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	815	THEN
	816	tend(:,j,i) = 0.0
	817	IF ( ws_scheme_sca ) THEN
	818	CALL advec_s_ws( i, j, e, 'e', flux_s_e, &
	819	diss_s_e, flux_l_e, diss_l_e, i_omp_start, tn )
	820	ELSE
	821	CALL advec_s_pw( i, j, e )
	822	ENDIF
	823	ELSE
	824	IF ( scalar_advec /= 'ups-scheme' ) THEN
	825	tend(:,j,i) = 0.0
	826	CALL advec_s_up( i, j, e )
	827	ENDIF
	828	ENDIF
	829	ENDIF
	830	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	831	THEN
	832	IF ( .NOT. humidity ) THEN
	833	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	834	km_m, l_grid, pt_m, pt_reference, &
	835	rif_m, tend, zu, zw )
	836	ELSE
	837	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	838	km_m, l_grid, vpt_m, pt_reference, &
	839	rif_m, tend, zu, zw )
	840	ENDIF
	841	ELSE
	842	IF ( .NOT. humidity ) THEN
	843	IF ( ocean ) THEN
	844	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	845	km, l_grid, prho, prho_reference, &
	846	rif, tend, zu, zw )
	847	ELSE
	848	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	849	km, l_grid, pt, pt_reference, rif, &
	850	tend, zu, zw )
	851	ENDIF
	852	ELSE
	853	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	854	l_grid, vpt, pt_reference, rif, tend, &
	855	zu, zw )
	856	ENDIF
	857	ENDIF
	858	ENDIF
	859	CALL production_e( i, j )
	860
	861	!
	862	!-- Additional sink term for flows through plant canopies
	863	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
	864
	865	CALL user_actions( i, j, 'e-tendency' )
	866
	867	!
	868	!-- Prognostic equation for TKE.
	869	!-- Eliminate negative TKE values, which can occur due to numerical
	870	!-- reasons in the course of the integration. In such cases the old TKE
	871	!-- value is reduced by 90%.
	872	DO k = nzb_s_inner(j,i)+1, nzt
	873	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
	874	dt_3d * ( &
	875	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	876	)
	877	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	878	ENDDO
	879
	880	!
	881	!-- Calculate tendencies for the next Runge-Kutta step
	882	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	883	IF ( intermediate_timestep_count == 1 ) THEN
	884	DO k = nzb_s_inner(j,i)+1, nzt
	885	te_m(k,j,i) = tend(k,j,i)
	886	ENDDO
	887	ELSEIF ( intermediate_timestep_count < &
	888	intermediate_timestep_count_max ) THEN
	889	DO k = nzb_s_inner(j,i)+1, nzt
	890	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	891	ENDDO
	892	ENDIF
	893	ENDIF
	894
	895	ENDDO
	896	ENDDO
	897
	898	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	899
	900	ENDIF
	901
	902
	903	END SUBROUTINE prognostic_equations_noopt
	904
	905
	906	SUBROUTINE prognostic_equations_cache
	907
	908	!------------------------------------------------------------------------------!
	909	! Version with one optimized loop over all equations. It is only allowed to
	910	! be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
	911	!
	912	! Here the calls of most subroutines are embedded in two DO loops over i and j,
	913	! so communication between CPUs is not allowed (does not make sense) within
	914	! these loops.
	915	!
	916	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
	917	!------------------------------------------------------------------------------!
	918
	919	IMPLICIT NONE
	920
	921	CHARACTER (LEN=9) :: time_to_string
	922	INTEGER :: i, i_omp_start, j, k, omp_get_thread_num, tn = 0
	923	LOGICAL :: loop_start
	924
	925
	926	!
	927	!-- Time measurement can only be performed for the whole set of equations
	928	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
	929
	930
	931	!
	932	!-- Calculate those variables needed in the tendency terms which need
	933	!-- global communication
	934	CALL calc_mean_profile( pt, 4 )
	935	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	936	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
	937	IF ( .NOT. constant_diffusion ) CALL production_e_init
	938	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
	939	intermediate_timestep_count == 1 ) CALL ws_statistics
	940
	941	!
	942	!-- Loop over all prognostic equations
	943	!$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
	944
	945	!$ tn = omp_get_thread_num()
	946	loop_start = .TRUE.
	947	!$OMP DO
	948	DO i = nxl, nxr
	949
	950	!
	951	!-- Store the first loop index. It differs for each thread and is required
	952	!-- later in advec_ws
	953	IF ( loop_start ) THEN
	954	loop_start = .FALSE.
	955	i_omp_start = i
	956	ENDIF
	957
	958	DO j = nys, nyn
	959	!
	960	!-- Tendency terms for u-velocity component
	961	IF ( .NOT. outflow_l .OR. i > nxl ) THEN
	962
	963	tend(:,j,i) = 0.0
	964	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	965	IF ( ws_scheme_mom ) THEN
	966	IF ( outflow_l .AND. i_omp_start == nxl ) THEN
	967	! CALL local_diss( i, j, u) ! dissipation control
	968	CALL advec_u_ws( i, j, i_omp_start + 1, tn )
	969	ELSE
	970	CALL advec_u_ws( i, j, i_omp_start, tn )
	971	ENDIF
	972	ELSE
	973	CALL advec_u_pw( i, j )
	974	ENDIF
	975	ELSE
	976	CALL advec_u_up( i, j )
	977	ENDIF
	978	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	979	THEN
	980	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, &
	981	u_m, usws_m, uswst_m, v_m, w_m )
	982	ELSE
	983	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
	984	usws, uswst, v, w )
	985	ENDIF
	986	CALL coriolis( i, j, 1 )
	987	IF ( sloping_surface ) CALL buoyancy( i, j, pt, pt_reference, 1, &
	988	4 )
	989
	990	!
	991	!-- Drag by plant canopy
	992	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 1 )
	993
	994	!
	995	!-- External pressure gradient
	996	IF ( dp_external ) THEN
	997	DO k = dp_level_ind_b+1, nzt
	998	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	999	ENDDO
	1000	ENDIF
	1001
	1002	CALL user_actions( i, j, 'u-tendency' )
	1003
	1004	!
	1005	!-- Prognostic equation for u-velocity component
	1006	DO k = nzb_u_inner(j,i)+1, nzt
	1007	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
	1008	dt_3d * ( &
	1009	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
	1010	) - &
	1011	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	1012	ENDDO
	1013
	1014	!
	1015	!-- Calculate tendencies for the next Runge-Kutta step
	1016	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1017	IF ( intermediate_timestep_count == 1 ) THEN
	1018	DO k = nzb_u_inner(j,i)+1, nzt
	1019	tu_m(k,j,i) = tend(k,j,i)
	1020	ENDDO
	1021	ELSEIF ( intermediate_timestep_count < &
	1022	intermediate_timestep_count_max ) THEN
	1023	DO k = nzb_u_inner(j,i)+1, nzt
	1024	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	1025	ENDDO
	1026	ENDIF
	1027	ENDIF
	1028
	1029	ENDIF
	1030
	1031	!
	1032	!-- Tendency terms for v-velocity component
	1033	IF ( .NOT. outflow_s .OR. j > nys ) THEN
	1034
	1035	tend(:,j,i) = 0.0
	1036	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1037	IF ( ws_scheme_mom ) THEN
	1038	! CALL local_diss( i, j, v)
	1039	CALL advec_v_ws( i, j, i_omp_start, tn )
	1040	ELSE
	1041	CALL advec_v_pw( i, j )
	1042	ENDIF
	1043	ELSE
	1044	CALL advec_v_up( i, j )
	1045	ENDIF
	1046	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	1047	THEN
	1048	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, &
	1049	u_m, v_m, vsws_m, vswst_m, w_m )
	1050	ELSE
	1051	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
	1052	vsws, vswst, w )
	1053	ENDIF
	1054	CALL coriolis( i, j, 2 )
	1055
	1056	!
	1057	!-- Drag by plant canopy
	1058	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 2 )
	1059
	1060	!
	1061	!-- External pressure gradient
	1062	IF ( dp_external ) THEN
	1063	DO k = dp_level_ind_b+1, nzt
	1064	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1065	ENDDO
	1066	ENDIF
	1067
	1068	CALL user_actions( i, j, 'v-tendency' )
	1069
	1070	!
	1071	!-- Prognostic equation for v-velocity component
	1072	DO k = nzb_v_inner(j,i)+1, nzt
	1073	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
	1074	dt_3d * ( &
	1075	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
	1076	) - &
	1077	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	1078	ENDDO
	1079
	1080	!
	1081	!-- Calculate tendencies for the next Runge-Kutta step
	1082	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1083	IF ( intermediate_timestep_count == 1 ) THEN
	1084	DO k = nzb_v_inner(j,i)+1, nzt
	1085	tv_m(k,j,i) = tend(k,j,i)
	1086	ENDDO
	1087	ELSEIF ( intermediate_timestep_count < &
	1088	intermediate_timestep_count_max ) THEN
	1089	DO k = nzb_v_inner(j,i)+1, nzt
	1090	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	1091	ENDDO
	1092	ENDIF
	1093	ENDIF
	1094
	1095	ENDIF
	1096
	1097	!
	1098	!-- Tendency terms for w-velocity component
	1099	tend(:,j,i) = 0.0
	1100	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1101	IF ( ws_scheme_mom ) THEN
	1102	! CALL local_diss( i, j, w)
	1103	CALL advec_w_ws( i, j, i_omp_start, tn )
	1104	ELSE
	1105	CALL advec_w_pw( i, j )
	1106	END IF
	1107	ELSE
	1108	CALL advec_w_up( i, j )
	1109	ENDIF
	1110	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	1111	THEN
	1112	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, &
	1113	km_damp_y, tend, u_m, v_m, w_m )
	1114	ELSE
	1115	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
	1116	tend, u, v, w )
	1117	ENDIF
	1118	CALL coriolis( i, j, 3 )
	1119	IF ( ocean ) THEN
	1120	CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
	1121	ELSE
	1122	IF ( .NOT. humidity ) THEN
	1123	CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
	1124	ELSE
	1125	CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
	1126	ENDIF
	1127	ENDIF
	1128
	1129	!
	1130	!-- Drag by plant canopy
	1131	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 3 )
	1132
	1133	CALL user_actions( i, j, 'w-tendency' )
	1134
	1135	!
	1136	!-- Prognostic equation for w-velocity component
	1137	DO k = nzb_w_inner(j,i)+1, nzt-1
	1138	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
	1139	dt_3d * ( &
	1140	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
	1141	) - &
	1142	tsc(5) * rdf(k) * w(k,j,i)
	1143	ENDDO
	1144
	1145	!
	1146	!-- Calculate tendencies for the next Runge-Kutta step
	1147	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1148	IF ( intermediate_timestep_count == 1 ) THEN
	1149	DO k = nzb_w_inner(j,i)+1, nzt-1
	1150	tw_m(k,j,i) = tend(k,j,i)
	1151	ENDDO
	1152	ELSEIF ( intermediate_timestep_count < &
	1153	intermediate_timestep_count_max ) THEN
	1154	DO k = nzb_w_inner(j,i)+1, nzt-1
	1155	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1156	ENDDO
	1157	ENDIF
	1158	ENDIF
	1159
	1160	!
	1161	!-- Tendency terms for potential temperature
	1162	tend(:,j,i) = 0.0
	1163	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1164	IF ( ws_scheme_sca ) THEN
	1165	! CALL local_diss( i, j, pt )
	1166	CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
	1167	diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
	1168	ELSE
	1169	CALL advec_s_pw( i, j, pt )
	1170	ENDIF
	1171	ELSE
	1172	CALL advec_s_up( i, j, pt )
	1173	ENDIF
	1174	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	1175	THEN
	1176	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
	1177	tswst_m, wall_heatflux, tend )
	1178	ELSE
	1179	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
	1180	wall_heatflux, tend )
	1181	ENDIF
	1182
	1183	!
	1184	!-- If required compute heating/cooling due to long wave radiation
	1185	!-- processes
	1186	IF ( radiation ) THEN
	1187	CALL calc_radiation( i, j )
	1188	ENDIF
	1189
	1190	!
	1191	!-- If required compute impact of latent heat due to precipitation
	1192	IF ( precipitation ) THEN
	1193	CALL impact_of_latent_heat( i, j )
	1194	ENDIF
	1195
	1196	!
	1197	!-- Consideration of heat sources within the plant canopy
	1198	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	1199	CALL plant_canopy_model( i, j, 4 )
	1200	ENDIF
	1201
	1202
	1203	!-- If required compute influence of large-scale subsidence/ascent
	1204	IF ( large_scale_subsidence ) THEN
	1205	CALL subsidence ( i, j, tend, pt, pt_init )
	1206	ENDIF
	1207
	1208
	1209	CALL user_actions( i, j, 'pt-tendency' )
	1210
	1211	!
	1212	!-- Prognostic equation for potential temperature
	1213	DO k = nzb_s_inner(j,i)+1, nzt
	1214	pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
	1215	dt_3d * ( &
	1216	tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
	1217	) - &
	1218	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
	1219	ENDDO
	1220
	1221	!
	1222	!-- Calculate tendencies for the next Runge-Kutta step
	1223	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1224	IF ( intermediate_timestep_count == 1 ) THEN
	1225	DO k = nzb_s_inner(j,i)+1, nzt
	1226	tpt_m(k,j,i) = tend(k,j,i)
	1227	ENDDO
	1228	ELSEIF ( intermediate_timestep_count < &
	1229	intermediate_timestep_count_max ) THEN
	1230	DO k = nzb_s_inner(j,i)+1, nzt
	1231	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1232	5.3125 * tpt_m(k,j,i)
	1233	ENDDO
	1234	ENDIF
	1235	ENDIF
	1236
	1237	!
	1238	!-- If required, compute prognostic equation for salinity
	1239	IF ( ocean ) THEN
	1240
	1241	!
	1242	!-- Tendency-terms for salinity
	1243	tend(:,j,i) = 0.0
	1244	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	1245	THEN
	1246	IF ( ws_scheme_sca ) THEN
	1247	! CALL local_diss( i, j, sa )
	1248	CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, &
	1249	diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
	1250	ELSE
	1251	CALL advec_s_pw( i, j, sa )
	1252	ENDIF
	1253	ELSE
	1254	CALL advec_s_up( i, j, sa )
	1255	ENDIF
	1256	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	1257	wall_salinityflux, tend )
	1258
	1259	CALL user_actions( i, j, 'sa-tendency' )
	1260
	1261	!
	1262	!-- Prognostic equation for salinity
	1263	DO k = nzb_s_inner(j,i)+1, nzt
	1264	sa_p(k,j,i) = tsc(1) * sa(k,j,i) + &
	1265	dt_3d * ( &
	1266	tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	1267	) - &
	1268	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
	1269	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1270	ENDDO
	1271
	1272	!
	1273	!-- Calculate tendencies for the next Runge-Kutta step
	1274	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1275	IF ( intermediate_timestep_count == 1 ) THEN
	1276	DO k = nzb_s_inner(j,i)+1, nzt
	1277	tsa_m(k,j,i) = tend(k,j,i)
	1278	ENDDO
	1279	ELSEIF ( intermediate_timestep_count < &
	1280	intermediate_timestep_count_max ) THEN
	1281	DO k = nzb_s_inner(j,i)+1, nzt
	1282	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1283	5.3125 * tsa_m(k,j,i)
	1284	ENDDO
	1285	ENDIF
	1286	ENDIF
	1287
	1288	!
	1289	!-- Calculate density by the equation of state for seawater
	1290	CALL eqn_state_seawater( i, j )
	1291
	1292	ENDIF
	1293
	1294	!
	1295	!-- If required, compute prognostic equation for total water content /
	1296	!-- scalar
	1297	IF ( humidity .OR. passive_scalar ) THEN
	1298
	1299	!
	1300	!-- Tendency-terms for total water content / scalar
	1301	tend(:,j,i) = 0.0
	1302	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	1303	THEN
	1304	IF ( ws_scheme_sca ) THEN
	1305	! CALL local_diss( i, j, q )
	1306	CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
	1307	diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
	1308	ELSE
	1309	CALL advec_s_pw( i, j, q )
	1310	ENDIF
	1311	ELSE
	1312	CALL advec_s_up( i, j, q )
	1313	ENDIF
	1314	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	1315	THEN
	1316	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
	1317	qswst_m, wall_qflux, tend )
	1318	ELSE
	1319	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
	1320	wall_qflux, tend )
	1321	ENDIF
	1322
	1323	!
	1324	!-- If required compute decrease of total water content due to
	1325	!-- precipitation
	1326	IF ( precipitation ) THEN
	1327	CALL calc_precipitation( i, j )
	1328	ENDIF
	1329
	1330	!
	1331	!-- Sink or source of scalar concentration due to canopy elements
	1332	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
	1333
	1334	!-- If required compute influence of large-scale subsidence/ascent
	1335	IF ( large_scale_subsidence ) THEN
	1336	CALL subsidence ( i, j, tend, q, q_init )
	1337	ENDIF
	1338
	1339	CALL user_actions( i, j, 'q-tendency' )
	1340
	1341	!
	1342	!-- Prognostic equation for total water content / scalar
	1343	DO k = nzb_s_inner(j,i)+1, nzt
	1344	q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
	1345	dt_3d * ( &
	1346	tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	1347	) - &
	1348	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
	1349	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	1350	ENDDO
	1351
	1352	!
	1353	!-- Calculate tendencies for the next Runge-Kutta step
	1354	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1355	IF ( intermediate_timestep_count == 1 ) THEN
	1356	DO k = nzb_s_inner(j,i)+1, nzt
	1357	tq_m(k,j,i) = tend(k,j,i)
	1358	ENDDO
	1359	ELSEIF ( intermediate_timestep_count < &
	1360	intermediate_timestep_count_max ) THEN
	1361	DO k = nzb_s_inner(j,i)+1, nzt
	1362	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1363	5.3125 * tq_m(k,j,i)
	1364	ENDDO
	1365	ENDIF
	1366	ENDIF
	1367
	1368	ENDIF
	1369
	1370	!
	1371	!-- If required, compute prognostic equation for turbulent kinetic
	1372	!-- energy (TKE)
	1373	IF ( .NOT. constant_diffusion ) THEN
	1374
	1375	!
	1376	!-- Tendency-terms for TKE
	1377	tend(:,j,i) = 0.0
	1378	IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	1379	.AND. .NOT. use_upstream_for_tke ) THEN
	1380	IF ( ws_scheme_sca ) THEN
	1381	! CALL local_diss( i, j, e )
	1382	CALL advec_s_ws( i, j, e, 'e', flux_s_e, &
	1383	diss_s_e, flux_l_e, diss_l_e , i_omp_start, tn )
	1384	ELSE
	1385	CALL advec_s_pw( i, j, e )
	1386	ENDIF
	1387	ELSE
	1388	CALL advec_s_up( i, j, e )
	1389	ENDIF
	1390	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	1391	THEN
	1392	IF ( .NOT. humidity ) THEN
	1393	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	1394	km_m, l_grid, pt_m, pt_reference, &
	1395	rif_m, tend, zu, zw )
	1396	ELSE
	1397	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	1398	km_m, l_grid, vpt_m, pt_reference, &
	1399	rif_m, tend, zu, zw )
	1400	ENDIF
	1401	ELSE
	1402	IF ( .NOT. humidity ) THEN
	1403	IF ( ocean ) THEN
	1404	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	1405	km, l_grid, prho, prho_reference, &
	1406	rif, tend, zu, zw )
	1407	ELSE
	1408	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	1409	km, l_grid, pt, pt_reference, rif, &
	1410	tend, zu, zw )
	1411	ENDIF
	1412	ELSE
	1413	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	1414	l_grid, vpt, pt_reference, rif, tend, &
	1415	zu, zw )
	1416	ENDIF
	1417	ENDIF
	1418	CALL production_e( i, j )
	1419
	1420	!
	1421	!-- Additional sink term for flows through plant canopies
	1422	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
	1423
	1424	CALL user_actions( i, j, 'e-tendency' )
	1425
	1426	!
	1427	!-- Prognostic equation for TKE.
	1428	!-- Eliminate negative TKE values, which can occur due to numerical
	1429	!-- reasons in the course of the integration. In such cases the old
	1430	!-- TKE value is reduced by 90%.
	1431	DO k = nzb_s_inner(j,i)+1, nzt
	1432	e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
	1433	dt_3d * ( &
	1434	tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	1435	)
	1436	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	1437	ENDDO
	1438
	1439	!
	1440	!-- Calculate tendencies for the next Runge-Kutta step
	1441	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1442	IF ( intermediate_timestep_count == 1 ) THEN
	1443	DO k = nzb_s_inner(j,i)+1, nzt
	1444	te_m(k,j,i) = tend(k,j,i)
	1445	ENDDO
	1446	ELSEIF ( intermediate_timestep_count < &
	1447	intermediate_timestep_count_max ) THEN
	1448	DO k = nzb_s_inner(j,i)+1, nzt
	1449	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1450	5.3125 * te_m(k,j,i)
	1451	ENDDO
	1452	ENDIF
	1453	ENDIF
	1454
	1455	ENDIF ! TKE equation
	1456
	1457	ENDDO
	1458	ENDDO
	1459	!$OMP END PARALLEL
	1460
	1461	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
	1462
	1463
	1464	END SUBROUTINE prognostic_equations_cache
	1465
	1466
	1467	SUBROUTINE prognostic_equations_vector
	1468
	1469	!------------------------------------------------------------------------------!
	1470	! Version for vector machines
	1471	!------------------------------------------------------------------------------!
	1472
	1473	IMPLICIT NONE
	1474
	1475	CHARACTER (LEN=9) :: time_to_string
	1476	INTEGER :: i, j, k
	1477	REAL :: sat, sbt
	1478
	1479	!
	1480	!-- Calculate those variables needed in the tendency terms which need
	1481	!-- global communication
	1482	CALL calc_mean_profile( pt, 4 )
	1483	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	1484	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
	1485	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
	1486	intermediate_timestep_count == 1 ) CALL ws_statistics
	1487
	1488	!
	1489	!-- u-velocity component
	1490	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	1491
	1492	!
	1493	!-- u-tendency terms with communication
	1494	IF ( momentum_advec == 'ups-scheme' ) THEN
	1495	tend = 0.0
	1496	CALL advec_u_ups
	1497	ENDIF
	1498
	1499	!
	1500	!-- u-tendency terms with no communication
	1501	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1502	tend = 0.0
	1503	IF ( ws_scheme_mom ) THEN
	1504	CALL advec_u_ws
	1505	ELSE
	1506	CALL advec_u_pw
	1507	ENDIF
	1508	ELSE
	1509	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1510	tend = 0.0
	1511	CALL advec_u_up
	1512	ENDIF
	1513	ENDIF
	1514	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	1515	CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
	1516	uswst_m, v_m, w_m )
	1517	ELSE
	1518	CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
	1519	ENDIF
	1520	CALL coriolis( 1 )
	1521	IF ( sloping_surface ) CALL buoyancy( pt, pt_reference, 1, 4 )
	1522
	1523	!
	1524	!-- Drag by plant canopy
	1525	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
	1526
	1527	!
	1528	!-- External pressure gradient
	1529	IF ( dp_external ) THEN
	1530	DO i = nxlu, nxr
	1531	DO j = nys, nyn
	1532	DO k = dp_level_ind_b+1, nzt
	1533	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	1534	ENDDO
	1535	ENDDO
	1536	ENDDO
	1537	ENDIF
	1538
	1539	CALL user_actions( 'u-tendency' )
	1540
	1541	!
	1542	!-- Prognostic equation for u-velocity component
	1543	DO i = nxlu, nxr
	1544	DO j = nys, nyn
	1545	DO k = nzb_u_inner(j,i)+1, nzt
	1546	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
	1547	dt_3d * ( &
	1548	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
	1549	) - &
	1550	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	1551	ENDDO
	1552	ENDDO
	1553	ENDDO
	1554
	1555	!
	1556	!-- Calculate tendencies for the next Runge-Kutta step
	1557	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1558	IF ( intermediate_timestep_count == 1 ) THEN
	1559	DO i = nxlu, nxr
	1560	DO j = nys, nyn
	1561	DO k = nzb_u_inner(j,i)+1, nzt
	1562	tu_m(k,j,i) = tend(k,j,i)
	1563	ENDDO
	1564	ENDDO
	1565	ENDDO
	1566	ELSEIF ( intermediate_timestep_count < &
	1567	intermediate_timestep_count_max ) THEN
	1568	DO i = nxlu, nxr
	1569	DO j = nys, nyn
	1570	DO k = nzb_u_inner(j,i)+1, nzt
	1571	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	1572	ENDDO
	1573	ENDDO
	1574	ENDDO
	1575	ENDIF
	1576	ENDIF
	1577
	1578	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	1579
	1580	!
	1581	!-- v-velocity component
	1582	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	1583
	1584	!
	1585	!-- v-tendency terms with communication
	1586	IF ( momentum_advec == 'ups-scheme' ) THEN
	1587	tend = 0.0
	1588	CALL advec_v_ups
	1589	ENDIF
	1590
	1591	!
	1592	!-- v-tendency terms with no communication
	1593	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1594	tend = 0.0
	1595	IF ( ws_scheme_mom ) THEN
	1596	CALL advec_v_ws
	1597	ELSE
	1598	CALL advec_v_pw
	1599	END IF
	1600	ELSE
	1601	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1602	tend = 0.0
	1603	CALL advec_v_up
	1604	ENDIF
	1605	ENDIF
	1606	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	1607	CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
	1608	vswst_m, w_m )
	1609	ELSE
	1610	CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
	1611	ENDIF
	1612	CALL coriolis( 2 )
	1613
	1614	!
	1615	!-- Drag by plant canopy
	1616	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
	1617
	1618	!
	1619	!-- External pressure gradient
	1620	IF ( dp_external ) THEN
	1621	DO i = nxl, nxr
	1622	DO j = nysv, nyn
	1623	DO k = dp_level_ind_b+1, nzt
	1624	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1625	ENDDO
	1626	ENDDO
	1627	ENDDO
	1628	ENDIF
	1629
	1630	CALL user_actions( 'v-tendency' )
	1631
	1632	!
	1633	!-- Prognostic equation for v-velocity component
	1634	DO i = nxl, nxr
	1635	DO j = nysv, nyn
	1636	DO k = nzb_v_inner(j,i)+1, nzt
	1637	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
	1638	dt_3d * ( &
	1639	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
	1640	) - &
	1641	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	1642	ENDDO
	1643	ENDDO
	1644	ENDDO
	1645
	1646	!
	1647	!-- Calculate tendencies for the next Runge-Kutta step
	1648	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1649	IF ( intermediate_timestep_count == 1 ) THEN
	1650	DO i = nxl, nxr
	1651	DO j = nysv, nyn
	1652	DO k = nzb_v_inner(j,i)+1, nzt
	1653	tv_m(k,j,i) = tend(k,j,i)
	1654	ENDDO
	1655	ENDDO
	1656	ENDDO
	1657	ELSEIF ( intermediate_timestep_count < &
	1658	intermediate_timestep_count_max ) THEN
	1659	DO i = nxl, nxr
	1660	DO j = nysv, nyn
	1661	DO k = nzb_v_inner(j,i)+1, nzt
	1662	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	1663	ENDDO
	1664	ENDDO
	1665	ENDDO
	1666	ENDIF
	1667	ENDIF
	1668
	1669	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1670
	1671	!
	1672	!-- w-velocity component
	1673	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1674
	1675	!
	1676	!-- w-tendency terms with communication
	1677	IF ( momentum_advec == 'ups-scheme' ) THEN
	1678	tend = 0.0
	1679	CALL advec_w_ups
	1680	ENDIF
	1681
	1682	!
	1683	!-- w-tendency terms with no communication
	1684	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1685	tend = 0.0
	1686	IF ( ws_scheme_mom ) THEN
	1687	CALL advec_w_ws
	1688	ELSE
	1689	CALL advec_w_pw
	1690	ENDIF
	1691	ELSE
	1692	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1693	tend = 0.0
	1694	CALL advec_w_up
	1695	ENDIF
	1696	ENDIF
	1697	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	1698	CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
	1699	v_m, w_m )
	1700	ELSE
	1701	CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
	1702	ENDIF
	1703	CALL coriolis( 3 )
	1704	IF ( ocean ) THEN
	1705	CALL buoyancy( rho, rho_reference, 3, 64 )
	1706	ELSE
	1707	IF ( .NOT. humidity ) THEN
	1708	CALL buoyancy( pt, pt_reference, 3, 4 )
	1709	ELSE
	1710	CALL buoyancy( vpt, pt_reference, 3, 44 )
	1711	ENDIF
	1712	ENDIF
	1713
	1714	!
	1715	!-- Drag by plant canopy
	1716	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
	1717
	1718	CALL user_actions( 'w-tendency' )
	1719
	1720	!
	1721	!-- Prognostic equation for w-velocity component
	1722	DO i = nxl, nxr
	1723	DO j = nys, nyn
	1724	DO k = nzb_w_inner(j,i)+1, nzt-1
	1725	w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
	1726	dt_3d * ( &
	1727	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
	1728	) - &
	1729	tsc(5) * rdf(k) * w(k,j,i)
	1730	ENDDO
	1731	ENDDO
	1732	ENDDO
	1733
	1734	!
	1735	!-- Calculate tendencies for the next Runge-Kutta step
	1736	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1737	IF ( intermediate_timestep_count == 1 ) THEN
	1738	DO i = nxl, nxr
	1739	DO j = nys, nyn
	1740	DO k = nzb_w_inner(j,i)+1, nzt-1
	1741	tw_m(k,j,i) = tend(k,j,i)
	1742	ENDDO
	1743	ENDDO
	1744	ENDDO
	1745	ELSEIF ( intermediate_timestep_count < &
	1746	intermediate_timestep_count_max ) THEN
	1747	DO i = nxl, nxr
	1748	DO j = nys, nyn
	1749	DO k = nzb_w_inner(j,i)+1, nzt-1
	1750	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1751	ENDDO
	1752	ENDDO
	1753	ENDDO
	1754	ENDIF
	1755	ENDIF
	1756
	1757	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1758
	1759	!
	1760	!-- potential temperature
	1761	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1762
	1763	!
	1764	!-- pt-tendency terms with communication
	1765	sat = tsc(1)
	1766	sbt = tsc(2)
	1767	IF ( scalar_advec == 'bc-scheme' ) THEN
	1768
	1769	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1770	!
	1771	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1772	!-- switched on. Thus:
	1773	sat = 1.0
	1774	sbt = 1.0
	1775	ENDIF
	1776	tend = 0.0
	1777	CALL advec_s_bc( pt, 'pt' )
	1778	ELSE
	1779	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
	1780	tend = 0.0
	1781	CALL advec_s_ups( pt, 'pt' )
	1782	ENDIF
	1783	ENDIF
	1784
	1785	!
	1786	!-- pt-tendency terms with no communication
	1787	IF ( scalar_advec == 'bc-scheme' ) THEN
	1788	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
	1789	tend )
	1790	ELSE
	1791	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1792	tend = 0.0
	1793	IF ( ws_scheme_sca ) THEN
	1794	CALL advec_s_ws( pt, 'pt' )
	1795	ELSE
	1796	CALL advec_s_pw( pt )
	1797	ENDIF
	1798	ELSE
	1799	IF ( scalar_advec /= 'ups-scheme' ) THEN
	1800	tend = 0.0
	1801	CALL advec_s_up( pt )
	1802	ENDIF
	1803	ENDIF
	1804	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	1805	CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, &
	1806	wall_heatflux, tend )
	1807	ELSE
	1808	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
	1809	tend )
	1810	ENDIF
	1811	ENDIF
	1812
	1813	!
	1814	!-- If required compute heating/cooling due to long wave radiation
	1815	!-- processes
	1816	IF ( radiation ) THEN
	1817	CALL calc_radiation
	1818	ENDIF
	1819
	1820	!
	1821	!-- If required compute impact of latent heat due to precipitation
	1822	IF ( precipitation ) THEN
	1823	CALL impact_of_latent_heat
	1824	ENDIF
	1825
	1826	!
	1827	!-- Consideration of heat sources within the plant canopy
	1828	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	1829	CALL plant_canopy_model( 4 )
	1830	ENDIF
	1831
	1832	!--If required compute influence of large-scale subsidence/ascent
	1833	IF ( large_scale_subsidence ) THEN
	1834	CALL subsidence ( tend, pt, pt_init )
	1835	ENDIF
	1836
	1837	CALL user_actions( 'pt-tendency' )
	1838
	1839	!
	1840	!-- Prognostic equation for potential temperature
	1841	DO i = nxl, nxr
	1842	DO j = nys, nyn
	1843	DO k = nzb_s_inner(j,i)+1, nzt
	1844	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
	1845	dt_3d * ( &
	1846	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
	1847	) - &
	1848	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
	1849	ENDDO
	1850	ENDDO
	1851	ENDDO
	1852
	1853	!
	1854	!-- Calculate tendencies for the next Runge-Kutta step
	1855	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1856	IF ( intermediate_timestep_count == 1 ) THEN
	1857	DO i = nxl, nxr
	1858	DO j = nys, nyn
	1859	DO k = nzb_s_inner(j,i)+1, nzt
	1860	tpt_m(k,j,i) = tend(k,j,i)
	1861	ENDDO
	1862	ENDDO
	1863	ENDDO
	1864	ELSEIF ( intermediate_timestep_count < &
	1865	intermediate_timestep_count_max ) THEN
	1866	DO i = nxl, nxr
	1867	DO j = nys, nyn
	1868	DO k = nzb_s_inner(j,i)+1, nzt
	1869	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
	1870	ENDDO
	1871	ENDDO
	1872	ENDDO
	1873	ENDIF
	1874	ENDIF
	1875
	1876	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1877
	1878	!
	1879	!-- If required, compute prognostic equation for salinity
	1880	IF ( ocean ) THEN
	1881
	1882	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	1883
	1884	!
	1885	!-- sa-tendency terms with communication
	1886	sat = tsc(1)
	1887	sbt = tsc(2)
	1888	IF ( scalar_advec == 'bc-scheme' ) THEN
	1889
	1890	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1891	!
	1892	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1893	!-- switched on. Thus:
	1894	sat = 1.0
	1895	sbt = 1.0
	1896	ENDIF
	1897	tend = 0.0
	1898	CALL advec_s_bc( sa, 'sa' )
	1899	ELSE
	1900	IF ( tsc(2) /= 2.0 ) THEN
	1901	IF ( scalar_advec == 'ups-scheme' ) THEN
	1902	tend = 0.0
	1903	CALL advec_s_ups( sa, 'sa' )
	1904	ENDIF
	1905	ENDIF
	1906	ENDIF
	1907
	1908	!
	1909	!-- sa-tendency terms with no communication
	1910	IF ( scalar_advec == 'bc-scheme' ) THEN
	1911	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
	1912	wall_salinityflux, tend )
	1913	ELSE
	1914	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1915	tend = 0.0
	1916	IF ( ws_scheme_sca ) THEN
	1917	CALL advec_s_ws( sa, 'sa' )
	1918	ELSE
	1919	CALL advec_s_pw( sa )
	1920	ENDIF
	1921	ELSE
	1922	IF ( scalar_advec /= 'ups-scheme' ) THEN
	1923	tend = 0.0
	1924	CALL advec_s_up( sa )
	1925	ENDIF
	1926	ENDIF
	1927	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
	1928	wall_salinityflux, tend )
	1929	ENDIF
	1930
	1931	CALL user_actions( 'sa-tendency' )
	1932
	1933	!
	1934	!-- Prognostic equation for salinity
	1935	DO i = nxl, nxr
	1936	DO j = nys, nyn
	1937	DO k = nzb_s_inner(j,i)+1, nzt
	1938	sa_p(k,j,i) = sat * sa(k,j,i) + &
	1939	dt_3d * ( &
	1940	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	1941	) - &
	1942	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
	1943	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1944	ENDDO
	1945	ENDDO
	1946	ENDDO
	1947
	1948	!
	1949	!-- Calculate tendencies for the next Runge-Kutta step
	1950	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1951	IF ( intermediate_timestep_count == 1 ) THEN
	1952	DO i = nxl, nxr
	1953	DO j = nys, nyn
	1954	DO k = nzb_s_inner(j,i)+1, nzt
	1955	tsa_m(k,j,i) = tend(k,j,i)
	1956	ENDDO
	1957	ENDDO
	1958	ENDDO
	1959	ELSEIF ( intermediate_timestep_count < &
	1960	intermediate_timestep_count_max ) THEN
	1961	DO i = nxl, nxr
	1962	DO j = nys, nyn
	1963	DO k = nzb_s_inner(j,i)+1, nzt
	1964	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1965	5.3125 * tsa_m(k,j,i)
	1966	ENDDO
	1967	ENDDO
	1968	ENDDO
	1969	ENDIF
	1970	ENDIF
	1971
	1972	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	1973
	1974	!
	1975	!-- Calculate density by the equation of state for seawater
	1976	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	1977	CALL eqn_state_seawater
	1978	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	1979
	1980	ENDIF
	1981
	1982	!
	1983	!-- If required, compute prognostic equation for total water content / scalar
	1984	IF ( humidity .OR. passive_scalar ) THEN
	1985
	1986	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	1987
	1988	!
	1989	!-- Scalar/q-tendency terms with communication
	1990	sat = tsc(1)
	1991	sbt = tsc(2)
	1992	IF ( scalar_advec == 'bc-scheme' ) THEN
	1993
	1994	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1995	!
	1996	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1997	!-- switched on. Thus:
	1998	sat = 1.0
	1999	sbt = 1.0
	2000	ENDIF
	2001	tend = 0.0
	2002	CALL advec_s_bc( q, 'q' )
	2003	ELSE
	2004	IF ( tsc(2) /= 2.0 ) THEN
	2005	IF ( scalar_advec == 'ups-scheme' ) THEN
	2006	tend = 0.0
	2007	CALL advec_s_ups( q, 'q' )
	2008	ENDIF
	2009	ENDIF
	2010	ENDIF
	2011
	2012	!
	2013	!-- Scalar/q-tendency terms with no communication
	2014	IF ( scalar_advec == 'bc-scheme' ) THEN
	2015	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, wall_qflux, tend )
	2016	ELSE
	2017	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	2018	tend = 0.0
	2019	IF ( ws_scheme_sca ) THEN
	2020	CALL advec_s_ws( q, 'q' )
	2021	ELSE
	2022	CALL advec_s_pw( q )
	2023	ENDIF
	2024	ELSE
	2025	IF ( scalar_advec /= 'ups-scheme' ) THEN
	2026	tend = 0.0
	2027	CALL advec_s_up( q )
	2028	ENDIF
	2029	ENDIF
	2030	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	2031	CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, &
	2032	wall_qflux, tend )
	2033	ELSE
	2034	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, &
	2035	wall_qflux, tend )
	2036	ENDIF
	2037	ENDIF
	2038
	2039	!
	2040	!-- If required compute decrease of total water content due to
	2041	!-- precipitation
	2042	IF ( precipitation ) THEN
	2043	CALL calc_precipitation
	2044	ENDIF
	2045
	2046	!
	2047	!-- Sink or source of scalar concentration due to canopy elements
	2048	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
	2049
	2050	!
	2051	!-- If required compute influence of large-scale subsidence/ascent
	2052	IF ( large_scale_subsidence ) THEN
	2053	CALL subsidence ( tend, q, q_init )
	2054	ENDIF
	2055
	2056	CALL user_actions( 'q-tendency' )
	2057
	2058	!
	2059	!-- Prognostic equation for total water content / scalar
	2060	DO i = nxl, nxr
	2061	DO j = nys, nyn
	2062	DO k = nzb_s_inner(j,i)+1, nzt
	2063	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
	2064	dt_3d * ( &
	2065	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	2066	) - &
	2067	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
	2068	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	2069	ENDDO
	2070	ENDDO
	2071	ENDDO
	2072
	2073	!
	2074	!-- Calculate tendencies for the next Runge-Kutta step
	2075	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2076	IF ( intermediate_timestep_count == 1 ) THEN
	2077	DO i = nxl, nxr
	2078	DO j = nys, nyn
	2079	DO k = nzb_s_inner(j,i)+1, nzt
	2080	tq_m(k,j,i) = tend(k,j,i)
	2081	ENDDO
	2082	ENDDO
	2083	ENDDO
	2084	ELSEIF ( intermediate_timestep_count < &
	2085	intermediate_timestep_count_max ) THEN
	2086	DO i = nxl, nxr
	2087	DO j = nys, nyn
	2088	DO k = nzb_s_inner(j,i)+1, nzt
	2089	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	2090	ENDDO
	2091	ENDDO
	2092	ENDDO
	2093	ENDIF
	2094	ENDIF
	2095
	2096	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	2097
	2098	ENDIF
	2099
	2100	!
	2101	!-- If required, compute prognostic equation for turbulent kinetic
	2102	!-- energy (TKE)
	2103	IF ( .NOT. constant_diffusion ) THEN
	2104
	2105	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	2106
	2107	!
	2108	!-- TKE-tendency terms with communication
	2109	CALL production_e_init
	2110
	2111	sat = tsc(1)
	2112	sbt = tsc(2)
	2113	IF ( .NOT. use_upstream_for_tke ) THEN
	2114	IF ( scalar_advec == 'bc-scheme' ) THEN
	2115
	2116	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2117	!
	2118	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	2119	!-- switched on. Thus:
	2120	sat = 1.0
	2121	sbt = 1.0
	2122	ENDIF
	2123	tend = 0.0
	2124	CALL advec_s_bc( e, 'e' )
	2125	ELSE
	2126	IF ( tsc(2) /= 2.0 ) THEN
	2127	IF ( scalar_advec == 'ups-scheme' ) THEN
	2128	tend = 0.0
	2129	CALL advec_s_ups( e, 'e' )
	2130	ENDIF
	2131	ENDIF
	2132	ENDIF
	2133	ENDIF
	2134
	2135	!
	2136	!-- TKE-tendency terms with no communication
	2137	IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) &
	2138	THEN
	2139	IF ( .NOT. humidity ) THEN
	2140	IF ( ocean ) THEN
	2141	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	2142	prho, prho_reference, rif, tend, zu, zw )
	2143	ELSE
	2144	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
	2145	pt_reference, rif, tend, zu, zw )
	2146	ENDIF
	2147	ELSE
	2148	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
	2149	pt_reference, rif, tend, zu, zw )
	2150	ENDIF
	2151	ELSE
	2152	IF ( use_upstream_for_tke ) THEN
	2153	tend = 0.0
	2154	CALL advec_s_up( e )
	2155	ELSE
	2156	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	2157	tend = 0.0
	2158	IF ( ws_scheme_sca ) THEN
	2159	CALL advec_s_ws( e, 'e' )
	2160	ELSE
	2161	CALL advec_s_pw( e )
	2162	ENDIF
	2163	ELSE
	2164	IF ( scalar_advec /= 'ups-scheme' ) THEN
	2165	tend = 0.0
	2166	CALL advec_s_up( e )
	2167	ENDIF
	2168	ENDIF
	2169	ENDIF
	2170	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	2171	IF ( .NOT. humidity ) THEN
	2172	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
	2173	pt_m, pt_reference, rif_m, tend, zu, zw )
	2174	ELSE
	2175	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
	2176	vpt_m, pt_reference, rif_m, tend, zu, zw )
	2177	ENDIF
	2178	ELSE
	2179	IF ( .NOT. humidity ) THEN
	2180	IF ( ocean ) THEN
	2181	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	2182	prho, prho_reference, rif, tend, zu, zw )
	2183	ELSE
	2184	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	2185	pt, pt_reference, rif, tend, zu, zw )
	2186	ENDIF
	2187	ELSE
	2188	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
	2189	pt_reference, rif, tend, zu, zw )
	2190	ENDIF
	2191	ENDIF
	2192	ENDIF
	2193	CALL production_e
	2194
	2195	!
	2196	!-- Additional sink term for flows through plant canopies
	2197	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
	2198	CALL user_actions( 'e-tendency' )
	2199
	2200	!
	2201	!-- Prognostic equation for TKE.
	2202	!-- Eliminate negative TKE values, which can occur due to numerical
	2203	!-- reasons in the course of the integration. In such cases the old TKE
	2204	!-- value is reduced by 90%.
	2205	DO i = nxl, nxr
	2206	DO j = nys, nyn
	2207	DO k = nzb_s_inner(j,i)+1, nzt
	2208	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
	2209	dt_3d * ( &
	2210	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	2211	)
	2212	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	2213	ENDDO
	2214	ENDDO
	2215	ENDDO
	2216
	2217	!
	2218	!-- Calculate tendencies for the next Runge-Kutta step
	2219	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2220	IF ( intermediate_timestep_count == 1 ) THEN
	2221	DO i = nxl, nxr
	2222	DO j = nys, nyn
	2223	DO k = nzb_s_inner(j,i)+1, nzt
	2224	te_m(k,j,i) = tend(k,j,i)
	2225	ENDDO
	2226	ENDDO
	2227	ENDDO
	2228	ELSEIF ( intermediate_timestep_count < &
	2229	intermediate_timestep_count_max ) THEN
	2230	DO i = nxl, nxr
	2231	DO j = nys, nyn
	2232	DO k = nzb_s_inner(j,i)+1, nzt
	2233	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	2234	ENDDO
	2235	ENDDO
	2236	ENDDO
	2237	ENDIF
	2238	ENDIF
	2239
	2240	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	2241
	2242	ENDIF
	2243
	2244
	2245	END SUBROUTINE prognostic_equations_vector
	2246
	2247
	2248	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |