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

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

document last commit

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