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

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

extended by cloud_droplets option

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