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

Last change on this file since 3149 was 3149, checked in by maronga, 6 years ago

correct accidental commit

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