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

Last change on this file since 2007 was 2007, checked in by kanani, 8 years ago

changes in the course of urban surface model implementation

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