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

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

Separate balance equations for humidity and passive_scalar

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