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

Last change on this file since 3146 was 3146, checked in by maronga, 3 years ago

major bugfix in calculation of Obukhov length and subsequent handling of surface information

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