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

Last change on this file since 315 was 315, checked in by raasch, 12 years ago

bugfix: qsws was calculated in case of constant_heatflux = .FALSE.

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