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

Last change on this file since 1623 was 1552, checked in by maronga, 10 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 19.5 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 1552 2015-03-03 14:27:15Z heinze $
27!
28! 1551 2015-03-03 14:18:16Z maronga
29! Removed land surface model part. The surface fluxes are now always calculated
30! within prandtl_fluxes, based on the given surface temperature/humidity (which
31! is either provided by the land surface model, by large scale forcing data, or
32! directly prescribed by the user.
33!
34! 1496 2014-12-02 17:25:50Z maronga
35! Adapted for land surface model
36!
37! 1494 2014-11-21 17:14:03Z maronga
38! Bugfixes: qs is now calculated before calculation of Rif. Ccalculation of
39! buoyancy flux in Rif corrected (added missing humidity term), allow use of
40! topography for coupled runs (not tested)
41!
42! 1361 2014-04-16 15:17:48Z hoffmann
43! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
44! drop concentration (nrsws) added
45!
46! 1340 2014-03-25 19:45:13Z kanani
47! REAL constants defined as wp-kind
48!
49! 1320 2014-03-20 08:40:49Z raasch
50! ONLY-attribute added to USE-statements,
51! kind-parameters added to all INTEGER and REAL declaration statements,
52! kinds are defined in new module kinds,
53! old module precision_kind is removed,
54! revision history before 2012 removed,
55! comment fields (!:) to be used for variable explanations added to
56! all variable declaration statements
57!
58! 1276 2014-01-15 13:40:41Z heinze
59! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
60!
61! 1257 2013-11-08 15:18:40Z raasch
62! openACC "kernels do" replaced by "kernels loop", "loop independent" added
63!
64! 1036 2012-10-22 13:43:42Z raasch
65! code put under GPL (PALM 3.9)
66!
67! 1015 2012-09-27 09:23:24Z raasch
68! OpenACC statements added
69!
70! 978 2012-08-09 08:28:32Z fricke
71! roughness length for scalar quantities z0h added
72!
73! Revision 1.1  1998/01/23 10:06:06  raasch
74! Initial revision
75!
76!
77! Description:
78! ------------
79! Diagnostic computation of vertical fluxes in the Prandtl layer from the
80! values of the variables at grid point k=1
81!------------------------------------------------------------------------------!
82
83    USE arrays_3d,                                                             &
84        ONLY:  e, nr, nrs, nrsws, pt, q, qr, qrs, qrsws, qs, qsws, rif, shf,   &
85               ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
86
87    USE control_parameters,                                                    &
88        ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
89               coupling_mode, g, humidity, ibc_e_b, icloud_scheme, kappa,      &
90               large_scale_forcing, lsf_surf, passive_scalar, precipitation,   &
91               pt_surface, q_surface, rif_max, rif_min, run_coupled,           &
92               surface_pressure
93
94    USE indices,                                                               &
95        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner,        &
96               nzb_u_inner, nzb_v_inner
97
98    USE kinds
99
100    IMPLICIT NONE
101
102    INTEGER(iwp) ::  i            !:
103    INTEGER(iwp) ::  j            !:
104    INTEGER(iwp) ::  k            !:
105
106    LOGICAL      ::  coupled_run  !:
107
108    REAL(wp)     ::  a            !:
109    REAL(wp)     ::  b            !:
110    REAL(wp)     ::  e_q          !:
111    REAL(wp)     ::  rifm         !:
112    REAL(wp)     ::  uv_total     !:
113    REAL(wp)     ::  z_p          !:
114
115!
116!-- Data information for accelerators
117    !$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt )  &
118    !$acc      present( q, qs, qsws, qrsws, rif, shf, ts, u, us, usws, v )     &
119    !$acc      present( vpt, vsws, zu, zw, z0, z0h )
120!
121!-- Compute theta*
122    IF ( constant_heatflux )  THEN
123
124!
125!--    For a given heat flux in the Prandtl layer:
126!--    for u* use the value from the previous time step
127       !$OMP PARALLEL DO
128       !$acc kernels loop
129       DO  i = nxlg, nxrg
130          DO  j = nysg, nyng
131             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
132!
133!--          ts must be limited, because otherwise overflow may occur in case of
134!--          us=0 when computing rif further below
135             IF ( ts(j,i) < -1.05E5_wp )  ts(j,i) = -1.0E5_wp
136             IF ( ts(j,i) >   1.0E5_wp )  ts(j,i) =  1.0E5_wp
137          ENDDO
138       ENDDO
139
140    ELSE
141!
142!--    For a given surface temperature:
143!--    (the Richardson number is still the one from the previous time step)
144       IF ( large_scale_forcing .AND. lsf_surf )  THEN
145          !$OMP PARALLEL DO
146          !$acc kernels loop
147          DO  i = nxlg, nxrg
148             DO  j = nysg, nyng
149                k = nzb_s_inner(j,i)
150                pt(k,j,i) = pt_surface
151             ENDDO
152          ENDDO
153       ENDIF
154
155       !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
156       !$acc kernels loop
157       DO  i = nxlg, nxrg
158          DO  j = nysg, nyng
159
160             k   = nzb_s_inner(j,i)
161             z_p = zu(k+1) - zw(k)
162
163             IF ( rif(j,i) >= 0.0_wp )  THEN
164!
165!--             Stable stratification
166                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (           &
167                                  LOG( z_p / z0h(j,i) ) +                   &
168                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
169                                                                )
170             ELSE
171!
172!--             Unstable stratification
173                a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) )
174                b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p )
175
176                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (             &
177                          LOG( z_p / z0h(j,i) ) -                              &
178                          2.0_wp * LOG( ( 1.0_wp + a ) / ( 1.0_wp + b ) ) )
179             ENDIF
180
181          ENDDO
182       ENDDO
183    ENDIF
184
185!
186!-- If required compute q*
187    IF ( humidity  .OR.  passive_scalar )  THEN
188       IF ( constant_waterflux )  THEN
189!
190!--       For a given water flux in the Prandtl layer:
191          !$OMP PARALLEL DO
192          !$acc kernels loop
193          DO  i = nxlg, nxrg
194             DO  j = nysg, nyng
195                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
196             ENDDO
197          ENDDO
198         
199       ELSE
200          coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled )
201
202           IF ( large_scale_forcing .AND. lsf_surf )  THEN
203                !$OMP PARALLEL DO
204                !$acc kernels loop
205                 DO  i = nxlg, nxrg
206                    DO  j = nysg, nyng
207                       k = nzb_s_inner(j,i)
208                       q(k,j,i) = q_surface
209                    ENDDO
210                 ENDDO
211           ENDIF
212
213          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
214          !$acc kernels loop independent
215          DO  i = nxlg, nxrg
216             !$acc loop independent
217             DO  j = nysg, nyng
218
219                k   = nzb_s_inner(j,i)
220                z_p = zu(k+1) - zw(k)
221
222!
223!--             Assume saturation for atmosphere coupled to ocean (but not
224!--             in case of precursor runs)
225                IF ( coupled_run )  THEN
226                   e_q = 6.1_wp * &
227                           EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i)) - 273.15_wp ) )
228                   q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q )
229                ENDIF
230                IF ( rif(j,i) >= 0.0_wp )  THEN
231!
232!--                Stable stratification
233                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (             &
234                                  LOG( z_p / z0h(j,i) ) +                      &
235                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
236                                                                 )
237                ELSE
238!
239!--                Unstable stratification
240                   a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 
241                   b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 
242 
243                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (           &
244                             LOG( z_p / z0h(j,i) ) -                           &
245                              2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) )
246                ENDIF
247
248             ENDDO
249          ENDDO
250       ENDIF
251    ENDIF
252
253!
254!-- Compute z_p/L (corresponds to the Richardson-flux number)
255    IF ( .NOT. humidity )  THEN
256       !$OMP PARALLEL DO PRIVATE( k, z_p )
257       !$acc kernels loop
258       DO  i = nxlg, nxrg
259          DO  j = nysg, nyng
260             k   = nzb_s_inner(j,i)
261             z_p = zu(k+1) - zw(k)
262             rif(j,i) = z_p * kappa * g * ts(j,i) / &
263                        ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
264!
265!--          Limit the value range of the Richardson numbers.
266!--          This is necessary for very small velocities (u,v --> 0), because
267!--          the absolute value of rif can then become very large, which in
268!--          consequence would result in very large shear stresses and very
269!--          small momentum fluxes (both are generally unrealistic).
270             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
271             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
272          ENDDO
273       ENDDO
274    ELSE
275       !$OMP PARALLEL DO PRIVATE( k, z_p )
276       !$acc kernels loop
277       DO  i = nxlg, nxrg
278          DO  j = nysg, nyng
279             k   = nzb_s_inner(j,i)
280             z_p = zu(k+1) - zw(k)
281             rif(j,i) = z_p * kappa * g *                                      &
282                        ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) + 0.61_wp  &
283                          * q(k+1,j,i) * ts(j,i)) /                            &
284                        ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
285!
286!--          Limit the value range of the Richardson numbers.
287!--          This is necessary for very small velocities (u,v --> 0), because
288!--          the absolute value of rif can then become very large, which in
289!--          consequence would result in very large shear stresses and very
290!--          small momentum fluxes (both are generally unrealistic).
291             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
292             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
293          ENDDO
294       ENDDO       
295    ENDIF
296
297!
298!-- Compute u* at the scalars' grid points
299    !$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
300    !$acc kernels loop
301    DO  i = nxl, nxr
302       DO  j = nys, nyn
303
304          k   = nzb_s_inner(j,i)
305          z_p = zu(k+1) - zw(k)
306
307!
308!--       Compute the absolute value of the horizontal velocity
309!--       (relative to the surface)
310          uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)              &
311                                      - u(k,j,i)   - u(k,j,i+1) ) )**2 +       &
312                           ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)              &
313                                      - v(k,j,i)   - v(k,j+1,i) ) )**2 )   
314
315
316          IF ( rif(j,i) >= 0.0_wp )  THEN
317!
318!--          Stable stratification
319             us(j,i) = kappa * uv_total / (                                    &
320                                  LOG( z_p / z0(j,i) ) +                       &
321                                  5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p  &
322                                          )
323          ELSE
324!
325!--          Unstable stratification
326             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) )
327             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) / z_p * z0(j,i) ) )
328
329             us(j,i) = kappa * uv_total / (                                    &
330                       LOG( z_p / z0(j,i) ) -                                  &
331                       LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
332                            ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
333                               2.0_wp * ( ATAN( a ) - ATAN( b ) )              &
334                                           )
335          ENDIF
336       ENDDO
337    ENDDO
338
339!
340!-- Values of us at ghost point locations are needed for the evaluation of usws
341!-- and vsws.
342    !$acc update host( us )
343    CALL exchange_horiz_2d( us )
344    !$acc update device( us )
345
346!
347!-- Compute u'w' for the total model domain.
348!-- First compute the corresponding component of u* and square it.
349    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
350    !$acc kernels loop
351    DO  i = nxl, nxr
352       DO  j = nys, nyn
353
354          k   = nzb_u_inner(j,i)
355          z_p = zu(k+1) - zw(k)
356
357!
358!--       Compute Richardson-flux number for this point
359          rifm = 0.5_wp * ( rif(j,i-1) + rif(j,i) )
360          IF ( rifm >= 0.0_wp )  THEN
361!
362!--          Stable stratification
363             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ (                  &
364                                     LOG( z_p / z0(j,i) ) +                    &
365                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p   &
366                                                            )
367          ELSE
368!
369!--          Unstable stratification
370             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
371             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
372
373             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / (                 &
374                         LOG( z_p / z0(j,i) ) -                                &
375                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (         &
376                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +     &
377                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )            &
378                                                 )
379          ENDIF
380          usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
381       ENDDO
382    ENDDO
383
384!
385!-- Compute v'w' for the total model domain.
386!-- First compute the corresponding component of u* and square it.
387    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
388    !$acc kernels loop
389    DO  i = nxl, nxr
390       DO  j = nys, nyn
391
392          k   = nzb_v_inner(j,i)
393          z_p = zu(k+1) - zw(k)
394
395!
396!--       Compute Richardson-flux number for this point
397          rifm = 0.5_wp * ( rif(j-1,i) + rif(j,i) )
398          IF ( rifm >= 0.0_wp )  THEN
399!
400!--          Stable stratification
401             vsws(j,i) = kappa * ( v(k+1,j,i) -  v(k,j,i) ) / (                &
402                                     LOG( z_p / z0(j,i) ) +                    &
403                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p   &
404                                                              )
405          ELSE
406!
407!--          Unstable stratification
408             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
409             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
410
411             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / (                 &
412                         LOG( z_p / z0(j,i) ) -                                &
413                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (         &
414                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +     &
415                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )            &
416                                                 )
417          ENDIF
418          vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j-1,i) + us(j,i) )
419       ENDDO
420    ENDDO
421
422!
423!-- If required compute qr* and nr*
424    IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation )  THEN
425
426       !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
427       !$acc kernels loop independent
428       DO  i = nxlg, nxrg
429          !$acc loop independent
430          DO  j = nysg, nyng
431
432             k   = nzb_s_inner(j,i)
433             z_p = zu(k+1) - zw(k)
434
435             IF ( rif(j,i) >= 0.0 )  THEN
436!
437!--             Stable stratification
438                qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / (             &
439                               LOG( z_p / z0h(j,i) ) +                         &
440                               5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
441                nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / (             &
442                               LOG( z_p / z0h(j,i) ) +                         &
443                               5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
444
445             ELSE
446!
447!--             Unstable stratification
448                a = SQRT( 1.0 - 16.0 * rif(j,i) ) 
449                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 
450 
451                qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) /   (           &
452                          LOG( z_p / z0h(j,i) ) -                              &
453                          2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
454                nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) /   (           &
455                          LOG( z_p / z0h(j,i) ) -                              &
456                          2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
457
458             ENDIF
459
460          ENDDO
461       ENDDO
462
463    ENDIF
464
465!
466!-- Exchange the boundaries for the momentum fluxes (only for sake of
467!-- completeness)
468    !$acc update host( usws, vsws )
469    CALL exchange_horiz_2d( usws )
470    CALL exchange_horiz_2d( vsws )
471    !$acc update device( usws, vsws )
472    IF ( humidity  .OR.  passive_scalar )  THEN
473       !$acc update host( qsws )
474       CALL exchange_horiz_2d( qsws )
475       !$acc update device( qsws )
476       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
477            precipitation )  THEN
478          !$acc update host( qrsws, nrsws )
479          CALL exchange_horiz_2d( qrsws )
480          CALL exchange_horiz_2d( nrsws )
481          !$acc update device( qrsws, nrsws )
482       ENDIF
483    ENDIF
484
485!
486!-- Compute the vertical kinematic heat flux
487    IF ( .NOT. constant_heatflux  )  THEN
488       !$OMP PARALLEL DO
489       !$acc kernels loop independent
490       DO  i = nxlg, nxrg
491          !$acc loop independent
492          DO  j = nysg, nyng
493             shf(j,i) = -ts(j,i) * us(j,i)
494          ENDDO
495       ENDDO
496    ENDIF
497
498!
499!-- Compute the vertical water/scalar flux
500    IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar ) )  THEN
501       !$OMP PARALLEL DO
502       !$acc kernels loop independent
503       DO  i = nxlg, nxrg
504          !$acc loop independent
505          DO  j = nysg, nyng
506             qsws(j,i) = -qs(j,i) * us(j,i)
507          ENDDO
508       ENDDO
509    ENDIF
510
511!
512!-- Compute (turbulent) fluxes of rain water content and rain drop concentartion
513    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
514       !$OMP PARALLEL DO
515       !$acc kernels loop independent
516       DO  i = nxlg, nxrg
517          !$acc loop independent
518          DO  j = nysg, nyng
519             qrsws(j,i) = -qrs(j,i) * us(j,i)
520             nrsws(j,i) = -nrs(j,i) * us(j,i)
521          ENDDO
522       ENDDO
523    ENDIF
524
525!
526!-- Bottom boundary condition for the TKE
527    IF ( ibc_e_b == 2 )  THEN
528       !$OMP PARALLEL DO
529       !$acc kernels loop independent
530       DO  i = nxlg, nxrg
531          !$acc loop independent
532          DO  j = nysg, nyng
533             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2
534!
535!--          As a test: cm = 0.4
536!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2
537             e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
538          ENDDO
539       ENDDO
540    ENDIF
541
542    !$acc end data
543
544 END SUBROUTINE prandtl_fluxes
Note: See TracBrowser for help on using the repository browser.