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

Last change on this file since 2292 was 2292, checked in by schwenkel, 4 years ago

implementation of new bulk microphysics scheme

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