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

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

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

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