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

Last change on this file since 75 was 75, checked in by raasch, 14 years ago

preliminary update for changes concerning non-cyclic boundary conditions

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