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

Last change on this file since 3668 was 3668, checked in by maronga, 5 years ago

removed most_methods circular and lookup. added improved version of palm_csd

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