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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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