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

Last change on this file since 1 was 1, checked in by raasch, 17 years ago

Initial repository layout and content

File size: 15.7 KB
Line 
1 SUBROUTINE prandtl_fluxes
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: prandtl_fluxes.f90,v $
11! Revision 1.19  2006/04/26 12:24:35  raasch
12! +OpenMP directives and optimization (array assignments replaced by DO loops)
13!
14! Revision 1.18  2006/02/23 12:49:32  raasch
15! shf, ts, rif, us, usws, vsws, qs and qsws are now defined at the actual
16! surface of the model domain, i.e. either at the bottom or at the height
17! of the topography. Thus, zu(1) is replaced
18! by zp = ( zu(nzb_[variable]_inner(j,i)+1 - zw(nzb_[variable]_inner(j,i) ),
19! [variable](1,j,i) by [variable](nzb_[variable]_inner(j,i)+1,j,i) and
20! [variable](0,j,i) by [variable](nzb_[variable]_inner(j,i)  ,j,i).
21!
22! Revision 1.17  2004/04/30 12:43:50  raasch
23! rif_m replaced by rif (they are the same when used here)
24!
25! Revision 1.16  2003/11/20 15:13:26  raasch
26! Variable name (B) changed from capital to small
27!
28! Revision 1.15  2001/03/30 07:46:47  raasch
29! Translation of remaining German identifiers (variables, subroutines, etc.)
30!
31! Revision 1.14  2001/01/29 12:33:38  raasch
32! Passive scalar is considered
33!
34! Revision 1.13  2001/01/25 07:23:22  raasch
35! Range of ts is limited, since otherwise overflow may occur on REAL*4
36! machines
37!
38! Revision 1.11  2001/01/22 07:50:34  raasch
39! Module test_variables removed
40!
41! Revision 1.10  2000/04/13 13:44:41  schroeter
42! + implentation of the Prandtl layer for the total water content
43!
44! Revision 1.9  2000/01/26  14:48:48  14:48:48  letzel (Marcus Letzel)
45! correct wrong time level of rif in computation of theta*,
46! all comments translated into English
47!
48! Revision 1.8  1998/09/22 17:27:15  raasch
49! Testweise Randbedingung fuer TKE mit cm = 0.4, aber vorerst auskommentiert
50!
51! Revision 1.7  1998/08/05 06:54:23  raasch
52! Begrenzung von rif ist jetzt variabel (rif_max, rif_min)
53!
54! Revision 1.6  1998/07/06 12:30:01  raasch
55! + USE test_variables
56!
57! Revision 1.5  1998/04/15 11:22:58  raasch
58! Berechnung von theta* ueber Oberflaechentemperatur moeglich
59!
60! Revision 1.4  1998/03/11 11:53:22  raasch
61! Zusaetzliche untere Randbedingung fuer TKE ( (u*/0.1)**2 )
62!
63! Revision 1.3  1998/02/10 15:08:53  raasch
64! Grenzfall bei labiler Schichtung ( a=1, b=1 ) wird stabil gerechnet
65! Begrenzung von rif aktiviert
66!
67! Revision 1.2  1998/02/04 16:09:29  raasch
68! Berechnung der Waermefluesse
69!
70! Revision 1.1  1998/01/23 10:06:06  raasch
71! Initial revision
72!
73!
74! Description:
75! ------------
76! Diagnostic computation of vertical fluxes in the Prandtl layer from the
77! values of the variables at grid point k=1
78!------------------------------------------------------------------------------!
79
80    USE arrays_3d
81    USE control_parameters
82    USE grid_variables
83    USE indices
84
85    IMPLICIT NONE
86
87    INTEGER ::  i, j, k
88    REAL    ::  a, b, rifm, uv_total, z_p
89
90!
91!-- Compute theta*
92    IF ( constant_heatflux )  THEN
93!
94!--    For a given heat flux in the Prandtl layer:
95!--    for u* use the value from the previous time step
96       !$OMP PARALLEL DO
97       DO  i = nxl-1, nxr+1
98          DO  j = nys-1, nyn+1
99             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
100!
101!--          ts must be limited, because otherwise overflow may occur in case of
102!--          us=0 when computing rif further below
103             IF ( ts(j,i) < -1.05E5 )  ts = -1.0E5
104             IF ( ts(j,i) >   1.0E5 )  ts =  1.0E5
105          ENDDO
106       ENDDO
107
108    ELSE
109!
110!--    For a given surface temperature:
111!--    (the Richardson number is still the one from the previous time step)
112       !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
113       DO  i = nxl-1, nxr+1
114          DO  j = nys-1, nyn+1
115
116             k   = nzb_s_inner(j,i)
117             z_p = zu(k+1) - zw(k)
118
119             IF ( rif(j,i) >= 0.0 )  THEN
120!
121!--             Stable stratification
122                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
123                                  LOG( z_p / z0(j,i) ) +                   &
124                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
125                                                                )
126             ELSE
127!
128!--             Unstable stratification
129                a = SQRT( 1.0 - 16.0 * rif(j,i) )
130                b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
131!
132!--             If a borderline case occurs, the formula for stable
133!--             stratification must be used anyway, or else a zero division
134!--             would occur in the argument of the logarithm
135                IF ( a == 1.0  .OR.  b == 1.0 )  THEN
136                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
137                                     LOG( z_p / z0(j,i) ) +                   &
138                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
139                                                                   )
140                ELSE
141                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
142                                 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
143                                                                   )
144                ENDIF
145             ENDIF
146
147          ENDDO
148       ENDDO
149    ENDIF
150
151!
152!-- Compute z_p/L (corresponds to the Richardson-flux number)
153    IF ( .NOT. moisture ) THEN
154       !$OMP PARALLEL DO PRIVATE( k, z_p )
155       DO  i = nxl-1, nxr+1
156          DO  j = nys-1, nyn+1
157             k   = nzb_s_inner(j,i)
158             z_p = zu(k+1) - zw(k)
159             rif(j,i) = z_p * kappa * g * ts(j,i) / &
160                        ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
161!
162!--          Limit the value range of the Richardson numbers.
163!--          This is necessary for very small velocities (u,v --> 0), because
164!--          the absolute value of rif can then become very large, which in
165!--          consequence would result in very large shear stresses and very
166!--          small momentum fluxes (both are generally unrealistic).
167             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
168             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
169          ENDDO
170       ENDDO
171    ELSE
172       !$OMP PARALLEL DO PRIVATE( k, z_p )
173       DO  i = nxl-1, nxr+1
174          DO  j = nys-1, nyn+1
175             k   = nzb_s_inner(j,i)
176             z_p = zu(k+1) - zw(k)
177             rif(j,i) = z_p * kappa * g *                            &
178                        ( ts(j,i) + 0.61 * pt(k+1,j,i) * qs(j,i) ) / &
179                        ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
180!
181!--          Limit the value range of the Richardson numbers.
182!--          This is necessary for very small velocities (u,v --> 0), because
183!--          the absolute value of rif can then become very large, which in
184!--          consequence would result in very large shear stresses and very
185!--          small momentum fluxes (both are generally unrealistic).
186             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
187             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
188          ENDDO
189       ENDDO       
190    ENDIF
191
192!
193!-- Compute u* at the scalars' grid points
194    !$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
195    DO  i = nxl, nxr
196       DO  j = nys, nyn
197
198          k   = nzb_s_inner(j,i)
199          z_p = zu(k+1) - zw(k)
200
201!
202!--       Compute the absolute value of the horizontal velocity
203          uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) ) )**2 + &
204                           ( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) ) )**2   &
205                         )
206
207          IF ( rif(j,i) >= 0.0 )  THEN
208!
209!--          Stable stratification
210             us(j,i) = kappa * uv_total / (                                &
211                                  LOG( z_p / z0(j,i) ) +                   &
212                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
213                                          )
214          ELSE
215!
216!--          Unstable stratification
217             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
218             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
219!
220!--          If a borderline case occurs, the formula for stable stratification
221!--          must be used anyway, or else a zero division would occur in the
222!--          argument of the logarithm.
223             IF ( a == 1.0  .OR.  b == 1.0 )  THEN
224                us(j,i) = kappa * uv_total / (                                &
225                                     LOG( z_p / z0(j,i) ) +                   &
226                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
227                                             )
228             ELSE
229                us(j,i) = kappa * uv_total / (                               &
230                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
231                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
232                                             )
233             ENDIF
234          ENDIF
235       ENDDO
236    ENDDO
237
238!
239!-- Compute u'w' for the total model domain.
240!-- First compute the corresponding component of u* and square it.
241    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
242    DO  i = nxl, nxr
243       DO  j = nys, nyn
244
245          k   = nzb_u_inner(j,i)
246          z_p = zu(k+1) - zw(k)
247
248!
249!--       Compute Richardson-flux number for this point
250          rifm = 0.5 * ( rif(j,i-1) + rif(j,i) )
251          IF ( rifm >= 0.0 )  THEN
252!
253!--          Stable stratification
254             usws(j,i) = kappa * u(k+1,j,i) / (                           &
255                                     LOG( z_p / z0(j,i) ) +               &
256                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
257                                              )
258          ELSE
259!
260!--          Unstable stratification
261             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
262             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
263!
264!--          If a borderline case occurs, the formula for stable stratification
265!--          must be used anyway, or else a zero division would occur in the
266!--          argument of the logarithm.
267             IF ( a == 1.0  .OR.  B == 1.0 )  THEN
268                usws(j,i) = kappa * u(k+1,j,i) / (                           &
269                                        LOG( z_p / z0(j,i) ) +               &
270                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
271                                                 )
272             ELSE
273                usws(j,i) = kappa * u(k+1,j,i) / (                           &
274                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
275                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
276                                                 )
277             ENDIF
278          ENDIF
279          usws(j,i) = -usws(j,i) * ABS( usws(j,i) )
280       ENDDO
281    ENDDO
282
283!
284!-- Compute v'w' for the total model domain.
285!-- First compute the corresponding component of u* and square it.
286    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
287    DO  i = nxl, nxr
288       DO  j = nys, nyn
289
290          k   = nzb_v_inner(j,i)
291          z_p = zu(k+1) - zw(k)
292
293!
294!--       Compute Richardson-flux number for this point
295          rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )
296          IF ( rifm >= 0.0 )  THEN
297!
298!--          Stable stratification
299             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
300                                     LOG( z_p / z0(j,i) ) +               &
301                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
302                                              )
303          ELSE
304!
305!--          Unstable stratification
306             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
307             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
308!
309!--          If a borderline case occurs, the formula for stable stratification
310!--          must be used anyway, or else a zero division would occur in the
311!--          argument of the logarithm.
312             IF ( a == 1.0  .OR.  b == 1.0 )  THEN
313                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
314                                        LOG( z_p / z0(j,i) ) +               &
315                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
316                                                 )
317             ELSE
318                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
319                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
320                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
321                                                 )
322             ENDIF
323          ENDIF
324          vsws(j,i) = -vsws(j,i) * ABS( vsws(j,i) )
325       ENDDO
326    ENDDO
327
328!
329!-- If required compute q*
330    IF ( moisture  .OR.  passive_scalar )  THEN
331       IF ( constant_waterflux )  THEN
332!
333!--       For a given water flux in the Prandtl layer:
334          !$OMP PARALLEL DO
335          DO  i = nxl-1, nxr+1
336             DO  j = nys-1, nyn+1
337                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
338             ENDDO
339          ENDDO
340         
341       ELSE         
342          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
343          DO  i = nxl-1, nxr+1
344             DO  j = nys-1, nyn+1
345
346                k   = nzb_s_inner(j,i)
347                z_p = zu(k+1) - zw(k)
348
349                IF ( rif(j,i) >= 0.0 )  THEN
350!
351!--                Stable stratification
352                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
353                                  LOG( z_p / z0(j,i) ) +                   &
354                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
355                                                                 )
356                ELSE
357!
358!--                Unstable stratification
359                   a = SQRT( 1.0 - 16.0 * rif(j,i) )
360                   b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
361!
362!--                If a borderline case occurs, the formula for stable
363!--                stratification must be used anyway, or else a zero division
364!--                would occur in the argument of the logarithm.
365                   IF ( a == 1.0  .OR.  b == 1.0 )  THEN
366                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
367                                     LOG( z_p / z0(j,i) ) +                   &
368                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
369                                                                    )
370                   ELSE
371                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
372                                  LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
373                                                                    )
374                   ENDIF
375                ENDIF
376
377             ENDDO
378          ENDDO
379       ENDIF
380    ENDIF
381
382!
383!-- Exchange the boundaries for u* and the momentum fluxes (fluxes only for
384!-- completeness's sake).
385    CALL exchange_horiz_2d( us )
386    CALL exchange_horiz_2d( usws )
387    CALL exchange_horiz_2d( vsws )
388    IF ( moisture  .OR.  passive_scalar )  CALL exchange_horiz_2d( qsws )
389
390!
391!-- Compute the vertical kinematic heat flux
392    IF ( .NOT. constant_heatflux )  THEN
393       !$OMP PARALLEL DO
394       DO  i = nxl-1, nxr+1
395          DO  j = nys-1, nyn+1
396             shf(j,i) = -ts(j,i) * us(j,i)
397          ENDDO
398       ENDDO
399    ENDIF
400
401!
402!-- Compute the vertical water/scalar flux
403    IF ( .NOT. constant_heatflux .AND. ( moisture .OR. passive_scalar ) ) THEN
404       !$OMP PARALLEL DO
405       DO  i = nxl-1, nxr+1
406          DO  j = nys-1, nyn+1
407             qsws(j,i) = -qs(j,i) * us(j,i)
408          ENDDO
409       ENDDO
410    ENDIF
411
412!
413!-- Bottom boundary condition for the TKE
414    IF ( ibc_e_b == 2 )  THEN
415       !$OMP PARALLEL DO
416       DO  i = nxl-1, nxr+1
417          DO  j = nys-1, nyn+1
418             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
419!
420!--          As a test: cm = 0.4
421!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2
422             e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
423          ENDDO
424       ENDDO
425    ENDIF
426
427
428 END SUBROUTINE prandtl_fluxes
Note: See TracBrowser for help on using the repository browser.