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

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

last commit documented

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