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

Last change on this file since 4089 was 3987, checked in by kanani, 5 years ago

clean up location, debug and error messages

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