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

Last change on this file since 1551 was 1551, checked in by maronga, 7 years ago

land surface model released

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