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

Last change on this file since 2011 was 2011, checked in by kanani, 5 years ago

changes related to steering and formating of urban surface model

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