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

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

last commit documented

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