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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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