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

Last change on this file since 2101 was 2101, checked in by suehring, 4 years ago

last commit documented

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