source: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 @ 1920

Last change on this file since 1920 was 1920, checked in by suehring, 8 years ago

Avoid segmentation fault (see change in 1915) by different initialization of us instead of adding a very small number in the denominator

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