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

Last change on this file since 3046 was 3045, checked in by Giersch, 3 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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