Home

Context Navigation

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

Last change on this file since 268 was 240, checked in by letzel, 16 years ago
External pressure gradient (check_parameters, init_3d_model, header, modules, parin, prognostic_equations) New topography case 'single_street_canyon'
Property svn:keywords set to `Id`
File size: 70.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |