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

Last change on this file since 1993 was 1993, checked in by suehring, 8 years ago

last commit documented

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