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

Last change on this file since 2077 was 2077, checked in by raasch, 5 years ago

last commit documented

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