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

Last change on this file since 2233 was 2233, checked in by suehring, 4 years ago

last commit documented

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