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

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

Added missing OpenMP directives within "disturb_field", "surface_layer_fluxes" and "timestep"

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