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

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

typo in file headers removed

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