Home

Context Navigation

source: palm/trunk/SOURCE/prandtl_fluxes.f90 @ 1621

Last change on this file since 1621 was 1552, checked in by maronga, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 19.5 KB

Line
1	SUBROUTINE prandtl_fluxes
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: prandtl_fluxes.f90 1552 2015-03-03 14:27:15Z heinze $
27	!
28	! 1551 2015-03-03 14:18:16Z maronga
29	! Removed land surface model part. The surface fluxes are now always calculated
30	! within prandtl_fluxes, based on the given surface temperature/humidity (which
31	! is either provided by the land surface model, by large scale forcing data, or
32	! directly prescribed by the user.
33	!
34	! 1496 2014-12-02 17:25:50Z maronga
35	! Adapted for land surface model
36	!
37	! 1494 2014-11-21 17:14:03Z maronga
38	! Bugfixes: qs is now calculated before calculation of Rif. Ccalculation of
39	! buoyancy flux in Rif corrected (added missing humidity term), allow use of
40	! topography for coupled runs (not tested)
41	!
42	! 1361 2014-04-16 15:17:48Z hoffmann
43	! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
44	! drop concentration (nrsws) added
45	!
46	! 1340 2014-03-25 19:45:13Z kanani
47	! REAL constants defined as wp-kind
48	!
49	! 1320 2014-03-20 08:40:49Z raasch
50	! ONLY-attribute added to USE-statements,
51	! kind-parameters added to all INTEGER and REAL declaration statements,
52	! kinds are defined in new module kinds,
53	! old module precision_kind is removed,
54	! revision history before 2012 removed,
55	! comment fields (!:) to be used for variable explanations added to
56	! all variable declaration statements
57	!
58	! 1276 2014-01-15 13:40:41Z heinze
59	! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
60	!
61	! 1257 2013-11-08 15:18:40Z raasch
62	! openACC "kernels do" replaced by "kernels loop", "loop independent" added
63	!
64	! 1036 2012-10-22 13:43:42Z raasch
65	! code put under GPL (PALM 3.9)
66	!
67	! 1015 2012-09-27 09:23:24Z raasch
68	! OpenACC statements added
69	!
70	! 978 2012-08-09 08:28:32Z fricke
71	! roughness length for scalar quantities z0h added
72	!
73	! Revision 1.1 1998/01/23 10:06:06 raasch
74	! Initial revision
75	!
76	!
77	! Description:
78	! ------------
79	! Diagnostic computation of vertical fluxes in the Prandtl layer from the
80	! values of the variables at grid point k=1
81	!------------------------------------------------------------------------------!
82
83	USE arrays_3d, &
84	ONLY: e, nr, nrs, nrsws, pt, q, qr, qrs, qrsws, qs, qsws, rif, shf, &
85	ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
86
87	USE control_parameters, &
88	ONLY: cloud_physics, constant_heatflux, constant_waterflux, &
89	coupling_mode, g, humidity, ibc_e_b, icloud_scheme, kappa, &
90	large_scale_forcing, lsf_surf, passive_scalar, precipitation, &
91	pt_surface, q_surface, rif_max, rif_min, run_coupled, &
92	surface_pressure
93
94	USE indices, &
95	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner, &
96	nzb_u_inner, nzb_v_inner
97
98	USE kinds
99
100	IMPLICIT NONE
101
102	INTEGER(iwp) :: i !:
103	INTEGER(iwp) :: j !:
104	INTEGER(iwp) :: k !:
105
106	LOGICAL :: coupled_run !:
107
108	REAL(wp) :: a !:
109	REAL(wp) :: b !:
110	REAL(wp) :: e_q !:
111	REAL(wp) :: rifm !:
112	REAL(wp) :: uv_total !:
113	REAL(wp) :: z_p !:
114
115	!
116	!-- Data information for accelerators
117	!$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt ) &
118	!$acc present( q, qs, qsws, qrsws, rif, shf, ts, u, us, usws, v ) &
119	!$acc present( vpt, vsws, zu, zw, z0, z0h )
120	!
121	!-- Compute theta*
122	IF ( constant_heatflux ) THEN
123
124	!
125	!-- For a given heat flux in the Prandtl layer:
126	!-- for u* use the value from the previous time step
127	!$OMP PARALLEL DO
128	!$acc kernels loop
129	DO i = nxlg, nxrg
130	DO j = nysg, nyng
131	ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
132	!
133	!-- ts must be limited, because otherwise overflow may occur in case of
134	!-- us=0 when computing rif further below
135	IF ( ts(j,i) < -1.05E5_wp ) ts(j,i) = -1.0E5_wp
136	IF ( ts(j,i) > 1.0E5_wp ) ts(j,i) = 1.0E5_wp
137	ENDDO
138	ENDDO
139
140	ELSE
141	!
142	!-- For a given surface temperature:
143	!-- (the Richardson number is still the one from the previous time step)
144	IF ( large_scale_forcing .AND. lsf_surf ) THEN
145	!$OMP PARALLEL DO
146	!$acc kernels loop
147	DO i = nxlg, nxrg
148	DO j = nysg, nyng
149	k = nzb_s_inner(j,i)
150	pt(k,j,i) = pt_surface
151	ENDDO
152	ENDDO
153	ENDIF
154
155	!$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
156	!$acc kernels loop
157	DO i = nxlg, nxrg
158	DO j = nysg, nyng
159
160	k = nzb_s_inner(j,i)
161	z_p = zu(k+1) - zw(k)
162
163	IF ( rif(j,i) >= 0.0_wp ) THEN
164	!
165	!-- Stable stratification
166	ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &
167	LOG( z_p / z0h(j,i) ) + &
168	5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
169	)
170	ELSE
171	!
172	!-- Unstable stratification
173	a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) )
174	b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p )
175
176	ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &
177	LOG( z_p / z0h(j,i) ) - &
178	2.0_wp * LOG( ( 1.0_wp + a ) / ( 1.0_wp + b ) ) )
179	ENDIF
180
181	ENDDO
182	ENDDO
183	ENDIF
184
185	!
186	!-- If required compute q*
187	IF ( humidity .OR. passive_scalar ) THEN
188	IF ( constant_waterflux ) THEN
189	!
190	!-- For a given water flux in the Prandtl layer:
191	!$OMP PARALLEL DO
192	!$acc kernels loop
193	DO i = nxlg, nxrg
194	DO j = nysg, nyng
195	qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
196	ENDDO
197	ENDDO
198
199	ELSE
200	coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled )
201
202	IF ( large_scale_forcing .AND. lsf_surf ) THEN
203	!$OMP PARALLEL DO
204	!$acc kernels loop
205	DO i = nxlg, nxrg
206	DO j = nysg, nyng
207	k = nzb_s_inner(j,i)
208	q(k,j,i) = q_surface
209	ENDDO
210	ENDDO
211	ENDIF
212
213	!$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
214	!$acc kernels loop independent
215	DO i = nxlg, nxrg
216	!$acc loop independent
217	DO j = nysg, nyng
218
219	k = nzb_s_inner(j,i)
220	z_p = zu(k+1) - zw(k)
221
222	!
223	!-- Assume saturation for atmosphere coupled to ocean (but not
224	!-- in case of precursor runs)
225	IF ( coupled_run ) THEN
226	e_q = 6.1_wp * &
227	EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i)) - 273.15_wp ) )
228	q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q )
229	ENDIF
230	IF ( rif(j,i) >= 0.0_wp ) THEN
231	!
232	!-- Stable stratification
233	qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &
234	LOG( z_p / z0h(j,i) ) + &
235	5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
236	)
237	ELSE
238	!
239	!-- Unstable stratification
240	a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) )
241	b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p )
242
243	qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &
244	LOG( z_p / z0h(j,i) ) - &
245	2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) )
246	ENDIF
247
248	ENDDO
249	ENDDO
250	ENDIF
251	ENDIF
252
253	!
254	!-- Compute z_p/L (corresponds to the Richardson-flux number)
255	IF ( .NOT. humidity ) THEN
256	!$OMP PARALLEL DO PRIVATE( k, z_p )
257	!$acc kernels loop
258	DO i = nxlg, nxrg
259	DO j = nysg, nyng
260	k = nzb_s_inner(j,i)
261	z_p = zu(k+1) - zw(k)
262	rif(j,i) = z_p * kappa * g * ts(j,i) / &
263	( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
264	!
265	!-- Limit the value range of the Richardson numbers.
266	!-- This is necessary for very small velocities (u,v --> 0), because
267	!-- the absolute value of rif can then become very large, which in
268	!-- consequence would result in very large shear stresses and very
269	!-- small momentum fluxes (both are generally unrealistic).
270	IF ( rif(j,i) < rif_min ) rif(j,i) = rif_min
271	IF ( rif(j,i) > rif_max ) rif(j,i) = rif_max
272	ENDDO
273	ENDDO
274	ELSE
275	!$OMP PARALLEL DO PRIVATE( k, z_p )
276	!$acc kernels loop
277	DO i = nxlg, nxrg
278	DO j = nysg, nyng
279	k = nzb_s_inner(j,i)
280	z_p = zu(k+1) - zw(k)
281	rif(j,i) = z_p * kappa * g * &
282	( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) + 0.61_wp &
283	* q(k+1,j,i) * ts(j,i)) / &
284	( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
285	!
286	!-- Limit the value range of the Richardson numbers.
287	!-- This is necessary for very small velocities (u,v --> 0), because
288	!-- the absolute value of rif can then become very large, which in
289	!-- consequence would result in very large shear stresses and very
290	!-- small momentum fluxes (both are generally unrealistic).
291	IF ( rif(j,i) < rif_min ) rif(j,i) = rif_min
292	IF ( rif(j,i) > rif_max ) rif(j,i) = rif_max
293	ENDDO
294	ENDDO
295	ENDIF
296
297	!
298	!-- Compute u* at the scalars' grid points
299	!$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
300	!$acc kernels loop
301	DO i = nxl, nxr
302	DO j = nys, nyn
303
304	k = nzb_s_inner(j,i)
305	z_p = zu(k+1) - zw(k)
306
307	!
308	!-- Compute the absolute value of the horizontal velocity
309	!-- (relative to the surface)
310	uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) &
311	- u(k,j,i) - u(k,j,i+1) ) )**2 + &
312	( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) &
313	- v(k,j,i) - v(k,j+1,i) ) )**2 )
314
315
316	IF ( rif(j,i) >= 0.0_wp ) THEN
317	!
318	!-- Stable stratification
319	us(j,i) = kappa * uv_total / ( &
320	LOG( z_p / z0(j,i) ) + &
321	5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
322	)
323	ELSE
324	!
325	!-- Unstable stratification
326	a = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) )
327	b = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) / z_p * z0(j,i) ) )
328
329	us(j,i) = kappa * uv_total / ( &
330	LOG( z_p / z0(j,i) ) - &
331	LOG( ( 1.0_wp + a )*2 ( 1.0_wp + a**2 ) / ( &
332	( 1.0_wp + b )*2 ( 1.0_wp + b**2 ) ) ) + &
333	2.0_wp * ( ATAN( a ) - ATAN( b ) ) &
334	)
335	ENDIF
336	ENDDO
337	ENDDO
338
339	!
340	!-- Values of us at ghost point locations are needed for the evaluation of usws
341	!-- and vsws.
342	!$acc update host( us )
343	CALL exchange_horiz_2d( us )
344	!$acc update device( us )
345
346	!
347	!-- Compute u'w' for the total model domain.
348	!-- First compute the corresponding component of u* and square it.
349	!$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
350	!$acc kernels loop
351	DO i = nxl, nxr
352	DO j = nys, nyn
353
354	k = nzb_u_inner(j,i)
355	z_p = zu(k+1) - zw(k)
356
357	!
358	!-- Compute Richardson-flux number for this point
359	rifm = 0.5_wp * ( rif(j,i-1) + rif(j,i) )
360	IF ( rifm >= 0.0_wp ) THEN
361	!
362	!-- Stable stratification
363	usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ ( &
364	LOG( z_p / z0(j,i) ) + &
365	5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
366	)
367	ELSE
368	!
369	!-- Unstable stratification
370	a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
371	b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
372
373	usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / ( &
374	LOG( z_p / z0(j,i) ) - &
375	LOG( (1.0_wp + a )*2 ( 1.0_wp + a**2 ) / ( &
376	(1.0_wp + b )*2 ( 1.0_wp + b**2 ) ) ) + &
377	2.0_wp * ( ATAN( a ) - ATAN( b ) ) &
378	)
379	ENDIF
380	usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
381	ENDDO
382	ENDDO
383
384	!
385	!-- Compute v'w' for the total model domain.
386	!-- First compute the corresponding component of u* and square it.
387	!$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
388	!$acc kernels loop
389	DO i = nxl, nxr
390	DO j = nys, nyn
391
392	k = nzb_v_inner(j,i)
393	z_p = zu(k+1) - zw(k)
394
395	!
396	!-- Compute Richardson-flux number for this point
397	rifm = 0.5_wp * ( rif(j-1,i) + rif(j,i) )
398	IF ( rifm >= 0.0_wp ) THEN
399	!
400	!-- Stable stratification
401	vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / ( &
402	LOG( z_p / z0(j,i) ) + &
403	5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
404	)
405	ELSE
406	!
407	!-- Unstable stratification
408	a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
409	b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
410
411	vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / ( &
412	LOG( z_p / z0(j,i) ) - &
413	LOG( (1.0_wp + a )*2 ( 1.0_wp + a**2 ) / ( &
414	(1.0_wp + b )*2 ( 1.0_wp + b**2 ) ) ) + &
415	2.0_wp * ( ATAN( a ) - ATAN( b ) ) &
416	)
417	ENDIF
418	vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j-1,i) + us(j,i) )
419	ENDDO
420	ENDDO
421
422	!
423	!-- If required compute qr* and nr*
424	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN
425
426	!$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
427	!$acc kernels loop independent
428	DO i = nxlg, nxrg
429	!$acc loop independent
430	DO j = nysg, nyng
431
432	k = nzb_s_inner(j,i)
433	z_p = zu(k+1) - zw(k)
434
435	IF ( rif(j,i) >= 0.0 ) THEN
436	!
437	!-- Stable stratification
438	qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / ( &
439	LOG( z_p / z0h(j,i) ) + &
440	5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
441	nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / ( &
442	LOG( z_p / z0h(j,i) ) + &
443	5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
444
445	ELSE
446	!
447	!-- Unstable stratification
448	a = SQRT( 1.0 - 16.0 * rif(j,i) )
449	b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p )
450
451	qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / ( &
452	LOG( z_p / z0h(j,i) ) - &
453	2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
454	nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / ( &
455	LOG( z_p / z0h(j,i) ) - &
456	2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
457
458	ENDIF
459
460	ENDDO
461	ENDDO
462
463	ENDIF
464
465	!
466	!-- Exchange the boundaries for the momentum fluxes (only for sake of
467	!-- completeness)
468	!$acc update host( usws, vsws )
469	CALL exchange_horiz_2d( usws )
470	CALL exchange_horiz_2d( vsws )
471	!$acc update device( usws, vsws )
472	IF ( humidity .OR. passive_scalar ) THEN
473	!$acc update host( qsws )
474	CALL exchange_horiz_2d( qsws )
475	!$acc update device( qsws )
476	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
477	precipitation ) THEN
478	!$acc update host( qrsws, nrsws )
479	CALL exchange_horiz_2d( qrsws )
480	CALL exchange_horiz_2d( nrsws )
481	!$acc update device( qrsws, nrsws )
482	ENDIF
483	ENDIF
484
485	!
486	!-- Compute the vertical kinematic heat flux
487	IF ( .NOT. constant_heatflux ) THEN
488	!$OMP PARALLEL DO
489	!$acc kernels loop independent
490	DO i = nxlg, nxrg
491	!$acc loop independent
492	DO j = nysg, nyng
493	shf(j,i) = -ts(j,i) * us(j,i)
494	ENDDO
495	ENDDO
496	ENDIF
497
498	!
499	!-- Compute the vertical water/scalar flux
500	IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar ) ) THEN
501	!$OMP PARALLEL DO
502	!$acc kernels loop independent
503	DO i = nxlg, nxrg
504	!$acc loop independent
505	DO j = nysg, nyng
506	qsws(j,i) = -qs(j,i) * us(j,i)
507	ENDDO
508	ENDDO
509	ENDIF
510
511	!
512	!-- Compute (turbulent) fluxes of rain water content and rain drop concentartion
513	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN
514	!$OMP PARALLEL DO
515	!$acc kernels loop independent
516	DO i = nxlg, nxrg
517	!$acc loop independent
518	DO j = nysg, nyng
519	qrsws(j,i) = -qrs(j,i) * us(j,i)
520	nrsws(j,i) = -nrs(j,i) * us(j,i)
521	ENDDO
522	ENDDO
523	ENDIF
524
525	!
526	!-- Bottom boundary condition for the TKE
527	IF ( ibc_e_b == 2 ) THEN
528	!$OMP PARALLEL DO
529	!$acc kernels loop independent
530	DO i = nxlg, nxrg
531	!$acc loop independent
532	DO j = nysg, nyng
533	e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2
534	!
535	!-- As a test: cm = 0.4
536	! e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2
537	e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
538	ENDDO
539	ENDDO
540	ENDIF
541
542	!$acc end data
543
544	END SUBROUTINE prandtl_fluxes

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |