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

Last change on this file since 4198 was 4186, checked in by suehring, 5 years ago

Enable limitation of Obukhov length so that it does not become zero; to read input data from netcdf in init_3d_model use provided module variables instead of defining local ones; tests updated because changes in Obukhov lengths causes small differences during the initial phase of the run

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