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

Last change on this file since 2716 was 2716, checked in by kanani, 3 years ago

Correction of "Former revisions" section

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