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

Last change on this file since 3271 was 3271, checked in by suehring, 3 years ago

Several critical bugfixes: initialization of surface temperature at water bodies; consider special case in initialization of vertical natural surface elements

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