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

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

last commit documented

  • Property svn:keywords set to Id
File size: 85.2 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 3148 2018-07-19 05:45:25Z 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                              ( g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m)   &
988                              )**(2.0_wp/3.0_wp)                               &
989                                )
990
991                               
992
993       ENDDO
994
995    END SUBROUTINE calc_uvw_abs
996
997
998!------------------------------------------------------------------------------!
999! Description:
1000! ------------
1001!> Compute the absolute value of the horizontal velocity (relative to the
1002!> surface) for horizontal surface elements. This is required by all methods.
1003!------------------------------------------------------------------------------!
1004    SUBROUTINE calc_uvw_abs_v_ugrid
1005
1006       IMPLICIT NONE
1007
1008       INTEGER(iwp) ::  i             !< running index x direction
1009       INTEGER(iwp) ::  j             !< running index y direction
1010       INTEGER(iwp) ::  k             !< running index z direction
1011       INTEGER(iwp) ::  m             !< running index surface elements
1012
1013       REAL(wp)     ::  u_i
1014       REAL(wp)     ::  w_i
1015
1016
1017       DO  m = 1, surf%ns
1018          i   = surf%i(m)           
1019          j   = surf%j(m)
1020          k   = surf%k(m)
1021!
1022!--       Compute the absolute value of the surface parallel velocity on u-grid.
1023          u_i  = u(k,j,i)
1024          w_i  = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) +                       &
1025                             w(k,j,i-1)   + w(k,j,i) ) 
1026
1027          surf%uvw_abs(m) = SQRT( u_i**2 + w_i**2 ) 
1028
1029       ENDDO
1030
1031    END SUBROUTINE calc_uvw_abs_v_ugrid
1032
1033!------------------------------------------------------------------------------!
1034! Description:
1035! ------------
1036!> Compute the absolute value of the horizontal velocity (relative to the
1037!> surface) for horizontal surface elements. This is required by all methods.
1038!------------------------------------------------------------------------------!
1039    SUBROUTINE calc_uvw_abs_v_vgrid
1040
1041       IMPLICIT NONE
1042
1043       INTEGER(iwp) ::  i             !< running index x direction
1044       INTEGER(iwp) ::  j             !< running index y direction
1045       INTEGER(iwp) ::  k             !< running index z direction
1046       INTEGER(iwp) ::  m             !< running index surface elements
1047
1048       REAL(wp)     ::  v_i
1049       REAL(wp)     ::  w_i
1050
1051
1052       DO  m = 1, surf%ns
1053          i   = surf%i(m)           
1054          j   = surf%j(m)
1055          k   = surf%k(m)
1056
1057          v_i  = u(k,j,i)
1058          w_i  = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) +                       &
1059                             w(k,j-1,i)   + w(k,j,i) ) 
1060
1061          surf%uvw_abs(m) = SQRT( v_i**2 + w_i**2 ) 
1062
1063       ENDDO
1064
1065    END SUBROUTINE calc_uvw_abs_v_vgrid
1066
1067!------------------------------------------------------------------------------!
1068! Description:
1069! ------------
1070!> Compute the absolute value of the horizontal velocity (relative to the
1071!> surface) for horizontal surface elements. This is required by all methods.
1072!------------------------------------------------------------------------------!
1073    SUBROUTINE calc_uvw_abs_v_wgrid
1074
1075       IMPLICIT NONE
1076
1077       INTEGER(iwp) ::  i             !< running index x direction
1078       INTEGER(iwp) ::  j             !< running index y direction
1079       INTEGER(iwp) ::  k             !< running index z direction
1080       INTEGER(iwp) ::  m             !< running index surface elements
1081
1082       REAL(wp)     ::  u_i
1083       REAL(wp)     ::  v_i
1084       REAL(wp)     ::  w_i
1085!
1086!--    North- (l=0) and south-facing (l=1) surfaces
1087       IF ( l == 0  .OR.  l == 1 )  THEN
1088          DO  m = 1, surf%ns
1089             i   = surf%i(m)           
1090             j   = surf%j(m)
1091             k   = surf%k(m)
1092
1093             u_i  = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) +                    &
1094                                u(k,j,i+1)   + u(k,j,i) )
1095             v_i  = 0.0_wp
1096             w_i  = w(k,j,i)
1097
1098             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1099          ENDDO
1100!
1101!--    East- (l=2) and west-facing (l=3) surfaces
1102       ELSE
1103          DO  m = 1, surf%ns
1104             i   = surf%i(m)           
1105             j   = surf%j(m)
1106             k   = surf%k(m)
1107
1108             u_i  = 0.0_wp
1109             v_i  = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) +                    &
1110                                v(k,j+1,i)   + v(k,j,i) )
1111             w_i  = w(k,j,i)
1112
1113             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1114          ENDDO
1115       ENDIF           
1116
1117    END SUBROUTINE calc_uvw_abs_v_wgrid
1118
1119!------------------------------------------------------------------------------!
1120! Description:
1121! ------------
1122!> Compute the absolute value of the horizontal velocity (relative to the
1123!> surface) for horizontal surface elements. This is required by all methods.
1124!------------------------------------------------------------------------------!
1125    SUBROUTINE calc_uvw_abs_v_sgrid
1126
1127       IMPLICIT NONE
1128
1129       INTEGER(iwp) ::  i             !< running index x direction
1130       INTEGER(iwp) ::  j             !< running index y direction
1131       INTEGER(iwp) ::  k             !< running index z direction
1132       INTEGER(iwp) ::  m             !< running index surface elements
1133
1134       REAL(wp)     ::  u_i
1135       REAL(wp)     ::  v_i
1136       REAL(wp)     ::  w_i
1137
1138!
1139!--    North- (l=0) and south-facing (l=1) walls
1140       IF ( l == 0  .OR.  l == 1 )  THEN
1141          DO  m = 1, surf%ns
1142             i   = surf%i(m)           
1143             j   = surf%j(m)
1144             k   = surf%k(m)
1145
1146             u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
1147             v_i = 0.0_wp
1148             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1149
1150             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1151          ENDDO
1152!
1153!--    East- (l=2) and west-facing (l=3) walls
1154       ELSE
1155          DO  m = 1, surf%ns
1156             i   = surf%i(m)           
1157             j   = surf%j(m)
1158             k   = surf%k(m)
1159
1160             u_i = 0.0_wp
1161             v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
1162             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1163
1164             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1165          ENDDO
1166       ENDIF 
1167
1168    END SUBROUTINE calc_uvw_abs_v_sgrid
1169
1170
1171!------------------------------------------------------------------------------!
1172! Description:
1173! ------------
1174!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
1175!------------------------------------------------------------------------------!
1176    SUBROUTINE calc_ol
1177
1178       IMPLICIT NONE
1179
1180       INTEGER(iwp) ::  iter    !< Newton iteration step
1181       INTEGER(iwp) ::  li      !< look index
1182       INTEGER(iwp) ::  m       !< loop variable over all horizontal wall elements
1183
1184       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
1185                       f_d_ol, & !< Derivative of f
1186                       ol_l,   & !< Lower bound of L for Newton iteration
1187                       ol_m,   & !< Previous value of L for Newton iteration
1188                       ol_old, & !< Previous time step value of L
1189                       ol_u      !< Upper bound of L for Newton iteration
1190
1191       IF ( TRIM( most_method ) /= 'circular' )  THEN
1192!
1193!--       Evaluate bulk Richardson number (calculation depends on
1194!--       definition based on setting of boundary conditions
1195          IF ( ibc_pt_b /= 1 )  THEN
1196             IF ( humidity )  THEN
1197                !$OMP PARALLEL DO PRIVATE( z_mo )
1198                DO  m = 1, surf%ns
1199
1200                   z_mo = surf%z_mo(m)
1201
1202                   surf%rib(m) = g * z_mo                                      &
1203                                 * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
1204                                 / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
1205                                 + 1.0E-20_wp )
1206                ENDDO
1207             ELSE
1208                !$OMP PARALLEL DO PRIVATE( z_mo )
1209                DO  m = 1, surf%ns
1210
1211                   z_mo = surf%z_mo(m)
1212
1213                   surf%rib(m) = g * z_mo                                      &
1214                                 * ( surf%pt1(m) - surf%pt_surface(m) )        &
1215                                 / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
1216                                 + 1.0E-20_wp )
1217                ENDDO
1218             ENDIF
1219          ELSE
1220             IF ( humidity )  THEN
1221                !$OMP PARALLEL DO PRIVATE( k, z_mo )
1222                DO  m = 1, surf%ns
1223
1224                   k   = surf%k(m)
1225
1226                   z_mo = surf%z_mo(m)
1227
1228                   surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
1229                           * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
1230                           * surf%pt1(m) * surf%qsws(m) ) *                    &
1231                             drho_air_zw(k-1)                /                 &
1232                         ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
1233                           + 1.0E-20_wp )
1234                ENDDO
1235             ELSE
1236                !$OMP PARALLEL DO PRIVATE( k, z_mo )
1237                DO  m = 1, surf%ns
1238
1239                   k   = surf%k(m)
1240
1241                   z_mo = surf%z_mo(m)
1242
1243                   surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
1244                                        drho_air_zw(k-1)            /          &
1245                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
1246                           + 1.0E-20_wp )
1247                ENDDO
1248             ENDIF
1249          ENDIF
1250
1251       ENDIF
1252
1253
1254!
1255!--    Calculate the Obukhov length either using a Newton iteration
1256!--    method, via a lookup table, or using the old circular way
1257       IF ( TRIM( most_method ) == 'newton' )  THEN
1258
1259          DO  m = 1, surf%ns
1260
1261             i   = surf%i(m)           
1262             j   = surf%j(m)
1263
1264             z_mo = surf%z_mo(m)
1265
1266!
1267!--          Store current value in case the Newton iteration fails
1268             ol_old = surf%ol(m)
1269
1270!
1271!--          Ensure that the bulk Richardson number and the Obukhov
1272!--          length have the same sign
1273             IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
1274                  ABS( surf%ol(m) ) == ol_max )  THEN
1275                IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
1276                IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
1277             ENDIF
1278!
1279!--          Iteration to find Obukhov length
1280             iter = 0
1281             DO
1282                iter = iter + 1
1283!
1284!--             In case of divergence, use the value of the previous time step
1285                IF ( iter > 1000 )  THEN
1286                   surf%ol(m) = ol_old
1287                   EXIT
1288                ENDIF
1289
1290                ol_m = surf%ol(m)
1291                ol_l = ol_m - 0.001_wp * ol_m
1292                ol_u = ol_m + 0.001_wp * ol_m
1293
1294
1295                IF ( ibc_pt_b /= 1 )  THEN
1296!
1297!--                Calculate f = Ri - [...]/[...]^2 = 0
1298                   f = surf%rib(m) - ( z_mo / ol_m ) * (                       &
1299                                                 LOG( z_mo / surf%z0h(m) )     &
1300                                                 - psi_h( z_mo / ol_m )        &
1301                                                 + psi_h( surf%z0h(m) /        &
1302                                                          ol_m )               &
1303                                                               )               &
1304                                              / ( LOG( z_mo / surf%z0(m) )     &
1305                                                 - psi_m( z_mo / ol_m )        &
1306                                                 + psi_m( surf%z0(m) /         &
1307                                                          ol_m )               &
1308                                                 )**2
1309
1310!
1311!--                 Calculate df/dL
1312                    f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /               &
1313                                                             surf%z0h(m) )     &
1314                                            - psi_h( z_mo / ol_u )             &
1315                                            + psi_h( surf%z0h(m) / ol_u )      &
1316                                              )                                &
1317                                            / ( LOG( z_mo / surf%z0(m) )       &
1318                                            - psi_m( z_mo / ol_u )             &
1319                                            + psi_m( surf%z0(m) / ol_u )       &
1320                                              )**2                             &
1321                           + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )     &
1322                                            - psi_h( z_mo / ol_l )             &
1323                                            + psi_h( surf%z0h(m) / ol_l )      &
1324                                              )                                &
1325                                            / ( LOG( z_mo / surf%z0(m) )       &
1326                                            - psi_m( z_mo / ol_l )             &
1327                                            + psi_m( surf%z0(m) / ol_l )       &
1328                                              )**2                             &
1329                             ) / ( ol_u - ol_l )
1330                ELSE
1331!
1332!--                Calculate f = Ri - 1 /[...]^3 = 0
1333                   f = surf%rib(m) - ( z_mo / ol_m ) /                         &
1334                                                ( LOG( z_mo / surf%z0(m) )     &
1335                                            - psi_m( z_mo / ol_m )             &
1336                                            + psi_m( surf%z0(m) / ol_m )       &
1337                                                )**3
1338
1339!
1340!--                Calculate df/dL
1341                   f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo /                &
1342                                                            surf%z0(m) )       &
1343                                            - psi_m( z_mo / ol_u )             &
1344                                            + psi_m( surf%z0(m) / ol_u )       &
1345                                                     )**3                      &
1346                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
1347                                            - psi_m( z_mo / ol_l )             &
1348                                            + psi_m( surf%z0(m) / ol_l )       &
1349                                               )**3                            &
1350                                  ) / ( ol_u - ol_l )
1351                ENDIF
1352!
1353!--             Calculate new L
1354                surf%ol(m) = ol_m - f / f_d_ol
1355
1356!
1357!--             Ensure that the bulk Richardson number and the Obukhov
1358!--             length have the same sign and ensure convergence.
1359                IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
1360!
1361!--             If unrealistic value occurs, set L to the maximum
1362!--             value that is allowed
1363                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
1364                   surf%ol(m) = ol_max
1365                   EXIT
1366                ENDIF
1367!
1368!--             Check for convergence
1369                IF ( ABS( ( surf%ol(m) - ol_m ) /                              &
1370                     surf%ol(m) ) < 1.0E-4_wp )  THEN
1371                   EXIT
1372                ELSE
1373                   CYCLE
1374                ENDIF
1375
1376             ENDDO
1377          ENDDO
1378
1379       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
1380
1381          !$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
1382          DO  m = 1, surf%ns
1383
1384             i   = surf%i(m)           
1385             j   = surf%j(m)
1386!
1387!--          If the bulk Richardson number is outside the range of the lookup
1388!--          table, set it to the exceeding threshold value
1389             IF ( surf%rib(m) < rib_min )  surf%rib(m) = rib_min
1390             IF ( surf%rib(m) > rib_max )  surf%rib(m) = rib_max
1391
1392!
1393!--          Find the correct index bounds for linear interpolation. As the
1394!--          Richardson number will not differ very much from time step to
1395!--          time step , use the index from the last step and search in the
1396!--          correct direction
1397             li = li_bnd
1398             IF ( rib_tab(li) - surf%rib(m) > 0.0_wp )  THEN
1399                DO  WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp  .AND.  li > 0 )
1400                   li = li-1
1401                ENDDO
1402             ELSE
1403                DO  WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp                 &
1404                           .AND.  li < num_steps-1 )
1405                   li = li+1
1406                ENDDO
1407             ENDIF
1408             li_bnd = li
1409
1410!
1411!--          Linear interpolation to find the correct value of z/L
1412             surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) )       &
1413                         / (  rib_tab(li) - rib_tab(li-1) )                    &
1414                         * ( surf%rib(m) - rib_tab(li-1) ) )
1415
1416          ENDDO
1417
1418       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
1419
1420          IF ( .NOT. humidity )  THEN
1421             !$OMP PARALLEL DO PRIVATE( z_mo )
1422             DO  m = 1, surf%ns
1423
1424                z_mo = surf%z_mo(m)
1425
1426                surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
1427                                       ( kappa * g *                           &
1428                                         surf%ts(m) + 1E-30_wp )
1429!
1430!--             Limit the value range of the Obukhov length.
1431!--             This is necessary for very small velocities (u,v --> 1), because
1432!--             the absolute value of ol can then become very small, which in
1433!--             consequence would result in very large shear stresses and very
1434!--             small momentum fluxes (both are generally unrealistic).
1435                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1436                   surf%ol(m) = z_mo / zeta_min
1437                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1438                   surf%ol(m) = z_mo / zeta_max
1439
1440             ENDDO
1441          ELSEIF ( cloud_physics  .OR.  cloud_droplets )  THEN
1442             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1443             DO  m = 1, surf%ns
1444
1445                i   = surf%i(m)           
1446                j   = surf%j(m)
1447                k   = surf%k(m)
1448
1449                z_mo = surf%z_mo(m)
1450
1451
1452                surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
1453                    ( kappa * g * ( surf%ts(m) +                               &
1454                         0.61_wp * surf%pt1(m) * surf%us(m)                    &
1455                       + 0.61_wp * surf%qv1(m) * surf%ts(m) -                  &
1456                                   surf%ts(m)  * ql(k,j,i) ) + 1E-30_wp )
1457!
1458!--             Limit the value range of the Obukhov length.
1459!--             This is necessary for very small velocities (u,v --> 1), because
1460!--             the absolute value of ol can then become very small, which in
1461!--             consequence would result in very large shear stresses and very
1462!--             small momentum fluxes (both are generally unrealistic).
1463                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1464                   surf%ol(m) = z_mo / zeta_min
1465                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1466                   surf%ol(m) = z_mo / zeta_max
1467
1468             ENDDO
1469          ELSE
1470
1471             !$OMP PARALLEL DO PRIVATE( z_mo )
1472             DO  m = 1, surf%ns
1473
1474                z_mo = surf%z_mo(m)
1475
1476                surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
1477                    ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
1478                                    surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
1479                                    surf%ts(m) ) + 1E-30_wp )
1480
1481!
1482!--             Limit the value range of the Obukhov length.
1483!--             This is necessary for very small velocities (u,v --> 1), because
1484!--             the absolute value of ol can then become very small, which in
1485!--             consequence would result in very large shear stresses and very
1486!--             small momentum fluxes (both are generally unrealistic).
1487                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1488                   surf%ol(m) = z_mo / zeta_min
1489                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1490                   surf%ol(m) = z_mo / zeta_max
1491
1492             ENDDO
1493
1494          ENDIF
1495
1496       ENDIF
1497
1498    END SUBROUTINE calc_ol
1499
1500!
1501!-- Calculate friction velocity u*
1502    SUBROUTINE calc_us
1503
1504       IMPLICIT NONE
1505
1506       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1507
1508!
1509!--    Compute u* at horizontal surfaces at the scalars' grid points
1510       IF ( .NOT. surf_vertical )  THEN
1511!
1512!--       Compute u* at upward-facing surfaces
1513          IF ( .NOT. downward )  THEN
1514             !$OMP PARALLEL  DO PRIVATE( z_mo )
1515             DO  m = 1, surf%ns
1516
1517                z_mo = surf%z_mo(m)
1518!
1519!--             Compute u* at the scalars' grid points
1520                surf%us(m) = kappa * surf%uvw_abs(m) /                         &
1521                            ( LOG( z_mo / surf%z0(m) )                         &
1522                           - psi_m( z_mo / surf%ol(m) )                        &
1523                           + psi_m( surf%z0(m) / surf%ol(m) ) )
1524   
1525             ENDDO
1526!
1527!--       Compute u* at downward-facing surfaces. This case, do not consider
1528!--       any stability.
1529          ELSE
1530             !$OMP PARALLEL  DO PRIVATE( z_mo )
1531             DO  m = 1, surf%ns
1532
1533                z_mo = surf%z_mo(m)
1534!
1535!--             Compute u* at the scalars' grid points
1536                surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1537   
1538             ENDDO
1539          ENDIF
1540!
1541!--    Compute u* at vertical surfaces at the u/v/v grid, respectively.
1542!--    No stability is considered in this case.
1543       ELSE
1544          !$OMP PARALLEL DO PRIVATE( z_mo )
1545          DO  m = 1, surf%ns
1546             z_mo = surf%z_mo(m)
1547
1548             surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1549          ENDDO
1550       ENDIF
1551
1552    END SUBROUTINE calc_us
1553
1554!
1555!-- Calculate potential temperature, specific humidity, and virtual potential
1556!-- temperature at first grid level
1557!-- ( only for upward-facing surfs )
1558    SUBROUTINE calc_pt_q
1559
1560       IMPLICIT NONE
1561
1562       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1563
1564       !$OMP PARALLEL DO PRIVATE( i, j, k )
1565       DO  m = 1, surf%ns
1566
1567          i   = surf%i(m)           
1568          j   = surf%j(m)
1569          k   = surf%k(m)
1570
1571          IF ( cloud_physics ) THEN
1572             surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
1573             surf%qv1(m) = q(k,j,i) - ql(k,j,i)
1574          ELSEIF( cloud_droplets ) THEN
1575             surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
1576             surf%qv1(m) = q(k,j,i) 
1577          ELSE
1578             surf%pt1(m) = pt(k,j,i)
1579             IF ( humidity )  THEN
1580                surf%qv1(m) = q(k,j,i)
1581             ELSE
1582                surf%qv1(m) = 0.0_wp
1583             ENDIF
1584          ENDIF
1585
1586          IF ( humidity )  THEN
1587             surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
1588          ENDIF
1589         
1590       ENDDO
1591
1592    END SUBROUTINE calc_pt_q
1593
1594
1595!
1596!-- Calculate potential temperature at first grid level
1597!-- ( only for upward-facing surfs )
1598    SUBROUTINE calc_pt_surface
1599
1600       IMPLICIT NONE
1601
1602       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1603       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1604       
1605       k_off = surf%koff
1606       !$OMP PARALLEL DO PRIVATE( i, j, k )
1607       DO  m = 1, surf%ns
1608
1609          i   = surf%i(m)           
1610          j   = surf%j(m)
1611          k   = surf%k(m)
1612
1613          surf%pt_surface(m) = pt(k+k_off,j,i)
1614
1615       ENDDO
1616
1617    END SUBROUTINE calc_pt_surface
1618
1619!
1620!-- Calculate virtual potential temperature at first grid level
1621!-- ( only for upward-facing surfs )
1622    SUBROUTINE calc_vpt_surface
1623
1624       IMPLICIT NONE
1625
1626       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1627       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1628       
1629       k_off = surf%koff
1630       !$OMP PARALLEL DO PRIVATE( i, j, k )
1631       DO  m = 1, surf%ns
1632
1633          i   = surf%i(m)           
1634          j   = surf%j(m)
1635          k   = surf%k(m)
1636
1637          surf%vpt_surface(m) = vpt(k+k_off,j,i)
1638
1639       ENDDO
1640
1641    END SUBROUTINE calc_vpt_surface
1642   
1643!
1644!-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
1645    SUBROUTINE calc_scaling_parameters
1646
1647       IMPLICIT NONE
1648
1649
1650       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
1651       INTEGER(iwp)  ::  lsp     !< running index for chemical species
1652!
1653!--    Compute theta* at horizontal surfaces
1654       IF ( constant_heatflux  .AND.  .NOT. surf_vertical )  THEN
1655!
1656!--       For a given heat flux in the surface layer:
1657
1658          !$OMP PARALLEL DO PRIVATE( i, j, k )
1659          DO  m = 1, surf%ns
1660
1661             i   = surf%i(m)           
1662             j   = surf%j(m)
1663             k   = surf%k(m)
1664
1665             surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) /                    &
1666                                  ( surf%us(m) + 1E-30_wp )
1667
1668!
1669!--          ts must be limited, because otherwise overflow may occur in case
1670!--          of us=0 when computing ol further below
1671             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
1672             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
1673
1674          ENDDO
1675
1676       ELSEIF ( .NOT. surf_vertical ) THEN
1677!
1678!--       For a given surface temperature:
1679          IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1680
1681             !$OMP PARALLEL DO PRIVATE( i, j, k )
1682             DO  m = 1, surf%ns
1683                i   = surf%i(m)           
1684                j   = surf%j(m)
1685                k   = surf%k(m)
1686
1687                pt(k-1,j,i) = pt_surface
1688             ENDDO
1689          ENDIF
1690
1691          !$OMP PARALLEL DO PRIVATE( z_mo )
1692          DO  m = 1, surf%ns   
1693
1694             z_mo = surf%z_mo(m)
1695
1696             surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) )      &
1697                                  / ( LOG( z_mo / surf%z0h(m) )             &
1698                                      - psi_h( z_mo / surf%ol(m) )          &
1699                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
1700
1701          ENDDO
1702
1703       ENDIF
1704!
1705!--    Compute theta* at vertical surfaces. This is only required in case of
1706!--    land-surface model, in order to compute aerodynamical resistance.
1707       IF ( surf_vertical )  THEN
1708          !$OMP PARALLEL DO PRIVATE( i, j )
1709          DO  m = 1, surf%ns
1710
1711             i          =  surf%i(m)           
1712             j          =  surf%j(m)
1713             surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )
1714!
1715!--          ts must be limited, because otherwise overflow may occur in case
1716!--          of us=0 when computing ol further below
1717             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
1718             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
1719
1720          ENDDO
1721       ENDIF
1722
1723!
1724!--    If required compute q* at horizontal surfaces
1725       IF ( humidity )  THEN
1726          IF ( constant_waterflux  .AND.  .NOT. surf_vertical )  THEN
1727!
1728!--          For a given water flux in the surface layer
1729             !$OMP PARALLEL DO PRIVATE( i, j, k )
1730             DO  m = 1, surf%ns 
1731
1732                i   = surf%i(m)           
1733                j   = surf%j(m)
1734                k   = surf%k(m)
1735                surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) /                &
1736                                               ( surf%us(m) + 1E-30_wp )
1737
1738             ENDDO
1739
1740          ELSEIF ( .NOT. surf_vertical )  THEN
1741             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
1742                             run_coupled )
1743
1744             IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1745                !$OMP PARALLEL DO PRIVATE( i, j, k )
1746                DO  m = 1, surf%ns 
1747
1748                   i   = surf%i(m)           
1749                   j   = surf%j(m)
1750                   k   = surf%k(m)
1751                   q(k-1,j,i) = q_surface
1752                   
1753                ENDDO
1754             ENDIF
1755
1756!
1757!--          Assume saturation for atmosphere coupled to ocean (but not
1758!--          in case of precursor runs)
1759             IF ( coupled_run )  THEN
1760                !$OMP PARALLEL DO PRIVATE( i, j, k, e_s )
1761                DO  m = 1, surf%ns   
1762                   i   = surf%i(m)           
1763                   j   = surf%j(m)
1764                   k   = surf%k(m)
1765                   e_s = 6.1_wp * &
1766                              EXP( 0.07_wp * ( MIN(pt(k-1,j,i),pt(k,j,i))      &
1767                                               - 273.15_wp ) )
1768                   q(k-1,j,i) = 0.622_wp * e_s / ( surface_pressure - e_s )
1769                ENDDO
1770             ENDIF
1771
1772             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
1773               !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1774                DO  m = 1, surf%ns   
1775
1776                   i   = surf%i(m)           
1777                   j   = surf%j(m)
1778                   k   = surf%k(m)
1779   
1780                   z_mo = surf%z_mo(m)
1781
1782                   surf%qs(m) = kappa * ( surf%qv1(m) - q(k-1,j,i) )           &
1783                                        / ( LOG( z_mo / surf%z0q(m) )          &
1784                                            - psi_h( z_mo / surf%ol(m) )       &
1785                                            + psi_h( surf%z0q(m) /             &
1786                                                     surf%ol(m) ) )
1787                ENDDO
1788             ELSE
1789                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1790                DO  m = 1, surf%ns   
1791
1792                   i   = surf%i(m)           
1793                   j   = surf%j(m)
1794                   k   = surf%k(m)
1795   
1796                   z_mo = surf%z_mo(m)
1797
1798                   surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) )              &
1799                                        / ( LOG( z_mo / surf%z0q(m) )          &
1800                                            - psi_h( z_mo / surf%ol(m) )       &
1801                                            + psi_h( surf%z0q(m) /             &
1802                                                     surf%ol(m) ) )
1803                ENDDO
1804             ENDIF
1805          ENDIF
1806!
1807!--       Compute q* at vertical surfaces
1808          IF ( surf_vertical )  THEN
1809             !$OMP PARALLEL DO PRIVATE( i, j )
1810             DO  m = 1, surf%ns 
1811
1812                i          =  surf%i(m)           
1813                j          =  surf%j(m)
1814                surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )
1815
1816             ENDDO
1817          ENDIF
1818       ENDIF
1819       
1820!
1821!--    If required compute s*
1822       IF ( passive_scalar )  THEN
1823!
1824!--       At horizontal surfaces
1825          IF ( constant_scalarflux  .AND.  .NOT. surf_vertical )  THEN
1826!
1827!--          For a given scalar flux in the surface layer
1828             !$OMP PARALLEL DO PRIVATE( i, j )
1829             DO  m = 1, surf%ns 
1830                i   = surf%i(m)           
1831                j   = surf%j(m)
1832                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1833             ENDDO
1834          ENDIF
1835!
1836!--       At vertical surfaces
1837          IF ( surf_vertical )  THEN
1838             !$OMP PARALLEL DO PRIVATE( i, j )
1839             DO  m = 1, surf%ns 
1840                i          =  surf%i(m)           
1841                j          =  surf%j(m)
1842                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1843             ENDDO
1844          ENDIF
1845       ENDIF
1846
1847!
1848!--    If required compute cs* (chemical species)
1849       IF ( air_chemistry  )  THEN 
1850!
1851!--       At horizontal surfaces                             
1852          DO  lsp = 1, nvar
1853             IF ( constant_csflux(lsp)  .AND.  .NOT.  surf_vertical )  THEN
1854!--             For a given chemical species' flux in the surface layer
1855                !$OMP PARALLEL DO PRIVATE( i, j )
1856                DO  m = 1, surf%ns 
1857                   i   = surf%i(m)           
1858                   j   = surf%j(m)
1859                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1860                ENDDO
1861             ENDIF
1862          ENDDO
1863!
1864!--       At vertical surfaces
1865          IF ( surf_vertical )  THEN
1866             DO  lsp = 1, nvar
1867                !$OMP PARALLEL DO PRIVATE( i, j )
1868                DO  m = 1, surf%ns 
1869                   i   = surf%i(m)           
1870                   j   = surf%j(m)
1871                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1872                ENDDO
1873             ENDDO
1874          ENDIF
1875       ENDIF
1876
1877!
1878!--    If required compute qc* and nc*
1879       IF ( cloud_physics  .AND.  microphysics_morrison  .AND.                 &
1880            .NOT. surf_vertical )  THEN
1881          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1882          DO  m = 1, surf%ns   
1883
1884             i    = surf%i(m)           
1885             j    = surf%j(m)
1886             k    = surf%k(m)
1887
1888             z_mo = surf%z_mo(m)
1889
1890             surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) )                 &
1891                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1892                                 - psi_h( z_mo / surf%ol(m) )                  &
1893                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1894
1895             surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) )                 &
1896                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1897                                 - psi_h( z_mo / surf%ol(m) )                  &
1898                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1899          ENDDO
1900
1901       ENDIF
1902
1903!
1904!--    If required compute qr* and nr*
1905       IF ( cloud_physics  .AND.  microphysics_seifert  .AND.                  &
1906            .NOT. surf_vertical )  THEN
1907          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1908          DO  m = 1, surf%ns   
1909
1910             i    = surf%i(m)           
1911             j    = surf%j(m)
1912             k    = surf%k(m)
1913
1914             z_mo = surf%z_mo(m)
1915
1916             surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) )                 &
1917                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1918                                 - psi_h( z_mo / surf%ol(m) )                  &
1919                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1920
1921             surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) )                 &
1922                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1923                                 - psi_h( z_mo / surf%ol(m) )                  &
1924                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1925          ENDDO
1926
1927       ENDIF
1928
1929    END SUBROUTINE calc_scaling_parameters
1930
1931
1932
1933!
1934!-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
1935    SUBROUTINE calc_surface_fluxes
1936
1937       IMPLICIT NONE
1938
1939       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
1940       INTEGER(iwp)  ::  lsp     !< running index for chemical species
1941
1942       REAL(wp)                            ::  dum     !< dummy to precalculate logarithm
1943       REAL(wp)                            ::  flag_u  !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
1944       REAL(wp)                            ::  flag_v  !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
1945       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_i     !< u-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1946       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_i     !< v-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1947       REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_i     !< w-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1948
1949!
1950!--    Calcuate surface fluxes at horizontal walls
1951       IF ( .NOT. surf_vertical )  THEN
1952!
1953!--       Compute u'w' for the total model domain at upward-facing surfaces.
1954!--       First compute the corresponding component of u* and square it.
1955          IF ( .NOT. downward )  THEN
1956             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1957             DO  m = 1, surf%ns 
1958   
1959                i = surf%i(m)           
1960                j = surf%j(m)
1961                k = surf%k(m)
1962
1963                z_mo = surf%z_mo(m)
1964
1965                surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) )               &
1966                              / ( LOG( z_mo / surf%z0(m) )                     &
1967                                  - psi_m( z_mo / surf%ol(m) )                 &
1968                                  + psi_m( surf%z0(m) / surf%ol(m) ) )
1969!
1970!--             Please note, the computation of usws is not fully accurate. Actually
1971!--             a further interpolation of us onto the u-grid, where usws is defined,
1972!--             is required. However, this is not done as this would require several
1973!--             data transfers between 2D-grid and the surf-type.
1974!--             The impact of the missing interpolation is negligible as several
1975!--             tests had shown.
1976!--             Same also for ol. 
1977                surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)
1978
1979             ENDDO
1980!
1981!--       At downward-facing surfaces
1982          ELSE
1983             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1984             DO  m = 1, surf%ns 
1985   
1986                i = surf%i(m)           
1987                j = surf%j(m)
1988                k = surf%k(m)
1989
1990                z_mo = surf%z_mo(m)
1991
1992                surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
1993                surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
1994
1995             ENDDO     
1996          ENDIF
1997
1998!
1999!--       Compute v'w' for the total model domain.
2000!--       First compute the corresponding component of u* and square it.
2001!--       Upward-facing surfaces
2002          IF ( .NOT. downward )  THEN
2003             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2004             DO  m = 1, surf%ns 
2005                i = surf%i(m)           
2006                j = surf%j(m)
2007                k = surf%k(m)
2008
2009                z_mo = surf%z_mo(m)
2010
2011                surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) )               &
2012                           / ( LOG( z_mo / surf%z0(m) )                        &
2013                               - psi_m( z_mo / surf%ol(m) )                    &
2014                               + psi_m( surf%z0(m) / surf%ol(m) ) )
2015!
2016!--             Please note, the computation of vsws is not fully accurate. Actually
2017!--             a further interpolation of us onto the v-grid, where vsws is defined,
2018!--             is required. However, this is not done as this would require several
2019!--             data transfers between 2D-grid and the surf-type.
2020!--             The impact of the missing interpolation is negligible as several
2021!--             tests had shown.
2022!--             Same also for ol. 
2023                surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)
2024             ENDDO
2025!
2026!--       Downward-facing surfaces
2027          ELSE
2028             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2029             DO  m = 1, surf%ns 
2030                i = surf%i(m)           
2031                j = surf%j(m)
2032                k = surf%k(m)
2033
2034                z_mo = surf%z_mo(m)
2035
2036                surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
2037                surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
2038             ENDDO
2039          ENDIF
2040!
2041!--       Compute the vertical kinematic heat flux
2042          IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point&
2043               <=  skip_time_do_lsm  .AND. simulated_time > 0.0_wp ) .OR.      &
2044               .NOT.  land_surface )  .AND.  .NOT. urban_surface  .AND.        &
2045               .NOT. downward )  THEN
2046             !$OMP PARALLEL DO PRIVATE( i, j, k )
2047             DO  m = 1, surf%ns
2048                i    = surf%i(m)           
2049                j    = surf%j(m)
2050                k    = surf%k(m)
2051                surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
2052             ENDDO
2053          ENDIF
2054!
2055!--       Compute the vertical water flux
2056          IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
2057               ( ( time_since_reference_point <= skip_time_do_lsm  .AND.       &
2058               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface  )  .AND.   &
2059               .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
2060             !$OMP PARALLEL DO PRIVATE( i, j, k )
2061             DO  m = 1, surf%ns
2062                i    = surf%i(m)           
2063                j    = surf%j(m)
2064                k    = surf%k(m)
2065                surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
2066             ENDDO
2067          ENDIF
2068!
2069!--       Compute the vertical scalar flux
2070          IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.       &
2071                .NOT.  downward )  THEN
2072             !$OMP PARALLEL DO PRIVATE( i, j )
2073             DO  m = 1, surf%ns   
2074
2075                i    = surf%i(m)           
2076                j    = surf%j(m)
2077                surf%ssws(m) = -surf%ss(m) * surf%us(m)
2078
2079             ENDDO
2080          ENDIF   
2081!
2082!--       Compute the vertical chemical species' flux
2083          DO  lsp = 1, nvar
2084             IF (  .NOT.  constant_csflux(lsp)  .AND.  air_chemistry  .AND.    &
2085                   .NOT.  downward )  THEN
2086                !$OMP PARALLEL DO PRIVATE( i, j )
2087                DO  m = 1, surf%ns   
2088                   i     = surf%i(m)           
2089                   j     = surf%j(m)
2090                   surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
2091                ENDDO
2092             ENDIF   
2093          ENDDO
2094
2095!
2096!--       Compute (turbulent) fluxes of cloud water content and cloud drop conc.
2097          IF ( cloud_physics  .AND.  microphysics_morrison  .AND.              &
2098               .NOT. downward)  THEN
2099             !$OMP PARALLEL DO PRIVATE( i, j )
2100             DO  m = 1, surf%ns   
2101
2102                i    = surf%i(m)           
2103                j    = surf%j(m)
2104
2105                surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
2106                surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
2107             ENDDO
2108          ENDIF   
2109!
2110!--       Compute (turbulent) fluxes of rain water content and rain drop conc.
2111          IF ( cloud_physics  .AND.  microphysics_seifert  .AND.               &
2112               .NOT. downward)  THEN
2113             !$OMP PARALLEL DO PRIVATE( i, j )
2114             DO  m = 1, surf%ns   
2115
2116                i    = surf%i(m)           
2117                j    = surf%j(m)
2118
2119                surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
2120                surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
2121             ENDDO
2122          ENDIF
2123
2124!
2125!--       Bottom boundary condition for the TKE.
2126          IF ( ibc_e_b == 2 )  THEN
2127             !$OMP PARALLEL DO PRIVATE( i, j, k )
2128             DO  m = 1, surf%ns   
2129
2130                i    = surf%i(m)           
2131                j    = surf%j(m)
2132                k    = surf%k(m)
2133
2134                e(k,j,i) = ( surf%us(m) / 0.1_wp )**2
2135!
2136!--             As a test: cm = 0.4
2137!               e(k,j,i) = ( us(j,i) / 0.4_wp )**2
2138                e(k-1,j,i)   = e(k,j,i)
2139
2140             ENDDO
2141          ENDIF
2142!
2143!--    Calcuate surface fluxes at vertical surfaces. No stability is considered.
2144       ELSE
2145!
2146!--       Compute usvs l={0,1} and vsus l={2,3}
2147          IF ( mom_uv )  THEN
2148!
2149!--          Generalize computation by introducing flags. At north- and south-
2150!--          facing surfaces u-component is used, at east- and west-facing
2151!--          surfaces v-component is used.
2152             flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0  .OR.  l == 1 )   
2153             flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2  .OR.  l == 3 )   
2154             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
2155             DO  m = 1, surf%ns 
2156                i = surf%i(m)           
2157                j = surf%j(m)
2158                k = surf%k(m)
2159
2160                z_mo = surf%z_mo(m)
2161
2162                surf%mom_flux_uv(m) = kappa *                                  &
2163                                ( flag_u * u(k,j,i) + flag_v * v(k,j,i) )  /   &
2164                                                        LOG( z_mo / surf%z0(m) )
2165
2166               surf%mom_flux_uv(m) =                                           &
2167                                    - surf%mom_flux_uv(m) * surf%us(m)
2168             ENDDO
2169          ENDIF
2170!
2171!--       Compute wsus l={0,1} and wsvs l={2,3}
2172          IF ( mom_w )  THEN
2173             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
2174             DO  m = 1, surf%ns 
2175                i = surf%i(m)           
2176                j = surf%j(m)
2177                k = surf%k(m)
2178
2179                z_mo = surf%z_mo(m)
2180
2181                surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
2182
2183                surf%mom_flux_w(m) =                                           &
2184                                     - surf%mom_flux_w(m) * surf%us(m)
2185             ENDDO
2186          ENDIF
2187!
2188!--       Compute momentum fluxes used for subgrid-scale TKE production at
2189!--       vertical surfaces. In constrast to the calculated momentum fluxes at
2190!--       vertical surfaces before, which are defined on the u/v/w-grid,
2191!--       respectively), the TKE fluxes are defined at the scalar grid.
2192!--       
2193          IF ( mom_tke )  THEN
2194!
2195!--          Precalculate velocity components at scalar grid point.
2196             ALLOCATE( u_i(1:surf%ns) )
2197             ALLOCATE( v_i(1:surf%ns) )
2198             ALLOCATE( w_i(1:surf%ns) )
2199
2200             IF ( l == 0  .OR.  l == 1 )  THEN
2201                !$OMP PARALLEL DO PRIVATE( i, j, k )
2202                DO  m = 1, surf%ns 
2203                   i = surf%i(m)           
2204                   j = surf%j(m)
2205                   k = surf%k(m)
2206
2207                   u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
2208                   v_i(m) = 0.0_wp
2209                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
2210                ENDDO
2211             ELSE
2212                !$OMP PARALLEL DO PRIVATE( i, j, k )
2213                DO  m = 1, surf%ns 
2214                   i = surf%i(m)           
2215                   j = surf%j(m)
2216                   k = surf%k(m)
2217
2218                   u_i(m) = 0.0_wp
2219                   v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
2220                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
2221                ENDDO
2222             ENDIF
2223
2224             !$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )
2225             DO  m = 1, surf%ns 
2226                i = surf%i(m)           
2227                j = surf%j(m)
2228
2229                z_mo = surf%z_mo(m)
2230
2231                dum = kappa / LOG( z_mo / surf%z0(m) )
2232!
2233!--             usvs (l=0,1) and vsus (l=2,3)
2234                surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )
2235!
2236!--             wsvs (l=0,1) and wsus (l=2,3)
2237                surf%mom_flux_tke(1,m) = dum * w_i(m)
2238
2239                surf%mom_flux_tke(0:1,m) =                                     &
2240                               - surf%mom_flux_tke(0:1,m) * surf%us(m)
2241             ENDDO
2242!
2243!--          Deallocate temporary arrays
2244             DEALLOCATE( u_i )             
2245             DEALLOCATE( v_i ) 
2246             DEALLOCATE( w_i ) 
2247          ENDIF
2248       ENDIF
2249
2250    END SUBROUTINE calc_surface_fluxes
2251
2252
2253!
2254!-- Integrated stability function for momentum
2255    FUNCTION psi_m( zeta ) 
2256       
2257       USE kinds
2258
2259       IMPLICIT NONE
2260
2261       REAL(wp)            :: psi_m !< Integrated similarity function result
2262       REAL(wp)            :: zeta  !< Stability parameter z/L
2263       REAL(wp)            :: x     !< dummy variable
2264
2265       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
2266       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
2267       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
2268       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
2269       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
2270       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
2271
2272
2273       IF ( zeta < 0.0_wp )  THEN
2274          x = SQRT( SQRT( 1.0_wp  - 16.0_wp * zeta ) )
2275          psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
2276                  * ( 1.0_wp + x**2 ) * 0.125_wp )
2277       ELSE
2278
2279          psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta         &
2280                   - bc_d_d
2281!
2282!--       Old version for stable conditions (only valid for z/L < 0.5)
2283!--       psi_m = - 5.0_wp * zeta
2284
2285       ENDIF
2286
2287    END FUNCTION psi_m
2288
2289
2290!
2291!-- Integrated stability function for heat and moisture
2292    FUNCTION psi_h( zeta ) 
2293       
2294       USE kinds
2295
2296       IMPLICIT NONE
2297
2298       REAL(wp)            :: psi_h !< Integrated similarity function result
2299       REAL(wp)            :: zeta  !< Stability parameter z/L
2300       REAL(wp)            :: x     !< dummy variable
2301
2302       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
2303       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
2304       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
2305       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
2306       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
2307       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
2308
2309
2310       IF ( zeta < 0.0_wp )  THEN
2311          x = SQRT( 1.0_wp  - 16.0_wp * zeta )
2312          psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
2313       ELSE
2314          psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp          &
2315                  + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d             &
2316                  + 1.0_wp
2317!
2318!--       Old version for stable conditions (only valid for z/L < 0.5)
2319!--       psi_h = - 5.0_wp * zeta
2320       ENDIF
2321
2322    END FUNCTION psi_h
2323
2324
2325!------------------------------------------------------------------------------!
2326! Description:
2327! ------------
2328!> Calculates stability function for momentum
2329!>
2330!> @author Hauke Wurps
2331!------------------------------------------------------------------------------!
2332    FUNCTION phi_m( zeta ) 
2333   
2334       IMPLICIT NONE
2335   
2336       REAL(wp)            :: phi_m         !< Value of the function
2337       REAL(wp)            :: zeta          !< Stability parameter z/L
2338   
2339       REAL(wp), PARAMETER :: a = 16.0_wp   !< constant
2340       REAL(wp), PARAMETER :: c = 5.0_wp    !< constant
2341   
2342       IF ( zeta < 0.0_wp )  THEN
2343          phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) )
2344       ELSE
2345          phi_m = 1.0_wp + c * zeta
2346       ENDIF
2347   
2348    END FUNCTION phi_m
2349
2350 END MODULE surface_layer_fluxes_mod
Note: See TracBrowser for help on using the repository browser.