Home

Context Navigation

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

Last change on this file since 1001 was 1001, checked in by raasch, 12 years ago
leapfrog timestep scheme and upstream-spline advection scheme completely removed from the code, reading of dt_fixed from restart file removed
Property svn:keywords set to `Id`
File size: 62.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |