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

Last change on this file since 313 was 291, checked in by raasch, 15 years ago

changes for coupling with independent precursor runs; z_i calculation with Sullivan criterion

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