Home

Context Navigation

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

Last change on this file since 108 was 108, checked in by letzel, 17 years ago

Improved coupler: evaporation - salinity-flux coupling for humidity = .T.,

avoid MPI hangs when coupled runs terminate, add DOC/app/chapter_3.8;

Optional calculation of km and kh from initial TKE e_init;
Default initialization of km,kh = 0.00001 for ocean = .T.;
Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T.;
Bugfix: Rayleigh damping for ocean fixed.

Property svn:keywords set to Id

File size: 14.2 KB

Line
1	SUBROUTINE prandtl_fluxes
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: prandtl_fluxes.f90 108 2007-08-24 15:10:38Z letzel $
11	!
12	! 75 2007-03-22 09:54:05Z raasch
13	! moisture renamed humidity
14	!
15	! RCS Log replace by Id keyword, revision history cleaned up
16	!
17	! Revision 1.19 2006/04/26 12:24:35 raasch
18	! +OpenMP directives and optimization (array assignments replaced by DO loops)
19	!
20	! Revision 1.1 1998/01/23 10:06:06 raasch
21	! Initial revision
22	!
23	!
24	! Description:
25	! ------------
26	! Diagnostic computation of vertical fluxes in the Prandtl layer from the
27	! values of the variables at grid point k=1
28	!------------------------------------------------------------------------------!
29
30	USE arrays_3d
31	USE control_parameters
32	USE grid_variables
33	USE indices
34
35	IMPLICIT NONE
36
37	INTEGER :: i, j, k
38	REAL :: a, b, e_q, rifm, uv_total, z_p
39
40	!
41	!-- Compute theta*
42	IF ( constant_heatflux ) THEN
43	!
44	!-- For a given heat flux in the Prandtl layer:
45	!-- for u* use the value from the previous time step
46	!$OMP PARALLEL DO
47	DO i = nxl-1, nxr+1
48	DO j = nys-1, nyn+1
49	ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
50	!
51	!-- ts must be limited, because otherwise overflow may occur in case of
52	!-- us=0 when computing rif further below
53	IF ( ts(j,i) < -1.05E5 ) ts = -1.0E5
54	IF ( ts(j,i) > 1.0E5 ) ts = 1.0E5
55	ENDDO
56	ENDDO
57
58	ELSE
59	!
60	!-- For a given surface temperature:
61	!-- (the Richardson number is still the one from the previous time step)
62	!$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
63	DO i = nxl-1, nxr+1
64	DO j = nys-1, nyn+1
65
66	k = nzb_s_inner(j,i)
67	z_p = zu(k+1) - zw(k)
68
69	IF ( rif(j,i) >= 0.0 ) THEN
70	!
71	!-- Stable stratification
72	ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &
73	LOG( z_p / z0(j,i) ) + &
74	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
75	)
76	ELSE
77	!
78	!-- Unstable stratification
79	a = SQRT( 1.0 - 16.0 * rif(j,i) )
80	b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
81	!
82	!-- If a borderline case occurs, the formula for stable
83	!-- stratification must be used anyway, or else a zero division
84	!-- would occur in the argument of the logarithm
85	IF ( a == 1.0 .OR. b == 1.0 ) THEN
86	ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &
87	LOG( z_p / z0(j,i) ) + &
88	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
89	)
90	ELSE
91	ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &
92	LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
93	)
94	ENDIF
95	ENDIF
96
97	ENDDO
98	ENDDO
99	ENDIF
100
101	!
102	!-- Compute z_p/L (corresponds to the Richardson-flux number)
103	IF ( .NOT. humidity ) THEN
104	!$OMP PARALLEL DO PRIVATE( k, z_p )
105	DO i = nxl-1, nxr+1
106	DO j = nys-1, nyn+1
107	k = nzb_s_inner(j,i)
108	z_p = zu(k+1) - zw(k)
109	rif(j,i) = z_p * kappa * g * ts(j,i) / &
110	( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
111	!
112	!-- Limit the value range of the Richardson numbers.
113	!-- This is necessary for very small velocities (u,v --> 0), because
114	!-- the absolute value of rif can then become very large, which in
115	!-- consequence would result in very large shear stresses and very
116	!-- small momentum fluxes (both are generally unrealistic).
117	IF ( rif(j,i) < rif_min ) rif(j,i) = rif_min
118	IF ( rif(j,i) > rif_max ) rif(j,i) = rif_max
119	ENDDO
120	ENDDO
121	ELSE
122	!$OMP PARALLEL DO PRIVATE( k, z_p )
123	DO i = nxl-1, nxr+1
124	DO j = nys-1, nyn+1
125	k = nzb_s_inner(j,i)
126	z_p = zu(k+1) - zw(k)
127	rif(j,i) = z_p * kappa * g * &
128	( ts(j,i) + 0.61 * pt(k+1,j,i) * qs(j,i) ) / &
129	( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
130	!
131	!-- Limit the value range of the Richardson numbers.
132	!-- This is necessary for very small velocities (u,v --> 0), because
133	!-- the absolute value of rif can then become very large, which in
134	!-- consequence would result in very large shear stresses and very
135	!-- small momentum fluxes (both are generally unrealistic).
136	IF ( rif(j,i) < rif_min ) rif(j,i) = rif_min
137	IF ( rif(j,i) > rif_max ) rif(j,i) = rif_max
138	ENDDO
139	ENDDO
140	ENDIF
141
142	!
143	!-- Compute u* at the scalars' grid points
144	!$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
145	DO i = nxl, nxr
146	DO j = nys, nyn
147
148	k = nzb_s_inner(j,i)
149	z_p = zu(k+1) - zw(k)
150
151	!
152	!-- Compute the absolute value of the horizontal velocity
153	uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) ) )**2 + &
154	( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) ) )**2 &
155	)
156
157	IF ( rif(j,i) >= 0.0 ) THEN
158	!
159	!-- Stable stratification
160	us(j,i) = kappa * uv_total / ( &
161	LOG( z_p / z0(j,i) ) + &
162	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
163	)
164	ELSE
165	!
166	!-- Unstable stratification
167	a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
168	b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
169	!
170	!-- If a borderline case occurs, the formula for stable stratification
171	!-- must be used anyway, or else a zero division would occur in the
172	!-- argument of the logarithm.
173	IF ( a == 1.0 .OR. b == 1.0 ) THEN
174	us(j,i) = kappa * uv_total / ( &
175	LOG( z_p / z0(j,i) ) + &
176	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
177	)
178	ELSE
179	us(j,i) = kappa * uv_total / ( &
180	LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
181	2.0 * ( ATAN( b ) - ATAN( a ) ) &
182	)
183	ENDIF
184	ENDIF
185	ENDDO
186	ENDDO
187
188	!
189	!-- Compute u'w' for the total model domain.
190	!-- First compute the corresponding component of u* and square it.
191	!$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
192	DO i = nxl, nxr
193	DO j = nys, nyn
194
195	k = nzb_u_inner(j,i)
196	z_p = zu(k+1) - zw(k)
197
198	!
199	!-- Compute Richardson-flux number for this point
200	rifm = 0.5 * ( rif(j,i-1) + rif(j,i) )
201	IF ( rifm >= 0.0 ) THEN
202	!
203	!-- Stable stratification
204	usws(j,i) = kappa * u(k+1,j,i) / ( &
205	LOG( z_p / z0(j,i) ) + &
206	5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
207	)
208	ELSE
209	!
210	!-- Unstable stratification
211	a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
212	b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
213	!
214	!-- If a borderline case occurs, the formula for stable stratification
215	!-- must be used anyway, or else a zero division would occur in the
216	!-- argument of the logarithm.
217	IF ( a == 1.0 .OR. B == 1.0 ) THEN
218	usws(j,i) = kappa * u(k+1,j,i) / ( &
219	LOG( z_p / z0(j,i) ) + &
220	5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
221	)
222	ELSE
223	usws(j,i) = kappa * u(k+1,j,i) / ( &
224	LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
225	2.0 * ( ATAN( b ) - ATAN( a ) ) &
226	)
227	ENDIF
228	ENDIF
229	usws(j,i) = -usws(j,i) * ABS( usws(j,i) )
230	ENDDO
231	ENDDO
232
233	!
234	!-- Compute v'w' for the total model domain.
235	!-- First compute the corresponding component of u* and square it.
236	!$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
237	DO i = nxl, nxr
238	DO j = nys, nyn
239
240	k = nzb_v_inner(j,i)
241	z_p = zu(k+1) - zw(k)
242
243	!
244	!-- Compute Richardson-flux number for this point
245	rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )
246	IF ( rifm >= 0.0 ) THEN
247	!
248	!-- Stable stratification
249	vsws(j,i) = kappa * v(k+1,j,i) / ( &
250	LOG( z_p / z0(j,i) ) + &
251	5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
252	)
253	ELSE
254	!
255	!-- Unstable stratification
256	a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
257	b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
258	!
259	!-- If a borderline case occurs, the formula for stable stratification
260	!-- must be used anyway, or else a zero division would occur in the
261	!-- argument of the logarithm.
262	IF ( a == 1.0 .OR. b == 1.0 ) THEN
263	vsws(j,i) = kappa * v(k+1,j,i) / ( &
264	LOG( z_p / z0(j,i) ) + &
265	5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
266	)
267	ELSE
268	vsws(j,i) = kappa * v(k+1,j,i) / ( &
269	LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
270	2.0 * ( ATAN( b ) - ATAN( a ) ) &
271	)
272	ENDIF
273	ENDIF
274	vsws(j,i) = -vsws(j,i) * ABS( vsws(j,i) )
275	ENDDO
276	ENDDO
277
278	!
279	!-- If required compute q*
280	IF ( humidity .OR. passive_scalar ) THEN
281	IF ( constant_waterflux ) THEN
282	!
283	!-- For a given water flux in the Prandtl layer:
284	!$OMP PARALLEL DO
285	DO i = nxl-1, nxr+1
286	DO j = nys-1, nyn+1
287	qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
288	ENDDO
289	ENDDO
290
291	ELSE
292	!$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
293	DO i = nxl-1, nxr+1
294	DO j = nys-1, nyn+1
295
296	k = nzb_s_inner(j,i)
297	z_p = zu(k+1) - zw(k)
298
299	!
300	!-- assume saturation for atmosphere coupled to ocean
301	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
302	e_q = 6.1 * &
303	EXP( 0.07 * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15 ) )
304	q(k,j,i) = 0.622 * e_q / ( surface_pressure - e_q )
305	ENDIF
306	IF ( rif(j,i) >= 0.0 ) THEN
307	!
308	!-- Stable stratification
309	qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &
310	LOG( z_p / z0(j,i) ) + &
311	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
312	)
313	ELSE
314	!
315	!-- Unstable stratification
316	a = SQRT( 1.0 - 16.0 * rif(j,i) )
317	b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
318	!
319	!-- If a borderline case occurs, the formula for stable
320	!-- stratification must be used anyway, or else a zero division
321	!-- would occur in the argument of the logarithm.
322	IF ( a == 1.0 .OR. b == 1.0 ) THEN
323	qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &
324	LOG( z_p / z0(j,i) ) + &
325	5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
326	)
327	ELSE
328	qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &
329	LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
330	)
331	ENDIF
332	ENDIF
333
334	ENDDO
335	ENDDO
336	ENDIF
337	ENDIF
338
339	!
340	!-- Exchange the boundaries for u* and the momentum fluxes (fluxes only for
341	!-- completeness's sake).
342	CALL exchange_horiz_2d( us )
343	CALL exchange_horiz_2d( usws )
344	CALL exchange_horiz_2d( vsws )
345	IF ( humidity .OR. passive_scalar ) CALL exchange_horiz_2d( qsws )
346
347	!
348	!-- Compute the vertical kinematic heat flux
349	IF ( .NOT. constant_heatflux ) THEN
350	!$OMP PARALLEL DO
351	DO i = nxl-1, nxr+1
352	DO j = nys-1, nyn+1
353	shf(j,i) = -ts(j,i) * us(j,i)
354	ENDDO
355	ENDDO
356	ENDIF
357
358	!
359	!-- Compute the vertical water/scalar flux
360	IF ( .NOT. constant_heatflux .AND. ( humidity .OR. passive_scalar ) ) THEN
361	!$OMP PARALLEL DO
362	DO i = nxl-1, nxr+1
363	DO j = nys-1, nyn+1
364	qsws(j,i) = -qs(j,i) * us(j,i)
365	ENDDO
366	ENDDO
367	ENDIF
368
369	!
370	!-- Bottom boundary condition for the TKE
371	IF ( ibc_e_b == 2 ) THEN
372	!$OMP PARALLEL DO
373	DO i = nxl-1, nxr+1
374	DO j = nys-1, nyn+1
375	e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
376	!
377	!-- As a test: cm = 0.4
378	! e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2
379	e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
380	ENDDO
381	ENDDO
382	ENDIF
383
384
385	END SUBROUTINE prandtl_fluxes

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |