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

Last change on this file since 2232 was 2232, checked in by suehring, 7 years ago

Adjustments according new topography and surface-modelling concept implemented

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