Ignore:
Timestamp:
May 30, 2017 5:47:52 PM (7 years ago)
Author:
suehring
Message:

Adjustments according new topography and surface-modelling concept implemented

File:
1 edited

Legend:

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

    r2119 r2232  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Adjustments to new surface concept
     23! OpenMP bugfix
    2324!
    2425! Former revisions:
     
    171172
    172173    USE arrays_3d,                                                             &
    173         ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
    174                s, shf, ss, ssws, ts, u, us, usws, v, vpt, vsws, zu, zw, z0,    &
    175                z0h, z0q, drho_air_zw, rho_air_zw
     174        ONLY:  e, kh, nr, pt, q, ql, qr, s, u, v, vpt, w, zu, zw, drho_air_zw, &
     175               rho_air_zw
    176176
    177177    USE cloud_parameters,                                                      &
     
    187187               constant_waterflux, coupling_mode, g, humidity, ibc_e_b,        &
    188188               ibc_pt_b, initializing_actions, kappa,                          &
    189                intermediate_timestep_count,                                    &
    190                intermediate_timestep_count_max, large_scale_forcing, lsf_surf, &
     189               intermediate_timestep_count, intermediate_timestep_count_max,   &
     190               land_surface, large_scale_forcing, lsf_surf,                    &
    191191               message_string, microphysics_seifert, most_method, neutral,     &
    192192               passive_scalar, pt_surface, q_surface, run_coupled,             &
    193193               surface_pressure, simulated_time, terminate_run,                &
    194                urban_surface,                                                  &
    195                zeta_max, zeta_min
     194               urban_surface, zeta_max, zeta_min
     195
     196    USE grid_variables,                                                        &
     197        ONLY:  dx, dy 
    196198
    197199    USE indices,                                                               &
    198         ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner,        &
    199                nzb_u_inner, nzb_v_inner
     200        ONLY:  nxl, nxr, nys, nyn, nzb
    200201
    201202    USE kinds
     
    204205
    205206    USE land_surface_model_mod,                                                &
    206         ONLY:  land_surface, skip_time_do_lsm
    207 
     207        ONLY:  aero_resist_kray, skip_time_do_lsm
     208
     209    USE surface_mod,                                                           &
     210        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_type,     &
     211                surf_usm_h, surf_usm_v
    208212       
    209213
     
    213217    INTEGER(iwp) ::  j              !< loop index y direction
    214218    INTEGER(iwp) ::  k              !< loop index z direction
    215     INTEGER(iwp) ::  l_bnd  = 7500  !< Lookup table index of the last time step
    216 
    217     INTEGER(iwp), PARAMETER     :: num_steps = 15000  !< number of steps in the lookup table
    218 
    219     LOGICAL      ::  coupled_run  !< Flag for coupled atmosphere-ocean runs
    220 
    221     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt1,      & !< Potential temperature at first grid level (required for cloud_physics = .T.)
    222                                              qv1,      & !< Specific humidity at first grid level (required for cloud_physics = .T.)
    223                                              uv_total    !< Total velocity at first grid level
     219    INTEGER(iwp) ::  l              !< loop index for surf type
     220    INTEGER(iwp) ::  li_bnd  = 7500 !< Lookup table index of the last time step
     221
     222    INTEGER(iwp), PARAMETER ::  num_steps = 15000  !< number of steps in the lookup table
     223
     224    LOGICAL      ::  coupled_run       !< Flag for coupled atmosphere-ocean runs
     225    LOGICAL      ::  downward = .FALSE.!< Flag indicating downward-facing horizontal surface
     226    LOGICAL      ::  mom_uv  = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces
     227    LOGICAL      ::  mom_w   = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces
     228    LOGICAL      ::  mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production
     229    LOGICAL      ::  surf_vertical     !< Flag indicating vertical surfaces
    224230
    225231    REAL(wp), DIMENSION(0:num_steps-1) :: rib_tab,  & !< Lookup table bulk Richardson number
     
    232238                     z_mo                 !< Height of the constant flux layer where MOST is assumed
    233239
     240    TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
     241
    234242
    235243    SAVE
     
    237245    PRIVATE
    238246
    239     PUBLIC init_surface_layer_fluxes, pt1, qv1, surface_layer_fluxes, uv_total
     247    PUBLIC init_surface_layer_fluxes, surface_layer_fluxes
    240248
    241249    INTERFACE init_surface_layer_fluxes
     
    260268       IMPLICIT NONE
    261269
     270       surf_vertical = .FALSE.
     271       downward      = .FALSE.
    262272!
    263273!--    In case cloud physics is used, it is required to derive potential
     
    265275!--    and q
    266276       IF ( cloud_physics )  THEN
    267           CALL calc_pt_q
     277!
     278!--       First call for horizontal default-type surfaces (l=0 - upward facing,
     279!--       l=1 - downward facing)
     280          DO  l = 0, 1
     281             IF ( surf_def_h(l)%ns >= 1  )  THEN
     282                surf => surf_def_h(l)
     283                CALL calc_pt_q
     284             ENDIF
     285          ENDDO
     286!
     287!--       Call for natural-type surfaces
     288          IF ( surf_lsm_h%ns >= 1  )  THEN
     289             surf => surf_lsm_h
     290             CALL calc_pt_q
     291          ENDIF
     292!
     293!--       Call for urban-type surfaces
     294          IF ( surf_usm_h%ns >= 1  )  THEN
     295             surf => surf_usm_h
     296             CALL calc_pt_q
     297          ENDIF
    268298       ENDIF
    269299
     
    277307!
    278308!--    Depending on setting of most_method use the "old" routine
     309!--    Note, each routine is called for different surface types.
     310!--    First call for default-type horizontal surfaces, for natural- and
     311!--    urban-type surfaces. Note, at this place only upward-facing horizontal
     312!--    surfaces are treted.
    279313       IF ( most_method == 'circular' )  THEN
    280 
    281           CALL calc_scaling_parameters
    282 
    283           CALL calc_uv_total
    284 
    285           IF ( .NOT. neutral )  THEN
    286              CALL calc_ol
    287           ENDIF
    288 
    289           CALL calc_us
    290 
    291           CALL calc_surface_fluxes
    292 
     314!
     315!--       Default-type upward-facing horizontal surfaces
     316          IF ( surf_def_h(0)%ns >= 1  )  THEN
     317             surf => surf_def_h(0)
     318             CALL calc_scaling_parameters
     319             CALL calc_uvw_abs
     320             IF ( .NOT. neutral )  CALL calc_ol
     321             CALL calc_us
     322             CALL calc_surface_fluxes
     323          ENDIF
     324!
     325!--       Natural-type horizontal surfaces
     326          IF ( surf_lsm_h%ns >= 1  )  THEN
     327             surf => surf_lsm_h
     328             CALL calc_scaling_parameters
     329             CALL calc_uvw_abs
     330             IF ( .NOT. neutral )  CALL calc_ol
     331             CALL calc_us
     332             CALL calc_surface_fluxes
     333          ENDIF
     334!
     335!--       Urban-type horizontal surfaces
     336          IF ( surf_usm_h%ns >= 1  )  THEN
     337             surf => surf_usm_h
     338             CALL calc_scaling_parameters
     339             CALL calc_uvw_abs
     340             IF ( .NOT. neutral )  CALL calc_ol
     341             CALL calc_us
     342             CALL calc_surface_fluxes
     343          ENDIF
    293344!
    294345!--    Use either Newton iteration or a lookup table for the bulk Richardson
    295346!--    number to calculate the Obukhov length
    296347       ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
    297 
    298           CALL calc_uv_total
    299 
    300           IF ( .NOT. neutral )  THEN
    301              CALL calc_ol
    302           ENDIF
    303 
     348!
     349!--       Default-type upward-facing horizontal surfaces
     350          IF ( surf_def_h(0)%ns >= 1  )  THEN
     351             surf => surf_def_h(0)
     352             CALL calc_uvw_abs
     353             IF ( .NOT. neutral )  CALL calc_ol
     354             CALL calc_us
     355             CALL calc_scaling_parameters
     356             CALL calc_surface_fluxes
     357          ENDIF
     358!
     359!--       Natural-type horizontal surfaces
     360          IF ( surf_lsm_h%ns >= 1  )  THEN
     361             surf => surf_lsm_h
     362             CALL calc_uvw_abs
     363             IF ( .NOT. neutral )  CALL calc_ol
     364             CALL calc_us
     365             CALL calc_scaling_parameters
     366             CALL calc_surface_fluxes
     367          ENDIF
     368!
     369!--       Urban-type horizontal surfaces
     370          IF ( surf_usm_h%ns >= 1  )  THEN
     371             surf => surf_usm_h
     372             CALL calc_uvw_abs
     373             IF ( .NOT. neutral )  CALL calc_ol
     374             CALL calc_us
     375             CALL calc_scaling_parameters
     376             CALL calc_surface_fluxes
     377          ENDIF
     378
     379       ENDIF
     380!
     381!--    Treat downward-facing horizontal surfaces. Note, so far, these are
     382!--    always default type. Stratification is not considered
     383!--    in this case, hence, no further distinction between different
     384!--    most_method is required. 
     385       IF ( surf_def_h(1)%ns >= 1  )  THEN
     386          downward = .TRUE.
     387          surf => surf_def_h(1)
     388          CALL calc_uvw_abs
    304389          CALL calc_us
    305 
    306           CALL calc_scaling_parameters
    307 
    308390          CALL calc_surface_fluxes
    309 
     391          downward = .FALSE.
    310392       ENDIF
     393!
     394!--    Calculate surfaces fluxes at vertical surfaces for momentum
     395!--    and subgrid-scale TKE.
     396!--    No stability is considered. Therefore, scaling parameters and Obukhov-
     397!--    length do not need to be calculated and no distinction in 'circular',
     398!--    'Newton' or 'lookup' is necessary so far.
     399!--    Note, this will change if stability is once considered.
     400       surf_vertical = .TRUE.
     401!
     402!--    Calculate horizontal momentum fluxes at north- and south-facing
     403!--    surfaces(usvs).
     404!--    For default-type surfaces
     405       mom_uv = .TRUE.
     406       DO  l = 0, 1
     407          IF ( surf_def_v(l)%ns >= 1  )  THEN
     408             surf => surf_def_v(l)
     409!
     410!--          Compute surface-parallel velocity
     411             CALL calc_uvw_abs_v_ugrid
     412!
     413!--          Compute respective friction velocity on staggered grid
     414             CALL calc_us
     415!
     416!--          Compute respective surface fluxes for momentum and TKE
     417             CALL calc_surface_fluxes
     418          ENDIF
     419       ENDDO
     420!
     421!--    For natural-type surfaces. Please note, even though stability is not
     422!--    considered for the calculation of momentum fluxes at vertical surfaces,
     423!--    scaling parameters and Obukov length are calculated nevertheless in this
     424!--    case. This is due to the requirement of ts in parameterization of heat
     425!--    flux in land-surface model in case of aero_resist_kray is not true.
     426       IF ( .NOT. aero_resist_kray )  THEN
     427          IF ( most_method == 'circular' )  THEN
     428             DO  l = 0, 1
     429                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     430                   surf => surf_lsm_v(l)
     431!
     432!--                Compute scaling parameters
     433                   CALL calc_scaling_parameters
     434!
     435!--                Compute surface-parallel velocity
     436                   CALL calc_uvw_abs_v_ugrid
     437!
     438!--                Compute Obukhov length
     439                   IF ( .NOT. neutral )  CALL calc_ol
     440!
     441!--                Compute respective friction velocity on staggered grid
     442                   CALL calc_us
     443!
     444!--                Compute respective surface fluxes for momentum and TKE
     445                   CALL calc_surface_fluxes
     446                ENDIF
     447             ENDDO
     448          ELSE
     449             DO  l = 0, 1
     450                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     451                   surf => surf_lsm_v(l)
     452!
     453!--                Compute surface-parallel velocity
     454                   CALL calc_uvw_abs_v_ugrid
     455!
     456!--                Compute Obukhov length
     457                   IF ( .NOT. neutral )  CALL calc_ol
     458!
     459!--                Compute respective friction velocity on staggered grid
     460                   CALL calc_us
     461!
     462!--                Compute scaling parameters
     463                   CALL calc_scaling_parameters
     464!
     465!--                Compute respective surface fluxes for momentum and TKE
     466                   CALL calc_surface_fluxes
     467                ENDIF
     468             ENDDO
     469          ENDIF
     470!
     471!--    No ts is required, so scaling parameters and Obukhov length do not need
     472!--    to be computed.
     473       ELSE
     474          DO  l = 0, 1
     475             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     476                surf => surf_lsm_v(l)
     477!   
     478!--             Compute surface-parallel velocity
     479                CALL calc_uvw_abs_v_ugrid
     480!
     481!--             Compute respective friction velocity on staggered grid
     482                CALL calc_us
     483!
     484!--             Compute respective surface fluxes for momentum and TKE
     485                CALL calc_surface_fluxes
     486             ENDIF
     487          ENDDO
     488       ENDIF
     489!
     490!--    For urban-type surfaces
     491       DO  l = 0, 1
     492          IF ( surf_usm_v(l)%ns >= 1  )  THEN
     493             surf => surf_usm_v(l)
     494!
     495!--          Compute surface-parallel velocity
     496             CALL calc_uvw_abs_v_ugrid
     497!
     498!--          Compute respective friction velocity on staggered grid
     499             CALL calc_us
     500!
     501!--          Compute respective surface fluxes for momentum and TKE
     502             CALL calc_surface_fluxes
     503          ENDIF
     504       ENDDO
     505!
     506!--    Calculate horizontal momentum fluxes at east- and west-facing
     507!--    surfaces (vsus).
     508!--    For default-type surfaces
     509       DO  l = 2, 3
     510          IF ( surf_def_v(l)%ns >= 1  )  THEN
     511             surf => surf_def_v(l)
     512!
     513!--          Compute surface-parallel velocity
     514             CALL calc_uvw_abs_v_vgrid
     515!
     516!--          Compute respective friction velocity on staggered grid
     517             CALL calc_us
     518!
     519!--          Compute respective surface fluxes for momentum and TKE
     520             CALL calc_surface_fluxes
     521          ENDIF
     522       ENDDO
     523!
     524!--    For natural-type surfaces. Please note, even though stability is not
     525!--    considered for the calculation of momentum fluxes at vertical surfaces,
     526!--    scaling parameters and Obukov length are calculated nevertheless in this
     527!--    case. This is due to the requirement of ts in parameterization of heat
     528!--    flux in land-surface model in case of aero_resist_kray is not true.
     529       IF ( .NOT. aero_resist_kray )  THEN
     530          IF ( most_method == 'circular' )  THEN
     531             DO  l = 2, 3
     532                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     533                   surf => surf_lsm_v(l)
     534!
     535!--                Compute scaling parameters
     536                   CALL calc_scaling_parameters
     537!
     538!--                Compute surface-parallel velocity
     539                   CALL calc_uvw_abs_v_vgrid
     540!
     541!--                Compute Obukhov length
     542                   IF ( .NOT. neutral )  CALL calc_ol
     543!
     544!--                Compute respective friction velocity on staggered grid
     545                   CALL calc_us
     546!
     547!--                Compute respective surface fluxes for momentum and TKE
     548                   CALL calc_surface_fluxes
     549                ENDIF
     550             ENDDO
     551          ELSE
     552             DO  l = 2, 3
     553                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     554                   surf => surf_lsm_v(l)
     555!
     556!--                Compute surface-parallel velocity
     557                   CALL calc_uvw_abs_v_vgrid
     558!
     559!--                Compute Obukhov length
     560                   IF ( .NOT. neutral )  CALL calc_ol
     561!
     562!--                Compute respective friction velocity on staggered grid
     563                   CALL calc_us
     564!
     565!--                Compute scaling parameters
     566                   CALL calc_scaling_parameters
     567!
     568!--                Compute respective surface fluxes for momentum and TKE
     569                   CALL calc_surface_fluxes
     570                ENDIF
     571             ENDDO
     572          ENDIF
     573       ELSE
     574          DO  l = 2, 3
     575             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     576                surf => surf_lsm_v(l)
     577!   
     578!--             Compute surface-parallel velocity
     579                CALL calc_uvw_abs_v_vgrid
     580!
     581!--             Compute respective friction velocity on staggered grid
     582                CALL calc_us
     583!
     584!--             Compute respective surface fluxes for momentum and TKE
     585                CALL calc_surface_fluxes
     586             ENDIF
     587          ENDDO
     588       ENDIF
     589!
     590!--    For urban-type surfaces
     591       DO  l = 2, 3
     592          IF ( surf_usm_v(l)%ns >= 1  )  THEN
     593             surf => surf_usm_v(l)
     594!
     595!--          Compute surface-parallel velocity
     596             CALL calc_uvw_abs_v_vgrid
     597!
     598!--          Compute respective friction velocity on staggered grid
     599             CALL calc_us
     600!
     601!--          Compute respective surface fluxes for momentum and TKE
     602             CALL calc_surface_fluxes
     603          ENDIF
     604       ENDDO
     605       mom_uv = .FALSE.
     606!
     607!--    Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial
     608!--    surfaces.
     609       mom_w = .TRUE.
     610!
     611!--    Default-type surfaces
     612       DO  l = 0, 3
     613          IF ( surf_def_v(l)%ns >= 1  )  THEN
     614             surf => surf_def_v(l)
     615             CALL calc_uvw_abs_v_wgrid
     616             CALL calc_us
     617             CALL calc_surface_fluxes
     618          ENDIF
     619       ENDDO
     620!
     621!--    Natural-type surfaces
     622       DO  l = 0, 3
     623          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     624             surf => surf_lsm_v(l)
     625             CALL calc_uvw_abs_v_wgrid
     626             CALL calc_us
     627             CALL calc_surface_fluxes
     628          ENDIF
     629       ENDDO
     630!
     631!--    Urban-type surfaces
     632       DO  l = 0, 3
     633          IF ( surf_usm_v(l)%ns >= 1  )  THEN
     634             surf => surf_usm_v(l)
     635             CALL calc_uvw_abs_v_wgrid
     636             CALL calc_us
     637             CALL calc_surface_fluxes
     638          ENDIF
     639       ENDDO
     640       mom_w = .FALSE.   
     641!
     642!--    Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical
     643!--    surfaces for TKE production. Note, here, momentum fluxes are defined
     644!--    at grid center and are not staggered as before.
     645       mom_tke = .TRUE.
     646!
     647!--    Default-type surfaces
     648       DO  l = 0, 3
     649          IF ( surf_def_v(l)%ns >= 1  )  THEN
     650             surf => surf_def_v(l)
     651             CALL calc_uvw_abs_v_sgrid
     652             CALL calc_us
     653             CALL calc_surface_fluxes
     654          ENDIF
     655       ENDDO
     656!
     657!--    Natural-type surfaces
     658       DO  l = 0, 3
     659          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     660             surf => surf_lsm_v(l)
     661             CALL calc_uvw_abs_v_sgrid
     662             CALL calc_us
     663             CALL calc_surface_fluxes
     664          ENDIF
     665       ENDDO
     666!
     667!--    Urban-type surfaces
     668       DO  l = 0, 3
     669          IF ( surf_usm_v(l)%ns >= 1  )  THEN
     670             surf => surf_usm_v(l)
     671             CALL calc_uvw_abs_v_sgrid
     672             CALL calc_us
     673             CALL calc_surface_fluxes
     674          ENDIF
     675       ENDDO
     676       mom_tke = .FALSE.
     677 
    311678
    312679    END SUBROUTINE surface_layer_fluxes
     
    324691       IMPLICIT NONE
    325692
    326        INTEGER(iwp) :: l,          & !< Index for loop to create lookup table
     693       INTEGER(iwp) :: li,         & !< Index for loop to create lookup table
    327694                       num_steps_n   !< Number of non-stretched zeta steps
    328695
     
    339706                   z0h_min   = 0.0_wp,   & !< Minimum value of z0h to create table
    340707                   z0_min    = 0.0_wp      !< Minimum value of z0 to create table
    341 !
    342 !--    When cloud physics is used, arrays for storing potential temperature and
    343 !--    specific humidity at first grid level are required
    344        IF ( cloud_physics )  THEN
    345           ALLOCATE ( pt1(nysg:nyng,nxlg:nxrg) )
    346           ALLOCATE ( qv1(nysg:nyng,nxlg:nxrg) )
    347        ENDIF
    348 
    349 !
    350 !--    Allocate field for storing the horizontal velocity
    351        ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) )
     708
     709
    352710
    353711
     
    355713!--    In case of runs with neutral statification, set Obukhov length to a
    356714!--    large value
    357        IF ( neutral ) ol = 1.0E10_wp
     715       IF ( neutral )  THEN
     716          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%ol = 1.0E10_wp
     717          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
     718          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
     719       ENDIF
    358720
    359721       IF ( most_method == 'lookup' )  THEN
     
    361723!
    362724!--       Check for roughness heterogeneity. In that case terminate run and
    363 !--       inform user
    364           IF ( MINVAL( z0h ) /= MAXVAL( z0h )  .OR.                            &
    365                MINVAL( z0  ) /= MAXVAL( z0  ) )  THEN
    366              terminate_run_l = .TRUE.
    367           ENDIF
     725!--       inform user. Check for both, natural and non-natural walls.
     726          IF ( surf_def_h(0)%ns >= 1 )  THEN
     727             IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR. &
     728                  MINVAL( surf_def_h(0)%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
     729                terminate_run_l = .TRUE.
     730             ENDIF
     731          ENDIF
     732          IF ( surf_lsm_h%ns >= 1 )  THEN
     733             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
     734                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
     735                terminate_run_l = .TRUE.
     736             ENDIF
     737          ENDIF
     738          IF ( surf_usm_h%ns >= 1 )  THEN
     739             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h )  .OR.       &
     740                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_usm_h%z0  ) )  THEN
     741                terminate_run_l = .TRUE.
     742             ENDIF
     743          ENDIF
     744!
     745!--       Check roughness homogeneity between differt surface types.
     746          IF ( surf_lsm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
     747             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
     748                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
     749                terminate_run_l = .TRUE.
     750             ENDIF
     751          ENDIF
     752          IF ( surf_usm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
     753             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
     754                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
     755                terminate_run_l = .TRUE.
     756             ENDIF
     757          ENDIF
     758          IF ( surf_usm_h%ns >= 1  .AND.  surf_lsm_h%ns >= 1 )  THEN
     759             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
     760                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
     761                terminate_run_l = .TRUE.
     762             ENDIF
     763          ENDIF
     764
    368765
    369766#if defined( __parallel )
     
    389786!
    390787!--       Use the lowest possible value for z_mo
    391           k    = MINVAL(nzb_s_inner)
     788          k    = nzb
    392789          z_mo = zu(k+1) - zw(k)
    393790
     
    401798          zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
    402799
    403           DO l = 1, num_steps_n
    404              zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
     800          DO li = 1, num_steps_n
     801             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
    405802          ENDDO
    406803
     
    415812!
    416813!--       Calculate z/L range from zeta_min to zeta_stretch
    417           DO l = num_steps_n+1, num_steps-1
    418              zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
     814          DO li = num_steps_n+1, num_steps-1
     815             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
    419816             zeta_step = zeta_step * regr
    420817          ENDDO
     
    422819!
    423820!--       Save roughness lengths to temporary variables
    424           z0h_min = z0h(nys,nxl)
    425           z0_min  = z0(nys,nxl)
    426          
     821          IF ( surf_def_h(0)%ns >= 1  )  THEN
     822             z0h_min = surf_def_h(0)%z0h(1)
     823             z0_min  = surf_def_h(0)%z0(1)
     824          ELSEIF ( surf_lsm_h%ns >= 1 )  THEN
     825             z0h_min = surf_lsm_h%z0h(1)
     826             z0_min  = surf_lsm_h%z0(1)
     827          ELSEIF ( surf_usm_h%ns >= 1 )  THEN
     828             z0h_min = surf_usm_h%z0h(1)
     829             z0_min  = surf_usm_h%z0(1)
     830          ENDIF         
    427831!
    428832!--       Calculate lookup table for the Richardson number versus Obukhov length
     
    430834!--       boundary conditions for temperature
    431835          IF ( ibc_pt_b == 1 )  THEN
    432              DO l = 0, num_steps-1
    433                 ol_tab(l)  = - z_mo / zeta_tmp(num_steps-1-l)
    434                 rib_tab(l) = z_mo / ol_tab(l)  / ( LOG( z_mo / z0_min )        &
    435                                                 - psi_m( z_mo / ol_tab(l) )    &
    436                                                 + psi_m( z0_min / ol_tab(l) )  &
     836             DO li = 0, num_steps-1
     837                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
     838                rib_tab(li) = z_mo / ol_tab(li)  / ( LOG( z_mo / z0_min )       &
     839                                                - psi_m( z_mo / ol_tab(li) )    &
     840                                                + psi_m( z0_min / ol_tab(li) )  &
    437841                                                  )**3
    438842             ENDDO 
    439843          ELSE
    440              DO l = 0, num_steps-1
    441                 ol_tab(l)  = - z_mo / zeta_tmp(num_steps-1-l)
    442                 rib_tab(l) = z_mo / ol_tab(l)  * ( LOG( z_mo / z0h_min )       &
    443                                               - psi_h( z_mo / ol_tab(l) )      &
    444                                               + psi_h( z0h_min / ol_tab(l) )   &
     844             DO li = 0, num_steps-1
     845                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
     846                rib_tab(li) = z_mo / ol_tab(li)  * ( LOG( z_mo / z0h_min )     &
     847                                              - psi_h( z_mo / ol_tab(li) )     &
     848                                              + psi_h( z0h_min / ol_tab(li) )  &
    445849                                            )                                  &
    446850                                          / ( LOG( z_mo / z0_min )             &
    447                                               - psi_m( z_mo / ol_tab(l) )      &
    448                                               + psi_m( z0_min / ol_tab(l) )    &
     851                                              - psi_m( z_mo / ol_tab(li) )     &
     852                                              + psi_m( z0_min / ol_tab(li) )   &
    449853                                            )**2
    450854             ENDDO
     
    467871! ------------
    468872!> Compute the absolute value of the horizontal velocity (relative to the
    469 !> surface). This is required by all methods
     873!> surface) for horizontal surface elements. This is required by all methods.
    470874!------------------------------------------------------------------------------!
    471     SUBROUTINE calc_uv_total
     875    SUBROUTINE calc_uvw_abs
    472876
    473877       IMPLICIT NONE
    474878
    475 
    476        !$OMP PARALLEL DO PRIVATE( k )
    477        DO  i = nxl, nxr
    478           DO  j = nys, nyn
    479 
    480              k   = nzb_s_inner(j,i)
    481              uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
    482                                          - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
    483                               ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
    484                                          - v(k,j,i)   - v(k,j+1,i) ) )**2 )
    485 
    486 !
    487 !--          For too small values of the local wind, MOST does not work. A
    488 !--          threshold value is thus set if required
    489 !            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
    490 
     879       INTEGER(iwp) ::  i             !< running index x direction
     880       INTEGER(iwp) ::  ibit          !< flag to mask computation of relative velocity in case of downward-facing surfaces
     881       INTEGER(iwp) ::  j             !< running index y direction
     882       INTEGER(iwp) ::  k             !< running index z direction
     883       INTEGER(iwp) ::  m             !< running index surface elements
     884
     885!
     886!--    ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces.
     887       ibit = MERGE( 1, 0, .NOT. downward )
     888
     889       DO  m = 1, surf%ns
     890
     891          i   = surf%i(m)           
     892          j   = surf%j(m)
     893          k   = surf%k(m)
     894!
     895!--       Compute the absolute value of the horizontal velocity.
     896!--       (relative to the surface in case the lower surface is the ocean).
     897!--       Please note, in new surface modelling concept the index values changed,
     898!--       i.e. the reference grid point is not the surface-grid point itself but
     899!--       the first grid point outside of the topography.
     900!--       Note, in case of coupled ocean-atmosphere simulations relative velocity
     901!--       with respect to the ocean surface is used, hence, (k-1,j,i) values
     902!--       are used to calculate the absolute velocity.
     903!--       However, this do not apply for downward-facing walls. To mask this,
     904!--       use ibit, which checks for upward/downward-facing surfaces.
     905
     906          surf%uvw_abs(m) = SQRT(                                              &
     907                              ( 0.5_wp * (   u(k,j,i)   + u(k,j,i+1)           &
     908                                        -  ( u(k-1,j,i) + u(k-1,j,i+1)         &
     909                                           ) * ibit                            &
     910                                         )                                     &
     911                              )**2 +                                           &
     912                              ( 0.5_wp * (   v(k,j,i)   + v(k,j+1,i)           &
     913                                        -  ( v(k-1,j,i) + v(k-1,j+1,i)         &
     914                                           ) * ibit                            &
     915                                         )                                     &
     916                              )**2                                             &
     917                                )
     918
     919       ENDDO
     920
     921    END SUBROUTINE calc_uvw_abs
     922
     923
     924!------------------------------------------------------------------------------!
     925! Description:
     926! ------------
     927!> Compute the absolute value of the horizontal velocity (relative to the
     928!> surface) for horizontal surface elements. This is required by all methods.
     929!------------------------------------------------------------------------------!
     930    SUBROUTINE calc_uvw_abs_v_ugrid
     931
     932       IMPLICIT NONE
     933
     934       INTEGER(iwp) ::  i             !< running index x direction
     935       INTEGER(iwp) ::  j             !< running index y direction
     936       INTEGER(iwp) ::  k             !< running index z direction
     937       INTEGER(iwp) ::  m             !< running index surface elements
     938
     939       REAL(wp)     ::  u_i
     940       REAL(wp)     ::  w_i
     941
     942
     943       DO  m = 1, surf%ns
     944          i   = surf%i(m)           
     945          j   = surf%j(m)
     946          k   = surf%k(m)
     947!
     948!--       Compute the absolute value of the surface parallel velocity on u-grid.
     949          u_i  = u(k,j,i)
     950          w_i  = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) +                       &
     951                             w(k,j,i-1)   + w(k,j,i) )
     952
     953          surf%uvw_abs(m) = SQRT( u_i**2 + w_i**2 )
     954
     955       ENDDO
     956
     957    END SUBROUTINE calc_uvw_abs_v_ugrid
     958
     959!------------------------------------------------------------------------------!
     960! Description:
     961! ------------
     962!> Compute the absolute value of the horizontal velocity (relative to the
     963!> surface) for horizontal surface elements. This is required by all methods.
     964!------------------------------------------------------------------------------!
     965    SUBROUTINE calc_uvw_abs_v_vgrid
     966
     967       IMPLICIT NONE
     968
     969       INTEGER(iwp) ::  i             !< running index x direction
     970       INTEGER(iwp) ::  j             !< running index y direction
     971       INTEGER(iwp) ::  k             !< running index z direction
     972       INTEGER(iwp) ::  m             !< running index surface elements
     973
     974       REAL(wp)     ::  v_i
     975       REAL(wp)     ::  w_i
     976
     977
     978       DO  m = 1, surf%ns
     979          i   = surf%i(m)           
     980          j   = surf%j(m)
     981          k   = surf%k(m)
     982
     983          v_i  = u(k,j,i)
     984          w_i  = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) +                       &
     985                             w(k,j-1,i)   + w(k,j,i) )
     986
     987          surf%uvw_abs(m) = SQRT( v_i**2 + w_i**2 )
     988
     989       ENDDO
     990
     991    END SUBROUTINE calc_uvw_abs_v_vgrid
     992
     993!------------------------------------------------------------------------------!
     994! Description:
     995! ------------
     996!> Compute the absolute value of the horizontal velocity (relative to the
     997!> surface) for horizontal surface elements. This is required by all methods.
     998!------------------------------------------------------------------------------!
     999    SUBROUTINE calc_uvw_abs_v_wgrid
     1000
     1001       IMPLICIT NONE
     1002
     1003       INTEGER(iwp) ::  i             !< running index x direction
     1004       INTEGER(iwp) ::  j             !< running index y direction
     1005       INTEGER(iwp) ::  k             !< running index z direction
     1006       INTEGER(iwp) ::  m             !< running index surface elements
     1007
     1008       REAL(wp)     ::  u_i
     1009       REAL(wp)     ::  v_i
     1010       REAL(wp)     ::  w_i
     1011!
     1012!--    North- (l=0) and south-facing (l=1) surfaces
     1013       IF ( l == 0  .OR.  l == 1 )  THEN
     1014          DO  m = 1, surf%ns
     1015             i   = surf%i(m)           
     1016             j   = surf%j(m)
     1017             k   = surf%k(m)
     1018
     1019             u_i  = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) +                    &
     1020                                u(k,j,i+1)   + u(k,j,i) )
     1021             v_i  = 0.0_wp
     1022             w_i  = w(k,j,i)
     1023
     1024             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 )
    4911025          ENDDO
    492        ENDDO
    493 
    494 !
    495 !--    Values of uv_total need to be exchanged at the ghost boundaries
    496        CALL exchange_horiz_2d( uv_total )
    497 
    498     END SUBROUTINE calc_uv_total
     1026!
     1027!--    East- (l=2) and west-facing (l=3) surfaces
     1028       ELSE
     1029          DO  m = 1, surf%ns
     1030             i   = surf%i(m)           
     1031             j   = surf%j(m)
     1032             k   = surf%k(m)
     1033
     1034             u_i  = 0.0_wp
     1035             v_i  = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) +                    &
     1036                                v(k,j+1,i)   + v(k,j,i) )
     1037             w_i  = w(k,j,i)
     1038
     1039             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 )
     1040          ENDDO
     1041       ENDIF           
     1042
     1043    END SUBROUTINE calc_uvw_abs_v_wgrid
     1044
     1045!------------------------------------------------------------------------------!
     1046! Description:
     1047! ------------
     1048!> Compute the absolute value of the horizontal velocity (relative to the
     1049!> surface) for horizontal surface elements. This is required by all methods.
     1050!------------------------------------------------------------------------------!
     1051    SUBROUTINE calc_uvw_abs_v_sgrid
     1052
     1053       IMPLICIT NONE
     1054
     1055       INTEGER(iwp) ::  i             !< running index x direction
     1056       INTEGER(iwp) ::  j             !< running index y direction
     1057       INTEGER(iwp) ::  k             !< running index z direction
     1058       INTEGER(iwp) ::  m             !< running index surface elements
     1059
     1060       REAL(wp)     ::  u_i
     1061       REAL(wp)     ::  v_i
     1062       REAL(wp)     ::  w_i
     1063
     1064!
     1065!--    North- (l=0) and south-facing (l=1) walls
     1066       IF ( l == 0  .OR.  l == 1 )  THEN
     1067          DO  m = 1, surf%ns
     1068             i   = surf%i(m)           
     1069             j   = surf%j(m)
     1070             k   = surf%k(m)
     1071
     1072             u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
     1073             v_i = 0.0_wp
     1074             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     1075
     1076             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 )
     1077          ENDDO
     1078!
     1079!--    East- (l=2) and west-facing (l=3) walls
     1080       ELSE
     1081          DO  m = 1, surf%ns
     1082             i   = surf%i(m)           
     1083             j   = surf%j(m)
     1084             k   = surf%k(m)
     1085
     1086             u_i = 0.0_wp
     1087             v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
     1088             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     1089
     1090             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 )
     1091          ENDDO
     1092       ENDIF 
     1093
     1094    END SUBROUTINE calc_uvw_abs_v_sgrid
    4991095
    5001096
     
    5081104       IMPLICIT NONE
    5091105
    510        INTEGER(iwp) :: iter,  & !< Newton iteration step
    511                        l        !< look index
    512 
    513        REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number
     1106       INTEGER(iwp) ::  iter    !< Newton iteration step
     1107       INTEGER(iwp) ::  li      !< look index
     1108       INTEGER(iwp) ::  m       !< loop variable over all horizontal wall elements
    5141109
    5151110       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
     
    5211116
    5221117       IF ( TRIM( most_method ) /= 'circular' )  THEN
    523      
    524           !$OMP PARALLEL DO PRIVATE( k, z_mo )
    525           DO  i = nxl, nxr
    526              DO  j = nys, nyn
    527 
    528                 k   = nzb_s_inner(j,i)
    529                 z_mo = zu(k+1) - zw(k)
    530 
    531 !
    532 !--             Evaluate bulk Richardson number (calculation depends on
    533 !--             definition based on setting of boundary conditions
    534                 IF ( ibc_pt_b /= 1 )  THEN
    535                    IF ( humidity )  THEN
    536                       rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) )      &
    537                            / ( uv_total(j,i)**2 * vpt(k+1,j,i) + 1.0E-20_wp )
    538                    ELSE
    539                       rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) )        &
    540                            / ( uv_total(j,i)**2 * pt(k+1,j,i)  + 1.0E-20_wp )
    541                    ENDIF     
    542                 ELSE
    543 !
    544 !--                When using Neumann boundary conditions, the buoyancy flux
    545 !--                is required but cannot be calculated at the surface, as pt
    546 !--                and q are not known at the surface. Hence the values at
    547 !--                first grid level are used to estimate the buoyancy flux
    548                    IF ( humidity )  THEN
    549                       rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
    550                                  * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
    551                                  * pt(k+1,j,i) * qsws(j,i) ) * drho_air_zw(k)  &
    552                                  / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2&
    553                                  + 1.0E-20_wp)
    554                    ELSE
    555                       rib(j,i) = - g * z_mo * shf(j,i) * drho_air_zw(k)        &
    556                            / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2       &
     1118!
     1119!--       Evaluate bulk Richardson number (calculation depends on
     1120!--       definition based on setting of boundary conditions
     1121          IF ( ibc_pt_b /= 1 )  THEN
     1122             IF ( humidity )  THEN
     1123                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1124                DO  m = 1, surf%ns
     1125
     1126                   i   = surf%i(m)           
     1127                   j   = surf%j(m)
     1128                   k   = surf%k(m)
     1129
     1130                   z_mo = surf%z_mo(m)
     1131
     1132                   surf%rib(m) = g * z_mo *                                    &
     1133                                         ( vpt(k,j,i) - vpt(k-1,j,i) ) /       &
     1134                      ( surf%uvw_abs(m)**2 * vpt(k,j,i) + 1.0E-20_wp )
     1135                ENDDO
     1136             ELSE
     1137                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1138                DO  m = 1, surf%ns
     1139
     1140                   i   = surf%i(m)           
     1141                   j   = surf%j(m)
     1142                   k   = surf%k(m)
     1143
     1144                   z_mo = surf%z_mo(m)
     1145
     1146                   surf%rib(m) = g * z_mo *                                    &
     1147                                         ( pt(k,j,i) - pt(k-1,j,i)   ) /       &
     1148                      ( surf%uvw_abs(m)**2 * pt(k,j,i)  + 1.0E-20_wp )
     1149                ENDDO
     1150             ENDIF
     1151          ELSE
     1152             IF ( humidity )  THEN
     1153                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1154                DO  m = 1, surf%ns
     1155
     1156                   i   = surf%i(m)           
     1157                   j   = surf%j(m)
     1158                   k   = surf%k(m)
     1159
     1160                   z_mo = surf%z_mo(m)
     1161
     1162                   surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
     1163                           * q(k,j,i) ) * surf%shf(m) + 0.61_wp                &
     1164                           * pt(k,j,i) * surf%qsws(m) ) *                      &
     1165                             drho_air_zw(k-1)                /                 &
     1166                         ( surf%uvw_abs(m)**3 * vpt(k,j,i) * kappa**2          &
    5571167                           + 1.0E-20_wp )
    558                    ENDIF
    559                 ENDIF 
    560      
    561              ENDDO
    562           ENDDO
     1168                ENDDO
     1169             ELSE
     1170                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1171                DO  m = 1, surf%ns
     1172
     1173                   i   = surf%i(m)           
     1174                   j   = surf%j(m)
     1175                   k   = surf%k(m)
     1176
     1177                   z_mo = surf%z_mo(m)
     1178
     1179                   surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
     1180                                        drho_air_zw(k-1)            /          &
     1181                        ( surf%uvw_abs(m)**3 * pt(k,j,i) * kappa**2            &
     1182                           + 1.0E-20_wp )
     1183                ENDDO
     1184             ENDIF
     1185          ENDIF
    5631186
    5641187       ENDIF
     1188
    5651189
    5661190!
     
    5691193       IF ( TRIM( most_method ) == 'newton' )  THEN
    5701194
    571           !$OMP PARALLEL DO PRIVATE( k, z_mo )
    572           DO  i = nxl, nxr
    573              DO  j = nys, nyn
    574 
    575                 k   = nzb_s_inner(j,i)
    576                 z_mo = zu(k+1) - zw(k)
    577 
    578 !
    579 !--             Store current value in case the Newton iteration fails
    580                 ol_old = ol(j,i)
     1195          DO  m = 1, surf%ns
     1196
     1197             i   = surf%i(m)           
     1198             j   = surf%j(m)
     1199             k   = surf%k(m)
     1200
     1201             z_mo = surf%z_mo(m)
     1202
     1203!
     1204!--          Store current value in case the Newton iteration fails
     1205             ol_old = surf%ol(m)
     1206
     1207!
     1208!--          Ensure that the bulk Richardson number and the Obukhov
     1209!--          length have the same sign
     1210             IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
     1211                  ABS( surf%ol(m) ) == ol_max )  THEN
     1212                IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
     1213                IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
     1214             ENDIF
     1215!
     1216!--          Iteration to find Obukhov length
     1217             iter = 0
     1218             DO
     1219                iter = iter + 1
     1220!
     1221!--             In case of divergence, use the value of the previous time step
     1222                IF ( iter > 1000 )  THEN
     1223                   surf%ol(m) = ol_old
     1224                   EXIT
     1225                ENDIF
     1226
     1227                ol_m = surf%ol(m)
     1228                ol_l = ol_m - 0.001_wp * ol_m
     1229                ol_u = ol_m + 0.001_wp * ol_m
     1230
     1231
     1232                IF ( ibc_pt_b /= 1 )  THEN
     1233!
     1234!--                Calculate f = Ri - [...]/[...]^2 = 0
     1235                   f = surf%rib(m) - ( z_mo / ol_m ) * (                       &
     1236                                                 LOG( z_mo / surf%z0h(m) )     &
     1237                                                 - psi_h( z_mo / ol_m )        &
     1238                                                 + psi_h( surf%z0h(m) /        &
     1239                                                          ol_m )               &
     1240                                                               )               &
     1241                                              / ( LOG( z_mo / surf%z0(m) )     &
     1242                                                 - psi_m( z_mo / ol_m )        &
     1243                                                 + psi_m( surf%z0(m) /         &
     1244                                                          ol_m )               &
     1245                                                 )**2
     1246
     1247!
     1248!--                 Calculate df/dL
     1249                    f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /               &
     1250                                                             surf%z0h(m) )     &
     1251                                            - psi_h( z_mo / ol_u )             &
     1252                                            + psi_h( surf%z0h(m) / ol_u )      &
     1253                                              )                                &
     1254                                            / ( LOG( z_mo / surf%z0(m) )       &
     1255                                            - psi_m( z_mo / ol_u )             &
     1256                                            + psi_m( surf%z0(m) / ol_u )       &
     1257                                              )**2                             &
     1258                           + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )     &
     1259                                            - psi_h( z_mo / ol_l )             &
     1260                                            + psi_h( surf%z0h(m) / ol_l )      &
     1261                                              )                                &
     1262                                            / ( LOG( z_mo / surf%z0(m) )       &
     1263                                            - psi_m( z_mo / ol_l )             &
     1264                                            + psi_m( surf%z0(m) / ol_l )       &
     1265                                              )**2                             &
     1266                             ) / ( ol_u - ol_l )
     1267                ELSE
     1268!
     1269!--                Calculate f = Ri - 1 /[...]^3 = 0
     1270                   f = surf%rib(m) - ( z_mo / ol_m ) /                         &
     1271                                                ( LOG( z_mo / surf%z0(m) )     &
     1272                                            - psi_m( z_mo / ol_m )             &
     1273                                            + psi_m( surf%z0(m) / ol_m )       &
     1274                                                )**3
     1275
     1276!
     1277!--                Calculate df/dL
     1278                   f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo /                &
     1279                                                            surf%z0(m) )       &
     1280                                            - psi_m( z_mo / ol_u )             &
     1281                                            + psi_m( surf%z0(m) / ol_u )       &
     1282                                                     )**3                      &
     1283                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
     1284                                            - psi_m( z_mo / ol_l )             &
     1285                                            + psi_m( surf%z0(m) / ol_l )       &
     1286                                               )**3                            &
     1287                                  ) / ( ol_u - ol_l )
     1288                ENDIF
     1289!
     1290!--             Calculate new L
     1291                surf%ol(m) = ol_m - f / f_d_ol
    5811292
    5821293!
    5831294!--             Ensure that the bulk Richardson number and the Obukhov
    584 !--             lengtH have the same sign
    585                 IF ( rib(j,i) * ol(j,i) < 0.0_wp  .OR.                         &
    586                      ABS( ol(j,i) ) == ol_max )  THEN
    587                    IF ( rib(j,i) > 0.0_wp ) ol(j,i) =  0.01_wp
    588                    IF ( rib(j,i) < 0.0_wp ) ol(j,i) = -0.01_wp
     1295!--             length have the same sign and ensure convergence.
     1296                IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
     1297!
     1298!--             If unrealistic value occurs, set L to the maximum
     1299!--             value that is allowed
     1300                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
     1301                   surf%ol(m) = ol_max
     1302                   EXIT
    5891303                ENDIF
    5901304!
    591 !--             Iteration to find Obukhov length
    592                 iter = 0
    593                 DO
    594                    iter = iter + 1
    595 !
    596 !--                In case of divergence, use the value of the previous time step
    597                    IF ( iter > 1000 )  THEN
    598                       ol(j,i) = ol_old
    599                       EXIT
    600                    ENDIF
    601 
    602                    ol_m = ol(j,i)
    603                    ol_l = ol_m - 0.001_wp * ol_m
    604                    ol_u = ol_m + 0.001_wp * ol_m
    605 
    606 
    607                    IF ( ibc_pt_b /= 1 )  THEN
    608 !
    609 !--                   Calculate f = Ri - [...]/[...]^2 = 0
    610                       f = rib(j,i) - ( z_mo / ol_m ) * ( LOG( z_mo / z0h(j,i) )&
    611                                                     - psi_h( z_mo / ol_m )     &
    612                                                     + psi_h( z0h(j,i) / ol_m ) &
    613                                                    )                           &
    614                                                  / ( LOG( z_mo / z0(j,i) )     &
    615                                                     - psi_m( z_mo / ol_m )     &
    616                                                     + psi_m( z0(j,i) / ol_m )  &
    617                                                     )**2
    618 
    619 !
    620 !--                    Calculate df/dL
    621                        f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / z0h(j,i) ) &
    622                                                    - psi_h( z_mo / ol_u )      &
    623                                                    + psi_h( z0h(j,i) / ol_u )  &
    624                                                  )                             &
    625                                                / ( LOG( z_mo / z0(j,i) )       &
    626                                                    - psi_m( z_mo / ol_u )      &
    627                                                    + psi_m( z0(j,i) / ol_u )   &
    628                                                  )**2                          &
    629                               + ( z_mo / ol_l ) * ( LOG( z_mo / z0h(j,i) )     &
    630                                                    - psi_h( z_mo / ol_l )      &
    631                                                    + psi_h( z0h(j,i) / ol_l )  &
    632                                                  )                             &
    633                                                / ( LOG( z_mo / z0(j,i) )       &
    634                                                    - psi_m( z_mo / ol_l )      &
    635                                                    + psi_m( z0(j,i) / ol_l )   &
    636                                                  )**2                          &
    637                                 ) / ( ol_u - ol_l )
    638                    ELSE
    639 !
    640 !--                   Calculate f = Ri - 1 /[...]^3 = 0
    641                       f = rib(j,i) - ( z_mo / ol_m ) / ( LOG( z_mo / z0(j,i) )&
    642                                                     - psi_m( z_mo / ol_m )    &
    643                                                     + psi_m( z0(j,i) / ol_m ) &
    644                                                        )**3
    645 
    646 !
    647 !--                   Calculate df/dL
    648                       f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / z0(j,i) )  &
    649                                                    - psi_m( z_mo / ol_u )     &
    650                                                    + psi_m( z0(j,i) / ol_u )  &
    651                                                  )**3                         &
    652                               + ( z_mo / ol_l ) / ( LOG( z_mo / z0(j,i) )     &
    653                                                    - psi_m( z_mo / ol_l )     &
    654                                                    + psi_m( z0(j,i) / ol_l )  &
    655                                                  )**3                         &
    656                                      ) / ( ol_u - ol_l )
    657                    ENDIF
    658 !
    659 !--                Calculate new L
    660                    ol(j,i) = ol_m - f / f_d_ol
    661 
    662 !
    663 !--                Ensure that the bulk Richardson number and the Obukhov
    664 !--                length have the same sign and ensure convergence.
    665                    IF ( ol(j,i) * ol_m < 0.0_wp )  ol(j,i) = ol_m * 0.5_wp
    666 
    667 !
    668 !--                If unrealistic value occurs, set L to the maximum
    669 !--                value that is allowed
    670                    IF ( ABS( ol(j,i) ) > ol_max )  THEN
    671                       ol(j,i) = ol_max
    672                       EXIT
    673                    ENDIF
    674 !
    675 !--                Check for convergence
    676                    IF ( ABS( ( ol(j,i) - ol_m ) / ol(j,i) ) < 1.0E-4_wp )  THEN
    677                       EXIT
    678                    ELSE
    679                       CYCLE
    680                    ENDIF
    681 
     1305!--             Check for convergence
     1306                IF ( ABS( ( surf%ol(m) - ol_m ) /                              &
     1307                     surf%ol(m) ) < 1.0E-4_wp )  THEN
     1308                   EXIT
     1309                ELSE
     1310                   CYCLE
     1311                ENDIF
     1312
     1313             ENDDO
     1314          ENDDO
     1315
     1316       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
     1317
     1318          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
     1319          DO  m = 1, surf%ns
     1320
     1321             i   = surf%i(m)           
     1322             j   = surf%j(m)
     1323             k   = surf%k(m)
     1324!
     1325!--          If the bulk Richardson number is outside the range of the lookup
     1326!--          table, set it to the exceeding threshold value
     1327             IF ( surf%rib(m) < rib_min )  surf%rib(m) = rib_min
     1328             IF ( surf%rib(m) > rib_max )  surf%rib(m) = rib_max
     1329
     1330!
     1331!--          Find the correct index bounds for linear interpolation. As the
     1332!--          Richardson number will not differ very much from time step to
     1333!--          time step , use the index from the last step and search in the
     1334!--          correct direction
     1335             li = li_bnd
     1336             IF ( rib_tab(li) - surf%rib(m) > 1.0_wp )  THEN
     1337                DO  WHILE ( rib_tab(li-1) - surf%rib(m) > 1.0_wp  .AND.  li > 1 )
     1338                   li = li-1
    6821339                ENDDO
    683                        
    684              ENDDO
     1340             ELSE
     1341                DO  WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp                 &
     1342                           .AND.  li < num_steps-1 )
     1343                   li = li+1
     1344                ENDDO
     1345             ENDIF
     1346             li_bnd = li
     1347
     1348!
     1349!--          Linear interpolation to find the correct value of z/L
     1350             surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) )       &
     1351                         / (  rib_tab(li) - rib_tab(li-1) )                    &
     1352                         * ( surf%rib(m) - rib_tab(li-1) ) )
     1353
    6851354          ENDDO
    6861355
    687        ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
    688 
    689           !$OMP PARALLEL DO PRIVATE( k, l, z_mo ) FIRSTPRIVATE( l_bnd ) LASTPRIVATE( l_bnd )
    690           DO  i = nxl, nxr
    691              DO  j = nys, nyn
    692 
    693 !
    694 !--             If the bulk Richardson number is outside the range of the lookup
    695 !--             table, set it to the exceeding threshold value
    696                 IF ( rib(j,i) < rib_min )  rib(j,i) = rib_min
    697                 IF ( rib(j,i) > rib_max )  rib(j,i) = rib_max
    698 
    699 !
    700 !--             Find the correct index bounds for linear interpolation. As the
    701 !--             Richardson number will not differ very much from time step to
    702 !--             time step , use the index from the last step and search in the
    703 !--             correct direction
    704                 l = l_bnd
    705                 IF ( rib_tab(l) - rib(j,i) > 0.0_wp )  THEN
    706                    DO  WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp  .AND.  l > 0 )
    707                       l = l-1
    708                    ENDDO
    709                 ELSE
    710                    DO  WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                  &
    711                               .AND.  l < num_steps-1 )
    712                       l = l+1
    713                    ENDDO
    714                 ENDIF
    715                 l_bnd = l
    716 
    717 !
    718 !--             Linear interpolation to find the correct value of z/L
    719                 ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )          &
    720                             / (  rib_tab(l) - rib_tab(l-1) )                   &
    721                             * ( rib(j,i) - rib_tab(l-1) ) )
    722 
    723              ENDDO
    724           ENDDO
    725 
    7261356       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
    7271357
    728           !$OMP PARALLEL DO PRIVATE( k, z_mo )
    729           DO  i = nxl, nxr
    730              DO  j = nys, nyn
    731 
    732                 k   = nzb_s_inner(j,i)
    733                 z_mo = zu(k+1) - zw(k)
    734 
    735                 IF ( .NOT. humidity )  THEN
    736                    ol(j,i) =  ( pt(k+1,j,i) *  us(j,i)**2 ) / ( kappa * g      &
    737                               * ts(j,i) + 1E-30_wp )
    738                 ELSEIF ( cloud_physics )  THEN
    739 
    740                    ol(j,i) =  ( vpt(k+1,j,i) * us(j,i)**2 ) / ( kappa * g      &
    741                               * ( ts(j,i) + 0.61_wp * pt1(j,i) * qs(j,i)       &
    742                               + 0.61_wp * qv1(j,i) * ts(j,i) - ts(j,i)         &
    743                               * ql(k+1,j,i) ) + 1E-30_wp )
    744                 ELSE
    745                    ol(j,i) =  ( vpt(k+1,j,i) *  us(j,i)**2 ) / ( kappa * g     &
    746                               * ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i)    &
    747                                   + 0.61_wp * q(k+1,j,i) * ts(j,i) ) + 1E-30_wp )
    748                 ENDIF
     1358          IF ( .NOT. humidity )  THEN
     1359             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1360             DO  m = 1, surf%ns
     1361
     1362                i   = surf%i(m)           
     1363                j   = surf%j(m)
     1364                k   = surf%k(m)
     1365
     1366                z_mo = surf%z_mo(m)
     1367
     1368                surf%ol(m) =  ( pt(k,j,i) *  surf%us(m)**2 ) /                 &
     1369                                       ( kappa * g *                           &
     1370                                         surf%ts(m) + 1E-30_wp )
    7491371!
    7501372!--             Limit the value range of the Obukhov length.
    751 !--             This is necessary for very small velocities (u,v --> 0), because
     1373!--             This is necessary for very small velocities (u,v --> 1), because
    7521374!--             the absolute value of ol can then become very small, which in
    7531375!--             consequence would result in very large shear stresses and very
    7541376!--             small momentum fluxes (both are generally unrealistic).
    755                 IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) < zeta_min )            &
    756                    ol(j,i) = z_mo / zeta_min
    757                 IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) > zeta_max )            &
    758                    ol(j,i) = z_mo / zeta_max
    759              ENDDO
    760           ENDDO
     1377                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
     1378                   surf%ol(m) = z_mo / zeta_min
     1379                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
     1380                   surf%ol(m) = z_mo / zeta_max
     1381
     1382             ENDDO
     1383          ELSEIF ( cloud_physics )  THEN
     1384             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1385             DO  m = 1, surf%ns
     1386
     1387                i   = surf%i(m)           
     1388                j   = surf%j(m)
     1389                k   = surf%k(m)
     1390
     1391                z_mo = surf%z_mo(m)
     1392
     1393
     1394                surf%ol(m) =  ( vpt(k,j,i) * surf%us(m)**2 ) /                 &
     1395                    ( kappa * g * ( surf%ts(m) +                               &
     1396                         0.61_wp * surf%pt1(m) * surf%us(m)                    &
     1397                       + 0.61_wp * surf%qv1(m) * surf%ts(m) -                  &
     1398                                   surf%ts(m)  * ql(k,j,i) ) + 1E-30_wp )
     1399!
     1400!--             Limit the value range of the Obukhov length.
     1401!--             This is necessary for very small velocities (u,v --> 1), because
     1402!--             the absolute value of ol can then become very small, which in
     1403!--             consequence would result in very large shear stresses and very
     1404!--             small momentum fluxes (both are generally unrealistic).
     1405                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
     1406                   surf%ol(m) = z_mo / zeta_min
     1407                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
     1408                   surf%ol(m) = z_mo / zeta_max
     1409
     1410             ENDDO
     1411          ELSE
     1412
     1413             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1414             DO  m = 1, surf%ns
     1415
     1416                i   = surf%i(m)           
     1417                j   = surf%j(m)
     1418                k   = surf%k(m)
     1419
     1420                z_mo = surf%z_mo(m)
     1421
     1422                surf%ol(m) =  ( vpt(k,j,i) *  surf%us(m)**2 ) /                &
     1423                    ( kappa * g * ( surf%ts(m) + 0.61_wp * pt(k,j,i) *         &
     1424                                    surf%qs(m) + 0.61_wp * q(k,j,i)  *         &
     1425                                    surf%ts(m) ) + 1E-30_wp )
     1426
     1427!
     1428!--             Limit the value range of the Obukhov length.
     1429!--             This is necessary for very small velocities (u,v --> 1), because
     1430!--             the absolute value of ol can then become very small, which in
     1431!--             consequence would result in very large shear stresses and very
     1432!--             small momentum fluxes (both are generally unrealistic).
     1433                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
     1434                   surf%ol(m) = z_mo / zeta_min
     1435                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
     1436                   surf%ol(m) = z_mo / zeta_max
     1437
     1438             ENDDO
     1439
     1440          ENDIF
    7611441
    7621442       ENDIF
    763 
    764 !
    765 !--    Values of ol at ghost point locations are needed for the evaluation
    766 !--    of usws and vsws.
    767        CALL exchange_horiz_2d( ol )
    7681443
    7691444    END SUBROUTINE calc_ol
     
    7751450       IMPLICIT NONE
    7761451
    777        !$OMP PARALLEL DO PRIVATE( k, z_mo )
    778        DO  i = nxlg, nxrg
    779           DO  j = nysg, nyng
    780 
    781              k   = nzb_s_inner(j,i)+1
    782              z_mo = zu(k+1) - zw(k)
    783 
    784 !
    785 !--          Compute u* at the scalars' grid points
    786              us(j,i) = kappa * uv_total(j,i) / ( LOG( z_mo / z0(j,i) )         &
    787                                           - psi_m( z_mo / ol(j,i) )            &
    788                                           + psi_m( z0(j,i) / ol(j,i) ) )
     1452       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
     1453
     1454!
     1455!--    Compute u* at horizontal surfaces at the scalars' grid points
     1456       IF ( .NOT. surf_vertical )  THEN
     1457!
     1458!--       Compute u* at upward-facing surfaces
     1459          IF ( .NOT. downward )  THEN
     1460             !$OMP PARALLEL  DO PRIVATE( k, z_mo )
     1461             DO  m = 1, surf%ns
     1462
     1463                k    = surf%k(m)
     1464                z_mo = surf%z_mo(m)
     1465!
     1466!--             Compute u* at the scalars' grid points
     1467                surf%us(m) = kappa * surf%uvw_abs(m) /                         &
     1468                            ( LOG( z_mo / surf%z0(m) )                         &
     1469                           - psi_m( z_mo / surf%ol(m) )                        &
     1470                           + psi_m( surf%z0(m) / surf%ol(m) ) )
     1471   
     1472             ENDDO
     1473!
     1474!--       Compute u* at downward-facing surfaces. This case, do not consider
     1475!--       any stability.
     1476          ELSE
     1477             !$OMP PARALLEL  DO PRIVATE( k, z_mo )
     1478             DO  m = 1, surf%ns
     1479
     1480                k    = surf%k(m)
     1481                z_mo = surf%z_mo(m)
     1482!
     1483!--             Compute u* at the scalars' grid points
     1484                surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
     1485   
     1486             ENDDO
     1487          ENDIF
     1488!
     1489!--    Compute u* at vertical surfaces at the u/v/v grid, respectively.
     1490!--    No stability is considered in this case.
     1491       ELSE
     1492          !$OMP PARALLEL DO PRIVATE( z_mo )
     1493          DO  m = 1, surf%ns
     1494             z_mo = surf%z_mo(m)
     1495
     1496             surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
    7891497          ENDDO
    790        ENDDO
     1498       ENDIF
    7911499
    7921500    END SUBROUTINE calc_us
     
    7941502!
    7951503!-- Calculate potential temperature and specific humidity at first grid level
     1504!-- ( only for upward-facing surfs )
    7961505    SUBROUTINE calc_pt_q
    7971506
    7981507       IMPLICIT NONE
    7991508
    800        DO  i = nxlg, nxrg
    801           DO  j = nysg, nyng
    802              k   = nzb_s_inner(j,i)+1
    803              pt1(j,i) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
    804              qv1(j,i) = q(k,j,i) - ql(k,j,i)
    805           ENDDO
     1509       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
     1510
     1511       !$OMP PARALLEL DO PRIVATE( i, j, k )
     1512       DO  m = 1, surf%ns
     1513
     1514          i   = surf%i(m)           
     1515          j   = surf%j(m)
     1516          k   = surf%k(m)
     1517
     1518          surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
     1519          surf%qv1(m) = q(k,j,i) - ql(k,j,i)
     1520
    8061521       ENDDO
    8071522
     
    8141529       IMPLICIT NONE
    8151530
     1531
     1532       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
     1533
    8161534!
    817 !--    Compute theta*
    818        IF ( constant_heatflux )  THEN
    819 
     1535!--    Compute theta* at horizontal surfaces
     1536       IF ( constant_heatflux  .AND.  .NOT. surf_vertical )  THEN
    8201537!
    8211538!--       For a given heat flux in the surface layer:
    822           !$OMP PARALLEL DO
    823           DO  i = nxlg, nxrg
    824              DO  j = nysg, nyng
    825                 k   = nzb_s_inner(j,i)
    826                 ts(j,i) = -shf(j,i) * drho_air_zw(k) / ( us(j,i) + 1E-30_wp )
    827 !
    828 !--             ts must be limited, because otherwise overflow may occur in case
    829 !--             of us=0 when computing ol further below
    830                 IF ( ts(j,i) < -1.05E5_wp )  ts(j,i) = -1.0E5_wp
    831                 IF ( ts(j,i) >   1.0E5_wp )  ts(j,i) =  1.0E5_wp
    832              ENDDO
     1539
     1540          !$OMP PARALLEL DO PRIVATE( i, j, k )
     1541          DO  m = 1, surf%ns
     1542
     1543             i   = surf%i(m)           
     1544             j   = surf%j(m)
     1545             k   = surf%k(m)
     1546
     1547             surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) /                    &
     1548                                  ( surf%us(m) + 1E-30_wp )
     1549
     1550!
     1551!--          ts must be limited, because otherwise overflow may occur in case
     1552!--          of us=0 when computing ol further below
     1553             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
     1554             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
     1555
    8331556          ENDDO
    8341557
    835        ELSE
     1558       ELSEIF ( .NOT. surf_vertical ) THEN
    8361559!
    8371560!--       For a given surface temperature:
    8381561          IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
    839              !$OMP PARALLEL DO
    840              DO  i = nxlg, nxrg
    841                 DO  j = nysg, nyng
    842                    k = nzb_s_inner(j,i)
    843                    pt(k,j,i) = pt_surface
    844                 ENDDO
    845              ENDDO
    846           ENDIF
    847 
    848           !$OMP PARALLEL DO PRIVATE( k, z_mo )
    849           DO  i = nxlg, nxrg
    850              DO  j = nysg, nyng
    851 
    852                 k   = nzb_s_inner(j,i)
    853                 z_mo = zu(k+1) - zw(k)
    854 
    855                 IF ( cloud_physics )  THEN
    856                    ts(j,i) = kappa * ( pt1(j,i) - pt(k,j,i) )                  &
    857                                      / ( LOG( z_mo / z0h(j,i) )                &
    858                                          - psi_h( z_mo / ol(j,i) )             &
    859                                          + psi_h( z0h(j,i) / ol(j,i) ) )
    860                 ELSE
    861                    ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) )               &
    862                                      / ( LOG( z_mo / z0h(j,i) )                &
    863                                          - psi_h( z_mo / ol(j,i) )             &
    864                                          + psi_h( z0h(j,i) / ol(j,i) ) )
    865                 ENDIF
    866 
    867              ENDDO
     1562
     1563             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1564             DO  m = 1, surf%ns
     1565                i   = surf%i(m)           
     1566                j   = surf%j(m)
     1567                k   = surf%k(m)
     1568
     1569                pt(k-1,j,i) = pt_surface
     1570             ENDDO
     1571          ENDIF
     1572
     1573          IF ( cloud_physics )  THEN
     1574             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1575             DO  m = 1, surf%ns   
     1576
     1577                i   = surf%i(m)           
     1578                j   = surf%j(m)
     1579                k   = surf%k(m)
     1580
     1581                z_mo = surf%z_mo(m)
     1582
     1583                surf%ts(m) = kappa * ( surf%pt1(m) - pt(k-1,j,i) )             &
     1584                                     / ( LOG( z_mo / surf%z0h(m) )             &
     1585                                         - psi_h( z_mo / surf%ol(m) )          &
     1586                                         + psi_h( surf%z0h(m) / surf%ol(m) ) )
     1587
     1588             ENDDO
     1589          ELSE
     1590             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1591             DO  m = 1, surf%ns   
     1592
     1593                i   = surf%i(m)           
     1594                j   = surf%j(m)
     1595                k   = surf%k(m)
     1596
     1597                z_mo = surf%z_mo(m)
     1598
     1599                surf%ts(m) = kappa * ( pt(k,j,i) - pt(k-1,j,i) )               &
     1600                                     / ( LOG( z_mo / surf%z0h(m) )             &
     1601                                         - psi_h( z_mo / surf%ol(m) )          &
     1602                                         + psi_h( surf%z0h(m) / surf%ol(m) ) )
     1603             ENDDO
     1604          ENDIF
     1605       ENDIF
     1606!
     1607!--    Compute theta* at vertical surfaces. This is only required in case of
     1608!--    land-surface model, in order to compute aerodynamical resistance.
     1609       IF ( surf_vertical )  THEN
     1610          !$OMP PARALLEL DO PRIVATE( i, j, k )
     1611          DO  m = 1, surf%ns
     1612
     1613             i          =  surf%i(m)           
     1614             j          =  surf%j(m)
     1615             k          =  surf%k(m)
     1616             surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )
     1617!
     1618!--          ts must be limited, because otherwise overflow may occur in case
     1619!--          of us=0 when computing ol further below
     1620             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
     1621             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
     1622
    8681623          ENDDO
    8691624       ENDIF
    8701625
    8711626!
    872 !--    If required compute q*
     1627!--    If required compute q* at horizontal surfaces
    8731628       IF ( humidity )  THEN
    874           IF ( constant_waterflux )  THEN
     1629          IF ( constant_waterflux  .AND.  .NOT. surf_vertical )  THEN
    8751630!
    8761631!--          For a given water flux in the surface layer
    877              !$OMP PARALLEL DO
    878              DO  i = nxlg, nxrg
    879                 DO  j = nysg, nyng
    880                    k   = nzb_s_inner(j,i)
    881                    qs(j,i) = -qsws(j,i) * drho_air_zw(k) / ( us(j,i) + 1E-30_wp )
    882                 ENDDO
    883              ENDDO
    884 
    885           ELSE
     1632             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1633             DO  m = 1, surf%ns 
     1634
     1635                i   = surf%i(m)           
     1636                j   = surf%j(m)
     1637                k   = surf%k(m)
     1638                surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) /                &
     1639                                               ( surf%us(m) + 1E-30_wp )
     1640
     1641             ENDDO
     1642
     1643          ELSEIF ( .NOT. surf_vertical )  THEN
    8861644             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
    8871645                             run_coupled )
    8881646
    8891647             IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
    890                 !$OMP PARALLEL DO
    891                 DO  i = nxlg, nxrg
    892                    DO  j = nysg, nyng
    893                       k = nzb_s_inner(j,i)
    894                       q(k,j,i) = q_surface
    895                    ENDDO
     1648                !$OMP PARALLEL DO PRIVATE( i, j, k )
     1649                DO  m = 1, surf%ns 
     1650
     1651                   i   = surf%i(m)           
     1652                   j   = surf%j(m)
     1653                   k   = surf%k(m)
     1654                   q(k-1,j,i) = q_surface
     1655                   
    8961656                ENDDO
    8971657             ENDIF
    8981658
    899              !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
    900              DO  i = nxlg, nxrg
    901                 DO  j = nysg, nyng
    902 
    903                    k   = nzb_s_inner(j,i)
    904                    z_mo = zu(k+1) - zw(k)
    905 
    906 !
    907 !--                Assume saturation for atmosphere coupled to ocean (but not
    908 !--                in case of precursor runs)
    909                    IF ( coupled_run )  THEN
    910                       e_s = 6.1_wp * &
    911                               EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i))      &
     1659!
     1660!--          Assume saturation for atmosphere coupled to ocean (but not
     1661!--          in case of precursor runs)
     1662             IF ( coupled_run )  THEN
     1663                !$OMP PARALLEL DO PRIVATE( i, j, k, e_s )
     1664                DO  m = 1, surf%ns   
     1665                   i   = surf%i(m)           
     1666                   j   = surf%j(m)
     1667                   k   = surf%k(m)
     1668                   e_s = 6.1_wp * &
     1669                              EXP( 0.07_wp * ( MIN(pt(k-1,j,i),pt(k,j,i))      &
    9121670                                               - 273.15_wp ) )
    913                       q(k,j,i) = 0.622_wp * e_s / ( surface_pressure - e_s )
    914                    ENDIF
    915 
    916                    IF ( cloud_physics )  THEN
    917                       qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) )                &
    918                                         / ( LOG( z_mo / z0q(j,i) )             &
    919                                             - psi_h( z_mo / ol(j,i) )          &
    920                                             + psi_h( z0q(j,i) / ol(j,i) ) )
    921 
    922                    ELSE
    923                       qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) )              &
    924                                         / ( LOG( z_mo / z0q(j,i) )             &
    925                                             - psi_h( z_mo / ol(j,i) )          &
    926                                             + psi_h( z0q(j,i) / ol(j,i) ) )
    927                    ENDIF
    928 
     1671                   q(k-1,j,i) = 0.622_wp * e_s / ( surface_pressure - e_s )
    9291672                ENDDO
     1673             ENDIF
     1674
     1675             IF ( cloud_physics )  THEN
     1676               !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1677                DO  m = 1, surf%ns   
     1678
     1679                   i   = surf%i(m)           
     1680                   j   = surf%j(m)
     1681                   k   = surf%k(m)
     1682   
     1683                   z_mo = surf%z_mo(m)
     1684
     1685                   surf%qs(m) = kappa * ( surf%qv1(m) - q(k-1,j,i) )           &
     1686                                        / ( LOG( z_mo / surf%z0q(m) )          &
     1687                                            - psi_h( z_mo / surf%ol(m) )       &
     1688                                            + psi_h( surf%z0q(m) /             &
     1689                                                     surf%ol(m) ) )
     1690                ENDDO
     1691             ELSE
     1692                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1693                DO  m = 1, surf%ns   
     1694
     1695                   i   = surf%i(m)           
     1696                   j   = surf%j(m)
     1697                   k   = surf%k(m)
     1698   
     1699                   z_mo = surf%z_mo(m)
     1700
     1701                   surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) )              &
     1702                                        / ( LOG( z_mo / surf%z0q(m) )          &
     1703                                            - psi_h( z_mo / surf%ol(m) )       &
     1704                                            + psi_h( surf%z0q(m) /             &
     1705                                                     surf%ol(m) ) )
     1706                ENDDO
     1707             ENDIF
     1708          ENDIF
     1709!
     1710!--       Compute q* at vertical surfaces
     1711          IF ( surf_vertical )  THEN
     1712             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1713             DO  m = 1, surf%ns 
     1714
     1715                i          =  surf%i(m)           
     1716                j          =  surf%j(m)
     1717                k          =  surf%k(m)
     1718                surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )
     1719
    9301720             ENDDO
    9311721          ENDIF
     
    9351725!--    If required compute s*
    9361726       IF ( passive_scalar )  THEN
    937           IF ( constant_scalarflux )  THEN
    938 !
    939 !--          For a given water flux in the surface layer
    940              !$OMP PARALLEL DO
    941              DO  i = nxlg, nxrg
    942                 DO  j = nysg, nyng
    943                    ss(j,i) = -ssws(j,i) / ( us(j,i) + 1E-30_wp )
    944                 ENDDO
     1727!
     1728!--       At horizontal surfaces
     1729          IF ( constant_scalarflux  .AND.  .NOT. surf_vertical )  THEN
     1730!
     1731!--          For a given scalar flux in the surface layer
     1732             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1733             DO  m = 1, surf%ns 
     1734                i   = surf%i(m)           
     1735                j   = surf%j(m)
     1736                k   = surf%k(m)
     1737                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
     1738             ENDDO
     1739          ENDIF
     1740!
     1741!--       At vertical surfaces
     1742          IF ( surf_vertical )  THEN
     1743             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1744             DO  m = 1, surf%ns 
     1745                i          =  surf%i(m)           
     1746                j          =  surf%j(m)
     1747                k          =  surf%k(m)
     1748                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
    9451749             ENDDO
    9461750          ENDIF
     
    9501754!
    9511755!--    If required compute qr* and nr*
    952        IF ( cloud_physics  .AND.  microphysics_seifert )   &
    953        THEN
    954 
    955           !$OMP PARALLEL DO PRIVATE( k, z_mo )
    956           DO  i = nxlg, nxrg
    957              DO  j = nysg, nyng
    958 
    959                 k   = nzb_s_inner(j,i)
    960                 z_mo = zu(k+1) - zw(k)
    961 
    962                 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) )              &
    963                                  / ( LOG( z_mo / z0q(j,i) )                 &
    964                                      - psi_h( z_mo / ol(j,i) )              &
    965                                      + psi_h( z0q(j,i) / ol(j,i) ) )
    966 
    967                 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) )              &
    968                                  / ( LOG( z_mo / z0q(j,i) )                 &
    969                                      - psi_h( z_mo / ol(j,i) )              &
    970                                      + psi_h( z0q(j,i) / ol(j,i) ) )
    971              ENDDO
     1756       IF ( cloud_physics  .AND.  microphysics_seifert  .AND.                  &
     1757            .NOT. surf_vertical )  THEN
     1758          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1759          DO  m = 1, surf%ns   
     1760
     1761             i    = surf%i(m)           
     1762             j    = surf%j(m)
     1763             k    = surf%k(m)
     1764
     1765             z_mo = surf%z_mo(m)
     1766
     1767             surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) )                 &
     1768                                 / ( LOG( z_mo / surf%z0q(m) )                 &
     1769                                 - psi_h( z_mo / surf%ol(m) )                  &
     1770                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1771
     1772             surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) )                 &
     1773                                 / ( LOG( z_mo / surf%z0q(m) )                 &
     1774                                 - psi_h( z_mo / surf%ol(m) )                  &
     1775                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
    9721776          ENDDO
    9731777
     
    9841788       IMPLICIT NONE
    9851789
    986        REAL(wp) :: ol_mid !< Grid-interpolated L
    987 
    988 !
    989 !--    Compute u'w' for the total model domain.
    990 !--    First compute the corresponding component of u* and square it.
    991        !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
    992        DO  i = nxl, nxr
    993           DO  j = nys, nyn
    994 
    995              k   = nzb_u_inner(j,i)
    996              z_mo = zu(k+1) - zw(k)
    997 !
    998 !--          Compute bulk Obukhov length for this point
    999              ol_mid = 0.5_wp * ( ol(j,i-1) + ol(j,i) )
    1000 
    1001              IF ( ol_mid == 0.0_wp )  THEN
    1002                 ol_mid = MIN(ol(j,i-1), ol(j,i))
     1790       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
     1791
     1792       REAL(wp)                            ::  dum     !< dummy to precalculate logarithm
     1793       REAL(wp)                            ::  flag_u  !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
     1794       REAL(wp)                            ::  flag_v  !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
     1795       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_i     !< u-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
     1796       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_i     !< v-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
     1797       REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_i     !< w-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
     1798
     1799!
     1800!--    Calcuate surface fluxes at horizontal walls
     1801       IF ( .NOT. surf_vertical )  THEN
     1802!
     1803!--       Compute u'w' for the total model domain at upward-facing surfaces.
     1804!--       First compute the corresponding component of u* and square it.
     1805          IF ( .NOT. downward )  THEN
     1806             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1807             DO  m = 1, surf%ns 
     1808   
     1809                i = surf%i(m)           
     1810                j = surf%j(m)
     1811                k = surf%k(m)
     1812
     1813                z_mo = surf%z_mo(m)
     1814
     1815                surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) )               &
     1816                              / ( LOG( z_mo / surf%z0(m) )                     &
     1817                                  - psi_m( z_mo / surf%ol(m) )                 &
     1818                                  + psi_m( surf%z0(m) / surf%ol(m) ) )
     1819!
     1820!--             Please note, the computation of usws is not fully accurate. Actually
     1821!--             a further interpolation of us onto the u-grid, where usws is defined,
     1822!--             is required. However, this is not done as this would require several
     1823!--             data transfers between 2D-grid and the surf-type.
     1824!--             The impact of the missing interpolation is negligible as several
     1825!--             tests had shown.
     1826!--             Same also for ol. 
     1827                surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)
     1828
     1829             ENDDO
     1830!
     1831!--       At downward-facing surfaces
     1832          ELSE
     1833             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1834             DO  m = 1, surf%ns 
     1835   
     1836                i = surf%i(m)           
     1837                j = surf%j(m)
     1838                k = surf%k(m)
     1839
     1840                z_mo = surf%z_mo(m)
     1841
     1842                surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
     1843!
     1844!--        To Do: Is the sign correct???
     1845                surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
     1846
     1847             ENDDO     
     1848          ENDIF
     1849
     1850!
     1851!--       Compute v'w' for the total model domain.
     1852!--       First compute the corresponding component of u* and square it.
     1853!--       Upward-facing surfaces
     1854          IF ( .NOT. downward )  THEN
     1855             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1856             DO  m = 1, surf%ns 
     1857                i = surf%i(m)           
     1858                j = surf%j(m)
     1859                k = surf%k(m)
     1860
     1861                z_mo = surf%z_mo(m)
     1862
     1863                surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) )               &
     1864                           / ( LOG( z_mo / surf%z0(m) )                        &
     1865                               - psi_m( z_mo / surf%ol(m) )                    &
     1866                               + psi_m( surf%z0(m) / surf%ol(m) ) )
     1867!
     1868!--             Please note, the computation of vsws is not fully accurate. Actually
     1869!--             a further interpolation of us onto the v-grid, where vsws is defined,
     1870!--             is required. However, this is not done as this would require several
     1871!--             data transfers between 2D-grid and the surf-type.
     1872!--             The impact of the missing interpolation is negligible as several
     1873!--             tests had shown.
     1874!--             Same also for ol. 
     1875                surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)
     1876             ENDDO
     1877!
     1878!--       Downward-facing surfaces
     1879          ELSE
     1880             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1881             DO  m = 1, surf%ns 
     1882                i = surf%i(m)           
     1883                j = surf%j(m)
     1884                k = surf%k(m)
     1885
     1886                z_mo = surf%z_mo(m)
     1887
     1888                surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
     1889 
     1890                surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
     1891             ENDDO
     1892          ENDIF
     1893!
     1894!--       Compute the vertical kinematic heat flux
     1895          IF (  .NOT.  constant_heatflux  .AND.  ( simulated_time <=           &
     1896               skip_time_do_lsm  .OR.  .NOT.  land_surface )  .AND.            &
     1897               .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
     1898             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1899             DO  m = 1, surf%ns
     1900                i    = surf%i(m)           
     1901                j    = surf%j(m)
     1902                k    = surf%k(m)
     1903                surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
     1904             ENDDO
     1905          ENDIF
     1906!
     1907!--       Compute the vertical water flux
     1908          IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
     1909                ( simulated_time <= skip_time_do_lsm                           &
     1910               .OR.  .NOT.  land_surface )  .AND.  .NOT. downward )  THEN
     1911             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1912             DO  m = 1, surf%ns
     1913                i    = surf%i(m)           
     1914                j    = surf%j(m)
     1915                k    = surf%k(m)
     1916                surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
     1917             ENDDO
     1918          ENDIF
     1919!
     1920!--       Compute the vertical scalar flux
     1921          IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.       &
     1922                .NOT.  downward )  THEN
     1923             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1924             DO  m = 1, surf%ns   
     1925
     1926                i    = surf%i(m)           
     1927                j    = surf%j(m)
     1928                k    = surf%k(m)
     1929                surf%ssws(m) = -surf%ss(m) * surf%us(m)
     1930
     1931             ENDDO
     1932          ENDIF       
     1933!
     1934!--       Compute (turbulent) fluxes of rain water content and rain drop conc.
     1935          IF ( cloud_physics  .AND.  microphysics_seifert  .AND.              &
     1936               .NOT. downward)  THEN
     1937             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1938             DO  m = 1, surf%ns   
     1939
     1940                i    = surf%i(m)           
     1941                j    = surf%j(m)
     1942                k    = surf%k(m)
     1943
     1944                surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
     1945                surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
     1946             ENDDO
     1947          ENDIF
     1948
     1949!
     1950!--       Bottom boundary condition for the TKE.
     1951          IF ( ibc_e_b == 2 )  THEN
     1952             !$OMP PARALLEL DO PRIVATE( i, j, k )
     1953             DO  m = 1, surf%ns   
     1954
     1955                i    = surf%i(m)           
     1956                j    = surf%j(m)
     1957                k    = surf%k(m)
     1958
     1959                e(k,j,i) = ( surf%us(m) / 0.1_wp )**2
     1960!
     1961!--             As a test: cm = 0.4
     1962!               e(k,j,i) = ( us(j,i) / 0.4_wp )**2
     1963                e(k-1,j,i)   = e(k,j,i)
     1964
     1965             ENDDO
     1966          ENDIF
     1967!
     1968!--    Calcuate surface fluxes at vertical surfaces. No stability is considered.
     1969       ELSE
     1970!
     1971!--       Compute usvs l={0,1} and vsus l={2,3}
     1972          IF ( mom_uv )  THEN
     1973!
     1974!--          Generalize computation by introducing flags. At north- and south-
     1975!--          facing surfaces u-component is used, at east- and west-facing
     1976!--          surfaces v-component is used.
     1977             flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0  .OR.  l == 1 )   
     1978             flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2  .OR.  l == 3 )   
     1979             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
     1980             DO  m = 1, surf%ns 
     1981                i = surf%i(m)           
     1982                j = surf%j(m)
     1983                k = surf%k(m)
     1984
     1985                z_mo = surf%z_mo(m)
     1986
     1987                surf%mom_flux_uv(m) = kappa *                                  &
     1988                                ( flag_u * u(k,j,i) + flag_v * v(k,j,i) )  /   &
     1989                                                        LOG( z_mo / surf%z0(m) )
     1990
     1991               surf%mom_flux_uv(m) =                                           &
     1992                                    - surf%mom_flux_uv(m) * surf%us(m)
     1993             ENDDO
     1994          ENDIF
     1995!
     1996!--       Compute wsus l={0,1} and wsvs l={2,3}
     1997          IF ( mom_w )  THEN
     1998             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
     1999             DO  m = 1, surf%ns 
     2000                i = surf%i(m)           
     2001                j = surf%j(m)
     2002                k = surf%k(m)
     2003
     2004                z_mo = surf%z_mo(m)
     2005
     2006                surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
     2007
     2008                surf%mom_flux_w(m) =                                           &
     2009                                     - surf%mom_flux_w(m) * surf%us(m)
     2010             ENDDO
     2011          ENDIF
     2012!
     2013!--       Compute momentum fluxes used for subgrid-scale TKE production at
     2014!--       vertical surfaces. In constrast to the calculated momentum fluxes at
     2015!--       vertical surfaces before, which are defined on the u/v/w-grid,
     2016!--       respectively), the TKE fluxes are defined at the scalar grid.
     2017!--       
     2018          IF ( mom_tke )  THEN
     2019!
     2020!--          Precalculate velocity components at scalar grid point.
     2021             ALLOCATE( u_i(1:surf%ns) )
     2022             ALLOCATE( v_i(1:surf%ns) )
     2023             ALLOCATE( w_i(1:surf%ns) )
     2024
     2025             IF ( l == 0  .OR.  l == 1 )  THEN
     2026                !$OMP PARALLEL DO PRIVATE( i, j, k )
     2027                DO  m = 1, surf%ns 
     2028                   i = surf%i(m)           
     2029                   j = surf%j(m)
     2030                   k = surf%k(m)
     2031
     2032                   u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
     2033                   v_i(m) = 0.0_wp
     2034                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     2035                ENDDO
     2036             ELSE
     2037                !$OMP PARALLEL DO PRIVATE( i, j, k )
     2038                DO  m = 1, surf%ns 
     2039                   i = surf%i(m)           
     2040                   j = surf%j(m)
     2041                   k = surf%k(m)
     2042
     2043                   u_i(m) = 0.0_wp
     2044                   v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
     2045                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     2046                ENDDO
    10032047             ENDIF
    10042048
    1005              usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )                     &
    1006                                  / ( LOG( z_mo / z0(j,i) )                     &
    1007                                      - psi_m( z_mo / ol_mid )                  &
    1008                                      + psi_m( z0(j,i) / ol_mid ) )
    1009 
    1010              usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )         &
    1011                                     * rho_air_zw(k)
    1012           ENDDO
    1013        ENDDO
    1014 
    1015 !
    1016 !--    Compute v'w' for the total model domain.
    1017 !--    First compute the corresponding component of u* and square it.
    1018        !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
    1019        DO  i = nxl, nxr
    1020           DO  j = nys, nyn
    1021 
    1022              k   = nzb_v_inner(j,i)
    1023              z_mo = zu(k+1) - zw(k)
    1024 !
    1025 !--          Compute bulk Obukhov length for this point
    1026              ol_mid = 0.5_wp * ( ol(j-1,i) + ol(j,i) )
    1027 
    1028              IF ( ol_mid == 0.0_wp )  THEN
    1029                 ol_mid = MIN(ol(j-1,i), ol(j-1,i))
    1030              ENDIF
    1031 
    1032              vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) )                     &
    1033                                  / ( LOG( z_mo / z0(j,i) )                     &
    1034                                      - psi_m( z_mo / ol_mid )                  &
    1035                                      + psi_m( z0(j,i) / ol_mid ) )
    1036 
    1037              vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j-1,i) + us(j,i) )         &
    1038                                     * rho_air_zw(k)
    1039 
    1040           ENDDO
    1041        ENDDO
    1042 
    1043 !
    1044 !--    Exchange the boundaries for the momentum fluxes (is this still required?)
    1045        CALL exchange_horiz_2d( usws )
    1046        CALL exchange_horiz_2d( vsws )
    1047 
    1048 !
    1049 !--    Compute the vertical kinematic heat flux
    1050        IF (  .NOT.  constant_heatflux  .AND.  ( simulated_time <=            &
    1051             skip_time_do_lsm  .OR.  .NOT.  land_surface )  .AND.             &
    1052             .NOT.  urban_surface )  THEN
    1053           !$OMP PARALLEL DO
    1054           DO  i = nxlg, nxrg
    1055              DO  j = nysg, nyng
    1056                 k   = nzb_s_inner(j,i)
    1057                 shf(j,i) = -ts(j,i) * us(j,i) * rho_air_zw(k)
    1058              ENDDO
    1059           ENDDO
    1060 
    1061        ENDIF
    1062 
    1063 !
    1064 !--    Compute the vertical water flux
    1065        IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.                 &
    1066              ( simulated_time <= skip_time_do_lsm                              &
    1067             .OR.  .NOT.  land_surface ) )  THEN
    1068           !$OMP PARALLEL DO
    1069           DO  i = nxlg, nxrg
    1070              DO  j = nysg, nyng
    1071                 k   = nzb_s_inner(j,i)
    1072                 qsws(j,i) = -qs(j,i) * us(j,i) * rho_air_zw(k)
    1073              ENDDO
    1074           ENDDO
    1075 
    1076        ENDIF
    1077        
    1078 !
    1079 !--    Compute the vertical scalar flux
    1080        IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.          &
    1081              ( simulated_time <= skip_time_do_lsm                              &
    1082             .OR.  .NOT.  land_surface ) )  THEN
    1083           !$OMP PARALLEL DO
    1084           DO  i = nxlg, nxrg
    1085              DO  j = nysg, nyng
    1086                 ssws(j,i) = -ss(j,i) * us(j,i)
    1087              ENDDO
    1088           ENDDO
    1089 
    1090        ENDIF       
    1091 
    1092 !
    1093 !--    Compute (turbulent) fluxes of rain water content and rain drop conc.
    1094        IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    1095           !$OMP PARALLEL DO
    1096           DO  i = nxlg, nxrg
    1097              DO  j = nysg, nyng
    1098                 qrsws(j,i) = -qrs(j,i) * us(j,i)
    1099                 nrsws(j,i) = -nrs(j,i) * us(j,i)
    1100              ENDDO
    1101           ENDDO
    1102        ENDIF
    1103 
    1104 !
    1105 !--    Bottom boundary condition for the TKE
    1106        IF ( ibc_e_b == 2 )  THEN
    1107           !$OMP PARALLEL DO
    1108           DO  i = nxlg, nxrg
    1109              DO  j = nysg, nyng
    1110                 k = nzb_s_inner(j,i)
    1111                 e(k+1,j,i) = ( us(j,i) / 0.1_wp )**2
    1112 !
    1113 !--             As a test: cm = 0.4
    1114 !               e(k+1,j,i) = ( us(j,i) / 0.4_wp )**2
    1115                 e(k,j,i)   = e(k+1,j,i)
    1116              ENDDO
    1117           ENDDO
     2049             !$OMP PARALLEL DO PRIVATE( i, j, k, dum, z_mo )
     2050             DO  m = 1, surf%ns 
     2051                i = surf%i(m)           
     2052                j = surf%j(m)
     2053                k = surf%k(m)
     2054
     2055                z_mo = surf%z_mo(m)
     2056
     2057                dum = kappa / LOG( z_mo / surf%z0(m) )
     2058!
     2059!--             usvs (l=0,1) and vsus (l=2,3)
     2060                surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )
     2061!
     2062!--             wsvs (l=0,1) and wsus (l=2,3)
     2063                surf%mom_flux_tke(1,m) = dum * w_i(m)
     2064
     2065                surf%mom_flux_tke(0:1,m) =                                     &
     2066                               - surf%mom_flux_tke(0:1,m) * surf%us(m)
     2067             ENDDO
     2068!
     2069!--          Deallocate temporary arrays
     2070             DEALLOCATE( u_i )             
     2071             DEALLOCATE( v_i ) 
     2072             DEALLOCATE( w_i ) 
     2073          ENDIF
    11182074       ENDIF
    11192075
Note: See TracChangeset for help on using the changeset viewer.