Home

Context Navigation

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

Last change on this file since 172 was 153, checked in by steinfeld, 17 years ago
Update for the plant canopy model
Property svn:keywords set to `Id`
File size: 68.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |