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

Last change on this file since 2037 was 2037, checked in by knoop, 7 years ago

Anelastic approximation implemented

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