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

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

last commit documented

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