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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago

Forced header and separation lines into 80 columns

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