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

Last change on this file since 3634 was 3634, checked in by knoop, 5 years ago

OpenACC port for SPEC

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