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

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

Bugfix: avoid segmentation fault in case of most_method = 'circular' at first timstep

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