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

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

openacc adjustments and bugfixes

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