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

Last change on this file since 2321 was 2321, checked in by schwenkel, 7 years ago

Bugfix for calculate ol via lookup table

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