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

Last change on this file since 3581 was 3547, checked in by suehring, 6 years ago

variable description added some routines

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