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

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

Removal of chem directive, plus minor changes

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