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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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