source: palm/trunk/SOURCE/surface_layer_fluxes.f90 @ 1749

Last change on this file since 1749 was 1749, checked in by raasch, 6 years ago

further OpenACC adjustments and one bugfix

  • Property svn:keywords set to Id
File size: 42.4 KB
Line 
1!> @file surface_layer_fluxes.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2015 Leibniz Universitaet Hannover
17!
18!--------------------------------------------------------------------------------!
19! Current revisions:
20! ------------------
21! further OpenACC adjustments
22!
23! Former revisions:
24! -----------------
25! $Id: surface_layer_fluxes.f90 1749 2016-02-09 12:19:56Z raasch $
26!
27! 1747 2016-02-08 12:25:53Z raasch
28! adjustments for OpenACC usage
29!
30! 1709 2015-11-04 14:47:01Z maronga
31! Bugfix: division by zero could occur when calculating rib at low wind speeds
32! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
33! neutral = .T.
34!
35! 1705 2015-11-02 14:28:56Z maronga
36! Typo removed
37!
38! 1697 2015-10-28 17:14:10Z raasch
39! FORTRAN and OpenMP errors removed
40!
41! 1696 2015-10-27 10:03:34Z maronga
42! Modularized and completely re-written version of prandtl_fluxes.f90. In the
43! course of the re-writing two additional methods have been implemented. See
44! updated description.
45!
46! 1551 2015-03-03 14:18:16Z maronga
47! Removed land surface model part. The surface fluxes are now always calculated
48! within prandtl_fluxes, based on the given surface temperature/humidity (which
49! is either provided by the land surface model, by large scale forcing data, or
50! directly prescribed by the user.
51!
52! 1496 2014-12-02 17:25:50Z maronga
53! Adapted for land surface model
54!
55! 1494 2014-11-21 17:14:03Z maronga
56! Bugfixes: qs is now calculated before calculation of Rif. calculation of
57! buoyancy flux in Rif corrected (added missing humidity term), allow use of
58! topography for coupled runs (not tested)
59!
60! 1361 2014-04-16 15:17:48Z hoffmann
61! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
62! drop concentration (nrsws) added
63!
64! 1340 2014-03-25 19:45:13Z kanani
65! REAL constants defined as wp-kind
66!
67! 1320 2014-03-20 08:40:49Z raasch
68! ONLY-attribute added to USE-statements,
69! kind-parameters added to all INTEGER and REAL declaration statements,
70! kinds are defined in new module kinds,
71! old module precision_kind is removed,
72! revision history before 2012 removed,
73! comment fields (!:) to be used for variable explanations added to
74! all variable declaration statements
75!
76! 1276 2014-01-15 13:40:41Z heinze
77! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
78!
79! 1257 2013-11-08 15:18:40Z raasch
80! openACC "kernels do" replaced by "kernels loop", "loop independent" added
81!
82! 1036 2012-10-22 13:43:42Z raasch
83! code put under GPL (PALM 3.9)
84!
85! 1015 2012-09-27 09:23:24Z raasch
86! OpenACC statements added
87!
88! 978 2012-08-09 08:28:32Z fricke
89! roughness length for scalar quantities z0h added
90!
91! Revision 1.1  1998/01/23 10:06:06  raasch
92! Initial revision
93!
94!
95! Description:
96! ------------
97!> Diagnostic computation of vertical fluxes in the constant flux layer from the
98!> values of the variables at grid point k=1. Three different methods are
99!> available:
100!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
101!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
102!>    slower
103!> 3) a method using a lookup table which is fast and accurate. Note, however,
104!>    that this method cannot be used in case of roughness heterogeneity
105!>
106!> @todo (re)move large_scale_forcing actions
107!> @todo check/optimize OpenMP and OpenACC directives
108!------------------------------------------------------------------------------!
109 MODULE surface_layer_fluxes_mod
110
111    USE arrays_3d,                                                             &
112        ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
113               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
114
115    USE cloud_parameters,                                                      &
116        ONLY:  l_d_cp, pt_d_t
117
118    USE constants,                                                             &
119        ONLY:  pi
120
121    USE cpulog
122
123    USE control_parameters,                                                    &
124        ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
125               coupling_mode, g, humidity, ibc_e_b, ibc_pt_b, icloud_scheme,   &
126               initializing_actions, kappa, intermediate_timestep_count,       &
127               intermediate_timestep_count_max, large_scale_forcing, lsf_surf, &
128               message_string, most_method, neutral, passive_scalar,           &
129               precipitation, pt_surface, q_surface, run_coupled,              &
130               surface_pressure, simulated_time, terminate_run, zeta_max,      &
131               zeta_min
132
133    USE indices,                                                               &
134        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner,        &
135               nzb_u_inner, nzb_v_inner
136
137    USE kinds
138
139    USE pegrid
140
141    USE land_surface_model_mod,                                                &
142        ONLY:  land_surface, skip_time_do_lsm
143
144    IMPLICIT NONE
145
146    INTEGER(iwp) ::  i            !< loop index x direction
147    INTEGER(iwp) ::  j            !< loop index y direction
148    INTEGER(iwp) ::  k            !< loop index z direction
149
150    INTEGER(iwp), PARAMETER     :: num_steps = 15000  !< number of steps in the lookup table
151
152    LOGICAL      ::  coupled_run  !< Flag for coupled atmosphere-ocean runs
153
154    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt1,      & !< Potential temperature at first grid level (required for cloud_physics = .T.)
155                                             qv1,      & !< Specific humidity at first grid level (required for cloud_physics = .T.)
156                                             uv_total    !< Total velocity at first grid level
157
158    REAL(wp), DIMENSION(0:num_steps-1) :: rib_tab,  & !< Lookup table bulk Richardson number
159                                          ol_tab      !< Lookup table values of L
160
161    REAL(wp)     ::  e_s,               & !< Saturation water vapor pressure
162                     l_bnd  = 7500,     & !< Lookup table index of the last time step
163                     ol_max = 1.0E6_wp, & !< Maximum Obukhov length
164                     rib_max,           & !< Maximum Richardson number in lookup table
165                     rib_min,           & !< Minimum Richardson number in lookup table
166                     z_mo                 !< Height of the constant flux layer where MOST is assumed
167
168
169    SAVE
170
171    PRIVATE
172
173    PUBLIC init_surface_layer_fluxes, pt1, qv1, surface_layer_fluxes, uv_total
174
175    INTERFACE init_surface_layer_fluxes
176       MODULE PROCEDURE init_surface_layer_fluxes
177    END INTERFACE init_surface_layer_fluxes
178
179    INTERFACE surface_layer_fluxes
180       MODULE PROCEDURE surface_layer_fluxes
181    END INTERFACE surface_layer_fluxes
182
183
184 CONTAINS
185
186
187!------------------------------------------------------------------------------!
188! Description:
189! ------------
190!> Main routine to compute the surface fluxes
191!------------------------------------------------------------------------------!
192    SUBROUTINE surface_layer_fluxes
193
194       IMPLICIT NONE
195
196!
197!--    In case cloud physics is used, it is required to derive potential
198!--    temperature and specific humidity at first grid level from the fields pt
199!--    and q
200       IF ( cloud_physics )  THEN
201          CALL calc_pt_q
202       ENDIF
203
204!
205!--    First, calculate the new Obukhov length, then new friction velocity,
206!--    followed by the new scaling parameters (th*, q*, etc.), and the new
207!--    surface fluxes if required. The old routine ("circular") requires a
208!--    different order of calls as the scaling parameters from the previous time
209!--    steps are used to calculate the Obukhov length
210
211!
212!--    Depending on setting of most_method use the "old" routine
213       IF ( most_method == 'circular' )  THEN
214
215          CALL calc_scaling_parameters
216
217          CALL calc_uv_total
218
219          IF ( .NOT. neutral )  THEN
220             CALL calc_ol
221          ENDIF
222
223          CALL calc_us
224
225          CALL calc_surface_fluxes
226
227!
228!--    Use either Newton iteration or a lookup table for the bulk Richardson
229!--    number to calculate the Obukhov length
230       ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' )  THEN
231
232          CALL calc_uv_total
233
234          IF ( .NOT. neutral )  THEN
235             CALL calc_ol
236          ENDIF
237
238          CALL calc_us
239
240          CALL calc_scaling_parameters
241
242          CALL calc_surface_fluxes
243
244       ENDIF
245
246    END SUBROUTINE surface_layer_fluxes
247
248
249!------------------------------------------------------------------------------!
250! Description:
251! ------------
252!> Initializing actions for the surface layer routine. Basically, this involves
253!> the preparation of a lookup table for the the bulk Richardson number vs
254!> Obukhov length L when using the lookup table method.
255!------------------------------------------------------------------------------!
256    SUBROUTINE init_surface_layer_fluxes
257
258       IMPLICIT NONE
259
260       INTEGER(iwp) :: l,          & !< Index for loop to create lookup table
261                       num_steps_n   !< Number of non-stretched zeta steps
262
263       LOGICAL :: terminate_run_l = .FALSE.    !< Flag to terminate run (global)
264
265       REAL(wp), PARAMETER ::  zeta_stretch = -10.0_wp !< Start of stretching in the free convection limit
266                               
267       REAL(wp), DIMENSION(:), ALLOCATABLE :: zeta_tmp
268
269
270       REAL(wp) :: zeta_step,            & !< Increment of zeta
271                   regr      = 1.01_wp,  & !< Stretching factor of zeta_step in the free convection limit
272                   regr_old  = 1.0E9_wp, & !< Stretching factor of last iteration step
273                   z0h_min   = 0.0_wp,   & !< Minimum value of z0h to create table
274                   z0_min    = 0.0_wp      !< Minimum value of z0 to create table
275!
276!--    When cloud physics is used, arrays for storing potential temperature and
277!--    specific humidity at first grid level are required
278       IF ( cloud_physics )  THEN
279          ALLOCATE ( pt1(nysg:nyng,nxlg:nxrg) )
280          ALLOCATE ( qv1(nysg:nyng,nxlg:nxrg) )
281       ENDIF
282
283!
284!--    Allocate field for storing the horizontal velocity
285       ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) )
286
287
288!
289!--    In case of runs with neutral statification, set Obukhov length to a
290!--    large value
291       IF ( neutral ) ol = 1.0E10_wp
292
293       IF ( most_method == 'lookup' )  THEN
294
295!
296!--       Check for roughness heterogeneity. In that case terminate run and
297!--       inform user
298          IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR.                             &
299               MINVAL( z0  ) /= MAXVAL( z0  ) )  THEN
300             terminate_run_l = .TRUE.
301          ENDIF
302
303#if defined( __parallel )
304!
305!--       Make a logical OR for all processes. Force termiation of model if result
306!--       is TRUE
307          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
308          CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,  &
309                              MPI_LOR, comm2d, ierr )
310#else
311          terminate_run = terminate_run_l
312#endif
313
314          IF ( terminate_run )  THEN
315             message_string = 'most_method = "lookup" cannot be used in ' //   &
316                              'combination with a prescribed roughness '  //   &
317                              'heterogeneity'
318             CALL message( 'surface_layer_fluxes', 'PA0417', 1, 2, 0, 6, 0 )
319          ENDIF
320
321          ALLOCATE(  zeta_tmp(0:num_steps-1) )
322
323!
324!--       Use the lowest possible value for z_mo
325          k    = MINVAL(nzb_s_inner)
326          z_mo = zu(k+1) - zw(k)
327
328!
329!--       Calculate z/L range from zeta_stretch to zeta_max using 90% of the
330!--       available steps (num_steps). The calculation is done with negative
331!--       values of zeta in order to simplify the stretching in the free
332!--       convection limit for the remaining 10% of steps.
333          zeta_tmp(0) = - zeta_max
334          num_steps_n = ( num_steps * 9 / 10 ) - 1
335          zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
336
337          DO l = 1, num_steps_n
338             zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
339          ENDDO
340
341!
342!--       Calculate stretching factor for the free convection range
343          DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
344             regr_old = regr
345             regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
346                    )**( 10.0_wp / REAL(num_steps) )
347          ENDDO
348
349!
350!--       Calculate z/L range from zeta_min to zeta_stretch
351          DO l = num_steps_n+1, num_steps-1
352             zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
353             zeta_step = zeta_step * regr
354          ENDDO
355
356!
357!--       Invert array and switch sign, then calculate Obukhov length and bulk
358!--       Richardson number
359          z0h_min = MINVAL(z0h)
360          z0_min  = MINVAL(z0)
361         
362!
363!--       Calculate lookup table for the Richardson number versus Obukhov length
364!--       The Richardson number (rib) is defined depending on the choice of
365!--       boundary conditions for temperature
366          IF ( ibc_pt_b == 1 )  THEN
367             DO l = 0, num_steps-1
368                ol_tab(l)  = - z_mo / zeta_tmp(num_steps-1-l)
369                rib_tab(l) = z_mo / ol_tab(l)  / ( LOG( z_mo / z0_min )        &
370                                                - psi_m( z_mo / ol_tab(l) )    &
371                                                + psi_m( z0_min / ol_tab(l) )  &
372                                                  )**3
373             ENDDO 
374          ELSE
375             DO l = 0, num_steps-1
376                ol_tab(l)  = - z_mo / zeta_tmp(num_steps-1-l)
377                rib_tab(l) = z_mo / ol_tab(l)  * ( LOG( z_mo / z0h_min )       &
378                                              - psi_h( z_mo / ol_tab(l) )      &
379                                              + psi_h( z0h_min / ol_tab(l) )   &
380                                            )                                  &
381                                          / ( LOG( z_mo / z0_min )             &
382                                              - psi_m( z_mo / ol_tab(l) )      &
383                                              + psi_m( z0_min / ol_tab(l) )    &
384                                            )**2
385             ENDDO
386          ENDIF
387
388!
389!--       Determine minimum values of rib in the lookup table. Set upper limit
390!--       to critical Richardson number (0.25)
391          rib_min  = MINVAL(rib_tab)
392          rib_max  = 0.25 !MAXVAL(rib_tab)
393
394          DEALLOCATE( zeta_tmp )
395       ENDIF
396
397    END SUBROUTINE init_surface_layer_fluxes
398
399
400!------------------------------------------------------------------------------!
401! Description:
402! ------------
403!> Compute the absolute value of the horizontal velocity (relative to the
404!> surface). This is required by all methods
405!------------------------------------------------------------------------------!
406    SUBROUTINE calc_uv_total
407
408       IMPLICIT NONE
409
410
411       !$OMP PARALLEL DO PRIVATE( k )
412       !$acc kernels loop present( nzb_s_inner, u, uv_total, v ) private( j, k )
413       DO  i = nxl, nxr
414          DO  j = nys, nyn
415
416             k   = nzb_s_inner(j,i)
417             uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
418                                         - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
419                              ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
420                                         - v(k,j,i)   - v(k,j+1,i) ) )**2 )
421
422!
423!--          For too small values of the local wind, MOST does not work. A
424!--          threshold value is thus set if required
425!            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
426
427          ENDDO
428       ENDDO
429
430!
431!--    Values of uv_total need to be exchanged at the ghost boundaries
432       !$acc update host( uv_total )
433       CALL exchange_horiz_2d( uv_total )
434       !$acc update device( uv_total )
435
436    END SUBROUTINE calc_uv_total
437
438
439!------------------------------------------------------------------------------!
440! Description:
441! ------------
442!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
443!------------------------------------------------------------------------------!
444    SUBROUTINE calc_ol
445
446       IMPLICIT NONE
447
448       INTEGER(iwp) :: iter,  & !< Newton iteration step
449                       l        !< look index
450
451       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number
452
453       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
454                       f_d_ol, & !< Derivative of f
455                       ol_l,   & !< Lower bound of L for Newton iteration
456                       ol_m,   & !< Previous value of L for Newton iteration
457                       ol_old, & !< Previous time step value of L
458                       ol_u      !< Upper bound of L for Newton iteration
459
460       IF ( TRIM( most_method ) /= 'circular' )  THEN
461     
462          !$acc data present( nzb_s_inner, pt, q, qsws, rib, shf, uv_total, vpt, zu, zw )
463
464          !$OMP PARALLEL DO PRIVATE( k, z_mo )
465          !$acc kernels loop private( j, k, z_mo )
466          DO  i = nxl, nxr
467             DO  j = nys, nyn
468
469                k   = nzb_s_inner(j,i)
470                z_mo = zu(k+1) - zw(k)
471
472!
473!--             Evaluate bulk Richardson number (calculation depends on
474!--             definition based on setting of boundary conditions
475                IF ( ibc_pt_b /= 1 )  THEN
476                   IF ( humidity )  THEN
477                      rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) )      &
478                           / ( uv_total(j,i)**2 * vpt(k+1,j,i) + 1.0E-20_wp )
479                   ELSE
480                      rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) )        &
481                           / ( uv_total(j,i)**2 * pt(k+1,j,i)  + 1.0E-20_wp )
482                   ENDIF     
483                ELSE
484!
485!--                When using Neumann boundary conditions, the buoyancy flux
486!--                is required but cannot be calculated at the surface, as pt
487!--                and q are not known at the surface. Hence the values at
488!--                first grid level are used to estimate the buoyancy flux
489                   IF ( humidity )  THEN
490                      rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
491                                 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
492                                 * pt(k+1,j,i) * qsws(j,i) )   &
493                                 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2&
494                                 + 1.0E-20_wp)
495                   ELSE
496                      rib(j,i) = - g * z_mo * shf(j,i)                         &
497                           / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2       &
498                           + 1.0E-20_wp )
499                   ENDIF
500                ENDIF 
501     
502             ENDDO
503          ENDDO
504          !$acc end data
505
506       ENDIF
507
508!
509!--    Calculate the Obukhov length either using a Newton iteration
510!--    method, via a lookup table, or using the old circular way
511       IF ( TRIM( most_method ) == 'newton' )  THEN
512
513          !$OMP PARALLEL DO PRIVATE( k, z_mo )
514          !# WARNING: does not work on GPU so far because of DO-loop with
515          !#          undetermined iterations
516          !!!!!!$acc kernels loop
517          DO  i = nxl, nxr
518             DO  j = nys, nyn
519
520                k   = nzb_s_inner(j,i)
521                z_mo = zu(k+1) - zw(k)
522
523!
524!--             Store current value in case the Newton iteration fails
525                ol_old = ol(j,i)
526
527!
528!--             Ensure that the bulk Richardson number and the Obukhov
529!--             lengtH have the same sign
530                IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR.                        &
531                     ABS( ol(j,i) ) == ol_max )  THEN
532                   IF ( rib(j,i) > 0.0_wp ) ol(j,i) =  0.01_wp
533                   IF ( rib(j,i) < 0.0_wp ) ol(j,i) = -0.01_wp
534                ENDIF
535!
536!--             Iteration to find Obukhov length
537                iter = 0
538                DO
539                   iter = iter + 1
540!
541!--                In case of divergence, use the value of the previous time step
542                   IF ( iter > 1000 )  THEN
543                      ol(j,i) = ol_old
544                      EXIT
545                   ENDIF
546
547                   ol_m = ol(j,i)
548                   ol_l = ol_m - 0.001_wp * ol_m
549                   ol_u = ol_m + 0.001_wp * ol_m
550
551
552                   IF ( ibc_pt_b /= 1 )  THEN
553!
554!--                   Calculate f = Ri - [...]/[...]^2 = 0
555                      f = rib(j,i) - ( z_mo / ol_m ) * ( LOG( z_mo / z0h(j,i) )&
556                                                    - psi_h( z_mo / ol_m )     &
557                                                    + psi_h( z0h(j,i) / ol_m ) &
558                                                   )                           &
559                                                 / ( LOG( z_mo / z0(j,i) )     &
560                                                    - psi_m( z_mo / ol_m )     &
561                                                    + psi_m( z0(j,i) / ol_m )  &
562                                                    )**2
563
564!
565!--                    Calculate df/dL
566                       f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / z0h(j,i) ) &
567                                                   - psi_h( z_mo / ol_u )      &
568                                                   + psi_h( z0h(j,i) / ol_u )  &
569                                                 )                             &
570                                               / ( LOG( z_mo / z0(j,i) )       &
571                                                   - psi_m( z_mo / ol_u )      &
572                                                   + psi_m( z0(j,i) / ol_u )   &
573                                                 )**2                          &
574                              + ( z_mo / ol_l ) * ( LOG( z_mo / z0h(j,i) )     &
575                                                   - psi_h( z_mo / ol_l )      &
576                                                   + psi_h( z0h(j,i) / ol_l )  &
577                                                 )                             &
578                                               / ( LOG( z_mo / z0(j,i) )       &
579                                                   - psi_m( z_mo / ol_l )      &
580                                                   + psi_m( z0(j,i) / ol_l )   &
581                                                 )**2                          &
582                                ) / ( ol_u - ol_l )
583                   ELSE
584!
585!--                   Calculate f = Ri - 1 /[...]^3 = 0
586                      f = rib(j,i) - ( z_mo / ol_m ) / ( LOG( z_mo / z0(j,i) )&
587                                                    - psi_m( z_mo / ol_m )    &
588                                                    + psi_m( z0(j,i) / ol_m ) &
589                                                       )**3
590
591!
592!--                   Calculate df/dL
593                      f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / z0(j,i) )  &
594                                                   - psi_m( z_mo / ol_u )     &
595                                                   + psi_m( z0(j,i) / ol_u )  &
596                                                 )**3                         &
597                              + ( z_mo / ol_l ) / ( LOG( z_mo / z0(j,i) )     &
598                                                   - psi_m( z_mo / ol_l )     &
599                                                   + psi_m( z0(j,i) / ol_l )  &
600                                                 )**3                         &
601                                     ) / ( ol_u - ol_l )
602                   ENDIF
603!
604!--                Calculate new L
605                   ol(j,i) = ol_m - f / f_d_ol
606
607!
608!--                Ensure that the bulk Richardson number and the Obukhov
609!--                length have the same sign and ensure convergence.
610                   IF ( ol(j,i) * ol_m < 0.0_wp )  ol(j,i) = ol_m * 0.5_wp
611
612!
613!--                If unrealistic value occurs, set L to the maximum
614!--                value that is allowed
615                   IF ( ABS( ol(j,i) ) > ol_max )  THEN
616                      ol(j,i) = ol_max
617                      EXIT
618                   ENDIF
619!
620!--                Check for convergence
621                   IF ( ABS( ( ol(j,i) - ol_m ) / ol(j,i) ) < 1.0E-4_wp )  THEN
622                      EXIT
623                   ELSE
624                      CYCLE
625                   ENDIF
626
627                ENDDO
628                       
629             ENDDO
630          ENDDO
631
632       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
633
634          !$OMP PARALLEL DO PRIVATE( k, z_mo )
635          !# WARNING: does not work on GPU so far because of DO WHILE construct
636          !!!!!!$acc kernels loop
637          DO  i = nxl, nxr
638             DO  j = nys, nyn
639
640                k   = nzb_s_inner(j,i)
641
642!
643!--             If the bulk Richardson number is outside the range of the lookup
644!--             table, set it to the exceeding threshold value
645                IF ( rib(j,i) < rib_min )  rib(j,i) = rib_min
646                IF ( rib(j,i) > rib_max )  rib(j,i) = rib_max
647
648!
649!--             Find the correct index bounds for linear interpolation. As the
650!--             Richardson number will not differ very much from time step to
651!--             time step , use the index from the last step and search in the
652!--             correct direction
653                l = l_bnd
654                IF ( rib_tab(l) - rib(j,i) > 0.0_wp )  THEN
655                   DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND. l > 0 )
656                      l = l-1
657                   ENDDO
658                ELSE
659                   DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                &
660                              .AND. l < num_steps-1 )
661                      l = l+1
662                   ENDDO
663                ENDIF
664                l_bnd = l
665
666!
667!--             Linear interpolation to find the correct value of z/L
668                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )       &
669                            / (  rib_tab(l) - rib_tab(l-1) )                &
670                            * ( rib(j,i) - rib_tab(l-1) ) )
671
672             ENDDO
673          ENDDO
674
675       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
676
677          !$OMP PARALLEL DO PRIVATE( k, z_mo )
678          !$acc kernels loop present( nzb_s_inner, ol, pt, pt1, q, ql, qs, qv1, ts, us, vpt, zu, zw ) private( j, k, z_mo )
679          DO  i = nxl, nxr
680             DO  j = nys, nyn
681
682                k   = nzb_s_inner(j,i)
683                z_mo = zu(k+1) - zw(k)
684
685                IF ( .NOT. humidity )  THEN
686                   ol(j,i) =  ( pt(k+1,j,i) *  us(j,i)**2 ) / ( kappa * g      &
687                              * ts(j,i) + 1E-30_wp )
688                ELSEIF ( cloud_physics )  THEN
689
690                   ol(j,i) =  ( vpt(k+1,j,i) * us(j,i)**2 ) / ( kappa * g      &
691                              * ( ts(j,i) + 0.61_wp * pt1(j,i) * qs(j,i)       &
692                              + 0.61_wp * qv1(j,i) * ts(j,i) - ts(j,i)         &
693                              * ql(k+1,j,i) ) + 1E-30_wp )
694                ELSE
695                   ol(j,i) =  ( vpt(k+1,j,i) *  us(j,i)**2 ) / ( kappa * g     &
696                              * ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i)    &
697                                  + 0.61_wp * q(k+1,j,i) * ts(j,i) ) + 1E-30_wp )
698                ENDIF
699!
700!--             Limit the value range of the Obukhov length.
701!--             This is necessary for very small velocities (u,v --> 0), because
702!--             the absolute value of ol can then become very small, which in
703!--             consequence would result in very large shear stresses and very
704!--             small momentum fluxes (both are generally unrealistic).
705                IF ( ( z_mo / ol(j,i) ) < zeta_min )  ol(j,i) = z_mo / zeta_min
706                IF ( ( z_mo / ol(j,i) ) > zeta_max )  ol(j,i) = z_mo / zeta_max
707
708             ENDDO
709          ENDDO
710
711       ENDIF
712
713!
714!--    Values of ol at ghost point locations are needed for the evaluation
715!--    of usws and vsws.
716       !$acc update host( ol )
717       CALL exchange_horiz_2d( ol )
718       !$acc update device( ol )
719
720    END SUBROUTINE calc_ol
721
722!
723!-- Calculate friction velocity u*
724    SUBROUTINE calc_us
725
726       IMPLICIT NONE
727
728       !$OMP PARALLEL DO PRIVATE( k, z_mo )
729       !$acc kernels loop present( nzb_s_inner, ol, us, uv_total, zu, zw, z0 ) private( j, k, z_mo )
730       DO  i = nxlg, nxrg
731          DO  j = nysg, nyng
732
733             k   = nzb_s_inner(j,i)+1
734             z_mo = zu(k+1) - zw(k)
735
736!
737!--          Compute u* at the scalars' grid points
738             us(j,i) = kappa * uv_total(j,i) / ( LOG( z_mo / z0(j,i) )         &
739                                          - psi_m( z_mo / ol(j,i) )            &
740                                          + psi_m( z0(j,i) / ol(j,i) ) )
741          ENDDO
742       ENDDO
743
744    END SUBROUTINE calc_us
745
746!
747!-- Calculate potential temperature and specific humidity at first grid level
748    SUBROUTINE calc_pt_q
749
750       IMPLICIT NONE
751
752       !$acc kernels loop present( nzb_s_inner, pt, pt1, pt_d_t, q, ql, qv1 ) private( j, k )
753       DO  i = nxlg, nxrg
754          DO  j = nysg, nyng
755             k   = nzb_s_inner(j,i)+1
756             pt1(j,i) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
757             qv1(j,i) = q(k,j,i) - ql(k,j,i)
758          ENDDO
759       ENDDO
760
761    END SUBROUTINE calc_pt_q
762
763!
764!-- Calculate the other MOST scaling parameters theta*, q*, (qr*, nr*)
765    SUBROUTINE calc_scaling_parameters
766
767       IMPLICIT NONE
768
769!
770!--    Data information for accelerators
771       !$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt )  &
772       !$acc      present( q, qs, qsws, qrsws, shf, ts, u, us, usws, v )     &
773       !$acc      present( vpt, vsws, zu, zw, z0, z0h )
774!
775!--    Compute theta*
776       IF ( constant_heatflux )  THEN
777
778!
779!--       For a given heat flux in the surface layer:
780          !$OMP PARALLEL DO
781          !$acc kernels loop private( j )
782          DO  i = nxlg, nxrg
783             DO  j = nysg, nyng
784                ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
785!
786!--             ts must be limited, because otherwise overflow may occur in case
787!--             of us=0 when computing ol further below
788                IF ( ts(j,i) < -1.05E5_wp )  ts(j,i) = -1.0E5_wp
789                IF ( ts(j,i) >   1.0E5_wp )  ts(j,i) =  1.0E5_wp
790             ENDDO
791          ENDDO
792
793       ELSE
794!
795!--       For a given surface temperature:
796          IF ( large_scale_forcing .AND. lsf_surf )  THEN
797             !$OMP PARALLEL DO
798             !$acc kernels loop private( j, k )
799             DO  i = nxlg, nxrg
800                DO  j = nysg, nyng
801                   k = nzb_s_inner(j,i)
802                   pt(k,j,i) = pt_surface
803                ENDDO
804             ENDDO
805          ENDIF
806
807          !$OMP PARALLEL DO PRIVATE( k, z_mo )
808          !$acc kernels loop present( nzb_s_inner, ol, pt, pt1, ts, zu, zw, z0h ) private( j, k, z_mo )
809          DO  i = nxlg, nxrg
810             DO  j = nysg, nyng
811
812                k   = nzb_s_inner(j,i)
813                z_mo = zu(k+1) - zw(k)
814
815                IF ( cloud_physics )  THEN
816                   ts(j,i) = kappa * ( pt1(j,i) - pt(k,j,i) )                  &
817                                     / ( LOG( z_mo / z0h(j,i) )                &
818                                         - psi_h( z_mo / ol(j,i) )             &
819                                         + psi_h( z0h(j,i) / ol(j,i) ) )
820                ELSE
821                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) )               &
822                                     / ( LOG( z_mo / z0h(j,i) )                &
823                                         - psi_h( z_mo / ol(j,i) )             &
824                                         + psi_h( z0h(j,i) / ol(j,i) ) )
825                ENDIF
826
827             ENDDO
828          ENDDO
829       ENDIF
830
831!
832!--    If required compute q*
833       IF ( humidity  .OR.  passive_scalar )  THEN
834          IF ( constant_waterflux )  THEN
835!
836!--          For a given water flux in the Prandtl layer:
837             !$OMP PARALLEL DO
838             !$acc kernels loop private( j )
839             DO  i = nxlg, nxrg
840                DO  j = nysg, nyng
841                   qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
842                ENDDO
843             ENDDO
844
845          ELSE
846             coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND.      &
847                             run_coupled )
848
849             IF ( large_scale_forcing .AND. lsf_surf )  THEN
850                !$OMP PARALLEL DO
851                !$acc kernels loop private( j, k )
852                DO  i = nxlg, nxrg
853                   DO  j = nysg, nyng
854                      k = nzb_s_inner(j,i)
855                      q(k,j,i) = q_surface
856                   ENDDO
857                ENDDO
858             ENDIF
859
860             !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
861             !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0h ) private( e_s, j, k, z_mo )
862             DO  i = nxlg, nxrg
863                !$acc loop independent
864                DO  j = nysg, nyng
865
866                   k   = nzb_s_inner(j,i)
867                   z_mo = zu(k+1) - zw(k)
868
869!
870!--                Assume saturation for atmosphere coupled to ocean (but not
871!--                in case of precursor runs)
872                   IF ( coupled_run )  THEN
873                      e_s = 6.1_wp * &
874                              EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i))      &
875                                               - 273.15_wp ) )
876                      q(k,j,i) = 0.622_wp * e_s / ( surface_pressure - e_s )
877                   ENDIF
878
879                   IF ( cloud_physics )  THEN
880                      qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) )                &
881                                        / ( LOG( z_mo / z0h(j,i) )             &
882                                            - psi_h( z_mo / ol(j,i) )          &
883                                            + psi_h( z0h(j,i) / ol(j,i) ) )
884
885                   ELSE
886                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) )              &
887                                        / ( LOG( z_mo / z0h(j,i) )             &
888                                            - psi_h( z_mo / ol(j,i) )          &
889                                            + psi_h( z0h(j,i) / ol(j,i) ) )
890                   ENDIF
891
892                ENDDO
893             ENDDO
894          ENDIF
895       ENDIF
896
897
898!
899!--    If required compute qr* and nr*
900       IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation )  THEN
901
902          !$OMP PARALLEL DO PRIVATE( k, z_mo )
903          !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0h ) private( j, k, z_mo )
904          DO  i = nxlg, nxrg
905             !$acc loop independent
906             DO  j = nysg, nyng
907
908                k   = nzb_s_inner(j,i)
909                z_mo = zu(k+1) - zw(k)
910
911                qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) )              &
912                                 / ( LOG( z_mo / z0h(j,i) )                 &
913                                     - psi_h( z_mo / ol(j,i) )              &
914                                     + psi_h( z0h(j,i) / ol(j,i) ) )
915
916                nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) )              &
917                                 / ( LOG( z_mo / z0h(j,i) )                 &
918                                     - psi_h( z_mo / ol(j,i) )              &
919                                     + psi_h( z0h(j,i) / ol(j,i) ) )
920             ENDDO
921          ENDDO
922
923       ENDIF
924       !$acc end data
925
926    END SUBROUTINE calc_scaling_parameters
927
928
929
930!
931!-- Calculate surface fluxes usws, vsws, shf, qsws, (qrsws, nrsws)
932    SUBROUTINE calc_surface_fluxes
933
934       IMPLICIT NONE
935
936       REAL(wp) :: ol_mid !< Grid-interpolated L
937
938!
939!--    Compute u'w' for the total model domain.
940!--    First compute the corresponding component of u* and square it.
941       !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
942       !$acc kernels loop present( nzb_u_inner, ol, u, us, usws, zu, zw, z0 ) private( j, k, z_mo )
943       DO  i = nxl, nxr
944          DO  j = nys, nyn
945
946             k   = nzb_u_inner(j,i)
947             z_mo = zu(k+1) - zw(k)
948!
949!--          Compute bulk Obukhov length for this point
950             ol_mid = 0.5_wp * ( ol(j,i-1) + ol(j,i) )
951
952             IF ( ol_mid == 0.0_wp )  THEN
953                ol_mid = MIN(ol(j,i-1), ol(j,i))
954             ENDIF
955
956             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )                     &
957                                 / ( LOG( z_mo / z0(j,i) )                     &
958                                     - psi_m( z_mo / ol_mid )                  &
959                                     + psi_m( z0(j,i) / ol_mid ) )
960
961             usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
962          ENDDO
963       ENDDO
964
965!
966!--    Compute v'w' for the total model domain.
967!--    First compute the corresponding component of u* and square it.
968       !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
969       !$acc kernels loop present( nzb_v_inner, ol, v, us, vsws, zu, zw, z0 ) private( j, k, ol_mid, z_mo )
970       DO  i = nxl, nxr
971          DO  j = nys, nyn
972
973             k   = nzb_v_inner(j,i)
974             z_mo = zu(k+1) - zw(k)
975!
976!--          Compute bulk Obukhov length for this point
977             ol_mid = 0.5_wp * ( ol(j-1,i) + ol(j,i) )
978
979             IF ( ol_mid == 0.0_wp )  THEN
980                ol_mid = MIN(ol(j-1,i), ol(j-1,i))
981             ENDIF
982
983             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) )                     &
984                                 / ( LOG( z_mo / z0(j,i) )                     &
985                                     - psi_m( z_mo / ol_mid )                  &
986                                     + psi_m( z0(j,i) / ol_mid ) )
987
988             vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
989
990          ENDDO
991       ENDDO
992
993!
994!--    Exchange the boundaries for the momentum fluxes (is this still required?)
995       !$acc update host( usws, vsws )
996       CALL exchange_horiz_2d( usws )
997       CALL exchange_horiz_2d( vsws )
998       !$acc update device( usws, vsws )
999
1000!
1001!--    Compute the vertical kinematic heat flux
1002       IF ( .NOT. constant_heatflux .AND. ( simulated_time <=                &
1003            skip_time_do_lsm .OR. .NOT. land_surface ) )  THEN
1004          !$OMP PARALLEL DO
1005          !$acc kernels loop independent present( shf, ts, us )
1006          DO  i = nxlg, nxrg
1007             !$acc loop independent
1008             DO  j = nysg, nyng
1009                shf(j,i) = -ts(j,i) * us(j,i)
1010             ENDDO
1011          ENDDO
1012
1013       ENDIF
1014
1015!
1016!-- Compute the vertical water/scalar flux
1017       IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar )    &
1018            .AND. ( simulated_time <= skip_time_do_lsm .OR. .NOT.            &
1019            land_surface ) )  THEN
1020          !$OMP PARALLEL DO
1021          !$acc kernels loop independent present( qs, qsws, us )
1022          DO  i = nxlg, nxrg
1023             !$acc loop independent
1024             DO  j = nysg, nyng
1025                qsws(j,i) = -qs(j,i) * us(j,i)
1026             ENDDO
1027          ENDDO
1028
1029       ENDIF
1030
1031!
1032!--    Compute (turbulent) fluxes of rain water content and rain drop conc.
1033       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                    &
1034            precipitation )  THEN
1035          !$OMP PARALLEL DO
1036          !$acc kernels loop independent present( nrs, nrsws, qrs, qrsws, us )
1037          DO  i = nxlg, nxrg
1038             !$acc loop independent
1039             DO  j = nysg, nyng
1040                qrsws(j,i) = -qrs(j,i) * us(j,i)
1041                nrsws(j,i) = -nrs(j,i) * us(j,i)
1042             ENDDO
1043          ENDDO
1044       ENDIF
1045
1046!
1047!--    Bottom boundary condition for the TKE
1048       IF ( ibc_e_b == 2 )  THEN
1049          !$OMP PARALLEL DO
1050          !$acc kernels loop independent present( e, nzb_s_inner, us )
1051          DO  i = nxlg, nxrg
1052             !$acc loop independent
1053             DO  j = nysg, nyng
1054                e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2
1055!
1056!--             As a test: cm = 0.4
1057!               e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2
1058                e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
1059             ENDDO
1060          ENDDO
1061       ENDIF
1062
1063    END SUBROUTINE calc_surface_fluxes
1064
1065
1066!
1067!-- Integrated stability function for momentum
1068    FUNCTION psi_m( zeta ) 
1069       
1070       USE kinds
1071
1072       IMPLICIT NONE
1073
1074       REAL(wp)            :: psi_m !< Integrated similarity function result
1075       REAL(wp)            :: zeta  !< Stability parameter z/L
1076       REAL(wp)            :: x     !< dummy variable
1077
1078       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
1079       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
1080       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
1081       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
1082       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
1083       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
1084
1085
1086       IF ( zeta < 0.0_wp )  THEN
1087          x = SQRT( SQRT(1.0_wp  - 16.0_wp * zeta ) )
1088          psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
1089                  * ( 1.0_wp + x**2 ) * 0.125_wp )
1090       ELSE
1091
1092          psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta         &
1093                   - bc_d_d
1094!
1095!--       Old version for stable conditions (only valid for z/L < 0.5)
1096!--       psi_m = - 5.0_wp * zeta
1097
1098       ENDIF
1099
1100    END FUNCTION psi_m
1101
1102
1103!
1104!-- Integrated stability function for heat and moisture
1105    FUNCTION psi_h( zeta ) 
1106       
1107       USE kinds
1108
1109       IMPLICIT NONE
1110
1111       REAL(wp)            :: psi_h !< Integrated similarity function result
1112       REAL(wp)            :: zeta  !< Stability parameter z/L
1113       REAL(wp)            :: x     !< dummy variable
1114
1115       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
1116       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
1117       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
1118       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
1119       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
1120       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
1121
1122
1123       IF ( zeta < 0.0_wp )  THEN
1124          x = SQRT(1.0_wp  - 16.0_wp * zeta )
1125          psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
1126       ELSE
1127          psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp          &
1128                  + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d             &
1129                  + 1.0_wp
1130!
1131!--       Old version for stable conditions (only valid for z/L < 0.5)
1132!--       psi_h = - 5.0_wp * zeta
1133       ENDIF
1134
1135    END FUNCTION psi_h
1136
1137 END MODULE surface_layer_fluxes_mod
Note: See TracBrowser for help on using the repository browser.