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

Last change on this file since 2076 was 2076, checked in by raasch, 7 years ago

further openmp bugfix for lookup method

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