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

Last change on this file since 4745 was 4717, checked in by pavelkrc, 4 years ago

Fixes and optimizations of OpenMP parallelization, formatting of OpenMP directives

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