Home

Context Navigation

prognostic_equations.f90 @ 113

Last change on this file since 113 was 110, checked in by raasch, 17 years ago

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:

Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:

Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

Property svn:keywords set to Id

File size: 65.6 KB

Rev	Line
[1]	1	MODULE prognostic_equations_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[110]	6	!
[98]	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: prognostic_equations.f90 110 2007-10-05 05:13:14Z raasch $
	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
[19]	371	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
[1]	372	ELSE
	373	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	374	tend(:,j,i) = 0.0
	375	CALL advec_s_pw( i, j, pt )
	376	ELSE
	377	IF ( scalar_advec /= 'ups-scheme' ) THEN
	378	tend(:,j,i) = 0.0
	379	CALL advec_s_up( i, j, pt )
	380	ENDIF
	381	ENDIF
	382	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	383	THEN
[19]	384	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
	385	tswst_m, tend )
[1]	386	ELSE
[19]	387	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
[1]	388	ENDIF
	389	ENDIF
	390
	391	!
	392	!-- If required compute heating/cooling due to long wave radiation
	393	!-- processes
	394	IF ( radiation ) THEN
	395	CALL calc_radiation( i, j )
	396	ENDIF
	397
	398	!
	399	!-- If required compute impact of latent heat due to precipitation
	400	IF ( precipitation ) THEN
	401	CALL impact_of_latent_heat( i, j )
	402	ENDIF
	403	CALL user_actions( i, j, 'pt-tendency' )
	404
	405	!
	406	!-- Prognostic equation for potential temperature
[19]	407	DO k = nzb_s_inner(j,i)+1, nzt
[1]	408	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
	409	dt_3d * ( &
	410	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
	411	) - &
	412	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
	413	ENDDO
	414
	415	!
	416	!-- Calculate tendencies for the next Runge-Kutta step
	417	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	418	IF ( intermediate_timestep_count == 1 ) THEN
[19]	419	DO k = nzb_s_inner(j,i)+1, nzt
[1]	420	tpt_m(k,j,i) = tend(k,j,i)
	421	ENDDO
	422	ELSEIF ( intermediate_timestep_count < &
	423	intermediate_timestep_count_max ) THEN
[19]	424	DO k = nzb_s_inner(j,i)+1, nzt
[1]	425	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
	426	ENDDO
	427	ENDIF
	428	ENDIF
	429
	430	ENDDO
	431	ENDDO
	432
	433	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	434
	435	!
[95]	436	!-- If required, compute prognostic equation for salinity
	437	IF ( ocean ) THEN
	438
	439	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	440
	441	!
	442	!-- sa-tendency terms with communication
	443	sat = tsc(1)
	444	sbt = tsc(2)
	445	IF ( scalar_advec == 'bc-scheme' ) THEN
	446
	447	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	448	!
	449	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	450	!-- switched on. Thus:
	451	sat = 1.0
	452	sbt = 1.0
	453	ENDIF
	454	tend = 0.0
	455	CALL advec_s_bc( sa, 'sa' )
	456	ELSE
	457	IF ( tsc(2) /= 2.0 ) THEN
	458	IF ( scalar_advec == 'ups-scheme' ) THEN
	459	tend = 0.0
	460	CALL advec_s_ups( sa, 'sa' )
	461	ENDIF
	462	ENDIF
	463	ENDIF
	464
	465	!
	466	!-- sa terms with no communication
	467	DO i = nxl, nxr
	468	DO j = nys, nyn
	469	!
	470	!-- Tendency-terms
	471	IF ( scalar_advec == 'bc-scheme' ) THEN
	472	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	473	tend )
	474	ELSE
	475	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	476	tend(:,j,i) = 0.0
	477	CALL advec_s_pw( i, j, sa )
	478	ELSE
	479	IF ( scalar_advec /= 'ups-scheme' ) THEN
	480	tend(:,j,i) = 0.0
	481	CALL advec_s_up( i, j, sa )
	482	ENDIF
	483	ENDIF
	484	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	485	tend )
	486	ENDIF
	487
	488	CALL user_actions( i, j, 'sa-tendency' )
	489
	490	!
	491	!-- Prognostic equation for salinity
	492	DO k = nzb_s_inner(j,i)+1, nzt
	493	sa_p(k,j,i) = sat * sa(k,j,i) + &
	494	dt_3d * ( &
	495	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	496	) - &
	497	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
	498	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	499	ENDDO
	500
	501	!
	502	!-- Calculate tendencies for the next Runge-Kutta step
	503	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	504	IF ( intermediate_timestep_count == 1 ) THEN
	505	DO k = nzb_s_inner(j,i)+1, nzt
	506	tsa_m(k,j,i) = tend(k,j,i)
	507	ENDDO
	508	ELSEIF ( intermediate_timestep_count < &
	509	intermediate_timestep_count_max ) THEN
	510	DO k = nzb_s_inner(j,i)+1, nzt
	511	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	512	5.3125 * tsa_m(k,j,i)
	513	ENDDO
	514	ENDIF
	515	ENDIF
	516
[96]	517	!
	518	!-- Calculate density by the equation of state for seawater
	519	CALL eqn_state_seawater( i, j )
	520
[95]	521	ENDDO
	522	ENDDO
	523
	524	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	525
	526	ENDIF
	527
	528	!
[1]	529	!-- If required, compute prognostic equation for total water content / scalar
[75]	530	IF ( humidity .OR. passive_scalar ) THEN
[1]	531
	532	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	533
	534	!
	535	!-- Scalar/q-tendency terms with communication
[70]	536	sat = tsc(1)
	537	sbt = tsc(2)
[1]	538	IF ( scalar_advec == 'bc-scheme' ) THEN
[70]	539
	540	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[1]	541	!
[70]	542	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	543	!-- switched on. Thus:
	544	sat = 1.0
	545	sbt = 1.0
	546	ENDIF
[1]	547	tend = 0.0
	548	CALL advec_s_bc( q, 'q' )
	549	ELSE
	550	IF ( tsc(2) /= 2.0 ) THEN
	551	IF ( scalar_advec == 'ups-scheme' ) THEN
	552	tend = 0.0
	553	CALL advec_s_ups( q, 'q' )
	554	ENDIF
	555	ENDIF
	556	ENDIF
	557
	558	!
	559	!-- Scalar/q-tendency terms with no communication
	560	DO i = nxl, nxr
	561	DO j = nys, nyn
	562	!
	563	!-- Tendency-terms
	564	IF ( scalar_advec == 'bc-scheme' ) THEN
[19]	565	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, tend )
[1]	566	ELSE
	567	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	568	tend(:,j,i) = 0.0
	569	CALL advec_s_pw( i, j, q )
	570	ELSE
	571	IF ( scalar_advec /= 'ups-scheme' ) THEN
	572	tend(:,j,i) = 0.0
	573	CALL advec_s_up( i, j, q )
	574	ENDIF
	575	ENDIF
	576	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	577	THEN
	578	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
[19]	579	qswst_m, tend )
	580	ELSE
	581	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
[1]	582	tend )
	583	ENDIF
	584	ENDIF
	585
	586	!
	587	!-- If required compute decrease of total water content due to
	588	!-- precipitation
	589	IF ( precipitation ) THEN
	590	CALL calc_precipitation( i, j )
	591	ENDIF
	592	CALL user_actions( i, j, 'q-tendency' )
	593
	594	!
	595	!-- Prognostic equation for total water content / scalar
[19]	596	DO k = nzb_s_inner(j,i)+1, nzt
[1]	597	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
	598	dt_3d * ( &
	599	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	600	) - &
	601	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
[73]	602	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
[1]	603	ENDDO
	604
	605	!
	606	!-- Calculate tendencies for the next Runge-Kutta step
	607	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	608	IF ( intermediate_timestep_count == 1 ) THEN
[19]	609	DO k = nzb_s_inner(j,i)+1, nzt
[1]	610	tq_m(k,j,i) = tend(k,j,i)
	611	ENDDO
	612	ELSEIF ( intermediate_timestep_count < &
	613	intermediate_timestep_count_max ) THEN
[19]	614	DO k = nzb_s_inner(j,i)+1, nzt
[1]	615	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	616	ENDDO
	617	ENDIF
	618	ENDIF
	619
	620	ENDDO
	621	ENDDO
	622
	623	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	624
	625	ENDIF
	626
	627	!
	628	!-- If required, compute prognostic equation for turbulent kinetic
	629	!-- energy (TKE)
	630	IF ( .NOT. constant_diffusion ) THEN
	631
	632	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	633
	634	!
	635	!-- TKE-tendency terms with communication
	636	CALL production_e_init
[70]	637
	638	sat = tsc(1)
	639	sbt = tsc(2)
[1]	640	IF ( .NOT. use_upstream_for_tke ) THEN
	641	IF ( scalar_advec == 'bc-scheme' ) THEN
[70]	642
	643	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[1]	644	!
[70]	645	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	646	!-- switched on. Thus:
	647	sat = 1.0
	648	sbt = 1.0
	649	ENDIF
[1]	650	tend = 0.0
	651	CALL advec_s_bc( e, 'e' )
	652	ELSE
	653	IF ( tsc(2) /= 2.0 ) THEN
	654	IF ( scalar_advec == 'ups-scheme' ) THEN
	655	tend = 0.0
	656	CALL advec_s_ups( e, 'e' )
	657	ENDIF
	658	ENDIF
	659	ENDIF
	660	ENDIF
	661
	662	!
	663	!-- TKE-tendency terms with no communication
	664	DO i = nxl, nxr
	665	DO j = nys, nyn
	666	!
	667	!-- Tendency-terms
	668	IF ( scalar_advec == 'bc-scheme' .AND. &
	669	.NOT. use_upstream_for_tke ) THEN
[75]	670	IF ( .NOT. humidity ) THEN
[97]	671	IF ( ocean ) THEN
	672	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	673	l_grid, rho, prho_reference, rif, tend, &
	674	zu, zw )
	675	ELSE
	676	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	677	l_grid, pt, pt_reference, rif, tend, &
	678	zu, zw )
	679	ENDIF
[1]	680	ELSE
[97]	681	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	682	l_grid, vpt, pt_reference, rif, tend, zu, &
	683	zw )
[1]	684	ENDIF
	685	ELSE
	686	IF ( use_upstream_for_tke ) THEN
	687	tend(:,j,i) = 0.0
	688	CALL advec_s_up( i, j, e )
	689	ELSE
	690	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	691	THEN
	692	tend(:,j,i) = 0.0
	693	CALL advec_s_pw( i, j, e )
	694	ELSE
	695	IF ( scalar_advec /= 'ups-scheme' ) THEN
	696	tend(:,j,i) = 0.0
	697	CALL advec_s_up( i, j, e )
	698	ENDIF
	699	ENDIF
	700	ENDIF
	701	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	702	THEN
[75]	703	IF ( .NOT. humidity ) THEN
[1]	704	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
[97]	705	km_m, l_grid, pt_m, pt_reference, &
	706	rif_m, tend, zu, zw )
[1]	707	ELSE
	708	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
[97]	709	km_m, l_grid, vpt_m, pt_reference, &
	710	rif_m, tend, zu, zw )
[1]	711	ENDIF
	712	ELSE
[75]	713	IF ( .NOT. humidity ) THEN
[97]	714	IF ( ocean ) THEN
	715	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	716	km, l_grid, rho, prho_reference, &
	717	rif, tend, zu, zw )
	718	ELSE
	719	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	720	km, l_grid, pt, pt_reference, rif, &
	721	tend, zu, zw )
	722	ENDIF
[1]	723	ELSE
	724	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
[97]	725	l_grid, vpt, pt_reference, rif, tend, &
	726	zu, zw )
[1]	727	ENDIF
	728	ENDIF
	729	ENDIF
	730	CALL production_e( i, j )
	731	CALL user_actions( i, j, 'e-tendency' )
	732
	733	!
	734	!-- Prognostic equation for TKE.
	735	!-- Eliminate negative TKE values, which can occur due to numerical
	736	!-- reasons in the course of the integration. In such cases the old TKE
	737	!-- value is reduced by 90%.
[19]	738	DO k = nzb_s_inner(j,i)+1, nzt
[1]	739	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
	740	dt_3d * ( &
	741	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	742	)
	743	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	744	ENDDO
	745
	746	!
	747	!-- Calculate tendencies for the next Runge-Kutta step
	748	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	749	IF ( intermediate_timestep_count == 1 ) THEN
[19]	750	DO k = nzb_s_inner(j,i)+1, nzt
[1]	751	te_m(k,j,i) = tend(k,j,i)
	752	ENDDO
	753	ELSEIF ( intermediate_timestep_count < &
	754	intermediate_timestep_count_max ) THEN
[19]	755	DO k = nzb_s_inner(j,i)+1, nzt
[1]	756	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	757	ENDDO
	758	ENDIF
	759	ENDIF
	760
	761	ENDDO
	762	ENDDO
	763
	764	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	765
	766	ENDIF
	767
	768
[63]	769	END SUBROUTINE prognostic_equations_noopt
[1]	770
	771
[63]	772	SUBROUTINE prognostic_equations_cache
[1]	773
	774	!------------------------------------------------------------------------------!
	775	! Version with one optimized loop over all equations. It is only allowed to
	776	! be called for the standard Piascek-Williams advection scheme.
	777	!
[95]	778	! Here the calls of most subroutines are embedded in two DO loops over i and j,
	779	! so communication between CPUs is not allowed (does not make sense) within
	780	! these loops.
[1]	781	!
	782	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
	783	!------------------------------------------------------------------------------!
	784
	785	IMPLICIT NONE
	786
	787	CHARACTER (LEN=9) :: time_to_string
	788	INTEGER :: i, j, k
	789
	790
	791	!
	792	!-- Time measurement can only be performed for the whole set of equations
	793	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
	794
	795
	796	!
	797	!-- Calculate those variables needed in the tendency terms which need
	798	!-- global communication
[96]	799	CALL calc_mean_profile( pt, 4 )
	800	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	801	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
[1]	802	IF ( .NOT. constant_diffusion ) CALL production_e_init
	803
	804
	805	!
	806	!-- Loop over all prognostic equations
	807	!$OMP PARALLEL private (i,j,k)
	808	!$OMP DO
[75]	809	DO i = nxl, nxr
	810	DO j = nys, nyn
[1]	811	!
	812	!-- Tendency terms for u-velocity component
[106]	813	IF ( .NOT. outflow_l .OR. i > nxl ) THEN
[1]	814
	815	tend(:,j,i) = 0.0
	816	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	817	CALL advec_u_pw( i, j )
	818	ELSE
	819	CALL advec_u_up( i, j )
	820	ENDIF
	821	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	822	THEN
	823	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, &
[102]	824	u_m, usws_m, uswst_m, v_m, w_m )
[1]	825	ELSE
	826	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
[102]	827	usws, uswst, v, w )
[1]	828	ENDIF
	829	CALL coriolis( i, j, 1 )
[97]	830	IF ( sloping_surface ) CALL buoyancy( i, j, pt, pt_reference, 1, &
	831	4 )
[1]	832	CALL user_actions( i, j, 'u-tendency' )
	833
	834	!
	835	!-- Prognostic equation for u-velocity component
	836	DO k = nzb_u_inner(j,i)+1, nzt
	837	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
	838	dt_3d * ( &
	839	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
	840	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
	841	) - &
	842	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	843	ENDDO
	844
	845	!
	846	!-- Calculate tendencies for the next Runge-Kutta step
	847	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	848	IF ( intermediate_timestep_count == 1 ) THEN
	849	DO k = nzb_u_inner(j,i)+1, nzt
	850	tu_m(k,j,i) = tend(k,j,i)
	851	ENDDO
	852	ELSEIF ( intermediate_timestep_count < &
	853	intermediate_timestep_count_max ) THEN
	854	DO k = nzb_u_inner(j,i)+1, nzt
	855	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	856	ENDDO
	857	ENDIF
	858	ENDIF
	859
	860	ENDIF
	861
	862	!
	863	!-- Tendency terms for v-velocity component
[106]	864	IF ( .NOT. outflow_s .OR. j > nys ) THEN
[1]	865
	866	tend(:,j,i) = 0.0
	867	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	868	CALL advec_v_pw( i, j )
	869	ELSE
	870	CALL advec_v_up( i, j )
	871	ENDIF
	872	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	873	THEN
	874	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, &
[102]	875	u_m, v_m, vsws_m, vswst_m, w_m )
[1]	876	ELSE
	877	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
[102]	878	vsws, vswst, w )
[1]	879	ENDIF
	880	CALL coriolis( i, j, 2 )
	881	CALL user_actions( i, j, 'v-tendency' )
	882
	883	!
	884	!-- Prognostic equation for v-velocity component
	885	DO k = nzb_v_inner(j,i)+1, nzt
	886	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
	887	dt_3d * ( &
	888	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
	889	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
	890	) - &
	891	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	892	ENDDO
	893
	894	!
	895	!-- Calculate tendencies for the next Runge-Kutta step
	896	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	897	IF ( intermediate_timestep_count == 1 ) THEN
	898	DO k = nzb_v_inner(j,i)+1, nzt
	899	tv_m(k,j,i) = tend(k,j,i)
	900	ENDDO
	901	ELSEIF ( intermediate_timestep_count < &
	902	intermediate_timestep_count_max ) THEN
	903	DO k = nzb_v_inner(j,i)+1, nzt
	904	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	905	ENDDO
	906	ENDIF
	907	ENDIF
	908
	909	ENDIF
	910
	911	!
	912	!-- Tendency terms for w-velocity component
[106]	913	tend(:,j,i) = 0.0
	914	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	915	CALL advec_w_pw( i, j )
	916	ELSE
	917	CALL advec_w_up( i, j )
	918	ENDIF
	919	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	920	THEN
	921	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, &
	922	km_damp_y, tend, u_m, v_m, w_m )
	923	ELSE
	924	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
	925	tend, u, v, w )
	926	ENDIF
	927	CALL coriolis( i, j, 3 )
	928	IF ( ocean ) THEN
	929	CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
	930	ELSE
	931	IF ( .NOT. humidity ) THEN
	932	CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
	933	ELSE
	934	CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
	935	ENDIF
	936	ENDIF
	937	CALL user_actions( i, j, 'w-tendency' )
[1]	938
[106]	939	!
	940	!-- Prognostic equation for w-velocity component
	941	DO k = nzb_w_inner(j,i)+1, nzt-1
	942	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
	943	dt_3d * ( &
	944	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
	945	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
	946	) - &
	947	tsc(5) * rdf(k) * w(k,j,i)
	948	ENDDO
	949
	950	!
	951	!-- Calculate tendencies for the next Runge-Kutta step
	952	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	953	IF ( intermediate_timestep_count == 1 ) THEN
	954	DO k = nzb_w_inner(j,i)+1, nzt-1
	955	tw_m(k,j,i) = tend(k,j,i)
	956	ENDDO
	957	ELSEIF ( intermediate_timestep_count < &
	958	intermediate_timestep_count_max ) THEN
	959	DO k = nzb_w_inner(j,i)+1, nzt-1
	960	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	961	ENDDO
	962	ENDIF
	963	ENDIF
	964
	965	!
	966	!-- Tendency terms for potential temperature
	967	tend(:,j,i) = 0.0
	968	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	969	CALL advec_s_pw( i, j, pt )
	970	ELSE
	971	CALL advec_s_up( i, j, pt )
	972	ENDIF
	973	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
	974	THEN
	975	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
	976	tswst_m, tend )
	977	ELSE
	978	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
	979	ENDIF
	980
	981	!
	982	!-- If required compute heating/cooling due to long wave radiation
	983	!-- processes
	984	IF ( radiation ) THEN
	985	CALL calc_radiation( i, j )
	986	ENDIF
	987
	988	!
	989	!-- If required compute impact of latent heat due to precipitation
	990	IF ( precipitation ) THEN
	991	CALL impact_of_latent_heat( i, j )
	992	ENDIF
	993	CALL user_actions( i, j, 'pt-tendency' )
	994
	995	!
	996	!-- Prognostic equation for potential temperature
	997	DO k = nzb_s_inner(j,i)+1, nzt
	998	pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
	999	dt_3d * ( &
	1000	tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
	1001	) - &
	1002	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
	1003	ENDDO
	1004
	1005	!
	1006	!-- Calculate tendencies for the next Runge-Kutta step
	1007	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1008	IF ( intermediate_timestep_count == 1 ) THEN
	1009	DO k = nzb_s_inner(j,i)+1, nzt
	1010	tpt_m(k,j,i) = tend(k,j,i)
	1011	ENDDO
	1012	ELSEIF ( intermediate_timestep_count < &
	1013	intermediate_timestep_count_max ) THEN
	1014	DO k = nzb_s_inner(j,i)+1, nzt
	1015	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1016	5.3125 * tpt_m(k,j,i)
	1017	ENDDO
	1018	ENDIF
	1019	ENDIF
	1020
	1021	!
	1022	!-- If required, compute prognostic equation for salinity
	1023	IF ( ocean ) THEN
	1024
	1025	!
	1026	!-- Tendency-terms for salinity
[1]	1027	tend(:,j,i) = 0.0
[106]	1028	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
[1]	1029	THEN
[106]	1030	CALL advec_s_pw( i, j, sa )
[1]	1031	ELSE
[106]	1032	CALL advec_s_up( i, j, sa )
[1]	1033	ENDIF
[106]	1034	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
	1035	tend )
[1]	1036
[106]	1037	CALL user_actions( i, j, 'sa-tendency' )
	1038
[1]	1039	!
[106]	1040	!-- Prognostic equation for salinity
	1041	DO k = nzb_s_inner(j,i)+1, nzt
	1042	sa_p(k,j,i) = tsc(1) * sa(k,j,i) + &
	1043	dt_3d * ( &
	1044	tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	1045	) - &
	1046	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
	1047	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
[1]	1048	ENDDO
	1049
	1050	!
	1051	!-- Calculate tendencies for the next Runge-Kutta step
	1052	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1053	IF ( intermediate_timestep_count == 1 ) THEN
[106]	1054	DO k = nzb_s_inner(j,i)+1, nzt
	1055	tsa_m(k,j,i) = tend(k,j,i)
[1]	1056	ENDDO
	1057	ELSEIF ( intermediate_timestep_count < &
	1058	intermediate_timestep_count_max ) THEN
[106]	1059	DO k = nzb_s_inner(j,i)+1, nzt
	1060	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1061	5.3125 * tsa_m(k,j,i)
[1]	1062	ENDDO
	1063	ENDIF
	1064	ENDIF
	1065
	1066	!
[106]	1067	!-- Calculate density by the equation of state for seawater
	1068	CALL eqn_state_seawater( i, j )
	1069
	1070	ENDIF
	1071
	1072	!
	1073	!-- If required, compute prognostic equation for total water content /
	1074	!-- scalar
	1075	IF ( humidity .OR. passive_scalar ) THEN
	1076
	1077	!
	1078	!-- Tendency-terms for total water content / scalar
[1]	1079	tend(:,j,i) = 0.0
[106]	1080	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	1081	THEN
	1082	CALL advec_s_pw( i, j, q )
[1]	1083	ELSE
[106]	1084	CALL advec_s_up( i, j, q )
[1]	1085	ENDIF
[106]	1086	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
[1]	1087	THEN
[106]	1088	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
	1089	qswst_m, tend )
[1]	1090	ELSE
[106]	1091	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
	1092	tend )
[1]	1093	ENDIF
[106]	1094
[1]	1095	!
[106]	1096	!-- If required compute decrease of total water content due to
	1097	!-- precipitation
[1]	1098	IF ( precipitation ) THEN
[106]	1099	CALL calc_precipitation( i, j )
[1]	1100	ENDIF
[106]	1101	CALL user_actions( i, j, 'q-tendency' )
[1]	1102
	1103	!
[106]	1104	!-- Prognostic equation for total water content / scalar
[19]	1105	DO k = nzb_s_inner(j,i)+1, nzt
[106]	1106	q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
	1107	dt_3d * ( &
	1108	tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	1109	) - &
	1110	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
	1111	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
[1]	1112	ENDDO
	1113
	1114	!
	1115	!-- Calculate tendencies for the next Runge-Kutta step
	1116	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1117	IF ( intermediate_timestep_count == 1 ) THEN
[19]	1118	DO k = nzb_s_inner(j,i)+1, nzt
[106]	1119	tq_m(k,j,i) = tend(k,j,i)
[1]	1120	ENDDO
	1121	ELSEIF ( intermediate_timestep_count < &
	1122	intermediate_timestep_count_max ) THEN
[19]	1123	DO k = nzb_s_inner(j,i)+1, nzt
[106]	1124	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1125	5.3125 * tq_m(k,j,i)
[1]	1126	ENDDO
	1127	ENDIF
	1128	ENDIF
	1129
[106]	1130	ENDIF
[95]	1131
	1132	!
[106]	1133	!-- If required, compute prognostic equation for turbulent kinetic
	1134	!-- energy (TKE)
	1135	IF ( .NOT. constant_diffusion ) THEN
[95]	1136
	1137	!
[106]	1138	!-- Tendency-terms for TKE
	1139	tend(:,j,i) = 0.0
	1140	IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
	1141	.AND. .NOT. use_upstream_for_tke ) THEN
	1142	CALL advec_s_pw( i, j, e )
	1143	ELSE
	1144	CALL advec_s_up( i, j, e )
[95]	1145	ENDIF
[106]	1146	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
	1147	THEN
	1148	IF ( .NOT. humidity ) THEN
	1149	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	1150	km_m, l_grid, pt_m, pt_reference, &
	1151	rif_m, tend, zu, zw )
[1]	1152	ELSE
[106]	1153	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
	1154	km_m, l_grid, vpt_m, pt_reference, &
	1155	rif_m, tend, zu, zw )
[1]	1156	ENDIF
[106]	1157	ELSE
	1158	IF ( .NOT. humidity ) THEN
	1159	IF ( ocean ) THEN
	1160	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	1161	km, l_grid, rho, prho_reference, &
	1162	rif, tend, zu, zw )
[1]	1163	ELSE
[106]	1164	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
	1165	km, l_grid, pt, pt_reference, rif, &
	1166	tend, zu, zw )
[1]	1167	ENDIF
	1168	ELSE
[106]	1169	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
	1170	l_grid, vpt, pt_reference, rif, tend, &
	1171	zu, zw )
[1]	1172	ENDIF
[106]	1173	ENDIF
	1174	CALL production_e( i, j )
	1175	CALL user_actions( i, j, 'e-tendency' )
[1]	1176
	1177	!
[106]	1178	!-- Prognostic equation for TKE.
	1179	!-- Eliminate negative TKE values, which can occur due to numerical
	1180	!-- reasons in the course of the integration. In such cases the old
	1181	!-- TKE value is reduced by 90%.
	1182	DO k = nzb_s_inner(j,i)+1, nzt
	1183	e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
	1184	dt_3d * ( &
	1185	tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	1186	)
	1187	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	1188	ENDDO
[1]	1189
	1190	!
[106]	1191	!-- Calculate tendencies for the next Runge-Kutta step
	1192	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1193	IF ( intermediate_timestep_count == 1 ) THEN
	1194	DO k = nzb_s_inner(j,i)+1, nzt
	1195	te_m(k,j,i) = tend(k,j,i)
	1196	ENDDO
	1197	ELSEIF ( intermediate_timestep_count < &
	1198	intermediate_timestep_count_max ) THEN
	1199	DO k = nzb_s_inner(j,i)+1, nzt
	1200	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1201	5.3125 * te_m(k,j,i)
	1202	ENDDO
[1]	1203	ENDIF
[106]	1204	ENDIF
[1]	1205
[106]	1206	ENDIF ! TKE equation
[1]	1207
	1208	ENDDO
	1209	ENDDO
	1210	!$OMP END PARALLEL
	1211
	1212	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
	1213
	1214
[63]	1215	END SUBROUTINE prognostic_equations_cache
[1]	1216
	1217
[63]	1218	SUBROUTINE prognostic_equations_vector
[1]	1219
	1220	!------------------------------------------------------------------------------!
	1221	! Version for vector machines
	1222	!------------------------------------------------------------------------------!
	1223
	1224	IMPLICIT NONE
	1225
	1226	CHARACTER (LEN=9) :: time_to_string
	1227	INTEGER :: i, j, k
	1228	REAL :: sat, sbt
	1229
	1230	!
	1231	!-- Calculate those variables needed in the tendency terms which need
	1232	!-- global communication
[96]	1233	CALL calc_mean_profile( pt, 4 )
	1234	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	1235	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
[1]	1236
	1237	!
	1238	!-- u-velocity component
	1239	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	1240
	1241	!
	1242	!-- u-tendency terms with communication
	1243	IF ( momentum_advec == 'ups-scheme' ) THEN
	1244	tend = 0.0
	1245	CALL advec_u_ups
	1246	ENDIF
	1247
	1248	!
	1249	!-- u-tendency terms with no communication
	1250	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1251	tend = 0.0
	1252	CALL advec_u_pw
	1253	ELSE
	1254	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1255	tend = 0.0
	1256	CALL advec_u_up
	1257	ENDIF
	1258	ENDIF
	1259	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[102]	1260	CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
	1261	uswst_m, v_m, w_m )
[1]	1262	ELSE
[102]	1263	CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
[1]	1264	ENDIF
	1265	CALL coriolis( 1 )
[97]	1266	IF ( sloping_surface ) CALL buoyancy( pt, pt_reference, 1, 4 )
[1]	1267	CALL user_actions( 'u-tendency' )
	1268
	1269	!
	1270	!-- Prognostic equation for u-velocity component
[106]	1271	DO i = nxlu, nxr
[1]	1272	DO j = nys, nyn
	1273	DO k = nzb_u_inner(j,i)+1, nzt
	1274	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
	1275	dt_3d * ( &
	1276	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
	1277	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
	1278	) - &
	1279	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	1280	ENDDO
	1281	ENDDO
	1282	ENDDO
	1283
	1284	!
	1285	!-- Calculate tendencies for the next Runge-Kutta step
	1286	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1287	IF ( intermediate_timestep_count == 1 ) THEN
[106]	1288	DO i = nxlu, nxr
[1]	1289	DO j = nys, nyn
	1290	DO k = nzb_u_inner(j,i)+1, nzt
	1291	tu_m(k,j,i) = tend(k,j,i)
	1292	ENDDO
	1293	ENDDO
	1294	ENDDO
	1295	ELSEIF ( intermediate_timestep_count < &
	1296	intermediate_timestep_count_max ) THEN
[106]	1297	DO i = nxlu, nxr
[1]	1298	DO j = nys, nyn
	1299	DO k = nzb_u_inner(j,i)+1, nzt
	1300	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	1301	ENDDO
	1302	ENDDO
	1303	ENDDO
	1304	ENDIF
	1305	ENDIF
	1306
	1307	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	1308
	1309	!
	1310	!-- v-velocity component
	1311	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	1312
	1313	!
	1314	!-- v-tendency terms with communication
	1315	IF ( momentum_advec == 'ups-scheme' ) THEN
	1316	tend = 0.0
	1317	CALL advec_v_ups
	1318	ENDIF
	1319
	1320	!
	1321	!-- v-tendency terms with no communication
	1322	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1323	tend = 0.0
	1324	CALL advec_v_pw
	1325	ELSE
	1326	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1327	tend = 0.0
	1328	CALL advec_v_up
	1329	ENDIF
	1330	ENDIF
	1331	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	1332	CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
[102]	1333	vswst_m, w_m )
[1]	1334	ELSE
[102]	1335	CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
[1]	1336	ENDIF
	1337	CALL coriolis( 2 )
	1338	CALL user_actions( 'v-tendency' )
	1339
	1340	!
	1341	!-- Prognostic equation for v-velocity component
	1342	DO i = nxl, nxr
[106]	1343	DO j = nysv, nyn
[1]	1344	DO k = nzb_v_inner(j,i)+1, nzt
	1345	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
	1346	dt_3d * ( &
	1347	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
	1348	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
	1349	) - &
	1350	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	1351	ENDDO
	1352	ENDDO
	1353	ENDDO
	1354
	1355	!
	1356	!-- Calculate tendencies for the next Runge-Kutta step
	1357	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1358	IF ( intermediate_timestep_count == 1 ) THEN
	1359	DO i = nxl, nxr
[106]	1360	DO j = nysv, nyn
[1]	1361	DO k = nzb_v_inner(j,i)+1, nzt
	1362	tv_m(k,j,i) = tend(k,j,i)
	1363	ENDDO
	1364	ENDDO
	1365	ENDDO
	1366	ELSEIF ( intermediate_timestep_count < &
	1367	intermediate_timestep_count_max ) THEN
	1368	DO i = nxl, nxr
[106]	1369	DO j = nysv, nyn
[1]	1370	DO k = nzb_v_inner(j,i)+1, nzt
	1371	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	1372	ENDDO
	1373	ENDDO
	1374	ENDDO
	1375	ENDIF
	1376	ENDIF
	1377
	1378	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1379
	1380	!
	1381	!-- w-velocity component
	1382	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1383
	1384	!
	1385	!-- w-tendency terms with communication
	1386	IF ( momentum_advec == 'ups-scheme' ) THEN
	1387	tend = 0.0
	1388	CALL advec_w_ups
	1389	ENDIF
	1390
	1391	!
	1392	!-- w-tendency terms with no communication
	1393	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1394	tend = 0.0
	1395	CALL advec_w_pw
	1396	ELSE
	1397	IF ( momentum_advec /= 'ups-scheme' ) THEN
	1398	tend = 0.0
	1399	CALL advec_w_up
	1400	ENDIF
	1401	ENDIF
	1402	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
	1403	CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
[57]	1404	v_m, w_m )
[1]	1405	ELSE
[57]	1406	CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
[1]	1407	ENDIF
	1408	CALL coriolis( 3 )
[97]	1409	IF ( ocean ) THEN
	1410	CALL buoyancy( rho, prho_reference, 3, 64 )
[1]	1411	ELSE
[97]	1412	IF ( .NOT. humidity ) THEN
	1413	CALL buoyancy( pt, pt_reference, 3, 4 )
	1414	ELSE
	1415	CALL buoyancy( vpt, pt_reference, 3, 44 )
	1416	ENDIF
[1]	1417	ENDIF
	1418	CALL user_actions( 'w-tendency' )
	1419
	1420	!
	1421	!-- Prognostic equation for w-velocity component
	1422	DO i = nxl, nxr
	1423	DO j = nys, nyn
	1424	DO k = nzb_w_inner(j,i)+1, nzt-1
	1425	w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
	1426	dt_3d * ( &
	1427	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
	1428	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
	1429	) - &
	1430	tsc(5) * rdf(k) * w(k,j,i)
	1431	ENDDO
	1432	ENDDO
	1433	ENDDO
	1434
	1435	!
	1436	!-- Calculate tendencies for the next Runge-Kutta step
	1437	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1438	IF ( intermediate_timestep_count == 1 ) THEN
	1439	DO i = nxl, nxr
	1440	DO j = nys, nyn
	1441	DO k = nzb_w_inner(j,i)+1, nzt-1
	1442	tw_m(k,j,i) = tend(k,j,i)
	1443	ENDDO
	1444	ENDDO
	1445	ENDDO
	1446	ELSEIF ( intermediate_timestep_count < &
	1447	intermediate_timestep_count_max ) THEN
	1448	DO i = nxl, nxr
	1449	DO j = nys, nyn
	1450	DO k = nzb_w_inner(j,i)+1, nzt-1
	1451	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1452	ENDDO
	1453	ENDDO
	1454	ENDDO
	1455	ENDIF
	1456	ENDIF
	1457
	1458	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1459
	1460	!
	1461	!-- potential temperature
	1462	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1463
	1464	!
	1465	!-- pt-tendency terms with communication
[63]	1466	sat = tsc(1)
	1467	sbt = tsc(2)
[1]	1468	IF ( scalar_advec == 'bc-scheme' ) THEN
[63]	1469
	1470	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[1]	1471	!
[63]	1472	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1473	!-- switched on. Thus:
	1474	sat = 1.0
	1475	sbt = 1.0
	1476	ENDIF
[1]	1477	tend = 0.0
	1478	CALL advec_s_bc( pt, 'pt' )
	1479	ELSE
	1480	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
	1481	tend = 0.0
	1482	CALL advec_s_ups( pt, 'pt' )
	1483	ENDIF
	1484	ENDIF
	1485
	1486	!
	1487	!-- pt-tendency terms with no communication
	1488	IF ( scalar_advec == 'bc-scheme' ) THEN
[19]	1489	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
[1]	1490	ELSE
	1491	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1492	tend = 0.0
	1493	CALL advec_s_pw( pt )
	1494	ELSE
	1495	IF ( scalar_advec /= 'ups-scheme' ) THEN
	1496	tend = 0.0
	1497	CALL advec_s_up( pt )
	1498	ENDIF
	1499	ENDIF
	1500	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[19]	1501	CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, tend )
[1]	1502	ELSE
[19]	1503	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
[1]	1504	ENDIF
	1505	ENDIF
	1506
	1507	!
	1508	!-- If required compute heating/cooling due to long wave radiation
	1509	!-- processes
	1510	IF ( radiation ) THEN
	1511	CALL calc_radiation
	1512	ENDIF
	1513
	1514	!
	1515	!-- If required compute impact of latent heat due to precipitation
	1516	IF ( precipitation ) THEN
	1517	CALL impact_of_latent_heat
	1518	ENDIF
	1519	CALL user_actions( 'pt-tendency' )
	1520
	1521	!
	1522	!-- Prognostic equation for potential temperature
	1523	DO i = nxl, nxr
	1524	DO j = nys, nyn
[19]	1525	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1526	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
	1527	dt_3d * ( &
	1528	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
	1529	) - &
	1530	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
	1531	ENDDO
	1532	ENDDO
	1533	ENDDO
	1534
	1535	!
	1536	!-- Calculate tendencies for the next Runge-Kutta step
	1537	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1538	IF ( intermediate_timestep_count == 1 ) THEN
	1539	DO i = nxl, nxr
	1540	DO j = nys, nyn
[19]	1541	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1542	tpt_m(k,j,i) = tend(k,j,i)
	1543	ENDDO
	1544	ENDDO
	1545	ENDDO
	1546	ELSEIF ( intermediate_timestep_count < &
	1547	intermediate_timestep_count_max ) THEN
	1548	DO i = nxl, nxr
	1549	DO j = nys, nyn
[19]	1550	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1551	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
	1552	ENDDO
	1553	ENDDO
	1554	ENDDO
	1555	ENDIF
	1556	ENDIF
	1557
	1558	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1559
	1560	!
[95]	1561	!-- If required, compute prognostic equation for salinity
	1562	IF ( ocean ) THEN
	1563
	1564	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	1565
	1566	!
	1567	!-- sa-tendency terms with communication
	1568	sat = tsc(1)
	1569	sbt = tsc(2)
	1570	IF ( scalar_advec == 'bc-scheme' ) THEN
	1571
	1572	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1573	!
	1574	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1575	!-- switched on. Thus:
	1576	sat = 1.0
	1577	sbt = 1.0
	1578	ENDIF
	1579	tend = 0.0
	1580	CALL advec_s_bc( sa, 'sa' )
	1581	ELSE
	1582	IF ( tsc(2) /= 2.0 ) THEN
	1583	IF ( scalar_advec == 'ups-scheme' ) THEN
	1584	tend = 0.0
	1585	CALL advec_s_ups( sa, 'sa' )
	1586	ENDIF
	1587	ENDIF
	1588	ENDIF
	1589
	1590	!
	1591	!-- Scalar/q-tendency terms with no communication
	1592	IF ( scalar_advec == 'bc-scheme' ) THEN
	1593	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, tend )
	1594	ELSE
	1595	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1596	tend = 0.0
	1597	CALL advec_s_pw( sa )
	1598	ELSE
	1599	IF ( scalar_advec /= 'ups-scheme' ) THEN
	1600	tend = 0.0
	1601	CALL advec_s_up( sa )
	1602	ENDIF
	1603	ENDIF
	1604	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, tend )
	1605	ENDIF
	1606
	1607	CALL user_actions( 'sa-tendency' )
	1608
	1609	!
	1610	!-- Prognostic equation for salinity
	1611	DO i = nxl, nxr
	1612	DO j = nys, nyn
	1613	DO k = nzb_s_inner(j,i)+1, nzt
	1614	sa_p(k,j,i) = sat * sa(k,j,i) + &
	1615	dt_3d * ( &
	1616	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
	1617	) - &
	1618	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
	1619	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1620	ENDDO
	1621	ENDDO
	1622	ENDDO
	1623
	1624	!
	1625	!-- Calculate tendencies for the next Runge-Kutta step
	1626	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1627	IF ( intermediate_timestep_count == 1 ) THEN
	1628	DO i = nxl, nxr
	1629	DO j = nys, nyn
	1630	DO k = nzb_s_inner(j,i)+1, nzt
	1631	tsa_m(k,j,i) = tend(k,j,i)
	1632	ENDDO
	1633	ENDDO
	1634	ENDDO
	1635	ELSEIF ( intermediate_timestep_count < &
	1636	intermediate_timestep_count_max ) 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) = -9.5625 * tend(k,j,i) + &
	1641	5.3125 * tsa_m(k,j,i)
	1642	ENDDO
	1643	ENDDO
	1644	ENDDO
	1645	ENDIF
	1646	ENDIF
	1647
	1648	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	1649
[96]	1650	!
	1651	!-- Calculate density by the equation of state for seawater
	1652	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	1653	CALL eqn_state_seawater
	1654	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	1655
[95]	1656	ENDIF
	1657
	1658	!
[1]	1659	!-- If required, compute prognostic equation for total water content / scalar
[75]	1660	IF ( humidity .OR. passive_scalar ) THEN
[1]	1661
	1662	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	1663
	1664	!
	1665	!-- Scalar/q-tendency terms with communication
[63]	1666	sat = tsc(1)
	1667	sbt = tsc(2)
[1]	1668	IF ( scalar_advec == 'bc-scheme' ) THEN
[63]	1669
	1670	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[1]	1671	!
[63]	1672	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1673	!-- switched on. Thus:
	1674	sat = 1.0
	1675	sbt = 1.0
	1676	ENDIF
[1]	1677	tend = 0.0
	1678	CALL advec_s_bc( q, 'q' )
	1679	ELSE
	1680	IF ( tsc(2) /= 2.0 ) THEN
	1681	IF ( scalar_advec == 'ups-scheme' ) THEN
	1682	tend = 0.0
	1683	CALL advec_s_ups( q, 'q' )
	1684	ENDIF
	1685	ENDIF
	1686	ENDIF
	1687
	1688	!
	1689	!-- Scalar/q-tendency terms with no communication
	1690	IF ( scalar_advec == 'bc-scheme' ) THEN
[19]	1691	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
[1]	1692	ELSE
	1693	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1694	tend = 0.0
	1695	CALL advec_s_pw( q )
	1696	ELSE
	1697	IF ( scalar_advec /= 'ups-scheme' ) THEN
	1698	tend = 0.0
	1699	CALL advec_s_up( q )
	1700	ENDIF
	1701	ENDIF
	1702	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[19]	1703	CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, tend )
[1]	1704	ELSE
[19]	1705	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
[1]	1706	ENDIF
	1707	ENDIF
	1708
	1709	!
	1710	!-- If required compute decrease of total water content due to
	1711	!-- precipitation
	1712	IF ( precipitation ) THEN
	1713	CALL calc_precipitation
	1714	ENDIF
	1715	CALL user_actions( 'q-tendency' )
	1716
	1717	!
	1718	!-- Prognostic equation for total water content / scalar
	1719	DO i = nxl, nxr
	1720	DO j = nys, nyn
[19]	1721	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1722	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
	1723	dt_3d * ( &
	1724	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
	1725	) - &
	1726	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
[73]	1727	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
[1]	1728	ENDDO
	1729	ENDDO
	1730	ENDDO
	1731
	1732	!
	1733	!-- Calculate tendencies for the next Runge-Kutta step
	1734	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1735	IF ( intermediate_timestep_count == 1 ) THEN
	1736	DO i = nxl, nxr
	1737	DO j = nys, nyn
[19]	1738	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1739	tq_m(k,j,i) = tend(k,j,i)
	1740	ENDDO
	1741	ENDDO
	1742	ENDDO
	1743	ELSEIF ( intermediate_timestep_count < &
	1744	intermediate_timestep_count_max ) THEN
	1745	DO i = nxl, nxr
	1746	DO j = nys, nyn
[19]	1747	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1748	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	1749	ENDDO
	1750	ENDDO
	1751	ENDDO
	1752	ENDIF
	1753	ENDIF
	1754
	1755	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	1756
	1757	ENDIF
	1758
	1759	!
	1760	!-- If required, compute prognostic equation for turbulent kinetic
	1761	!-- energy (TKE)
	1762	IF ( .NOT. constant_diffusion ) THEN
	1763
	1764	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	1765
	1766	!
	1767	!-- TKE-tendency terms with communication
	1768	CALL production_e_init
[63]	1769
	1770	sat = tsc(1)
	1771	sbt = tsc(2)
[1]	1772	IF ( .NOT. use_upstream_for_tke ) THEN
	1773	IF ( scalar_advec == 'bc-scheme' ) THEN
[63]	1774
	1775	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[1]	1776	!
[63]	1777	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
	1778	!-- switched on. Thus:
	1779	sat = 1.0
	1780	sbt = 1.0
	1781	ENDIF
[1]	1782	tend = 0.0
	1783	CALL advec_s_bc( e, 'e' )
	1784	ELSE
	1785	IF ( tsc(2) /= 2.0 ) THEN
	1786	IF ( scalar_advec == 'ups-scheme' ) THEN
	1787	tend = 0.0
	1788	CALL advec_s_ups( e, 'e' )
	1789	ENDIF
	1790	ENDIF
	1791	ENDIF
	1792	ENDIF
	1793
	1794	!
	1795	!-- TKE-tendency terms with no communication
	1796	IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) &
	1797	THEN
[75]	1798	IF ( .NOT. humidity ) THEN
[97]	1799	IF ( ocean ) THEN
	1800	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, rho, &
	1801	prho_reference, rif, tend, zu, zw )
	1802	ELSE
	1803	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
	1804	pt_reference, rif, tend, zu, zw )
	1805	ENDIF
[1]	1806	ELSE
	1807	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
[97]	1808	pt_reference, rif, tend, zu, zw )
[1]	1809	ENDIF
	1810	ELSE
	1811	IF ( use_upstream_for_tke ) THEN
	1812	tend = 0.0
	1813	CALL advec_s_up( e )
	1814	ELSE
	1815	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
	1816	tend = 0.0
	1817	CALL advec_s_pw( e )
	1818	ELSE
	1819	IF ( scalar_advec /= 'ups-scheme' ) THEN
	1820	tend = 0.0
	1821	CALL advec_s_up( e )
	1822	ENDIF
	1823	ENDIF
	1824	ENDIF
	1825	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
[75]	1826	IF ( .NOT. humidity ) THEN
[1]	1827	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
[97]	1828	pt_m, pt_reference, rif_m, tend, zu, zw )
[1]	1829	ELSE
	1830	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
[97]	1831	vpt_m, pt_reference, rif_m, tend, zu, zw )
[1]	1832	ENDIF
	1833	ELSE
[75]	1834	IF ( .NOT. humidity ) THEN
[97]	1835	IF ( ocean ) THEN
	1836	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	1837	rho, prho_reference, rif, tend, zu, zw )
	1838	ELSE
	1839	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	1840	pt, pt_reference, rif, tend, zu, zw )
	1841	ENDIF
[1]	1842	ELSE
	1843	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
[97]	1844	pt_reference, rif, tend, zu, zw )
[1]	1845	ENDIF
	1846	ENDIF
	1847	ENDIF
	1848	CALL production_e
	1849	CALL user_actions( 'e-tendency' )
	1850
	1851	!
	1852	!-- Prognostic equation for TKE.
	1853	!-- Eliminate negative TKE values, which can occur due to numerical
	1854	!-- reasons in the course of the integration. In such cases the old TKE
	1855	!-- value is reduced by 90%.
	1856	DO i = nxl, nxr
	1857	DO j = nys, nyn
[19]	1858	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1859	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
	1860	dt_3d * ( &
	1861	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
	1862	)
	1863	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	1864	ENDDO
	1865	ENDDO
	1866	ENDDO
	1867
	1868	!
	1869	!-- Calculate tendencies for the next Runge-Kutta step
	1870	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1871	IF ( intermediate_timestep_count == 1 ) THEN
	1872	DO i = nxl, nxr
	1873	DO j = nys, nyn
[19]	1874	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1875	te_m(k,j,i) = tend(k,j,i)
	1876	ENDDO
	1877	ENDDO
	1878	ENDDO
	1879	ELSEIF ( intermediate_timestep_count < &
	1880	intermediate_timestep_count_max ) THEN
	1881	DO i = nxl, nxr
	1882	DO j = nys, nyn
[19]	1883	DO k = nzb_s_inner(j,i)+1, nzt
[1]	1884	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	1885	ENDDO
	1886	ENDDO
	1887	ENDDO
	1888	ENDIF
	1889	ENDIF
	1890
	1891	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	1892
	1893	ENDIF
	1894
	1895
[63]	1896	END SUBROUTINE prognostic_equations_vector
[1]	1897
	1898
	1899	END MODULE prognostic_equations_mod

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

Context Navigation

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

Download in other formats: