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

Last change on this file since 2001 was 2001, checked in by knoop, 5 years ago

last commit documented

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