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

Last change on this file since 3130 was 3130, checked in by gronemeier, 3 years ago

calculate km according to MOST within the surface layer

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