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

Last change on this file since 3881 was 3881, checked in by suehring, 5 years ago

Make level 3 initialization of urban-surfaces consistent to input data standard; revise flagging of ground-floor level surfaces at buidlings; bugfixes in level 3 initialization of albedo; bugfix in OpenMP directive; check for zero output timestep in surface output; assure that Obukhov lenght does not become zero

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