Changeset 4562 for palm


Ignore:
Timestamp:
Jun 12, 2020 8:38:47 AM (6 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/check_parameters.f90

    r4559 r4562  
    2424! -----------------
    2525! $Id$
     26! bugfix: revised error message for exceeding allow number of time series
     27!
     28! 4559 2020-06-11 08:51:48Z raasch
    2629! file re-formatted to follow the PALM coding standard
    2730!
     
    671674       WRITE( message_string, * ) 'number of time series quantities exceeds',                      &
    672675                                  ' its maximum of dots_max = ', dots_max,                         &
    673                                   '&Please increase dots_max in modules.f90.'
    674        CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )
     676                                  '&Please increase dots_max in netcdf_interface_mod.f90.'
     677       CALL message( 'check_parameters', 'PA0194', 1, 2, 0, 6, 0 )
    675678    ENDIF
    676679
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

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

    r4559 r4562  
    2525! -----------------
    2626! $Id$
     27! Parts of r4559 re-formatted
     28!
     29! 4559 2020-06-11 08:51:48Z raasch
    2730! File re-formatted to follow the PALM coding standard
    2831!
     
    13751378!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
    13761379!--             Additional factors are added to improve the variance of v and w
    1377                 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) * MERGE( 1.0_wp, 0.0_wp,       &
    1378                                  BTEST( wall_flags_total_0(k,j,i), 1 ) )
     1380                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) *                              &
     1381                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
    13791382             ENDDO
    13801383          ENDDO
     
    13931396                dist_yz(k,j,3) = MIN( ( SQRT( a33(k) / MAXVAL( a33 ) ) * 1.3_wp ) *                &
    13941397                                      ( a31(k) * fu_yz(k,j) + a32(k) * fv_yz(k,j) + a33(k)         &
    1395                                         * fw_yz(k,j) ), 3.0_wp ) *  MERGE( 1.0_wp, 0.0_wp,         &
    1396                                         BTEST( wall_flags_total_0(k,j,i), 3 ) )
     1398                                        * fw_yz(k,j) ), 3.0_wp ) *                                 &
     1399                                      MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
    13971400             ENDDO
    13981401          ENDDO
     
    14121415          IF ( myidx == id_stg_right )  i = nxr+1
    14131416
    1414           mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1) * MERGE( 1.0_wp, 0.0_wp,                &
    1415                            BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 1 ) ) )
     1417          mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1) *                                       &
     1418                       MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 1 ) ) )
    14161419
    14171420          IF ( myidx == id_stg_left  )  i = nxl-1
    14181421          IF ( myidx == id_stg_right )  i = nxr+1
    14191422
    1420           mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2) * MERGE( 1.0_wp, 0.0_wp,               &
    1421                            BTEST( wall_flags_total_0(nzb:nzt,nysv:nyn,i), 2 ) ) )
    1422           mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3) * MERGE( 1.0_wp, 0.0_wp,                &
    1423                            BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 3 ) ) )
     1423          mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2) *                                      &
     1424                       MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,nysv:nyn,i), 2 ) ) )
     1425          mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3) *                                       &
     1426                       MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 3 ) ) )
    14241427
    14251428#if defined( __parallel )
     
    14351438          IF ( myidx == id_stg_right )  i = nxr+1
    14361439
    1437           dist_yz(:,nys:nyn,1) = ( dist_yz(:,nys:nyn,1) - mc_factor(1) ) * MERGE( 1.0_wp, 0.0_wp,  &
    1438                                 BTEST( wall_flags_total_0(:,nys:nyn,i), 1 ) )
     1440          dist_yz(:,nys:nyn,1) = ( dist_yz(:,nys:nyn,1) - mc_factor(1) ) *                         &
     1441                                MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,nys:nyn,i), 1 ) )
    14391442
    14401443
     
    14421445          IF ( myidx == id_stg_right )  i = nxr+1
    14431446
    1444           dist_yz(:,nys:nyn,2) = ( dist_yz(:,nys:nyn,2) - mc_factor(2) ) * MERGE( 1.0_wp, 0.0_wp,  &
    1445                                 BTEST( wall_flags_total_0(:,nys:nyn,i), 2 ) )
    1446 
    1447           dist_yz(:,nys:nyn,3) = ( dist_yz(:,nys:nyn,3) - mc_factor(3) ) * MERGE( 1.0_wp, 0.0_wp,  &
    1448                                 BTEST( wall_flags_total_0(:,nys:nyn,i), 3 ) )
     1447          dist_yz(:,nys:nyn,2) = ( dist_yz(:,nys:nyn,2) - mc_factor(2) ) *                         &
     1448                                MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,nys:nyn,i), 2 ) )
     1449
     1450          dist_yz(:,nys:nyn,3) = ( dist_yz(:,nys:nyn,3) - mc_factor(3) ) *                         &
     1451                                MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,nys:nyn,i), 3 ) )
    14491452!
    14501453!--       Add disturbances
     
    14591462                DO  j = nys, nyn
    14601463                   DO  k = nzb, nzt+1
    1461                       u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )            &
    1462                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,0), 1 ) )
    1463                       v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + dist_yz(k,j,2) )             &
    1464                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 2 ) )
    1465                       w(k,j,-nbgp:-1) =  dist_yz(k,j,3) * MERGE( 1.0_wp, 0.0_wp,                   &
    1466                                          BTEST( wall_flags_total_0(k,j,-1), 3 ) )
     1464                      u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) ) *          &
     1465                                    MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,0), 1 ) )
     1466                      v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) + dist_yz(k,j,2) ) *          &
     1467                                     MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 2 ) )
     1468                      w(k,j,-nbgp:-1)  = dist_yz(k,j,3) *                                          &
     1469                                     MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 3 ) )
    14671470                   ENDDO
    14681471                ENDDO
     
    14711474                DO  j = nys, nyn
    14721475                   DO  k = nzb+1, nzt
    1473                       u(k,j,0)   = ( u(k,j,0) + dist_yz(k,j,1) ) * MERGE( 1.0_wp, 0.0_wp,          &
    1474                                    BTEST( wall_flags_total_0(k,j,0), 1 ) )
     1476                      u(k,j,0)   = ( u(k,j,0)  + dist_yz(k,j,1) ) *                                &
     1477                                   MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,0), 1 ) )
    14751478                      u(k,j,-1)  = u(k,j,0)
    1476                       v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2) ) * MERGE( 1.0_wp, 0.0_wp,        &
    1477                                    BTEST( wall_flags_total_0(k,j,-1), 2 ) )
    1478                       w(k,j,-1)  = ( w(k,j,-1)  + dist_yz(k,j,3) ) * MERGE( 1.0_wp, 0.0_wp,        &
    1479                                    BTEST( wall_flags_total_0(k,j,-1), 3 ) )
     1479                      v(k,j,-1)  = ( v(k,j,-1) + dist_yz(k,j,2) ) *                                &
     1480                                   MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 2 ) )
     1481                      w(k,j,-1)  = ( w(k,j,-1) + dist_yz(k,j,3) ) *                                &
     1482                                   MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 3 ) )
    14801483                   ENDDO
    14811484                ENDDO
     
    14851488             DO  j = nys, nyn
    14861489                DO  k = nzb+1, nzt
    1487                    u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) ) * MERGE( 1.0_wp, 0.0_wp,       &
    1488                                   BTEST( wall_flags_total_0(k,j,nxr+1), 1 ) )
    1489                    v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) ) * MERGE( 1.0_wp, 0.0_wp,       &
    1490                                     BTEST( wall_flags_total_0(k,j,nxr+1), 2 ) )
    1491                    w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) ) * MERGE( 1.0_wp, 0.0_wp,       &
    1492                                   BTEST( wall_flags_total_0(k,j,nxr+1), 3 ) )
     1490                   u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) ) *                              &
     1491                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 1 ) )
     1492                   v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) ) *                              &
     1493                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 2 ) )
     1494                   w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) ) *                              &
     1495                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 3 ) )
    14931496                ENDDO
    14941497             ENDDO
     
    15561559!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
    15571560!--             Additional factors are added to improve the variance of v and w
    1558                 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) * MERGE( 1.0_wp, 0.0_wp,       &
    1559                                  BTEST( wall_flags_total_0(k,j,i), 1 ) )
     1561                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) *                              &
     1562                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
    15601563
    15611564                dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL( a33 ) ) * 1.3_wp )                   &
    15621565                                 * ( a31(k) * fu_xz(k,i) + a32(k) * fv_xz(k,i) + a33(k)            &
    1563                                  * fw_xz(k,i) ), 3.0_wp ) * MERGE( 1.0_wp, 0.0_wp,                 &
    1564                                  BTEST( wall_flags_total_0(k,j,i), 3 ) )
     1566                                 * fw_xz(k,i) ), 3.0_wp ) *                                        &
     1567                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
    15651568             ENDDO
    15661569          ENDDO
     
    15791582          IF ( myidy == id_stg_north )  j = nyn+1
    15801583
    1581           mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2) * MERGE( 1.0_wp, 0.0_wp,                &
    1582                            BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 2 ) ) )
     1584          mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2) *                                       &
     1585                       MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 2 ) ) )
    15831586
    15841587          IF ( myidy == id_stg_south ) j = nys-1
    15851588          IF ( myidy == id_stg_north ) j = nyn+1
    15861589
    1587           mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1) * MERGE( 1.0_wp, 0.0_wp,               &
    1588                            BTEST( wall_flags_total_0(nzb:nzt,j,nxlu:nxr), 1 ) ) )
    1589           mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3) * MERGE( 1.0_wp, 0.0_wp,                &
    1590                            BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 3 ) ) )
     1590          mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1) *                                      &
     1591                       MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,j,nxlu:nxr), 1 ) ) )
     1592          mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3) *                                       &
     1593                       MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 3 ) ) )
    15911594
    15921595#if defined( __parallel )
     
    16011604          IF ( myidy == id_stg_north ) j = nyn+1
    16021605
    1603           dist_xz(:,nxl:nxr,2) = ( dist_xz(:,nxl:nxr,2) - mc_factor(2) ) * MERGE( 1.0_wp, 0.0_wp,  &
    1604                                 BTEST( wall_flags_total_0(:,j,nxl:nxr), 2 ) )
     1606          dist_xz(:,nxl:nxr,2) = ( dist_xz(:,nxl:nxr,2) - mc_factor(2) ) *                         &
     1607                                MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,j,nxl:nxr), 2 ) )
    16051608
    16061609
     
    16081611          IF ( myidy == id_stg_north ) j = nyn+1
    16091612
    1610           dist_xz(:,nxl:nxr,1) = ( dist_xz(:,nxl:nxr,1) - mc_factor(1) ) * MERGE( 1.0_wp, 0.0_wp,  &
    1611                                 BTEST( wall_flags_total_0(:,j,nxl:nxr), 1 ) )
    1612 
    1613           dist_xz(:,nxl:nxr,3) = ( dist_xz(:,nxl:nxr,3) - mc_factor(3) ) * MERGE( 1.0_wp, 0.0_wp,  &
    1614                                 BTEST( wall_flags_total_0(:,j,nxl:nxr), 3 ) )
     1613          dist_xz(:,nxl:nxr,1) = ( dist_xz(:,nxl:nxr,1) - mc_factor(1) ) *                         &
     1614                                MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,j,nxl:nxr), 1 ) )
     1615
     1616          dist_xz(:,nxl:nxr,3) = ( dist_xz(:,nxl:nxr,3) - mc_factor(3) ) *                         &
     1617                                MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,j,nxl:nxr), 3 ) )
    16151618!
    16161619!--       Add disturbances
     
    16321635             DO  i = nxl, nxr
    16331636                DO  k = nzb+1, nzt
    1634                    u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) * MERGE( 1.0_wp, 0.0_wp,       &
    1635                                   BTEST( wall_flags_total_0(k,nyn+1,i), 1 ) )
    1636                    v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) ) * MERGE( 1.0_wp, 0.0_wp,       &
    1637                                   BTEST( wall_flags_total_0(k,nyn+1,i), 2 ) )
    1638                    w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) ) * MERGE( 1.0_wp, 0.0_wp,       &
    1639                                   BTEST( wall_flags_total_0(k,nyn+1,i), 3 ) )
     1637                   u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) *                              &
     1638                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 1 ) )
     1639                   v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) ) *                              &
     1640                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 2 ) )
     1641                   w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) ) *                              &
     1642                                  MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 3 ) )
    16401643                ENDDO
    16411644             ENDDO
     
    20222025!--    u'u' and v'v'. Assume isotropy. Note, add a small negative number to the denominator, else
    20232026!--    the mergpe-function can crash if scale_l is zero.
    2024        r11(k) = scale_us**2 * ( MERGE( 0.35_wp                                                     &
    2025                 * ABS( - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) )**( 2.0_wp / 3.0_wp ),        &
    2026                 0.0_wp, scale_l < 0.0_wp ) + 5.0_wp - 4.0_wp * zzi ) * blend
     2027       r11(k) = scale_us**2 * ( MERGE( 0.35_wp *                                                   &
     2028                         ABS( - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) )**( 2.0_wp / 3.0_wp ), &
     2029                         0.0_wp, scale_l < 0.0_wp ) + 5.0_wp - 4.0_wp * zzi                                 &
     2030                              ) * blend
    20272031
    20282032       r22(k) = r11(k)
Note: See TracChangeset for help on using the changeset viewer.