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

Last change on this file since 3834 was 3834, checked in by forkel, 2 years ago

removed forgotten ,nvar from USE chem_modules

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