Home

Context Navigation

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

Last change on this file since 1035 was 1020, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 59.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |