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

Last change on this file since 3352 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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