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

Last change on this file since 1496 was 1496, checked in by maronga, 9 years ago

added beta version of a land surface model and a simple radiation model for clear sky conditions

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