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

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

Clean-up unnecessary index access to surface type

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