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

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

last commit documented

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