Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |