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

Last change on this file since 3744 was 3744, checked in by suehring, 2 years ago

Coupling of indoor model to atmosphere; output of indoor temperatures and waste heat; enable restarts with indoor model; bugfix plant transpiration; bugfix - missing calculation of 10cm temperature

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