Home

Context Navigation

prognostic_equations.f90 @ 131

Last change on this file since 131 was 129, checked in by letzel, 17 years ago

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

Bugfix: assignment of fluxes at walls

Updates to documentation:
chapter_4.2.html#mode_dvrp
chapter_3.5.4.html#time_series

Default for mrun options -q and -n is "sla3" for lctit. Queues bes1 and bes2
removed. DOC/misc/Tsubame.html updated.

Modified default paths (/work/...) for lctit in .mrun.config.default

Property svn:keywords set to Id

File size: 66.3 KB

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

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

Context Navigation

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

Download in other formats: