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

Last change on this file since 3157 was 3157, checked in by maronga, 3 years ago

added local free convection velocity scale in calculation of horizontal wind at first grid level

  • Property svn:keywords set to Id
File size: 86.7 KB
Line 
1!> @file surface_layer_fluxes_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
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-2018 Leibniz Universitaet Hannover
18!
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! ------------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: surface_layer_fluxes_mod.f90 3157 2018-07-19 21:08:49Z maronga $
28! Added local free convection velocity scale w_lfc in calculation of uvw_abs.
29! This can be switche on/off by the user via the flag namelist parameter
30! use_free_convection_scaling
31!
32! 3152 2018-07-19 13:26:52Z suehring
33! q_surface is now part of surface structure
34!
35! 3149 2018-07-19 05:48:09Z maronga
36! Correct accidental last commit
37!
38! 3148 2018-07-19 05:45:25Z maronga
39! Major bugfix in calculation of Obukhov length
40!
41! 3146 2018-07-18 22:36:19Z maronga
42! move phi_m from turbulence_closure_mod
43!
44! 3045 2018-05-28 07:55:41Z Giersch
45! Error message revised
46!
47! 2766 2018-01-22 17:17:47Z kanani
48! Removed preprocessor directive __chem
49!
50! 2718 2018-01-02 08:49:38Z maronga
51! Corrected "Former revisions" section
52!
53! 2696 2017-12-14 17:12:51Z kanani
54! - Change in file header (GPL part)
55! - Implementation of chemistry module (FK)
56! - Added calculation of pt1 and qv1 for all surface types. Added calculation of
57!   pt_surface for default-type surfaces (BM)
58! - Add flag to disable computation of qsws in case of urban surface (MS)
59!
60! 2547 2017-10-16 12:41:56Z schwenkel
61! extended by cloud_droplets option
62!
63! 2321 2017-07-24 15:57:07Z schwenkel
64! Bugfix: Correct index in lookup table for Obukhov length
65!
66! 2299 2017-06-29 10:14:38Z suehring
67! Adjusted for allow separate spinups of LSM and atmosphere code
68!
69! 2292 2017-06-20 09:51:42Z schwenkel
70! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
71! includes two more prognostic equations for cloud drop concentration (nc) 
72! and cloud water content (qc).
73!
74! 2281 2017-06-13 11:34:50Z suehring
75! Clean-up unnecessary index access to surface type
76!
77! 2233 2017-05-30 18:08:54Z suehring
78!
79! 2232 2017-05-30 17:47:52Z suehring
80! Adjustments to new surface concept
81! OpenMP bugfix
82!
83! 2118 2017-01-17 16:38:49Z raasch
84! OpenACC directives and related code removed
85!
86! 2091 2016-12-21 16:38:18Z suehring
87! Bugfix in calculation of vsws ( incorrect linear interpolation of us )
88!
89! 2076 2016-12-02 13:54:20Z raasch
90! further openmp bugfix for lookup method
91!
92! 2073 2016-11-30 14:34:05Z raasch
93! openmp bugfix for lookup method
94!
95! 2037 2016-10-26 11:15:40Z knoop
96! Anelastic approximation implemented
97!
98! 2011 2016-09-19 17:29:57Z kanani
99! Flag urban_surface is now defined in module control_parameters.
100!
101! 2007 2016-08-24 15:47:17Z kanani
102! Account for urban surface model in computation of vertical kinematic heatflux
103!
104! 2000 2016-08-20 18:09:15Z knoop
105! Forced header and separation lines into 80 columns
106!
107! 1992 2016-08-12 15:14:59Z suehring
108! Minor bug, declaration of look-up index as INTEGER
109!
110! 1960 2016-07-12 16:34:24Z suehring
111! Treat humidity and passive scalar separately
112!
113! 1929 2016-06-09 16:25:25Z suehring
114! Bugfix: avoid segmentation fault in case one grid point is horizontally
115! completely surrounded by topography
116!
117! 1920 2016-05-30 10:50:15Z suehring
118! Avoid segmentation fault (see change in 1915) by different initialization of
119! us instead of adding a very small number in the denominator
120!
121! 1915 2016-05-27 11:05:02Z suehring
122! Bugfix: avoid segmentation fault in case of most_method = 'circular' at first
123! timestep
124!
125! 1850 2016-04-08 13:29:27Z maronga
126! Module renamed
127!
128!
129! 1822 2016-04-07 07:49:42Z hoffmann
130! icloud_scheme replaced by microphysics_*
131!
132! 1788 2016-03-10 11:01:04Z maronga
133! Added parameter z0q which replaces z0h in the similarity functions for
134! humidity.
135! Syntax layout improved.
136!
137! 1757 2016-02-22 15:49:32Z maronga
138! Minor fixes.
139!
140! 1749 2016-02-09 12:19:56Z raasch
141! further OpenACC adjustments
142!
143! 1747 2016-02-08 12:25:53Z raasch
144! adjustments for OpenACC usage
145!
146! 1709 2015-11-04 14:47:01Z maronga
147! Bugfix: division by zero could occur when calculating rib at low wind speeds
148! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
149! neutral = .T.
150!
151! 1705 2015-11-02 14:28:56Z maronga
152! Typo removed
153!
154! 1697 2015-10-28 17:14:10Z raasch
155! FORTRAN and OpenMP errors removed
156!
157! 1696 2015-10-27 10:03:34Z maronga
158! Modularized and completely re-written version of prandtl_fluxes.f90. In the
159! course of the re-writing two additional methods have been implemented. See
160! updated description.
161!
162! 1551 2015-03-03 14:18:16Z maronga
163! Removed land surface model part. The surface fluxes are now always calculated
164! within prandtl_fluxes, based on the given surface temperature/humidity (which
165! is either provided by the land surface model, by large scale forcing data, or
166! directly prescribed by the user.
167!
168! 1496 2014-12-02 17:25:50Z maronga
169! Adapted for land surface model
170!
171! 1494 2014-11-21 17:14:03Z maronga
172! Bugfixes: qs is now calculated before calculation of Rif. calculation of
173! buoyancy flux in Rif corrected (added missing humidity term), allow use of
174! topography for coupled runs (not tested)
175!
176! 1361 2014-04-16 15:17:48Z hoffmann
177! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
178! drop concentration (nrsws) added
179!
180! 1340 2014-03-25 19:45:13Z kanani
181! REAL constants defined as wp-kind
182!
183! 1320 2014-03-20 08:40:49Z raasch
184! ONLY-attribute added to USE-statements,
185! kind-parameters added to all INTEGER and REAL declaration statements,
186! kinds are defined in new module kinds,
187! old module precision_kind is removed,
188! revision history before 2012 removed,
189! comment fields (!:) to be used for variable explanations added to
190! all variable declaration statements
191!
192! 1276 2014-01-15 13:40:41Z heinze
193! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
194!
195! 1257 2013-11-08 15:18:40Z raasch
196! openACC "kernels do" replaced by "kernels loop", "loop independent" added
197!
198! 1036 2012-10-22 13:43:42Z raasch
199! code put under GPL (PALM 3.9)
200!
201! 1015 2012-09-27 09:23:24Z raasch
202! OpenACC statements added
203!
204! 978 2012-08-09 08:28:32Z fricke
205! roughness length for scalar quantities z0h added
206!
207! Revision 1.1  1998/01/23 10:06:06  raasch
208! Initial revision
209!
210!
211! Description:
212! ------------
213!> Diagnostic computation of vertical fluxes in the constant flux layer from the
214!> values of the variables at grid point k=1. Three different methods are
215!> available:
216!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
217!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
218!>    slower
219!> 3) a method using a lookup table which is fast and accurate. Note, however,
220!>    that this method cannot be used in case of roughness heterogeneity
221!>
222!> @todo (re)move large_scale_forcing actions
223!> @todo check/optimize OpenMP directives
224!> @todo simplify if conditions (which flux need to be computed in which case)
225!------------------------------------------------------------------------------!
226 MODULE surface_layer_fluxes_mod
227
228    USE arrays_3d,                                                             &
229        ONLY:  e, kh, nc, nr, pt, q, ql, qc, qr, s, u, v, vpt, w, zu, zw,      &
230               drho_air_zw, rho_air_zw
231
232    USE chem_modules,                                                          &
233        ONLY:  constant_csflux, nvar
234
235    USE cloud_parameters,                                                      &
236        ONLY:  l_d_cp, pt_d_t
237
238    USE constants,                                                             &
239        ONLY:  pi
240
241    USE cpulog
242
243    USE control_parameters,                                                    &
244        ONLY:  air_chemistry, cloud_droplets, cloud_physics,                   &
245               constant_heatflux, constant_scalarflux,                         &     
246               constant_waterflux, coupling_mode, g, humidity, ibc_e_b,        &
247               ibc_pt_b, initializing_actions, kappa,                          &
248               intermediate_timestep_count, intermediate_timestep_count_max,   &
249               land_surface, large_scale_forcing, lsf_surf,                    &
250               message_string, microphysics_morrison, microphysics_seifert,    &
251               most_method, neutral, passive_scalar, pt_surface, q_surface,    &
252               run_coupled, surface_pressure, simulated_time, terminate_run,   &
253               time_since_reference_point, urban_surface,                      &
254               use_free_convection_scaling, zeta_max, zeta_min
255
256    USE grid_variables,                                                        &
257        ONLY:  dx, dy 
258
259    USE indices,                                                               &
260        ONLY:  nxl, nxr, nys, nyn, nzb
261
262    USE kinds
263
264    USE pegrid
265
266    USE land_surface_model_mod,                                                &
267        ONLY:  aero_resist_kray, skip_time_do_lsm
268
269    USE surface_mod,                                                           &
270        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_type,     &
271                surf_usm_h, surf_usm_v
272       
273
274    IMPLICIT NONE
275
276    INTEGER(iwp) ::  i              !< loop index x direction
277    INTEGER(iwp) ::  j              !< loop index y direction
278    INTEGER(iwp) ::  k              !< loop index z direction
279    INTEGER(iwp) ::  l              !< loop index for surf type
280    INTEGER(iwp) ::  li_bnd  = 7500 !< Lookup table index of the last time step
281
282    INTEGER(iwp), PARAMETER ::  num_steps = 15000  !< number of steps in the lookup table
283
284    LOGICAL      ::  coupled_run       !< Flag for coupled atmosphere-ocean runs
285    LOGICAL      ::  downward = .FALSE.!< Flag indicating downward-facing horizontal surface
286    LOGICAL      ::  mom_uv  = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces
287    LOGICAL      ::  mom_w   = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces
288    LOGICAL      ::  mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production
289    LOGICAL      ::  surf_vertical     !< Flag indicating vertical surfaces
290
291    REAL(wp), DIMENSION(0:num_steps-1) :: rib_tab,  & !< Lookup table bulk Richardson number
292                                          ol_tab      !< Lookup table values of L
293
294    REAL(wp)     ::  e_s,               & !< Saturation water vapor pressure
295                     ol_max = 1.0E6_wp, & !< Maximum Obukhov length
296                     rib_max,           & !< Maximum Richardson number in lookup table
297                     rib_min,           & !< Minimum Richardson number in lookup table
298                     z_mo                 !< Height of the constant flux layer where MOST is assumed
299
300    TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
301
302
303    SAVE
304
305    PRIVATE
306
307    PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes
308
309    INTERFACE init_surface_layer_fluxes
310       MODULE PROCEDURE init_surface_layer_fluxes
311    END INTERFACE init_surface_layer_fluxes
312
313    INTERFACE phi_m
314       MODULE PROCEDURE phi_m
315    END INTERFACE phi_m
316
317    INTERFACE surface_layer_fluxes
318       MODULE PROCEDURE surface_layer_fluxes
319    END INTERFACE surface_layer_fluxes
320
321
322 CONTAINS
323
324
325!------------------------------------------------------------------------------!
326! Description:
327! ------------
328!> Main routine to compute the surface fluxes
329!------------------------------------------------------------------------------!
330    SUBROUTINE surface_layer_fluxes
331
332       IMPLICIT NONE
333
334       surf_vertical = .FALSE.
335       downward      = .FALSE.
336!
337!--    Derive potential temperature and specific humidity at first grid level
338!--    from the fields pt and q
339!
340!--    First call for horizontal default-type surfaces (l=0 - upward facing,
341!--    l=1 - downward facing)
342       DO  l = 0, 1
343          IF ( surf_def_h(l)%ns >= 1  )  THEN
344             surf => surf_def_h(l)
345             CALL calc_pt_q
346             IF ( .NOT. neutral )  THEN
347                CALL calc_pt_surface
348                IF ( humidity )  THEN
349                   CALL calc_q_surface
350                   CALL calc_vpt_surface
351                ENDIF
352             ENDIF
353          ENDIF
354       ENDDO
355!
356!--    Call for natural-type horizontal surfaces
357       IF ( surf_lsm_h%ns >= 1  )  THEN
358          surf => surf_lsm_h
359          CALL calc_pt_q
360       ENDIF
361
362!
363!--    Call for urban-type horizontal surfaces
364       IF ( surf_usm_h%ns >= 1  )  THEN
365          surf => surf_usm_h
366          CALL calc_pt_q
367       ENDIF
368
369!
370!--    Call for natural-type vertical surfaces
371       DO  l = 0, 3
372          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
373             surf => surf_lsm_v(l)
374             CALL calc_pt_q
375          ENDIF
376
377!--       Call for urban-type vertical surfaces
378          IF ( surf_usm_v(l)%ns >= 1  )  THEN
379             surf => surf_usm_v(l)
380             CALL calc_pt_q
381          ENDIF
382       ENDDO
383
384!
385!--    First, calculate the new Obukhov length, then new friction velocity,
386!--    followed by the new scaling parameters (th*, q*, etc.), and the new
387!--    surface fluxes if required. The old routine ("circular") requires a
388!--    different order of calls as the scaling parameters from the previous time
389!--    steps are used to calculate the Obukhov length
390
391!
392!--    Depending on setting of most_method use the "old" routine
393!--    Note, each routine is called for different surface types.
394!--    First call for default-type horizontal surfaces, for natural- and
395!--    urban-type surfaces. Note, at this place only upward-facing horizontal
396!--    surfaces are treted.
397       IF ( most_method == 'circular' )  THEN
398!
399!--       Default-type upward-facing horizontal surfaces
400          IF ( surf_def_h(0)%ns >= 1  )  THEN
401             surf => surf_def_h(0)
402             CALL calc_scaling_parameters
403             CALL calc_uvw_abs
404             IF ( .NOT. neutral )  CALL calc_ol
405             CALL calc_us
406             CALL calc_surface_fluxes
407          ENDIF
408!
409!--       Natural-type horizontal surfaces
410          IF ( surf_lsm_h%ns >= 1  )  THEN
411             surf => surf_lsm_h
412             CALL calc_scaling_parameters
413             CALL calc_uvw_abs
414             IF ( .NOT. neutral )  CALL calc_ol
415             CALL calc_us
416             CALL calc_surface_fluxes
417          ENDIF
418!
419!--       Urban-type horizontal surfaces
420          IF ( surf_usm_h%ns >= 1  )  THEN
421             surf => surf_usm_h
422             CALL calc_scaling_parameters
423             CALL calc_uvw_abs
424             IF ( .NOT. neutral )  CALL calc_ol
425             CALL calc_us
426             CALL calc_surface_fluxes
427          ENDIF
428!
429!--    Use either Newton iteration or a lookup table for the bulk Richardson
430!--    number to calculate the Obukhov length
431       ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
432!
433!--       Default-type upward-facing horizontal surfaces
434          IF ( surf_def_h(0)%ns >= 1  )  THEN
435             surf => surf_def_h(0)
436             CALL calc_uvw_abs
437             IF ( .NOT. neutral )  CALL calc_ol
438             CALL calc_us
439             CALL calc_scaling_parameters
440             CALL calc_surface_fluxes
441          ENDIF
442!
443!--       Natural-type horizontal surfaces
444          IF ( surf_lsm_h%ns >= 1  )  THEN
445             surf => surf_lsm_h
446             CALL calc_uvw_abs
447             IF ( .NOT. neutral )  CALL calc_ol
448             CALL calc_us
449             CALL calc_scaling_parameters
450             CALL calc_surface_fluxes
451          ENDIF
452!
453!--       Urban-type horizontal surfaces
454          IF ( surf_usm_h%ns >= 1  )  THEN
455             surf => surf_usm_h
456             CALL calc_uvw_abs
457             IF ( .NOT. neutral )  CALL calc_ol
458             CALL calc_us
459             CALL calc_scaling_parameters
460             CALL calc_surface_fluxes
461          ENDIF
462
463       ENDIF
464!
465!--    Treat downward-facing horizontal surfaces. Note, so far, these are
466!--    always default type. Stratification is not considered
467!--    in this case, hence, no further distinction between different
468!--    most_method is required. 
469       IF ( surf_def_h(1)%ns >= 1  )  THEN
470          downward = .TRUE.
471          surf => surf_def_h(1)
472          CALL calc_uvw_abs
473          CALL calc_us
474          CALL calc_surface_fluxes
475          downward = .FALSE.
476       ENDIF
477!
478!--    Calculate surfaces fluxes at vertical surfaces for momentum
479!--    and subgrid-scale TKE.
480!--    No stability is considered. Therefore, scaling parameters and Obukhov-
481!--    length do not need to be calculated and no distinction in 'circular',
482!--    'Newton' or 'lookup' is necessary so far.
483!--    Note, this will change if stability is once considered.
484       surf_vertical = .TRUE.
485!
486!--    Calculate horizontal momentum fluxes at north- and south-facing
487!--    surfaces(usvs).
488!--    For default-type surfaces
489       mom_uv = .TRUE. 
490       DO  l = 0, 1
491          IF ( surf_def_v(l)%ns >= 1  )  THEN
492             surf => surf_def_v(l)
493!
494!--          Compute surface-parallel velocity
495             CALL calc_uvw_abs_v_ugrid
496!
497!--          Compute respective friction velocity on staggered grid
498             CALL calc_us
499!
500!--          Compute respective surface fluxes for momentum and TKE
501             CALL calc_surface_fluxes
502          ENDIF
503       ENDDO
504!
505!--    For natural-type surfaces. Please note, even though stability is not
506!--    considered for the calculation of momentum fluxes at vertical surfaces,
507!--    scaling parameters and Obukov length are calculated nevertheless in this
508!--    case. This is due to the requirement of ts in parameterization of heat
509!--    flux in land-surface model in case of aero_resist_kray is not true.
510       IF ( .NOT. aero_resist_kray )  THEN
511          IF ( most_method == 'circular' )  THEN
512             DO  l = 0, 1
513                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
514                   surf => surf_lsm_v(l)
515!
516!--                Compute scaling parameters
517                   CALL calc_scaling_parameters
518!
519!--                Compute surface-parallel velocity
520                   CALL calc_uvw_abs_v_ugrid
521!
522!--                Compute Obukhov length
523                   IF ( .NOT. neutral )  CALL calc_ol
524!
525!--                Compute respective friction velocity on staggered grid
526                   CALL calc_us
527!
528!--                Compute respective surface fluxes for momentum and TKE
529                   CALL calc_surface_fluxes
530                ENDIF
531             ENDDO
532          ELSE
533             DO  l = 0, 1
534                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
535                   surf => surf_lsm_v(l)
536!
537!--                Compute surface-parallel velocity
538                   CALL calc_uvw_abs_v_ugrid
539!
540!--                Compute Obukhov length
541                   IF ( .NOT. neutral )  CALL calc_ol
542!
543!--                Compute respective friction velocity on staggered grid
544                   CALL calc_us
545!
546!--                Compute scaling parameters
547                   CALL calc_scaling_parameters
548!
549!--                Compute respective surface fluxes for momentum and TKE
550                   CALL calc_surface_fluxes
551                ENDIF
552             ENDDO
553          ENDIF
554!
555!--    No ts is required, so scaling parameters and Obukhov length do not need
556!--    to be computed.
557       ELSE
558          DO  l = 0, 1
559             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
560                surf => surf_lsm_v(l)
561!   
562!--             Compute surface-parallel velocity
563                CALL calc_uvw_abs_v_ugrid
564!
565!--             Compute respective friction velocity on staggered grid
566                CALL calc_us
567!
568!--             Compute respective surface fluxes for momentum and TKE
569                CALL calc_surface_fluxes
570             ENDIF
571          ENDDO
572       ENDIF
573!
574!--    For urban-type surfaces
575       DO  l = 0, 1
576          IF ( surf_usm_v(l)%ns >= 1  )  THEN
577             surf => surf_usm_v(l)
578!
579!--          Compute surface-parallel velocity
580             CALL calc_uvw_abs_v_ugrid
581!
582!--          Compute respective friction velocity on staggered grid
583             CALL calc_us
584!
585!--          Compute respective surface fluxes for momentum and TKE
586             CALL calc_surface_fluxes
587          ENDIF
588       ENDDO
589!
590!--    Calculate horizontal momentum fluxes at east- and west-facing
591!--    surfaces (vsus).
592!--    For default-type surfaces
593       DO  l = 2, 3
594          IF ( surf_def_v(l)%ns >= 1  )  THEN
595             surf => surf_def_v(l)
596!
597!--          Compute surface-parallel velocity
598             CALL calc_uvw_abs_v_vgrid
599!
600!--          Compute respective friction velocity on staggered grid
601             CALL calc_us
602!
603!--          Compute respective surface fluxes for momentum and TKE
604             CALL calc_surface_fluxes
605          ENDIF
606       ENDDO
607!
608!--    For natural-type surfaces. Please note, even though stability is not
609!--    considered for the calculation of momentum fluxes at vertical surfaces,
610!--    scaling parameters and Obukov length are calculated nevertheless in this
611!--    case. This is due to the requirement of ts in parameterization of heat
612!--    flux in land-surface model in case of aero_resist_kray is not true.
613       IF ( .NOT. aero_resist_kray )  THEN
614          IF ( most_method == 'circular' )  THEN
615             DO  l = 2, 3
616                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
617                   surf => surf_lsm_v(l)
618!
619!--                Compute scaling parameters
620                   CALL calc_scaling_parameters
621!
622!--                Compute surface-parallel velocity
623                   CALL calc_uvw_abs_v_vgrid
624!
625!--                Compute Obukhov length
626                   IF ( .NOT. neutral )  CALL calc_ol
627!
628!--                Compute respective friction velocity on staggered grid
629                   CALL calc_us
630!
631!--                Compute respective surface fluxes for momentum and TKE
632                   CALL calc_surface_fluxes
633                ENDIF
634             ENDDO
635          ELSE
636             DO  l = 2, 3
637                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
638                   surf => surf_lsm_v(l)
639!
640!--                Compute surface-parallel velocity
641                   CALL calc_uvw_abs_v_vgrid
642!
643!--                Compute Obukhov length
644                   IF ( .NOT. neutral )  CALL calc_ol
645!
646!--                Compute respective friction velocity on staggered grid
647                   CALL calc_us
648!
649!--                Compute scaling parameters
650                   CALL calc_scaling_parameters
651!
652!--                Compute respective surface fluxes for momentum and TKE
653                   CALL calc_surface_fluxes
654                ENDIF
655             ENDDO
656          ENDIF
657       ELSE
658          DO  l = 2, 3
659             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
660                surf => surf_lsm_v(l)
661!   
662!--             Compute surface-parallel velocity
663                CALL calc_uvw_abs_v_vgrid
664!
665!--             Compute respective friction velocity on staggered grid
666                CALL calc_us
667!
668!--             Compute respective surface fluxes for momentum and TKE
669                CALL calc_surface_fluxes
670             ENDIF
671          ENDDO
672       ENDIF
673!
674!--    For urban-type surfaces
675       DO  l = 2, 3
676          IF ( surf_usm_v(l)%ns >= 1  )  THEN
677             surf => surf_usm_v(l)
678!
679!--          Compute surface-parallel velocity
680             CALL calc_uvw_abs_v_vgrid
681!
682!--          Compute respective friction velocity on staggered grid
683             CALL calc_us
684!
685!--          Compute respective surface fluxes for momentum and TKE
686             CALL calc_surface_fluxes
687          ENDIF
688       ENDDO
689       mom_uv = .FALSE.
690!
691!--    Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial
692!--    surfaces.
693       mom_w = .TRUE.
694!
695!--    Default-type surfaces
696       DO  l = 0, 3
697          IF ( surf_def_v(l)%ns >= 1  )  THEN
698             surf => surf_def_v(l)
699             CALL calc_uvw_abs_v_wgrid
700             CALL calc_us
701             CALL calc_surface_fluxes
702          ENDIF
703       ENDDO
704!
705!--    Natural-type surfaces
706       DO  l = 0, 3
707          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
708             surf => surf_lsm_v(l)
709             CALL calc_uvw_abs_v_wgrid
710             CALL calc_us
711             CALL calc_surface_fluxes
712          ENDIF
713       ENDDO
714!
715!--    Urban-type surfaces
716       DO  l = 0, 3
717          IF ( surf_usm_v(l)%ns >= 1  )  THEN
718             surf => surf_usm_v(l)
719             CALL calc_uvw_abs_v_wgrid
720             CALL calc_us
721             CALL calc_surface_fluxes
722          ENDIF
723       ENDDO
724       mom_w = .FALSE.   
725!
726!--    Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical
727!--    surfaces for TKE production. Note, here, momentum fluxes are defined
728!--    at grid center and are not staggered as before.
729       mom_tke = .TRUE.
730!
731!--    Default-type surfaces
732       DO  l = 0, 3
733          IF ( surf_def_v(l)%ns >= 1  )  THEN
734             surf => surf_def_v(l)
735             CALL calc_uvw_abs_v_sgrid
736             CALL calc_us
737             CALL calc_surface_fluxes
738          ENDIF
739       ENDDO
740!
741!--    Natural-type surfaces
742       DO  l = 0, 3
743          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
744             surf => surf_lsm_v(l)
745             CALL calc_uvw_abs_v_sgrid
746             CALL calc_us
747             CALL calc_surface_fluxes
748          ENDIF
749       ENDDO
750!
751!--    Urban-type surfaces
752       DO  l = 0, 3
753          IF ( surf_usm_v(l)%ns >= 1  )  THEN
754             surf => surf_usm_v(l)
755             CALL calc_uvw_abs_v_sgrid
756             CALL calc_us
757             CALL calc_surface_fluxes
758          ENDIF
759       ENDDO
760       mom_tke = .FALSE.
761 
762
763    END SUBROUTINE surface_layer_fluxes
764
765
766!------------------------------------------------------------------------------!
767! Description:
768! ------------
769!> Initializing actions for the surface layer routine. Basically, this involves
770!> the preparation of a lookup table for the the bulk Richardson number vs
771!> Obukhov length L when using the lookup table method.
772!------------------------------------------------------------------------------!
773    SUBROUTINE init_surface_layer_fluxes
774
775       IMPLICIT NONE
776
777       INTEGER(iwp) :: li,         & !< Index for loop to create lookup table
778                       num_steps_n   !< Number of non-stretched zeta steps
779
780       LOGICAL :: terminate_run_l = .FALSE.    !< Flag to terminate run (global)
781
782       REAL(wp), PARAMETER ::  zeta_stretch = -10.0_wp !< Start of stretching in the free convection limit
783                               
784       REAL(wp), DIMENSION(:), ALLOCATABLE :: zeta_tmp
785
786
787       REAL(wp) :: zeta_step,            & !< Increment of zeta
788                   regr      = 1.01_wp,  & !< Stretching factor of zeta_step in the free convection limit
789                   regr_old  = 1.0E9_wp, & !< Stretching factor of last iteration step
790                   z0h_min   = 0.0_wp,   & !< Minimum value of z0h to create table
791                   z0_min    = 0.0_wp      !< Minimum value of z0 to create table
792
793
794
795
796!
797!--    In case of runs with neutral statification, set Obukhov length to a
798!--    large value
799       IF ( neutral )  THEN
800          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%ol = 1.0E10_wp
801          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
802          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
803       ENDIF
804
805       IF ( most_method == 'lookup' )  THEN
806
807!
808!--       Check for roughness heterogeneity. In that case terminate run and
809!--       inform user. Check for both, natural and non-natural walls.
810          IF ( surf_def_h(0)%ns >= 1 )  THEN
811             IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR. &
812                  MINVAL( surf_def_h(0)%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
813                terminate_run_l = .TRUE.
814             ENDIF
815          ENDIF
816          IF ( surf_lsm_h%ns >= 1 )  THEN
817             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
818                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
819                terminate_run_l = .TRUE.
820             ENDIF
821          ENDIF
822          IF ( surf_usm_h%ns >= 1 )  THEN
823             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h )  .OR.       &
824                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_usm_h%z0  ) )  THEN
825                terminate_run_l = .TRUE.
826             ENDIF
827          ENDIF
828!
829!--       Check roughness homogeneity between differt surface types.
830          IF ( surf_lsm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
831             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
832                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
833                terminate_run_l = .TRUE.
834             ENDIF
835          ENDIF
836          IF ( surf_usm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
837             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
838                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
839                terminate_run_l = .TRUE.
840             ENDIF
841          ENDIF
842          IF ( surf_usm_h%ns >= 1  .AND.  surf_lsm_h%ns >= 1 )  THEN
843             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
844                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
845                terminate_run_l = .TRUE.
846             ENDIF
847          ENDIF
848
849
850#if defined( __parallel )
851!
852!--       Make a logical OR for all processes. Force termiation of model if result
853!--       is TRUE
854          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
855          CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,  &
856                              MPI_LOR, comm2d, ierr )
857#else
858          terminate_run = terminate_run_l
859#endif
860
861          IF ( terminate_run )  THEN
862             message_string = 'most_method = "lookup" cannot be used in ' //   &
863                              'combination with a prescribed roughness '  //   &
864                              'heterogeneity'
865             CALL message( 'surface_layer_fluxes', 'PA0116', 1, 2, 0, 6, 0 )
866          ENDIF
867
868          ALLOCATE(  zeta_tmp(0:num_steps-1) )
869
870!
871!--       Use the lowest possible value for z_mo
872          k    = nzb
873          z_mo = zu(k+1) - zw(k)
874
875!
876!--       Calculate z/L range from zeta_stretch to zeta_max using 90% of the
877!--       available steps (num_steps). The calculation is done with negative
878!--       values of zeta in order to simplify the stretching in the free
879!--       convection limit for the remaining 10% of steps.
880          zeta_tmp(0) = - zeta_max
881          num_steps_n = ( num_steps * 9 / 10 ) - 1
882          zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
883
884          DO li = 1, num_steps_n
885             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
886          ENDDO
887
888!
889!--       Calculate stretching factor for the free convection range
890          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
891             regr_old = regr
892             regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
893                    )**( 10.0_wp / REAL(num_steps) )
894          ENDDO
895
896!
897!--       Calculate z/L range from zeta_min to zeta_stretch
898          DO li = num_steps_n+1, num_steps-1
899             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
900             zeta_step = zeta_step * regr
901          ENDDO
902
903!
904!--       Save roughness lengths to temporary variables
905          IF ( surf_def_h(0)%ns >= 1  )  THEN
906             z0h_min = surf_def_h(0)%z0h(1)
907             z0_min  = surf_def_h(0)%z0(1)
908          ELSEIF ( surf_lsm_h%ns >= 1 )  THEN
909             z0h_min = surf_lsm_h%z0h(1)
910             z0_min  = surf_lsm_h%z0(1)
911          ELSEIF ( surf_usm_h%ns >= 1 )  THEN
912             z0h_min = surf_usm_h%z0h(1)
913             z0_min  = surf_usm_h%z0(1)
914          ENDIF         
915!
916!--       Calculate lookup table for the Richardson number versus Obukhov length
917!--       The Richardson number (rib) is defined depending on the choice of
918!--       boundary conditions for temperature
919          IF ( ibc_pt_b == 1 )  THEN
920             DO li = 0, num_steps-1
921                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
922                rib_tab(li) = z_mo / ol_tab(li)  / ( LOG( z_mo / z0_min )       &
923                                                - psi_m( z_mo / ol_tab(li) )    &
924                                                + psi_m( z0_min / ol_tab(li) )  &
925                                                  )**3
926             ENDDO 
927          ELSE
928             DO li = 0, num_steps-1
929                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
930                rib_tab(li) = z_mo / ol_tab(li)  * ( LOG( z_mo / z0h_min )     &
931                                              - psi_h( z_mo / ol_tab(li) )     &
932                                              + psi_h( z0h_min / ol_tab(li) )  &
933                                            )                                  &
934                                          / ( LOG( z_mo / z0_min )             &
935                                              - psi_m( z_mo / ol_tab(li) )     &
936                                              + psi_m( z0_min / ol_tab(li) )   &
937                                            )**2
938             ENDDO
939          ENDIF
940
941!
942!--       Determine minimum values of rib in the lookup table. Set upper limit
943!--       to critical Richardson number (0.25)
944          rib_min  = MINVAL(rib_tab)
945          rib_max  = 0.25 !MAXVAL(rib_tab)
946
947          DEALLOCATE( zeta_tmp )
948       ENDIF
949
950    END SUBROUTINE init_surface_layer_fluxes
951
952
953!------------------------------------------------------------------------------!
954! Description:
955! ------------
956!> Compute the absolute value of the horizontal velocity (relative to the
957!> surface) for horizontal surface elements. This is required by all methods.
958!------------------------------------------------------------------------------!
959    SUBROUTINE calc_uvw_abs
960   
961       IMPLICIT NONE
962
963       INTEGER(iwp) ::  i             !< running index x direction
964       INTEGER(iwp) ::  ibit          !< flag to mask computation of relative velocity in case of downward-facing surfaces
965       INTEGER(iwp) ::  j             !< running index y direction
966       INTEGER(iwp) ::  k             !< running index z direction
967       INTEGER(iwp) ::  m             !< running index surface elements
968
969       REAL(wp)     :: w_lfc          !< local free convection velocity scale
970!
971!--    ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces.
972       ibit = MERGE( 1, 0, .NOT. downward )
973
974       DO  m = 1, surf%ns
975
976          i   = surf%i(m)           
977          j   = surf%j(m)
978          k   = surf%k(m)
979       
980!
981!--       Calculate free convection velocity scale w_lfc is
982!--       use_free_convection_scaling = .T.. This will maintain a horizontal
983!--       velocity even for very weak wind convective conditions. SIGN is used
984!--       to set w_lfc to zero under stable conditions.
985          IF ( use_free_convection_scaling )  THEN
986             w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m))
987             w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp)
988          ELSE
989             w_lfc = 0.0_wp
990          ENDIF
991
992!
993!--       Compute the absolute value of the horizontal velocity.
994!--       (relative to the surface in case the lower surface is the ocean).
995!--       Please note, in new surface modelling concept the index values changed,
996!--       i.e. the reference grid point is not the surface-grid point itself but
997!--       the first grid point outside of the topography.
998!--       Note, in case of coupled ocean-atmosphere simulations relative velocity
999!--       with respect to the ocean surface is used, hence, (k-1,j,i) values
1000!--       are used to calculate the absolute velocity.
1001!--       However, this do not apply for downward-facing walls. To mask this,
1002!--       use ibit, which checks for upward/downward-facing surfaces.
1003          surf%uvw_abs(m) = SQRT(                                              &
1004                              ( 0.5_wp * (   u(k,j,i)   + u(k,j,i+1)           &
1005                                        -  ( u(k-1,j,i) + u(k-1,j,i+1)         &
1006                                           ) * ibit                            &
1007                                         )                                     &
1008                              )**2 +                                           &
1009                              ( 0.5_wp * (   v(k,j,i)   + v(k,j+1,i)           &
1010                                        -  ( v(k-1,j,i) + v(k-1,j+1,i)         &
1011                                           ) * ibit                            &
1012                                         )                                     &
1013                              )**2  + w_lfc**2                                 &
1014                                )
1015
1016                               
1017
1018       ENDDO
1019
1020    END SUBROUTINE calc_uvw_abs
1021
1022
1023!------------------------------------------------------------------------------!
1024! Description:
1025! ------------
1026!> Compute the absolute value of the horizontal velocity (relative to the
1027!> surface) for horizontal surface elements. This is required by all methods.
1028!------------------------------------------------------------------------------!
1029    SUBROUTINE calc_uvw_abs_v_ugrid
1030
1031       IMPLICIT NONE
1032
1033       INTEGER(iwp) ::  i             !< running index x direction
1034       INTEGER(iwp) ::  j             !< running index y direction
1035       INTEGER(iwp) ::  k             !< running index z direction
1036       INTEGER(iwp) ::  m             !< running index surface elements
1037
1038       REAL(wp)     ::  u_i
1039       REAL(wp)     ::  w_i
1040
1041
1042       DO  m = 1, surf%ns
1043          i   = surf%i(m)           
1044          j   = surf%j(m)
1045          k   = surf%k(m)
1046!
1047!--       Compute the absolute value of the surface parallel velocity on u-grid.
1048          u_i  = u(k,j,i)
1049          w_i  = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) +                       &
1050                             w(k,j,i-1)   + w(k,j,i) ) 
1051
1052          surf%uvw_abs(m) = SQRT( u_i**2 + w_i**2 ) 
1053
1054       ENDDO
1055
1056    END SUBROUTINE calc_uvw_abs_v_ugrid
1057
1058!------------------------------------------------------------------------------!
1059! Description:
1060! ------------
1061!> Compute the absolute value of the horizontal velocity (relative to the
1062!> surface) for horizontal surface elements. This is required by all methods.
1063!------------------------------------------------------------------------------!
1064    SUBROUTINE calc_uvw_abs_v_vgrid
1065
1066       IMPLICIT NONE
1067
1068       INTEGER(iwp) ::  i             !< running index x direction
1069       INTEGER(iwp) ::  j             !< running index y direction
1070       INTEGER(iwp) ::  k             !< running index z direction
1071       INTEGER(iwp) ::  m             !< running index surface elements
1072
1073       REAL(wp)     ::  v_i
1074       REAL(wp)     ::  w_i
1075
1076
1077       DO  m = 1, surf%ns
1078          i   = surf%i(m)           
1079          j   = surf%j(m)
1080          k   = surf%k(m)
1081
1082          v_i  = u(k,j,i)
1083          w_i  = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) +                       &
1084                             w(k,j-1,i)   + w(k,j,i) ) 
1085
1086          surf%uvw_abs(m) = SQRT( v_i**2 + w_i**2 ) 
1087
1088       ENDDO
1089
1090    END SUBROUTINE calc_uvw_abs_v_vgrid
1091
1092!------------------------------------------------------------------------------!
1093! Description:
1094! ------------
1095!> Compute the absolute value of the horizontal velocity (relative to the
1096!> surface) for horizontal surface elements. This is required by all methods.
1097!------------------------------------------------------------------------------!
1098    SUBROUTINE calc_uvw_abs_v_wgrid
1099
1100       IMPLICIT NONE
1101
1102       INTEGER(iwp) ::  i             !< running index x direction
1103       INTEGER(iwp) ::  j             !< running index y direction
1104       INTEGER(iwp) ::  k             !< running index z direction
1105       INTEGER(iwp) ::  m             !< running index surface elements
1106
1107       REAL(wp)     ::  u_i
1108       REAL(wp)     ::  v_i
1109       REAL(wp)     ::  w_i
1110!
1111!--    North- (l=0) and south-facing (l=1) surfaces
1112       IF ( l == 0  .OR.  l == 1 )  THEN
1113          DO  m = 1, surf%ns
1114             i   = surf%i(m)           
1115             j   = surf%j(m)
1116             k   = surf%k(m)
1117
1118             u_i  = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) +                    &
1119                                u(k,j,i+1)   + u(k,j,i) )
1120             v_i  = 0.0_wp
1121             w_i  = w(k,j,i)
1122
1123             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1124          ENDDO
1125!
1126!--    East- (l=2) and west-facing (l=3) surfaces
1127       ELSE
1128          DO  m = 1, surf%ns
1129             i   = surf%i(m)           
1130             j   = surf%j(m)
1131             k   = surf%k(m)
1132
1133             u_i  = 0.0_wp
1134             v_i  = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) +                    &
1135                                v(k,j+1,i)   + v(k,j,i) )
1136             w_i  = w(k,j,i)
1137
1138             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1139          ENDDO
1140       ENDIF           
1141
1142    END SUBROUTINE calc_uvw_abs_v_wgrid
1143
1144!------------------------------------------------------------------------------!
1145! Description:
1146! ------------
1147!> Compute the absolute value of the horizontal velocity (relative to the
1148!> surface) for horizontal surface elements. This is required by all methods.
1149!------------------------------------------------------------------------------!
1150    SUBROUTINE calc_uvw_abs_v_sgrid
1151
1152       IMPLICIT NONE
1153
1154       INTEGER(iwp) ::  i             !< running index x direction
1155       INTEGER(iwp) ::  j             !< running index y direction
1156       INTEGER(iwp) ::  k             !< running index z direction
1157       INTEGER(iwp) ::  m             !< running index surface elements
1158
1159       REAL(wp)     ::  u_i
1160       REAL(wp)     ::  v_i
1161       REAL(wp)     ::  w_i
1162
1163!
1164!--    North- (l=0) and south-facing (l=1) walls
1165       IF ( l == 0  .OR.  l == 1 )  THEN
1166          DO  m = 1, surf%ns
1167             i   = surf%i(m)           
1168             j   = surf%j(m)
1169             k   = surf%k(m)
1170
1171             u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
1172             v_i = 0.0_wp
1173             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1174
1175             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1176          ENDDO
1177!
1178!--    East- (l=2) and west-facing (l=3) walls
1179       ELSE
1180          DO  m = 1, surf%ns
1181             i   = surf%i(m)           
1182             j   = surf%j(m)
1183             k   = surf%k(m)
1184
1185             u_i = 0.0_wp
1186             v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
1187             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1188
1189             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1190          ENDDO
1191       ENDIF 
1192
1193    END SUBROUTINE calc_uvw_abs_v_sgrid
1194
1195
1196!------------------------------------------------------------------------------!
1197! Description:
1198! ------------
1199!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
1200!------------------------------------------------------------------------------!
1201    SUBROUTINE calc_ol
1202
1203       IMPLICIT NONE
1204
1205       INTEGER(iwp) ::  iter    !< Newton iteration step
1206       INTEGER(iwp) ::  li      !< look index
1207       INTEGER(iwp) ::  m       !< loop variable over all horizontal wall elements
1208
1209       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
1210                       f_d_ol, & !< Derivative of f
1211                       ol_l,   & !< Lower bound of L for Newton iteration
1212                       ol_m,   & !< Previous value of L for Newton iteration
1213                       ol_old, & !< Previous time step value of L
1214                       ol_u      !< Upper bound of L for Newton iteration
1215
1216       IF ( TRIM( most_method ) /= 'circular' )  THEN
1217!
1218!--       Evaluate bulk Richardson number (calculation depends on
1219!--       definition based on setting of boundary conditions
1220          IF ( ibc_pt_b /= 1 )  THEN
1221             IF ( humidity )  THEN
1222                !$OMP PARALLEL DO PRIVATE( z_mo )
1223                DO  m = 1, surf%ns
1224
1225                   z_mo = surf%z_mo(m)
1226
1227                   surf%rib(m) = g * z_mo                                      &
1228                                 * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
1229                                 / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
1230                                 + 1.0E-20_wp )
1231                ENDDO
1232             ELSE
1233                !$OMP PARALLEL DO PRIVATE( z_mo )
1234                DO  m = 1, surf%ns
1235
1236                   z_mo = surf%z_mo(m)
1237
1238                   surf%rib(m) = g * z_mo                                      &
1239                                 * ( surf%pt1(m) - surf%pt_surface(m) )        &
1240                                 / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
1241                                 + 1.0E-20_wp )
1242                ENDDO
1243             ENDIF
1244          ELSE
1245             IF ( humidity )  THEN
1246                !$OMP PARALLEL DO PRIVATE( k, z_mo )
1247                DO  m = 1, surf%ns
1248
1249                   k   = surf%k(m)
1250
1251                   z_mo = surf%z_mo(m)
1252
1253                   surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
1254                           * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
1255                           * surf%pt1(m) * surf%qsws(m) ) *                    &
1256                             drho_air_zw(k-1)                /                 &
1257                         ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
1258                           + 1.0E-20_wp )
1259                ENDDO
1260             ELSE
1261                !$OMP PARALLEL DO PRIVATE( k, z_mo )
1262                DO  m = 1, surf%ns
1263
1264                   k   = surf%k(m)
1265
1266                   z_mo = surf%z_mo(m)
1267
1268                   surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
1269                                        drho_air_zw(k-1)            /          &
1270                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
1271                           + 1.0E-20_wp )
1272                ENDDO
1273             ENDIF
1274          ENDIF
1275
1276       ENDIF
1277
1278
1279!
1280!--    Calculate the Obukhov length either using a Newton iteration
1281!--    method, via a lookup table, or using the old circular way
1282       IF ( TRIM( most_method ) == 'newton' )  THEN
1283
1284          DO  m = 1, surf%ns
1285
1286             i   = surf%i(m)           
1287             j   = surf%j(m)
1288
1289             z_mo = surf%z_mo(m)
1290
1291!
1292!--          Store current value in case the Newton iteration fails
1293             ol_old = surf%ol(m)
1294
1295!
1296!--          Ensure that the bulk Richardson number and the Obukhov
1297!--          length have the same sign
1298             IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
1299                  ABS( surf%ol(m) ) == ol_max )  THEN
1300                IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
1301                IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
1302             ENDIF
1303!
1304!--          Iteration to find Obukhov length
1305             iter = 0
1306             DO
1307                iter = iter + 1
1308!
1309!--             In case of divergence, use the value of the previous time step
1310                IF ( iter > 1000 )  THEN
1311                   surf%ol(m) = ol_old
1312                   EXIT
1313                ENDIF
1314
1315                ol_m = surf%ol(m)
1316                ol_l = ol_m - 0.001_wp * ol_m
1317                ol_u = ol_m + 0.001_wp * ol_m
1318
1319
1320                IF ( ibc_pt_b /= 1 )  THEN
1321!
1322!--                Calculate f = Ri - [...]/[...]^2 = 0
1323                   f = surf%rib(m) - ( z_mo / ol_m ) * (                       &
1324                                                 LOG( z_mo / surf%z0h(m) )     &
1325                                                 - psi_h( z_mo / ol_m )        &
1326                                                 + psi_h( surf%z0h(m) /        &
1327                                                          ol_m )               &
1328                                                               )               &
1329                                              / ( LOG( z_mo / surf%z0(m) )     &
1330                                                 - psi_m( z_mo / ol_m )        &
1331                                                 + psi_m( surf%z0(m) /         &
1332                                                          ol_m )               &
1333                                                 )**2
1334
1335!
1336!--                 Calculate df/dL
1337                    f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /               &
1338                                                             surf%z0h(m) )     &
1339                                            - psi_h( z_mo / ol_u )             &
1340                                            + psi_h( surf%z0h(m) / ol_u )      &
1341                                              )                                &
1342                                            / ( LOG( z_mo / surf%z0(m) )       &
1343                                            - psi_m( z_mo / ol_u )             &
1344                                            + psi_m( surf%z0(m) / ol_u )       &
1345                                              )**2                             &
1346                           + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )     &
1347                                            - psi_h( z_mo / ol_l )             &
1348                                            + psi_h( surf%z0h(m) / ol_l )      &
1349                                              )                                &
1350                                            / ( LOG( z_mo / surf%z0(m) )       &
1351                                            - psi_m( z_mo / ol_l )             &
1352                                            + psi_m( surf%z0(m) / ol_l )       &
1353                                              )**2                             &
1354                             ) / ( ol_u - ol_l )
1355                ELSE
1356!
1357!--                Calculate f = Ri - 1 /[...]^3 = 0
1358                   f = surf%rib(m) - ( z_mo / ol_m ) /                         &
1359                                                ( LOG( z_mo / surf%z0(m) )     &
1360                                            - psi_m( z_mo / ol_m )             &
1361                                            + psi_m( surf%z0(m) / ol_m )       &
1362                                                )**3
1363
1364!
1365!--                Calculate df/dL
1366                   f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo /                &
1367                                                            surf%z0(m) )       &
1368                                            - psi_m( z_mo / ol_u )             &
1369                                            + psi_m( surf%z0(m) / ol_u )       &
1370                                                     )**3                      &
1371                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
1372                                            - psi_m( z_mo / ol_l )             &
1373                                            + psi_m( surf%z0(m) / ol_l )       &
1374                                               )**3                            &
1375                                  ) / ( ol_u - ol_l )
1376                ENDIF
1377!
1378!--             Calculate new L
1379                surf%ol(m) = ol_m - f / f_d_ol
1380
1381!
1382!--             Ensure that the bulk Richardson number and the Obukhov
1383!--             length have the same sign and ensure convergence.
1384                IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
1385!
1386!--             If unrealistic value occurs, set L to the maximum
1387!--             value that is allowed
1388                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
1389                   surf%ol(m) = ol_max
1390                   EXIT
1391                ENDIF
1392!
1393!--             Check for convergence
1394                IF ( ABS( ( surf%ol(m) - ol_m ) /                              &
1395                     surf%ol(m) ) < 1.0E-4_wp )  THEN
1396                   EXIT
1397                ELSE
1398                   CYCLE
1399                ENDIF
1400
1401             ENDDO
1402          ENDDO
1403
1404       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
1405
1406          !$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
1407          DO  m = 1, surf%ns
1408
1409             i   = surf%i(m)           
1410             j   = surf%j(m)
1411!
1412!--          If the bulk Richardson number is outside the range of the lookup
1413!--          table, set it to the exceeding threshold value
1414             IF ( surf%rib(m) < rib_min )  surf%rib(m) = rib_min
1415             IF ( surf%rib(m) > rib_max )  surf%rib(m) = rib_max
1416
1417!
1418!--          Find the correct index bounds for linear interpolation. As the
1419!--          Richardson number will not differ very much from time step to
1420!--          time step , use the index from the last step and search in the
1421!--          correct direction
1422             li = li_bnd
1423             IF ( rib_tab(li) - surf%rib(m) > 0.0_wp )  THEN
1424                DO  WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp  .AND.  li > 0 )
1425                   li = li-1
1426                ENDDO
1427             ELSE
1428                DO  WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp                 &
1429                           .AND.  li < num_steps-1 )
1430                   li = li+1
1431                ENDDO
1432             ENDIF
1433             li_bnd = li
1434
1435!
1436!--          Linear interpolation to find the correct value of z/L
1437             surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) )       &
1438                         / (  rib_tab(li) - rib_tab(li-1) )                    &
1439                         * ( surf%rib(m) - rib_tab(li-1) ) )
1440
1441          ENDDO
1442
1443       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
1444
1445          IF ( .NOT. humidity )  THEN
1446             !$OMP PARALLEL DO PRIVATE( z_mo )
1447             DO  m = 1, surf%ns
1448
1449                z_mo = surf%z_mo(m)
1450
1451                surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
1452                                       ( kappa * g *                           &
1453                                         surf%ts(m) + 1E-30_wp )
1454!
1455!--             Limit the value range of the Obukhov length.
1456!--             This is necessary for very small velocities (u,v --> 1), because
1457!--             the absolute value of ol can then become very small, which in
1458!--             consequence would result in very large shear stresses and very
1459!--             small momentum fluxes (both are generally unrealistic).
1460                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1461                   surf%ol(m) = z_mo / zeta_min
1462                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1463                   surf%ol(m) = z_mo / zeta_max
1464
1465             ENDDO
1466          ELSEIF ( cloud_physics  .OR.  cloud_droplets )  THEN
1467             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1468             DO  m = 1, surf%ns
1469
1470                i   = surf%i(m)           
1471                j   = surf%j(m)
1472                k   = surf%k(m)
1473
1474                z_mo = surf%z_mo(m)
1475
1476
1477                surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
1478                    ( kappa * g * ( surf%ts(m) +                               &
1479                         0.61_wp * surf%pt1(m) * surf%us(m)                    &
1480                       + 0.61_wp * surf%qv1(m) * surf%ts(m) -                  &
1481                                   surf%ts(m)  * ql(k,j,i) ) + 1E-30_wp )
1482!
1483!--             Limit the value range of the Obukhov length.
1484!--             This is necessary for very small velocities (u,v --> 1), because
1485!--             the absolute value of ol can then become very small, which in
1486!--             consequence would result in very large shear stresses and very
1487!--             small momentum fluxes (both are generally unrealistic).
1488                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1489                   surf%ol(m) = z_mo / zeta_min
1490                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1491                   surf%ol(m) = z_mo / zeta_max
1492
1493             ENDDO
1494          ELSE
1495
1496             !$OMP PARALLEL DO PRIVATE( z_mo )
1497             DO  m = 1, surf%ns
1498
1499                z_mo = surf%z_mo(m)
1500
1501                surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
1502                    ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
1503                                    surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
1504                                    surf%ts(m) ) + 1E-30_wp )
1505
1506!
1507!--             Limit the value range of the Obukhov length.
1508!--             This is necessary for very small velocities (u,v --> 1), because
1509!--             the absolute value of ol can then become very small, which in
1510!--             consequence would result in very large shear stresses and very
1511!--             small momentum fluxes (both are generally unrealistic).
1512                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1513                   surf%ol(m) = z_mo / zeta_min
1514                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1515                   surf%ol(m) = z_mo / zeta_max
1516
1517             ENDDO
1518
1519          ENDIF
1520
1521       ENDIF
1522
1523    END SUBROUTINE calc_ol
1524
1525!
1526!-- Calculate friction velocity u*
1527    SUBROUTINE calc_us
1528
1529       IMPLICIT NONE
1530
1531       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1532
1533!
1534!--    Compute u* at horizontal surfaces at the scalars' grid points
1535       IF ( .NOT. surf_vertical )  THEN
1536!
1537!--       Compute u* at upward-facing surfaces
1538          IF ( .NOT. downward )  THEN
1539             !$OMP PARALLEL  DO PRIVATE( z_mo )
1540             DO  m = 1, surf%ns
1541
1542                z_mo = surf%z_mo(m)
1543!
1544!--             Compute u* at the scalars' grid points
1545                surf%us(m) = kappa * surf%uvw_abs(m) /                         &
1546                            ( LOG( z_mo / surf%z0(m) )                         &
1547                           - psi_m( z_mo / surf%ol(m) )                        &
1548                           + psi_m( surf%z0(m) / surf%ol(m) ) )
1549   
1550             ENDDO
1551!
1552!--       Compute u* at downward-facing surfaces. This case, do not consider
1553!--       any stability.
1554          ELSE
1555             !$OMP PARALLEL  DO PRIVATE( z_mo )
1556             DO  m = 1, surf%ns
1557
1558                z_mo = surf%z_mo(m)
1559!
1560!--             Compute u* at the scalars' grid points
1561                surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1562   
1563             ENDDO
1564          ENDIF
1565!
1566!--    Compute u* at vertical surfaces at the u/v/v grid, respectively.
1567!--    No stability is considered in this case.
1568       ELSE
1569          !$OMP PARALLEL DO PRIVATE( z_mo )
1570          DO  m = 1, surf%ns
1571             z_mo = surf%z_mo(m)
1572
1573             surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1574          ENDDO
1575       ENDIF
1576
1577    END SUBROUTINE calc_us
1578
1579!
1580!-- Calculate potential temperature, specific humidity, and virtual potential
1581!-- temperature at first grid level
1582!-- ( only for upward-facing surfs )
1583    SUBROUTINE calc_pt_q
1584
1585       IMPLICIT NONE
1586
1587       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1588
1589       !$OMP PARALLEL DO PRIVATE( i, j, k )
1590       DO  m = 1, surf%ns
1591
1592          i   = surf%i(m)           
1593          j   = surf%j(m)
1594          k   = surf%k(m)
1595
1596          IF ( cloud_physics ) THEN
1597             surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
1598             surf%qv1(m) = q(k,j,i) - ql(k,j,i)
1599          ELSEIF( cloud_droplets ) THEN
1600             surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
1601             surf%qv1(m) = q(k,j,i) 
1602          ELSE
1603             surf%pt1(m) = pt(k,j,i)
1604             IF ( humidity )  THEN
1605                surf%qv1(m) = q(k,j,i)
1606             ELSE
1607                surf%qv1(m) = 0.0_wp
1608             ENDIF
1609          ENDIF
1610
1611          IF ( humidity )  THEN
1612             surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
1613          ENDIF
1614         
1615       ENDDO
1616
1617    END SUBROUTINE calc_pt_q
1618
1619
1620!
1621!-- Set potential temperature at surface grid level.
1622!-- ( only for upward-facing surfs )
1623    SUBROUTINE calc_pt_surface
1624
1625       IMPLICIT NONE
1626
1627       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1628       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1629       
1630       k_off = surf%koff
1631       !$OMP PARALLEL DO PRIVATE( i, j, k )
1632       DO  m = 1, surf%ns
1633
1634          i   = surf%i(m)           
1635          j   = surf%j(m)
1636          k   = surf%k(m)
1637
1638          surf%pt_surface(m) = pt(k+k_off,j,i)
1639
1640       ENDDO
1641
1642    END SUBROUTINE calc_pt_surface
1643
1644!
1645!-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. )
1646    SUBROUTINE calc_q_surface
1647
1648       IMPLICIT NONE
1649
1650       INTEGER(iwp) ::  k_off   !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1651       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1652       
1653       k_off = surf%koff
1654       !$OMP PARALLEL DO PRIVATE( i, j, k )
1655       DO  m = 1, surf%ns
1656
1657          i   = surf%i(m)           
1658          j   = surf%j(m)
1659          k   = surf%k(m)
1660
1661          surf%q_surface(m) = q(k+k_off,j,i)
1662
1663       ENDDO
1664
1665    END SUBROUTINE calc_q_surface
1666   
1667!
1668!-- Set virtual potential temperature at surface grid level.
1669!-- ( only for upward-facing surfs )
1670    SUBROUTINE calc_vpt_surface
1671
1672       IMPLICIT NONE
1673
1674       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1675       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1676       
1677       k_off = surf%koff
1678       !$OMP PARALLEL DO PRIVATE( i, j, k )
1679       DO  m = 1, surf%ns
1680
1681          i   = surf%i(m)           
1682          j   = surf%j(m)
1683          k   = surf%k(m)
1684
1685          surf%vpt_surface(m) = vpt(k+k_off,j,i)
1686
1687       ENDDO
1688
1689    END SUBROUTINE calc_vpt_surface
1690   
1691!
1692!-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
1693    SUBROUTINE calc_scaling_parameters
1694
1695       IMPLICIT NONE
1696
1697
1698       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
1699       INTEGER(iwp)  ::  lsp     !< running index for chemical species
1700!
1701!--    Compute theta* at horizontal surfaces
1702       IF ( constant_heatflux  .AND.  .NOT. surf_vertical )  THEN
1703!
1704!--       For a given heat flux in the surface layer:
1705
1706          !$OMP PARALLEL DO PRIVATE( i, j, k )
1707          DO  m = 1, surf%ns
1708
1709             i   = surf%i(m)           
1710             j   = surf%j(m)
1711             k   = surf%k(m)
1712
1713             surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) /                    &
1714                                  ( surf%us(m) + 1E-30_wp )
1715
1716!
1717!--          ts must be limited, because otherwise overflow may occur in case
1718!--          of us=0 when computing ol further below
1719             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
1720             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
1721
1722          ENDDO
1723
1724       ELSEIF ( .NOT. surf_vertical ) THEN
1725!
1726!--       For a given surface temperature:
1727          IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1728
1729             !$OMP PARALLEL DO PRIVATE( i, j, k )
1730             DO  m = 1, surf%ns
1731                i   = surf%i(m)           
1732                j   = surf%j(m)
1733                k   = surf%k(m)
1734
1735                pt(k-1,j,i) = pt_surface
1736             ENDDO
1737          ENDIF
1738
1739          !$OMP PARALLEL DO PRIVATE( z_mo )
1740          DO  m = 1, surf%ns   
1741
1742             z_mo = surf%z_mo(m)
1743
1744             surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) )      &
1745                                  / ( LOG( z_mo / surf%z0h(m) )             &
1746                                      - psi_h( z_mo / surf%ol(m) )          &
1747                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
1748
1749          ENDDO
1750
1751       ENDIF
1752!
1753!--    Compute theta* at vertical surfaces. This is only required in case of
1754!--    land-surface model, in order to compute aerodynamical resistance.
1755       IF ( surf_vertical )  THEN
1756          !$OMP PARALLEL DO PRIVATE( i, j )
1757          DO  m = 1, surf%ns
1758
1759             i          =  surf%i(m)           
1760             j          =  surf%j(m)
1761             surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )
1762!
1763!--          ts must be limited, because otherwise overflow may occur in case
1764!--          of us=0 when computing ol further below
1765             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
1766             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
1767
1768          ENDDO
1769       ENDIF
1770
1771!
1772!--    If required compute q* at horizontal surfaces
1773       IF ( humidity )  THEN
1774          IF ( constant_waterflux  .AND.  .NOT. surf_vertical )  THEN
1775!
1776!--          For a given water flux in the surface layer
1777             !$OMP PARALLEL DO PRIVATE( i, j, k )
1778             DO  m = 1, surf%ns 
1779
1780                i   = surf%i(m)           
1781                j   = surf%j(m)
1782                k   = surf%k(m)
1783                surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) /                &
1784                                               ( surf%us(m) + 1E-30_wp )
1785
1786             ENDDO
1787
1788          ELSEIF ( .NOT. surf_vertical )  THEN
1789             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
1790                             run_coupled )
1791
1792             IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1793                !$OMP PARALLEL DO PRIVATE( i, j, k )
1794                DO  m = 1, surf%ns 
1795
1796                   i   = surf%i(m)           
1797                   j   = surf%j(m)
1798                   k   = surf%k(m)
1799                   q(k-1,j,i) = q_surface
1800                   
1801                ENDDO
1802             ENDIF
1803
1804!
1805!--          Assume saturation for atmosphere coupled to ocean (but not
1806!--          in case of precursor runs)
1807             IF ( coupled_run )  THEN
1808                !$OMP PARALLEL DO PRIVATE( i, j, k, e_s )
1809                DO  m = 1, surf%ns   
1810                   i   = surf%i(m)           
1811                   j   = surf%j(m)
1812                   k   = surf%k(m)
1813                   e_s = 6.1_wp * &
1814                              EXP( 0.07_wp * ( MIN(pt(k-1,j,i),pt(k,j,i))      &
1815                                               - 273.15_wp ) )
1816                   q(k-1,j,i) = 0.622_wp * e_s / ( surface_pressure - e_s )
1817                ENDDO
1818             ENDIF
1819
1820             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
1821               !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1822                DO  m = 1, surf%ns   
1823
1824                   i   = surf%i(m)           
1825                   j   = surf%j(m)
1826                   k   = surf%k(m)
1827   
1828                   z_mo = surf%z_mo(m)
1829
1830                   surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) )    &
1831                                        / ( LOG( z_mo / surf%z0q(m) )          &
1832                                            - psi_h( z_mo / surf%ol(m) )       &
1833                                            + psi_h( surf%z0q(m) /             &
1834                                                     surf%ol(m) ) )
1835                ENDDO
1836             ELSE
1837                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1838                DO  m = 1, surf%ns   
1839
1840                   i   = surf%i(m)           
1841                   j   = surf%j(m)
1842                   k   = surf%k(m)
1843   
1844                   z_mo = surf%z_mo(m)
1845
1846                   surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) )              &
1847                                        / ( LOG( z_mo / surf%z0q(m) )          &
1848                                            - psi_h( z_mo / surf%ol(m) )       &
1849                                            + psi_h( surf%z0q(m) /             &
1850                                                     surf%ol(m) ) )
1851                ENDDO
1852             ENDIF
1853          ENDIF
1854!
1855!--       Compute q* at vertical surfaces
1856          IF ( surf_vertical )  THEN
1857             !$OMP PARALLEL DO PRIVATE( i, j )
1858             DO  m = 1, surf%ns 
1859
1860                i          =  surf%i(m)           
1861                j          =  surf%j(m)
1862                surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )
1863
1864             ENDDO
1865          ENDIF
1866       ENDIF
1867       
1868!
1869!--    If required compute s*
1870       IF ( passive_scalar )  THEN
1871!
1872!--       At horizontal surfaces
1873          IF ( constant_scalarflux  .AND.  .NOT. surf_vertical )  THEN
1874!
1875!--          For a given scalar flux in the surface layer
1876             !$OMP PARALLEL DO PRIVATE( i, j )
1877             DO  m = 1, surf%ns 
1878                i   = surf%i(m)           
1879                j   = surf%j(m)
1880                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1881             ENDDO
1882          ENDIF
1883!
1884!--       At vertical surfaces
1885          IF ( surf_vertical )  THEN
1886             !$OMP PARALLEL DO PRIVATE( i, j )
1887             DO  m = 1, surf%ns 
1888                i          =  surf%i(m)           
1889                j          =  surf%j(m)
1890                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1891             ENDDO
1892          ENDIF
1893       ENDIF
1894
1895!
1896!--    If required compute cs* (chemical species)
1897       IF ( air_chemistry  )  THEN 
1898!
1899!--       At horizontal surfaces                             
1900          DO  lsp = 1, nvar
1901             IF ( constant_csflux(lsp)  .AND.  .NOT.  surf_vertical )  THEN
1902!--             For a given chemical species' flux in the surface layer
1903                !$OMP PARALLEL DO PRIVATE( i, j )
1904                DO  m = 1, surf%ns 
1905                   i   = surf%i(m)           
1906                   j   = surf%j(m)
1907                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1908                ENDDO
1909             ENDIF
1910          ENDDO
1911!
1912!--       At vertical surfaces
1913          IF ( surf_vertical )  THEN
1914             DO  lsp = 1, nvar
1915                !$OMP PARALLEL DO PRIVATE( i, j )
1916                DO  m = 1, surf%ns 
1917                   i   = surf%i(m)           
1918                   j   = surf%j(m)
1919                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1920                ENDDO
1921             ENDDO
1922          ENDIF
1923       ENDIF
1924
1925!
1926!--    If required compute qc* and nc*
1927       IF ( cloud_physics  .AND.  microphysics_morrison  .AND.                 &
1928            .NOT. surf_vertical )  THEN
1929          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1930          DO  m = 1, surf%ns   
1931
1932             i    = surf%i(m)           
1933             j    = surf%j(m)
1934             k    = surf%k(m)
1935
1936             z_mo = surf%z_mo(m)
1937
1938             surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) )                 &
1939                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1940                                 - psi_h( z_mo / surf%ol(m) )                  &
1941                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1942
1943             surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) )                 &
1944                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1945                                 - psi_h( z_mo / surf%ol(m) )                  &
1946                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1947          ENDDO
1948
1949       ENDIF
1950
1951!
1952!--    If required compute qr* and nr*
1953       IF ( cloud_physics  .AND.  microphysics_seifert  .AND.                  &
1954            .NOT. surf_vertical )  THEN
1955          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1956          DO  m = 1, surf%ns   
1957
1958             i    = surf%i(m)           
1959             j    = surf%j(m)
1960             k    = surf%k(m)
1961
1962             z_mo = surf%z_mo(m)
1963
1964             surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) )                 &
1965                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1966                                 - psi_h( z_mo / surf%ol(m) )                  &
1967                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1968
1969             surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) )                 &
1970                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1971                                 - psi_h( z_mo / surf%ol(m) )                  &
1972                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1973          ENDDO
1974
1975       ENDIF
1976
1977    END SUBROUTINE calc_scaling_parameters
1978
1979
1980
1981!
1982!-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
1983    SUBROUTINE calc_surface_fluxes
1984
1985       IMPLICIT NONE
1986
1987       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
1988       INTEGER(iwp)  ::  lsp     !< running index for chemical species
1989
1990       REAL(wp)                            ::  dum     !< dummy to precalculate logarithm
1991       REAL(wp)                            ::  flag_u  !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
1992       REAL(wp)                            ::  flag_v  !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
1993       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_i     !< u-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1994       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_i     !< v-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1995       REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_i     !< w-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1996
1997!
1998!--    Calcuate surface fluxes at horizontal walls
1999       IF ( .NOT. surf_vertical )  THEN
2000!
2001!--       Compute u'w' for the total model domain at upward-facing surfaces.
2002!--       First compute the corresponding component of u* and square it.
2003          IF ( .NOT. downward )  THEN
2004             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2005             DO  m = 1, surf%ns 
2006   
2007                i = surf%i(m)           
2008                j = surf%j(m)
2009                k = surf%k(m)
2010
2011                z_mo = surf%z_mo(m)
2012
2013                surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) )               &
2014                              / ( LOG( z_mo / surf%z0(m) )                     &
2015                                  - psi_m( z_mo / surf%ol(m) )                 &
2016                                  + psi_m( surf%z0(m) / surf%ol(m) ) )
2017!
2018!--             Please note, the computation of usws is not fully accurate. Actually
2019!--             a further interpolation of us onto the u-grid, where usws is defined,
2020!--             is required. However, this is not done as this would require several
2021!--             data transfers between 2D-grid and the surf-type.
2022!--             The impact of the missing interpolation is negligible as several
2023!--             tests had shown.
2024!--             Same also for ol. 
2025                surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)
2026
2027             ENDDO
2028!
2029!--       At downward-facing surfaces
2030          ELSE
2031             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2032             DO  m = 1, surf%ns 
2033   
2034                i = surf%i(m)           
2035                j = surf%j(m)
2036                k = surf%k(m)
2037
2038                z_mo = surf%z_mo(m)
2039
2040                surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
2041                surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
2042
2043             ENDDO     
2044          ENDIF
2045
2046!
2047!--       Compute v'w' for the total model domain.
2048!--       First compute the corresponding component of u* and square it.
2049!--       Upward-facing surfaces
2050          IF ( .NOT. downward )  THEN
2051             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2052             DO  m = 1, surf%ns 
2053                i = surf%i(m)           
2054                j = surf%j(m)
2055                k = surf%k(m)
2056
2057                z_mo = surf%z_mo(m)
2058
2059                surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) )               &
2060                           / ( LOG( z_mo / surf%z0(m) )                        &
2061                               - psi_m( z_mo / surf%ol(m) )                    &
2062                               + psi_m( surf%z0(m) / surf%ol(m) ) )
2063!
2064!--             Please note, the computation of vsws is not fully accurate. Actually
2065!--             a further interpolation of us onto the v-grid, where vsws is defined,
2066!--             is required. However, this is not done as this would require several
2067!--             data transfers between 2D-grid and the surf-type.
2068!--             The impact of the missing interpolation is negligible as several
2069!--             tests had shown.
2070!--             Same also for ol. 
2071                surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)
2072             ENDDO
2073!
2074!--       Downward-facing surfaces
2075          ELSE
2076             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2077             DO  m = 1, surf%ns 
2078                i = surf%i(m)           
2079                j = surf%j(m)
2080                k = surf%k(m)
2081
2082                z_mo = surf%z_mo(m)
2083
2084                surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
2085                surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
2086             ENDDO
2087          ENDIF
2088!
2089!--       Compute the vertical kinematic heat flux
2090          IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point&
2091               <=  skip_time_do_lsm  .AND. simulated_time > 0.0_wp ) .OR.      &
2092               .NOT.  land_surface )  .AND.  .NOT. urban_surface  .AND.        &
2093               .NOT. downward )  THEN
2094             !$OMP PARALLEL DO PRIVATE( i, j, k )
2095             DO  m = 1, surf%ns
2096                i    = surf%i(m)           
2097                j    = surf%j(m)
2098                k    = surf%k(m)
2099                surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
2100             ENDDO
2101          ENDIF
2102!
2103!--       Compute the vertical water flux
2104          IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
2105               ( ( time_since_reference_point <= skip_time_do_lsm  .AND.       &
2106               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface  )  .AND.   &
2107               .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
2108             !$OMP PARALLEL DO PRIVATE( i, j, k )
2109             DO  m = 1, surf%ns
2110                i    = surf%i(m)           
2111                j    = surf%j(m)
2112                k    = surf%k(m)
2113                surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
2114             ENDDO
2115          ENDIF
2116!
2117!--       Compute the vertical scalar flux
2118          IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.       &
2119                .NOT.  downward )  THEN
2120             !$OMP PARALLEL DO PRIVATE( i, j )
2121             DO  m = 1, surf%ns   
2122
2123                i    = surf%i(m)           
2124                j    = surf%j(m)
2125                surf%ssws(m) = -surf%ss(m) * surf%us(m)
2126
2127             ENDDO
2128          ENDIF   
2129!
2130!--       Compute the vertical chemical species' flux
2131          DO  lsp = 1, nvar
2132             IF (  .NOT.  constant_csflux(lsp)  .AND.  air_chemistry  .AND.    &
2133                   .NOT.  downward )  THEN
2134                !$OMP PARALLEL DO PRIVATE( i, j )
2135                DO  m = 1, surf%ns   
2136                   i     = surf%i(m)           
2137                   j     = surf%j(m)
2138                   surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
2139                ENDDO
2140             ENDIF   
2141          ENDDO
2142
2143!
2144!--       Compute (turbulent) fluxes of cloud water content and cloud drop conc.
2145          IF ( cloud_physics  .AND.  microphysics_morrison  .AND.              &
2146               .NOT. downward)  THEN
2147             !$OMP PARALLEL DO PRIVATE( i, j )
2148             DO  m = 1, surf%ns   
2149
2150                i    = surf%i(m)           
2151                j    = surf%j(m)
2152
2153                surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
2154                surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
2155             ENDDO
2156          ENDIF   
2157!
2158!--       Compute (turbulent) fluxes of rain water content and rain drop conc.
2159          IF ( cloud_physics  .AND.  microphysics_seifert  .AND.               &
2160               .NOT. downward)  THEN
2161             !$OMP PARALLEL DO PRIVATE( i, j )
2162             DO  m = 1, surf%ns   
2163
2164                i    = surf%i(m)           
2165                j    = surf%j(m)
2166
2167                surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
2168                surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
2169             ENDDO
2170          ENDIF
2171
2172!
2173!--       Bottom boundary condition for the TKE.
2174          IF ( ibc_e_b == 2 )  THEN
2175             !$OMP PARALLEL DO PRIVATE( i, j, k )
2176             DO  m = 1, surf%ns   
2177
2178                i    = surf%i(m)           
2179                j    = surf%j(m)
2180                k    = surf%k(m)
2181
2182                e(k,j,i) = ( surf%us(m) / 0.1_wp )**2
2183!
2184!--             As a test: cm = 0.4
2185!               e(k,j,i) = ( us(j,i) / 0.4_wp )**2
2186                e(k-1,j,i)   = e(k,j,i)
2187
2188             ENDDO
2189          ENDIF
2190!
2191!--    Calcuate surface fluxes at vertical surfaces. No stability is considered.
2192       ELSE
2193!
2194!--       Compute usvs l={0,1} and vsus l={2,3}
2195          IF ( mom_uv )  THEN
2196!
2197!--          Generalize computation by introducing flags. At north- and south-
2198!--          facing surfaces u-component is used, at east- and west-facing
2199!--          surfaces v-component is used.
2200             flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0  .OR.  l == 1 )   
2201             flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2  .OR.  l == 3 )   
2202             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
2203             DO  m = 1, surf%ns 
2204                i = surf%i(m)           
2205                j = surf%j(m)
2206                k = surf%k(m)
2207
2208                z_mo = surf%z_mo(m)
2209
2210                surf%mom_flux_uv(m) = kappa *                                  &
2211                                ( flag_u * u(k,j,i) + flag_v * v(k,j,i) )  /   &
2212                                                        LOG( z_mo / surf%z0(m) )
2213
2214               surf%mom_flux_uv(m) =                                           &
2215                                    - surf%mom_flux_uv(m) * surf%us(m)
2216             ENDDO
2217          ENDIF
2218!
2219!--       Compute wsus l={0,1} and wsvs l={2,3}
2220          IF ( mom_w )  THEN
2221             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
2222             DO  m = 1, surf%ns 
2223                i = surf%i(m)           
2224                j = surf%j(m)
2225                k = surf%k(m)
2226
2227                z_mo = surf%z_mo(m)
2228
2229                surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
2230
2231                surf%mom_flux_w(m) =                                           &
2232                                     - surf%mom_flux_w(m) * surf%us(m)
2233             ENDDO
2234          ENDIF
2235!
2236!--       Compute momentum fluxes used for subgrid-scale TKE production at
2237!--       vertical surfaces. In constrast to the calculated momentum fluxes at
2238!--       vertical surfaces before, which are defined on the u/v/w-grid,
2239!--       respectively), the TKE fluxes are defined at the scalar grid.
2240!--       
2241          IF ( mom_tke )  THEN
2242!
2243!--          Precalculate velocity components at scalar grid point.
2244             ALLOCATE( u_i(1:surf%ns) )
2245             ALLOCATE( v_i(1:surf%ns) )
2246             ALLOCATE( w_i(1:surf%ns) )
2247
2248             IF ( l == 0  .OR.  l == 1 )  THEN
2249                !$OMP PARALLEL DO PRIVATE( i, j, k )
2250                DO  m = 1, surf%ns 
2251                   i = surf%i(m)           
2252                   j = surf%j(m)
2253                   k = surf%k(m)
2254
2255                   u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
2256                   v_i(m) = 0.0_wp
2257                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
2258                ENDDO
2259             ELSE
2260                !$OMP PARALLEL DO PRIVATE( i, j, k )
2261                DO  m = 1, surf%ns 
2262                   i = surf%i(m)           
2263                   j = surf%j(m)
2264                   k = surf%k(m)
2265
2266                   u_i(m) = 0.0_wp
2267                   v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
2268                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
2269                ENDDO
2270             ENDIF
2271
2272             !$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )
2273             DO  m = 1, surf%ns 
2274                i = surf%i(m)           
2275                j = surf%j(m)
2276
2277                z_mo = surf%z_mo(m)
2278
2279                dum = kappa / LOG( z_mo / surf%z0(m) )
2280!
2281!--             usvs (l=0,1) and vsus (l=2,3)
2282                surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )
2283!
2284!--             wsvs (l=0,1) and wsus (l=2,3)
2285                surf%mom_flux_tke(1,m) = dum * w_i(m)
2286
2287                surf%mom_flux_tke(0:1,m) =                                     &
2288                               - surf%mom_flux_tke(0:1,m) * surf%us(m)
2289             ENDDO
2290!
2291!--          Deallocate temporary arrays
2292             DEALLOCATE( u_i )             
2293             DEALLOCATE( v_i ) 
2294             DEALLOCATE( w_i ) 
2295          ENDIF
2296       ENDIF
2297
2298    END SUBROUTINE calc_surface_fluxes
2299
2300
2301!
2302!-- Integrated stability function for momentum
2303    FUNCTION psi_m( zeta ) 
2304       
2305       USE kinds
2306
2307       IMPLICIT NONE
2308
2309       REAL(wp)            :: psi_m !< Integrated similarity function result
2310       REAL(wp)            :: zeta  !< Stability parameter z/L
2311       REAL(wp)            :: x     !< dummy variable
2312
2313       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
2314       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
2315       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
2316       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
2317       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
2318       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
2319
2320
2321       IF ( zeta < 0.0_wp )  THEN
2322          x = SQRT( SQRT( 1.0_wp  - 16.0_wp * zeta ) )
2323          psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
2324                  * ( 1.0_wp + x**2 ) * 0.125_wp )
2325       ELSE
2326
2327          psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta         &
2328                   - bc_d_d
2329!
2330!--       Old version for stable conditions (only valid for z/L < 0.5)
2331!--       psi_m = - 5.0_wp * zeta
2332
2333       ENDIF
2334
2335    END FUNCTION psi_m
2336
2337
2338!
2339!-- Integrated stability function for heat and moisture
2340    FUNCTION psi_h( zeta ) 
2341       
2342       USE kinds
2343
2344       IMPLICIT NONE
2345
2346       REAL(wp)            :: psi_h !< Integrated similarity function result
2347       REAL(wp)            :: zeta  !< Stability parameter z/L
2348       REAL(wp)            :: x     !< dummy variable
2349
2350       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
2351       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
2352       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
2353       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
2354       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
2355       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
2356
2357
2358       IF ( zeta < 0.0_wp )  THEN
2359          x = SQRT( 1.0_wp  - 16.0_wp * zeta )
2360          psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
2361       ELSE
2362          psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp          &
2363                  + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d             &
2364                  + 1.0_wp
2365!
2366!--       Old version for stable conditions (only valid for z/L < 0.5)
2367!--       psi_h = - 5.0_wp * zeta
2368       ENDIF
2369
2370    END FUNCTION psi_h
2371
2372
2373!------------------------------------------------------------------------------!
2374! Description:
2375! ------------
2376!> Calculates stability function for momentum
2377!>
2378!> @author Hauke Wurps
2379!------------------------------------------------------------------------------!
2380    FUNCTION phi_m( zeta ) 
2381   
2382       IMPLICIT NONE
2383   
2384       REAL(wp)            :: phi_m         !< Value of the function
2385       REAL(wp)            :: zeta          !< Stability parameter z/L
2386   
2387       REAL(wp), PARAMETER :: a = 16.0_wp   !< constant
2388       REAL(wp), PARAMETER :: c = 5.0_wp    !< constant
2389   
2390       IF ( zeta < 0.0_wp )  THEN
2391          phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) )
2392       ELSE
2393          phi_m = 1.0_wp + c * zeta
2394       ENDIF
2395   
2396    END FUNCTION phi_m
2397
2398 END MODULE surface_layer_fluxes_mod
Note: See TracBrowser for help on using the repository browser.