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

Last change on this file since 3147 was 3147, checked in by maronga, 3 years ago

last commit documented

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