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

Last change on this file since 1353 was 1341, checked in by kanani, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 15.7 KB
Line 
1 SUBROUTINE prandtl_fluxes
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: prandtl_fluxes.f90 1341 2014-03-25 19:48:09Z heinze $
27!
28! 1340 2014-03-25 19:45:13Z kanani
29! REAL constants defined as wp-kind
30!
31! 1320 2014-03-20 08:40:49Z raasch
32! ONLY-attribute added to USE-statements,
33! kind-parameters added to all INTEGER and REAL declaration statements,
34! kinds are defined in new module kinds,
35! old module precision_kind is removed,
36! revision history before 2012 removed,
37! comment fields (!:) to be used for variable explanations added to
38! all variable declaration statements
39!
40! 1276 2014-01-15 13:40:41Z heinze
41! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
42!
43! 1257 2013-11-08 15:18:40Z raasch
44! openACC "kernels do" replaced by "kernels loop", "loop independent" added
45!
46! 1036 2012-10-22 13:43:42Z raasch
47! code put under GPL (PALM 3.9)
48!
49! 1015 2012-09-27 09:23:24Z raasch
50! OpenACC statements added
51!
52! 978 2012-08-09 08:28:32Z fricke
53! roughness length for scalar quantities z0h added
54!
55! Revision 1.1  1998/01/23 10:06:06  raasch
56! Initial revision
57!
58!
59! Description:
60! ------------
61! Diagnostic computation of vertical fluxes in the Prandtl layer from the
62! values of the variables at grid point k=1
63!------------------------------------------------------------------------------!
64
65    USE arrays_3d,                                                             &
66        ONLY:  e, pt, q, qs, qsws, rif, shf, ts, u, us, usws, v, vpt, vsws,    &
67               zu, zw, z0, z0h
68
69    USE control_parameters,                                                    &
70        ONLY:  constant_heatflux, constant_waterflux, coupling_mode, g,        &
71               humidity, ibc_e_b, kappa, large_scale_forcing, lsf_surf,        &
72               passive_scalar, pt_surface, q_surface, rif_max, rif_min,        &
73               run_coupled, surface_pressure
74
75    USE indices,                                                               &
76        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner,        &
77               nzb_u_inner, nzb_v_inner
78
79    USE kinds
80
81    IMPLICIT NONE
82
83    INTEGER(iwp) ::  i            !:
84    INTEGER(iwp) ::  j            !:
85    INTEGER(iwp) ::  k            !:
86
87    LOGICAL      ::  coupled_run  !:
88
89    REAL(wp)     ::  a            !:
90    REAL(wp)     ::  b            !:
91    REAL(wp)     ::  e_q          !:
92    REAL(wp)     ::  rifm         !:
93    REAL(wp)     ::  uv_total     !:
94    REAL(wp)     ::  z_p          !:
95
96!
97!-- Data information for accelerators
98    !$acc data present( e, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt, q, qs ) &
99    !$acc      present( qsws, rif, shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h )
100!
101!-- Compute theta*
102    IF ( constant_heatflux )  THEN
103!
104!--    For a given heat flux in the Prandtl layer:
105!--    for u* use the value from the previous time step
106       !$OMP PARALLEL DO
107       !$acc kernels loop
108       DO  i = nxlg, nxrg
109          DO  j = nysg, nyng
110             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
111!
112!--          ts must be limited, because otherwise overflow may occur in case of
113!--          us=0 when computing rif further below
114             IF ( ts(j,i) < -1.05E5_wp )  ts(j,i) = -1.0E5_wp
115             IF ( ts(j,i) >   1.0E5_wp )  ts(j,i) =  1.0E5_wp
116          ENDDO
117       ENDDO
118
119    ELSE
120!
121!--    For a given surface temperature:
122!--    (the Richardson number is still the one from the previous time step)
123   
124       IF ( large_scale_forcing .AND. lsf_surf )  THEN
125          pt(0,:,:) = pt_surface
126       ENDIF
127
128       !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
129       !$acc kernels loop
130       DO  i = nxlg, nxrg
131          DO  j = nysg, nyng
132
133             k   = nzb_s_inner(j,i)
134             z_p = zu(k+1) - zw(k)
135
136             IF ( rif(j,i) >= 0.0_wp )  THEN
137!
138!--             Stable stratification
139                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (           &
140                                  LOG( z_p / z0h(j,i) ) +                   &
141                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
142                                                                )
143             ELSE
144!
145!--             Unstable stratification
146                a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) )
147                b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p )
148
149                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (          &
150                          LOG( z_p / z0h(j,i) ) -                           &
151                          2.0_wp * LOG( ( 1.0_wp + a ) / ( 1.0_wp + b ) ) )
152             ENDIF
153
154          ENDDO
155       ENDDO
156    ENDIF
157
158!
159!-- Compute z_p/L (corresponds to the Richardson-flux number)
160    IF ( .NOT. humidity )  THEN
161       !$OMP PARALLEL DO PRIVATE( k, z_p )
162       !$acc kernels loop
163       DO  i = nxlg, nxrg
164          DO  j = nysg, nyng
165             k   = nzb_s_inner(j,i)
166             z_p = zu(k+1) - zw(k)
167             rif(j,i) = z_p * kappa * g * ts(j,i) / &
168                        ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
169!
170!--          Limit the value range of the Richardson numbers.
171!--          This is necessary for very small velocities (u,v --> 0), because
172!--          the absolute value of rif can then become very large, which in
173!--          consequence would result in very large shear stresses and very
174!--          small momentum fluxes (both are generally unrealistic).
175             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
176             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
177          ENDDO
178       ENDDO
179    ELSE
180       !$OMP PARALLEL DO PRIVATE( k, z_p )
181       !$acc kernels loop
182       DO  i = nxlg, nxrg
183          DO  j = nysg, nyng
184             k   = nzb_s_inner(j,i)
185             z_p = zu(k+1) - zw(k)
186             rif(j,i) = z_p * kappa * g *                            &
187                        ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) ) / &
188                        ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
189!
190!--          Limit the value range of the Richardson numbers.
191!--          This is necessary for very small velocities (u,v --> 0), because
192!--          the absolute value of rif can then become very large, which in
193!--          consequence would result in very large shear stresses and very
194!--          small momentum fluxes (both are generally unrealistic).
195             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
196             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
197          ENDDO
198       ENDDO       
199    ENDIF
200
201!
202!-- Compute u* at the scalars' grid points
203    !$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
204    !$acc kernels loop
205    DO  i = nxl, nxr
206       DO  j = nys, nyn
207
208          k   = nzb_s_inner(j,i)
209          z_p = zu(k+1) - zw(k)
210
211!
212!--       Compute the absolute value of the horizontal velocity
213!--       (relative to the surface)
214          uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)        &
215                                      - u(k,j,i)   - u(k,j,i+1) ) )**2 + &
216                           ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)        &
217                                      - v(k,j,i)   - v(k,j+1,i) ) )**2 )   
218
219
220          IF ( rif(j,i) >= 0.0_wp )  THEN
221!
222!--          Stable stratification
223             us(j,i) = kappa * uv_total / (                                &
224                                  LOG( z_p / z0(j,i) ) +                   &
225                                  5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
226                                          )
227          ELSE
228!
229!--          Unstable stratification
230             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) )
231             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) / z_p * z0(j,i) ) )
232
233             us(j,i) = kappa * uv_total / (                                  &
234                       LOG( z_p / z0(j,i) ) -                                &
235                       LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (              &
236                            ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +          &
237                               2.0_wp * ( ATAN( a ) - ATAN( b ) )                  &
238                                           )
239          ENDIF
240       ENDDO
241    ENDDO
242
243!
244!-- Values of us at ghost point locations are needed for the evaluation of usws
245!-- and vsws.
246    !$acc update host( us )
247    CALL exchange_horiz_2d( us )
248    !$acc update device( us )
249
250!
251!-- Compute u'w' for the total model domain.
252!-- First compute the corresponding component of u* and square it.
253    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
254    !$acc kernels loop
255    DO  i = nxl, nxr
256       DO  j = nys, nyn
257
258          k   = nzb_u_inner(j,i)
259          z_p = zu(k+1) - zw(k)
260
261!
262!--       Compute Richardson-flux number for this point
263          rifm = 0.5_wp * ( rif(j,i-1) + rif(j,i) )
264          IF ( rifm >= 0.0_wp )  THEN
265!
266!--          Stable stratification
267             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ (              &
268                                     LOG( z_p / z0(j,i) ) +               &
269                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
270                                                            )
271          ELSE
272!
273!--          Unstable stratification
274             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
275             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
276
277             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / (            &
278                         LOG( z_p / z0(j,i) ) -                           &
279                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
280                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
281                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )             &
282                                                 )
283          ENDIF
284          usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
285       ENDDO
286    ENDDO
287
288!
289!-- Compute v'w' for the total model domain.
290!-- First compute the corresponding component of u* and square it.
291    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
292    !$acc kernels loop
293    DO  i = nxl, nxr
294       DO  j = nys, nyn
295
296          k   = nzb_v_inner(j,i)
297          z_p = zu(k+1) - zw(k)
298
299!
300!--       Compute Richardson-flux number for this point
301          rifm = 0.5_wp * ( rif(j-1,i) + rif(j,i) )
302          IF ( rifm >= 0.0_wp )  THEN
303!
304!--          Stable stratification
305             vsws(j,i) = kappa * ( v(k+1,j,i) -  v(k,j,i) ) / (           &
306                                     LOG( z_p / z0(j,i) ) +               &
307                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
308                                                              )
309          ELSE
310!
311!--          Unstable stratification
312             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
313             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
314
315             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / (            &
316                         LOG( z_p / z0(j,i) ) -                           &
317                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
318                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
319                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )             &
320                                                 )
321          ENDIF
322          vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j-1,i) + us(j,i) )
323       ENDDO
324    ENDDO
325
326!
327!-- If required compute q*
328    IF ( humidity  .OR.  passive_scalar )  THEN
329       IF ( constant_waterflux )  THEN
330!
331!--       For a given water flux in the Prandtl layer:
332          !$OMP PARALLEL DO
333          !$acc kernels loop
334          DO  i = nxlg, nxrg
335             DO  j = nysg, nyng
336                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
337             ENDDO
338          ENDDO
339         
340       ELSE
341          coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled )
342
343           IF ( large_scale_forcing .AND. lsf_surf )  THEN
344              q(0,:,:) = q_surface
345           ENDIF
346
347          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
348          !$acc kernels loop independent
349          DO  i = nxlg, nxrg
350             !$acc loop independent
351             DO  j = nysg, nyng
352
353                k   = nzb_s_inner(j,i)
354                z_p = zu(k+1) - zw(k)
355
356!
357!--             Assume saturation for atmosphere coupled to ocean (but not
358!--             in case of precursor runs)
359                IF ( coupled_run )  THEN
360                   e_q = 6.1_wp * &
361                           EXP( 0.07_wp * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15_wp ) )
362                   q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q )
363                ENDIF
364                IF ( rif(j,i) >= 0.0_wp )  THEN
365!
366!--                Stable stratification
367                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
368                                  LOG( z_p / z0h(j,i) ) +                   &
369                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
370                                                                 )
371                ELSE
372!
373!--                Unstable stratification
374                   a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 
375                   b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 
376 
377                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (        &
378                             LOG( z_p / z0h(j,i) ) -                        &
379                              2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) )
380                ENDIF
381
382             ENDDO
383          ENDDO
384       ENDIF
385    ENDIF
386
387!
388!-- Exchange the boundaries for the momentum fluxes (only for sake of
389!-- completeness)
390    !$acc update host( usws, vsws )
391    CALL exchange_horiz_2d( usws )
392    CALL exchange_horiz_2d( vsws )
393    !$acc update device( usws, vsws )
394    IF ( humidity  .OR.  passive_scalar )  THEN
395       !$acc update host( qsws )
396       CALL exchange_horiz_2d( qsws )
397       !$acc update device( qsws )
398    ENDIF
399
400!
401!-- Compute the vertical kinematic heat flux
402    IF ( .NOT. constant_heatflux )  THEN
403       !$OMP PARALLEL DO
404       !$acc kernels loop independent
405       DO  i = nxlg, nxrg
406          !$acc loop independent
407          DO  j = nysg, nyng
408             shf(j,i) = -ts(j,i) * us(j,i)
409          ENDDO
410       ENDDO
411    ENDIF
412
413!
414!-- Compute the vertical water/scalar flux
415    IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar ) ) THEN
416       !$OMP PARALLEL DO
417       !$acc kernels loop independent
418       DO  i = nxlg, nxrg
419          !$acc loop independent
420          DO  j = nysg, nyng
421             qsws(j,i) = -qs(j,i) * us(j,i)
422          ENDDO
423       ENDDO
424    ENDIF
425
426!
427!-- Bottom boundary condition for the TKE
428    IF ( ibc_e_b == 2 )  THEN
429       !$OMP PARALLEL DO
430       !$acc kernels loop independent
431       DO  i = nxlg, nxrg
432          !$acc loop independent
433          DO  j = nysg, nyng
434             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2
435!
436!--          As a test: cm = 0.4
437!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2
438             e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
439          ENDDO
440       ENDDO
441    ENDIF
442
443    !$acc end data
444
445 END SUBROUTINE prandtl_fluxes
Note: See TracBrowser for help on using the repository browser.