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

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

Some interface calls moved to module_interface + cleanup

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