source: palm/tags/release-3.6/SOURCE/prandtl_fluxes.f90 @ 4003

Last change on this file since 4003 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1 SUBROUTINE prandtl_fluxes
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: prandtl_fluxes.f90 198 2008-09-17 08:55:28Z knoop $
11!
12! 187 2008-08-06 16:25:09Z letzel
13! Bugfix: modification of the calculation of the vertical turbulent momentum
14! fluxes u'w' and v'w'
15! Bugfix: change definition of us_wall from 1D to 2D
16! Change: modification of the integrated version of the profile function for
17! momentum for unstable stratification (does not effect results)
18!
19! 108 2007-08-24 15:10:38Z letzel
20! assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean
21!
22! 75 2007-03-22 09:54:05Z raasch
23! moisture renamed humidity
24!
25! RCS Log replace by Id keyword, revision history cleaned up
26!
27! Revision 1.19  2006/04/26 12:24:35  raasch
28! +OpenMP directives and optimization (array assignments replaced by DO loops)
29!
30! Revision 1.1  1998/01/23 10:06:06  raasch
31! Initial revision
32!
33!
34! Description:
35! ------------
36! Diagnostic computation of vertical fluxes in the Prandtl layer from the
37! values of the variables at grid point k=1
38!------------------------------------------------------------------------------!
39
40    USE arrays_3d
41    USE control_parameters
42    USE grid_variables
43    USE indices
44
45    IMPLICIT NONE
46
47    INTEGER ::  i, j, k
48    REAL    ::  a, b, e_q, rifm, uv_total, z_p
49
50!
51!-- Compute theta*
52    IF ( constant_heatflux )  THEN
53!
54!--    For a given heat flux in the Prandtl layer:
55!--    for u* use the value from the previous time step
56       !$OMP PARALLEL DO
57       DO  i = nxl-1, nxr+1
58          DO  j = nys-1, nyn+1
59             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
60!
61!--          ts must be limited, because otherwise overflow may occur in case of
62!--          us=0 when computing rif further below
63             IF ( ts(j,i) < -1.05E5 )  ts = -1.0E5
64             IF ( ts(j,i) >   1.0E5 )  ts =  1.0E5
65          ENDDO
66       ENDDO
67
68    ELSE
69!
70!--    For a given surface temperature:
71!--    (the Richardson number is still the one from the previous time step)
72       !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
73       DO  i = nxl-1, nxr+1
74          DO  j = nys-1, nyn+1
75
76             k   = nzb_s_inner(j,i)
77             z_p = zu(k+1) - zw(k)
78
79             IF ( rif(j,i) >= 0.0 )  THEN
80!
81!--             Stable stratification
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!
88!--             Unstable stratification
89                a = SQRT( 1.0 - 16.0 * rif(j,i) )
90                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
91
92                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (         &
93                          LOG( z_p / z0(j,i) ) -                           &
94                          2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
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 = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
168             b = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
169
170             us(j,i) = kappa * uv_total / (                                  &
171                       LOG( z_p / z0(j,i) ) -                                &
172                       LOG( ( 1.0 + a )**2 * ( 1.0 + a**2 ) / (              &
173                            ( 1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +          &
174                            2.0 * ( ATAN( a ) - ATAN( b ) )                  &
175                                           )
176          ENDIF
177       ENDDO
178    ENDDO
179
180!
181!-- Values of us at ghost point locations are needed for the evaluation of usws
182!-- and vsws.
183    CALL exchange_horiz_2d( us )
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 = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
208             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
209
210             usws(j,i) = kappa * u(k+1,j,i) / (                           &
211                         LOG( z_p / z0(j,i) ) -                           &
212                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
213                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
214                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
215                                                 )
216          ENDIF
217          usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) )
218       ENDDO
219    ENDDO
220
221!
222!-- Compute v'w' for the total model domain.
223!-- First compute the corresponding component of u* and square it.
224    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
225    DO  i = nxl, nxr
226       DO  j = nys, nyn
227
228          k   = nzb_v_inner(j,i)
229          z_p = zu(k+1) - zw(k)
230
231!
232!--       Compute Richardson-flux number for this point
233          rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )
234          IF ( rifm >= 0.0 )  THEN
235!
236!--          Stable stratification
237             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
238                                     LOG( z_p / z0(j,i) ) +               &
239                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
240                                              )
241          ELSE
242!
243!--          Unstable stratification
244             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
245             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
246
247             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
248                         LOG( z_p / z0(j,i) ) -                           &
249                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
250                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
251                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
252                                                 )
253          ENDIF
254          vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) )
255       ENDDO
256    ENDDO
257
258!
259!-- If required compute q*
260    IF ( humidity  .OR.  passive_scalar )  THEN
261       IF ( constant_waterflux )  THEN
262!
263!--       For a given water flux in the Prandtl layer:
264          !$OMP PARALLEL DO
265          DO  i = nxl-1, nxr+1
266             DO  j = nys-1, nyn+1
267                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
268             ENDDO
269          ENDDO
270         
271       ELSE         
272          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
273          DO  i = nxl-1, nxr+1
274             DO  j = nys-1, nyn+1
275
276                k   = nzb_s_inner(j,i)
277                z_p = zu(k+1) - zw(k)
278
279!
280!--             assume saturation for atmosphere coupled to ocean
281                IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
282                   e_q = 6.1 * &
283                        EXP( 0.07 * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15 ) )
284                   q(k,j,i) = 0.622 * e_q / ( surface_pressure - e_q )
285                ENDIF
286                IF ( rif(j,i) >= 0.0 )  THEN
287!
288!--                Stable stratification
289                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
290                                  LOG( z_p / z0(j,i) ) +                   &
291                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
292                                                                 )
293                ELSE
294!
295!--                Unstable stratification
296                   a = SQRT( 1.0 - 16.0 * rif(j,i) ) 
297                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p ) 
298 
299                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (   &
300                             LOG( z_p / z0(j,i) ) -                    &
301                              2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
302                ENDIF
303
304             ENDDO
305          ENDDO
306       ENDIF
307    ENDIF
308
309!
310!-- Exchange the boundaries for the momentum fluxes (only for sake of
311!-- completeness)
312    CALL exchange_horiz_2d( usws )
313    CALL exchange_horiz_2d( vsws )
314    IF ( humidity  .OR.  passive_scalar )  CALL exchange_horiz_2d( qsws )
315
316!
317!-- Compute the vertical kinematic heat flux
318    IF ( .NOT. constant_heatflux )  THEN
319       !$OMP PARALLEL DO
320       DO  i = nxl-1, nxr+1
321          DO  j = nys-1, nyn+1
322             shf(j,i) = -ts(j,i) * us(j,i)
323          ENDDO
324       ENDDO
325    ENDIF
326
327!
328!-- Compute the vertical water/scalar flux
329    IF ( .NOT. constant_heatflux .AND. ( humidity .OR. passive_scalar ) ) THEN
330       !$OMP PARALLEL DO
331       DO  i = nxl-1, nxr+1
332          DO  j = nys-1, nyn+1
333             qsws(j,i) = -qs(j,i) * us(j,i)
334          ENDDO
335       ENDDO
336    ENDIF
337
338!
339!-- Bottom boundary condition for the TKE
340    IF ( ibc_e_b == 2 )  THEN
341       !$OMP PARALLEL DO
342       DO  i = nxl-1, nxr+1
343          DO  j = nys-1, nyn+1
344             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
345!
346!--          As a test: cm = 0.4
347!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2
348             e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
349          ENDDO
350       ENDDO
351    ENDIF
352
353
354 END SUBROUTINE prandtl_fluxes
Note: See TracBrowser for help on using the repository browser.