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

Last change on this file since 2299 was 2299, checked in by maronga, 4 years ago

improvements for spinup mechanism

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