Changeset 3120 for palm/trunk


Ignore:
Timestamp:
Jul 11, 2018 6:30:57 PM (6 years ago)
Author:
gronemeier
Message:

implementation of dynamic sgs model

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

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

    r3083 r3120  
    2525! -----------------
    2626! $Id$
     27! +les_dynamic
     28!
     29! 3083 2018-06-19 14:03:12Z gronemeier
    2730! set dt_3d = 0.01
    2831!
     
    13061309    LOGICAL ::  large_scale_subsidence = .FALSE.                 !< namelist parameter
    13071310    LOGICAL ::  land_surface = .FALSE.                           !< use land surface model?
     1311    LOGICAL ::  les_dynamic = .FALSE.                            !< use dynamic subgrid model as turbulence closure for LES mode
    13081312    LOGICAL ::  les_mw = .FALSE.                                 !< use Moeng-Wyngaard turbulence closure for LES mode
    13091313    LOGICAL ::  lsf_exception = .FALSE.                          !< use of lsf with buildings (temporary)?
  • palm/trunk/SOURCE/timestep.f90

    r3084 r3120  
    2525! -----------------
    2626! $Id$
     27! Put ABS( km ) in computation of time step according to the diffusion criterion
     28!
     29! 3084 2018-06-19 15:30:55Z gronemeier
    2730! limit increase of dt_3d only in case of RANS mode
    2831!
     
    302305!
    303306!--    Compute time step according to the diffusion criterion.
    304 !--    First calculate minimum grid spacing which only depends on index k
     307!--    First calculate minimum grid spacing which only depends on index k.
     308!--    When using the dynamic subgrid model, negative km are possible.
    305309       dt_diff_l = 999999.0_wp
    306310
     
    315319             DO  k = nzb+1, nzt
    316320                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) /                     &
    317                                     ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
     321                                    ( MAX( kh(k,j,i), ABS( km(k,j,i) ) )       &
     322                                      + 1E-20_wp ) )
    318323             ENDDO
    319324          ENDDO
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r3086 r3120  
    2525! -----------------
    2626! $Id$
     27! - changed tcm_diffusivities to tcm_diffusivities_default
     28! - created subroutine tcm_diffusivities that calls tcm_diffusivities_default
     29!   and tcm_diffusivities_dynamic
     30!
     31! 3086 2018-06-25 09:08:04Z gronemeier
    2732! bugfix: set rans_const_sigma(1) = 1.3
    2833!
     
    95100! --------
    96101! @author Tobias Gronemeier
    97 !
     102! @author Hauke Wurps
    98103!
    99104! Description:
     
    126131        ONLY:  constant_diffusion, dt_3d, e_init, humidity, inflow_l,          &
    127132               initializing_actions, intermediate_timestep_count,              &
    128                intermediate_timestep_count_max, kappa, km_constant, les_mw,    &
    129                ocean, plant_canopy, prandtl_number, prho_reference,            &
    130                pt_reference, rans_mode, rans_tke_e, rans_tke_l, simulated_time,&
    131                timestep_scheme, turbulence_closure, turbulent_inflow,          &
    132                use_upstream_for_tke, vpt_reference, ws_scheme_sca
     133               intermediate_timestep_count_max, kappa, km_constant,            &
     134               les_dynamic, les_mw, ocean, plant_canopy, prandtl_number,       &
     135               prho_reference, pt_reference, rans_mode, rans_tke_e, rans_tke_l,&
     136               simulated_time,timestep_scheme, turbulence_closure,             &
     137               turbulent_inflow, use_upstream_for_tke, vpt_reference,          &
     138               ws_scheme_sca
    133139
    134140    USE advec_ws,                                                              &
     
    148154
    149155    USE indices,                                                               &
    150         ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg,               &
    151                nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt,   &
     156        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,     &
    152157               wall_flags_0
    153158
     
    296301
    297302!
    298 !-- Calculate diffusivities
     303!-- Call tcm_diffusivities_default and tcm_diffusivities_dynamic
    299304    INTERFACE tcm_diffusivities
    300305       MODULE PROCEDURE tcm_diffusivities
    301306    END INTERFACE tcm_diffusivities
     307
     308!
     309!-- Calculate diffusivities
     310    INTERFACE tcm_diffusivities_default
     311       MODULE PROCEDURE tcm_diffusivities_default
     312    END INTERFACE tcm_diffusivities_default
     313
     314!
     315!-- Calculate diffusivities according to dynamic sgs model
     316    INTERFACE tcm_diffusivities_dynamic
     317       MODULE PROCEDURE tcm_diffusivities_dynamic
     318    END INTERFACE tcm_diffusivities_dynamic
     319
     320!
     321!-- Box-filter method for dynamic sgs model
     322    INTERFACE tcm_box_filter_2d
     323       MODULE PROCEDURE tcm_box_filter_2d_single
     324       MODULE PROCEDURE tcm_box_filter_2d_array
     325    END INTERFACE tcm_box_filter_2d
    302326
    303327!
     
    390414          CASE ( 'Moeng_Wyngaard' )
    391415             les_mw = .TRUE.
     416
     417          CASE ( 'dynamic' )
     418             les_dynamic = .TRUE.
    392419
    393420          CASE DEFAULT
     
    19451972    IMPLICIT NONE
    19461973
    1947     INTEGER(iwp) ::  i   !< grid index x-direction
    1948     INTEGER(iwp) ::  j   !< grid index y-direction
    1949     INTEGER(iwp) ::  k   !< grid index z-direction
    1950     INTEGER(iwp) ::  m   !< running index surface elements
     1974    INTEGER(iwp) ::  i      !< grid index x-direction
     1975    INTEGER(iwp) ::  j      !< grid index y-direction
     1976    INTEGER(iwp) ::  k      !< grid index z-direction
     1977    INTEGER(iwp) ::  m      !< running index surface elements
    19511978
    19521979    IF ( constant_flux_layer )  THEN
     
    19732000!--       and km are not on the same grid. Actually, a further
    19742001!--       interpolation of km onto the u/v-grid is necessary. However, the
    1975 !--       effect of this error is negligible. 
     2002!--       effect of this error is negligible.
    19762003          surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *          &
    1977                                      drho_air_zw(k-1) *                        &
    1978                                      ( zu(k+1)    - zu(k-1)    )  /            &
    1979                                      ( km(k,j,i)  + 1.0E-20_wp ) 
     2004                                     drho_air_zw(k-1)               *          &
     2005                                     ( zu(k+1) - zu(k-1)    )       /          &
     2006                                     ( km(k,j,i)  + 1.0E-20_wp )
    19802007          surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *          &
    1981                                      drho_air_zw(k-1) *                        &
    1982                                      ( zu(k+1)    - zu(k-1)    )  /            &
    1983                                      ( km(k,j,i)  + 1.0E-20_wp ) 
     2008                                     drho_air_zw(k-1)               *          &
     2009                                     ( zu(k+1) - zu(k-1)    )       /          &
     2010                                     ( km(k,j,i)  + 1.0E-20_wp )   
    19842011
    19852012          IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
     
    20322059!--       interpolation of km onto the u/v-grid is necessary. However, the
    20332060!--       effect of this error is negligible.
    2034           surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)      *           &
    2035                                         drho_air_zw(k-1) *                     &
    2036                                         ( zu(k+1)   - zu(k-1)    )  /          &
    2037                                         ( km(k,j,i) + 1.0E-20_wp ) 
    2038           surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)      *           &
    2039                                         drho_air_zw(k-1) *                     &
    2040                                         ( zu(k+1)   - zu(k-1)    )  /          &
    2041                                         ( km(k,j,i) + 1.0E-20_wp )
     2061!--       effect of this error is negligible.
     2062          surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)    *             &
     2063                                        drho_air_zw(k-1)         *             &
     2064                                        ( zu(k+1) - zu(k-1)    ) /             &
     2065                                        ( km(k,j,i)  + 1.0E-20_wp ) 
     2066          surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)    *             &
     2067                                        drho_air_zw(k-1)         *             &
     2068                                        ( zu(k+1) - zu(k-1)    ) /             &
     2069                                        ( km(k,j,i)  + 1.0E-20_wp )
    20422070
    20432071          IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
     
    20632091!--       interpolation of km onto the u/v-grid is necessary. However, the
    20642092!--       effect of this error is negligible.
    2065           surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)      *           &
    2066                                         drho_air_zw(k-1) *                     &
    2067                                         ( zu(k+1)   - zu(k-1)    )  /          &
    2068                                         ( km(k,j,i) + 1.0E-20_wp ) 
    2069           surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)      *           &
    2070                                         drho_air_zw(k-1) *                     &
    2071                                         ( zu(k+1)   - zu(k-1)    )  /          &
    2072                                         ( km(k,j,i) + 1.0E-20_wp )
     2093          surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)               &
     2094                                        drho_air_zw(k-1)         *             &
     2095                                        ( zu(k+1) - zu(k-1)    ) /             &
     2096                                        ( km(k,j,i)  + 1.0E-20_wp )
     2097          surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)               &
     2098                                        drho_air_zw(k-1)         *             &
     2099                                        ( zu(k+1) - zu(k-1)    ) /             &
     2100                                        ( km(k,j,i)  + 1.0E-20_wp )
    20732101
    20742102          IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
     
    42304258!
    42314259!--          Calculate dissipation
    4232              IF ( les_mw )  THEN
     4260             IF ( les_dynamic .OR. les_mw )  THEN
    42334261
    42344262                CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
     
    43674395!--    Calculate dissipation...
    43684396!--    ...in case of LES
    4369        IF ( les_mw )  THEN
     4397       IF ( les_dynamic .OR. les_mw )  THEN
    43704398
    43714399          CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
     
    46784706! Description:
    46794707! ------------
     4708!> Computation of the turbulent diffusion coefficients for momentum and heat.
     4709!------------------------------------------------------------------------------!
     4710 SUBROUTINE tcm_diffusivities( var, var_reference )
     4711
     4712    REAL(wp) ::  var_reference       !< reference temperature
     4713
     4714#if defined( __nopointer )
     4715    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
     4716#else
     4717    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
     4718#endif
     4719!
     4720!-- Call default diffusivities routine. This is always used to calculate kh.
     4721    CALL tcm_diffusivities_default( var, var_reference )
     4722!
     4723!-- Call dynamic subgrid model to calculate km.
     4724    IF ( les_dynamic )  THEN
     4725       CALL tcm_diffusivities_dynamic
     4726    ENDIF
     4727
     4728 END SUBROUTINE tcm_diffusivities
     4729
     4730
     4731!------------------------------------------------------------------------------!
     4732! Description:
     4733! ------------
    46804734!> Computation of the turbulent diffusion coefficients for momentum and heat
    46814735!> according to Prandtl-Kolmogorov.
    46824736!> @todo consider non-default surfaces
    46834737!------------------------------------------------------------------------------!
    4684  SUBROUTINE tcm_diffusivities( var, var_reference )
     4738 SUBROUTINE tcm_diffusivities_default( var, var_reference )
    46854739 
    46864740
     
    47464800    ENDIF
    47474801
    4748     IF ( les_mw )  THEN
     4802    IF ( les_dynamic .OR. les_mw )  THEN
    47494803       !$OMP DO
    47504804       DO  i = nxlg, nxrg
     
    49264980    ENDIF
    49274981
    4928  END SUBROUTINE tcm_diffusivities
     4982 END SUBROUTINE tcm_diffusivities_default
     4983
     4984
     4985!------------------------------------------------------------------------------!
     4986! Description:
     4987! ------------
     4988!> Calculates the eddy viscosity dynamically using the linear dynamic model
     4989!> according to
     4990!> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via
     4991!> stochastic analysis."
     4992!> Monte Carlo Methods and Applications 14.4 (2008): 311-329.
     4993!>
     4994!> Furthermore dynamic bounds are used to limit the absolute value of c* as
     4995!> described in
     4996!> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation:
     4997!> Stability via realizability." Physics of Fluids 29.10 (2017): 105104.
     4998!>
     4999!> @author Hauke Wurps
     5000!------------------------------------------------------------------------------!
     5001 SUBROUTINE tcm_diffusivities_dynamic
     5002       
     5003    USE arrays_3d,                                                             &
     5004        ONLY:  ddzw, dzw, dd2zu, w, ug, vg
     5005
     5006    USE control_parameters,                                                    &
     5007        ONLY:  e_min, outflow_l, outflow_n, outflow_r, outflow_s
     5008   
     5009    USE grid_variables,                                                        &
     5010        ONLY : ddx, ddy, dx, dy
     5011
     5012    USE surface_mod,                                                           &
     5013        ONLY :  bc_h
     5014
     5015    IMPLICIT NONE
     5016
     5017    INTEGER(iwp) ::  i           !< running index x-direction
     5018    INTEGER(iwp) ::  j           !< running index y-direction
     5019    INTEGER(iwp) ::  k           !< running index z-direction
     5020    INTEGER(iwp) ::  l           !< running index
     5021    INTEGER(iwp) ::  m           !< running index
     5022    INTEGER(iwp) ::  n           !< running index
     5023
     5024    REAL(wp)     ::  dudx        !< Gradient of u-component in x-direction
     5025    REAL(wp)     ::  dudy        !< Gradient of u-component in y-direction
     5026    REAL(wp)     ::  dudz        !< Gradient of u-component in z-direction
     5027    REAL(wp)     ::  dvdx        !< Gradient of v-component in x-direction
     5028    REAL(wp)     ::  dvdy        !< Gradient of v-component in y-direction
     5029    REAL(wp)     ::  dvdz        !< Gradient of v-component in z-direction
     5030    REAL(wp)     ::  dwdx        !< Gradient of w-component in x-direction
     5031    REAL(wp)     ::  dwdy        !< Gradient of w-component in y-direction
     5032    REAL(wp)     ::  dwdz        !< Gradient of w-component in z-direction
     5033
     5034    REAL(wp)     ::  uc(-1:1,-1:1)  !< u on grid center
     5035    REAL(wp)     ::  vc(-1:1,-1:1)  !< v on grid center
     5036    REAL(wp)     ::  wc(-1:1,-1:1)  !< w on grid center
     5037    REAL(wp)     ::  u2(-1:1,-1:1)  !< u2 on grid center
     5038    REAL(wp)     ::  v2(-1:1,-1:1)  !< v2 on grid center
     5039    REAL(wp)     ::  w2(-1:1,-1:1)  !< w2 on grid center
     5040    REAL(wp)     ::  uv(-1:1,-1:1)  !< u*v on grid center
     5041    REAL(wp)     ::  uw(-1:1,-1:1)  !< u*w on grid center
     5042    REAL(wp)     ::  vw(-1:1,-1:1)  !< v*w on grid center
     5043
     5044    REAL(wp)     ::  ut(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered u
     5045    REAL(wp)     ::  vt(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered v
     5046    REAL(wp)     ::  wt(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered w
     5047
     5048    REAL(wp)     ::  uct         !< test filtered u on grid center
     5049    REAL(wp)     ::  vct         !< test filtered v on grid center
     5050    REAL(wp)     ::  wct         !< test filtered w on grid center
     5051    REAL(wp)     ::  u2t         !< test filtered u**2 on grid center
     5052    REAL(wp)     ::  v2t         !< test filtered v**2 on grid center
     5053    REAL(wp)     ::  w2t         !< test filtered w**2 on grid center
     5054    REAL(wp)     ::  uvt         !< test filtered u*v on grid center
     5055    REAL(wp)     ::  uwt         !< test filtered u*w on grid center
     5056    REAL(wp)     ::  vwt         !< test filtered v*w on grid center
     5057
     5058    REAL(wp)     ::  sd11        !< deviatoric shear tensor
     5059    REAL(wp)     ::  sd22        !< deviatoric shear tensor
     5060    REAL(wp)     ::  sd33        !< deviatoric shear tensor
     5061    REAL(wp)     ::  sd12        !< deviatoric shear tensor
     5062    REAL(wp)     ::  sd13        !< deviatoric shear tensor
     5063    REAL(wp)     ::  sd23        !< deviatoric shear tensor
     5064
     5065    REAL(wp)     ::  sd2         !< sum: sd_ij*sd_ij
     5066
     5067    REAL(wp)     ::  sdt11       !< filtered deviatoric shear tensor
     5068    REAL(wp)     ::  sdt22       !< filtered deviatoric shear tensor
     5069    REAL(wp)     ::  sdt33       !< filtered deviatoric shear tensor
     5070    REAL(wp)     ::  sdt12       !< filtered deviatoric shear tensor
     5071    REAL(wp)     ::  sdt13       !< filtered deviatoric shear tensor
     5072    REAL(wp)     ::  sdt23       !< filtered deviatoric shear tensor
     5073
     5074    REAL(wp)     ::  sdt2        !< sum: sdt_ij*sdt_ij
     5075
     5076    REAL(wp)     ::  ld11        !< deviatoric stress tensor
     5077    REAL(wp)     ::  ld22        !< deviatoric stress tensor
     5078    REAL(wp)     ::  ld33        !< deviatoric stress tensor
     5079    REAL(wp)     ::  ld12        !< deviatoric stress tensor
     5080    REAL(wp)     ::  ld13        !< deviatoric stress tensor
     5081    REAL(wp)     ::  ld23        !< deviatoric stress tensor
     5082
     5083    REAL(wp)     ::  lnn         !< sum ld_nn
     5084    REAL(wp)     ::  ldsd        !< sum: ld_ij*sd_ij
     5085
     5086    REAL(wp)     ::  delta       !< grid size
     5087    REAL(wp)     ::  cst         !< c*
     5088    REAL(wp)     ::  cstnust_t   !< product c*nu*
     5089    REAL(wp)     ::  cst_max     !< bounds of c*
     5090
     5091    REAL(wp), PARAMETER :: fac_cmax = 23.0_wp/(24.0_wp*sqrt(3.0_wp))  !< constant
     5092
     5093!
     5094!-- Introduce an optional minimum tke
     5095    IF ( e_min > 0.0_wp )  THEN
     5096       !$OMP DO
     5097       DO  i = nxlg, nxrg
     5098          DO  j = nysg, nyng
     5099             DO  k = nzb+1, nzt
     5100                e(k,j,i) = MAX( e(k,j,i), e_min ) *                            &
     5101                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     5102             ENDDO
     5103          ENDDO
     5104       ENDDO
     5105    ENDIF
     5106!
     5107!-- velocities on grid centers:
     5108    CALL tcm_box_filter_2d( u, ut )
     5109    CALL tcm_box_filter_2d( v, vt )
     5110    CALL tcm_box_filter_2d( w, wt )
     5111
     5112    DO  i = nxl, nxr
     5113       DO  j = nys, nyn
     5114          DO  k = nzb+1, nzt
     5115!
     5116!--          Compute the deviatoric shear tensor s_ij on grid centers:
     5117!--          s_ij =  0.5 * ( du_i/dx_j + du_j/dx_i )
     5118             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     5119             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     5120                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     5121             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     5122                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     5123
     5124             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     5125                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     5126             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     5127             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     5128                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     5129
     5130             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     5131                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     5132             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     5133                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     5134             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     5135
     5136             sd11 = dudx
     5137             sd22 = dvdy
     5138             sd33 = dwdz
     5139             sd12 = 0.5_wp * ( dudy + dvdx )
     5140             sd13 = 0.5_wp * ( dudz + dwdx )
     5141             sd23 = 0.5_wp * ( dvdz + dwdy )
     5142!
     5143!--          sum: sd_ij*sd_ij
     5144             sd2 = sd11**2 + sd22**2 + sd33**2                     &
     5145                   + 2.0_wp * ( sd12**2 + sd13**2 + sd23**2 )
     5146!
     5147!--          The filtered velocities are needed to calculate the filtered shear
     5148!--          tensor: sdt_ij = 0.5 * ( dut_i/dx_j + dut_j/dx_i )
     5149             dudx  =           ( ut(k,j,i+1) - ut(k,j,i)     ) * ddx
     5150             dudy  = 0.25_wp * ( ut(k,j+1,i) + ut(k,j+1,i+1) - &
     5151                                 ut(k,j-1,i) - ut(k,j-1,i+1) ) * ddy
     5152             dudz  = 0.5_wp  * ( ut(k+1,j,i) + ut(k+1,j,i+1) - &
     5153                                 ut(k-1,j,i) - ut(k-1,j,i+1) ) * dd2zu(k)
     5154
     5155             dvdx  = 0.25_wp * ( vt(k,j,i+1) + vt(k,j+1,i+1) - &
     5156                                 vt(k,j,i-1) - vt(k,j+1,i-1) ) * ddx
     5157             dvdy  =           ( vt(k,j+1,i) - vt(k,j,i)     ) * ddy
     5158             dvdz  = 0.5_wp  * ( vt(k+1,j,i) + vt(k+1,j+1,i) - &
     5159                                 vt(k-1,j,i) - vt(k-1,j+1,i) ) * dd2zu(k)
     5160
     5161             dwdx  = 0.25_wp * ( wt(k,j,i+1) + wt(k-1,j,i+1) - &
     5162                                 wt(k,j,i-1) - wt(k-1,j,i-1) ) * ddx
     5163             dwdy  = 0.25_wp * ( wt(k,j+1,i) + wt(k-1,j+1,i) - &
     5164                                 wt(k,j-1,i) - wt(k-1,j-1,i) ) * ddy
     5165             dwdz  =           ( wt(k,j,i)   - wt(k-1,j,i)   ) * ddzw(k)
     5166
     5167             sdt11 = dudx
     5168             sdt22 = dvdy
     5169             sdt33 = dwdz
     5170             sdt12 = 0.5_wp * ( dudy + dvdx )
     5171             sdt13 = 0.5_wp * ( dudz + dwdx )
     5172             sdt23 = 0.5_wp * ( dvdz + dwdy )
     5173!
     5174!--          sum: sd_ij*sd_ij
     5175             sdt2 = sdt11**2 + sdt22**2 + sdt33**2         &
     5176                   + 2.0_wp * ( sdt12**2 + sdt13**2 + sdt23**2 )
     5177!
     5178!--          Need filtered velocities and filtered squared velocities on grid
     5179!--          centers. Substraction of geostrophic velocity helps to avoid
     5180!--          numerical errors in the expression <u**2> - <u>*<u>, which can be
     5181!--          very small (<...> means filtered).
     5182             DO  l = -1, 1
     5183                DO  m = -1, 1
     5184                   uc(l,m) = 0.5_wp * ( u(k,j+l,i+m) + u(k,j+l,i+m+1) ) - ug(k)
     5185                   vc(l,m) = 0.5_wp * ( v(k,j+l,i+m) + v(k,j+l+1,i+m) ) - vg(k)
     5186                   wc(l,m) = 0.5_wp * ( w(k-1,j+l,i+m) + w(k,j+l,i+m) )
     5187                ENDDO
     5188             ENDDO
     5189
     5190             CALL tcm_box_filter_2d( uc, uct )
     5191             CALL tcm_box_filter_2d( vc, vct )
     5192             CALL tcm_box_filter_2d( wc, wct )
     5193             CALL tcm_box_filter_2d( uc**2, u2t )
     5194             CALL tcm_box_filter_2d( vc**2, v2t )
     5195             CALL tcm_box_filter_2d( wc**2, w2t )
     5196             CALL tcm_box_filter_2d( uc*vc, uvt )
     5197             CALL tcm_box_filter_2d( uc*wc, uwt )
     5198             CALL tcm_box_filter_2d( vc*wc, vwt )
     5199
     5200             ld11 = u2t - uct*uct
     5201             ld22 = v2t - vct*vct
     5202             ld33 = w2t - wct*wct
     5203             ld12 = uvt - uct*vct
     5204             ld13 = uwt - uct*wct
     5205             ld23 = vwt - vct*wct
     5206
     5207             lnn = ld11 + ld22 + ld33
     5208!
     5209!--          Substract trace to get deviatoric resolved stress
     5210             ld11 = ld11 - lnn / 3.0_wp
     5211             ld22 = ld22 - lnn / 3.0_wp
     5212             ld33 = ld33 - lnn / 3.0_wp
     5213
     5214             ldsd = ld11 * sdt11 + ld22 * sdt22 + ld33 * sdt33 + &
     5215                    2.0_wp * ( ld12 * sdt12 + ld13 * sdt13 + ld23 * sdt23 )
     5216!
     5217!--          c* nu*^T is SGS viscosity on test filter level:
     5218             cstnust_t = -ldsd / sdt2
     5219!
     5220!--          The model was only tested for an isotropic grid. The following
     5221!--          expression was a recommendation of Stefan Heinz.
     5222             delta = MAX( dx, dy, dzw(k) )
     5223             cst = cstnust_t / ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) )
     5224!
     5225!--          Calculate border according to Mokhtarpoor and Heinz (2017)
     5226             cst_max = fac_cmax * SQRT( e(k,j,i) ) / ( delta * SQRT( 2.0_wp * sd2 ) )
     5227
     5228             IF ( ABS( cst ) > cst_max )  THEN
     5229                cst = cst_max * cst / ABS( cst )
     5230             ENDIF
     5231
     5232             km(k,j,i) = cst * delta * SQRT( e(k,j,i) )
     5233
     5234          ENDDO
     5235       ENDDO
     5236    ENDDO
     5237!
     5238!-- Set vertical boundary values (Neumann conditions both at upward- and
     5239!-- downward facing walls. To set wall-boundary values, the surface data type
     5240!-- is applied.
     5241!-- Horizontal boundary conditions at vertical walls are not set because
     5242!-- so far vertical surfaces require usage of a constant-flux layer where the
     5243!-- boundary values of the diffusivities are not needed.
     5244
     5245!
     5246!-- Upward-facing
     5247    !$OMP PARALLEL DO PRIVATE( i, j, k )
     5248    DO  m = 1, bc_h(0)%ns
     5249        i = bc_h(0)%i(m)           
     5250        j = bc_h(0)%j(m)
     5251        k = bc_h(0)%k(m)
     5252        km(k-1,j,i) = km(k,j,i)
     5253    ENDDO
     5254!
     5255!-- Downward facing surfaces
     5256    !$OMP PARALLEL DO PRIVATE( i, j, k )
     5257    DO  m = 1, bc_h(1)%ns
     5258        i = bc_h(1)%i(m)           
     5259        j = bc_h(1)%j(m)
     5260        k = bc_h(1)%k(m)
     5261        km(k+1,j,i) = km(k,j,i)
     5262    ENDDO
     5263!
     5264!-- Model top
     5265    !$OMP PARALLEL DO
     5266    DO  i = nxlg, nxrg
     5267       DO  j = nysg, nyng
     5268          km(nzt+1,j,i) = km(nzt,j,i)
     5269       ENDDO
     5270    ENDDO
     5271!
     5272!-- Set Neumann boundary conditions at the outflow boundaries in case of
     5273!-- non-cyclic lateral boundaries
     5274    IF ( outflow_l )  THEN
     5275       km(:,:,nxl-1) = km(:,:,nxl)
     5276    ENDIF
     5277    IF ( outflow_r )  THEN
     5278       km(:,:,nxr+1) = km(:,:,nxr)
     5279    ENDIF
     5280    IF ( outflow_s )  THEN
     5281       km(:,nys-1,:) = km(:,nys,:)
     5282    ENDIF
     5283    IF ( outflow_n )  THEN
     5284       km(:,nyn+1,:) = km(:,nyn,:)
     5285    ENDIF
     5286
     5287 END SUBROUTINE tcm_diffusivities_dynamic
     5288
     5289
     5290!------------------------------------------------------------------------------!
     5291! Description:
     5292! ------------
     5293!> This subroutine acts as a box filter with filter width 2 * dx.
     5294!> Output is only one point.
     5295!------------------------------------------------------------------------------!
     5296 SUBROUTINE tcm_box_filter_2d_single( var, var_fil )
     5297 
     5298    IMPLICIT NONE
     5299
     5300    REAL(wp)     ::  var(-1:1,-1:1)      !< variable to be filtered
     5301    REAL(wp)     ::  var_fil             !< filtered variable
     5302!
     5303!-- It is assumed that a box with a side length of 2 * dx and centered at the
     5304!-- variable*s position contains one half of the four closest neigbours and one
     5305!-- forth of the diagonally closest neighbours.
     5306    var_fil = 0.25_wp * ( var(0,0) +                                           &
     5307                      0.5_wp  * ( var(0,1)  + var(1,0)   +                     &
     5308                                  var(0,-1) + var(-1,0)  ) +                   &
     5309                      0.25_wp * ( var(1,1)  + var(1,-1)  +                     &
     5310                                  var(-1,1) + var(-1,-1) ) )
     5311
     5312 END SUBROUTINE tcm_box_filter_2d_single
     5313
     5314!------------------------------------------------------------------------------!
     5315! Description:
     5316! ------------
     5317!> This subroutine acts as a box filter with filter width 2 * dx.
     5318!> The filtered variable var_fil is on the same grid as var.
     5319!------------------------------------------------------------------------------!
     5320 SUBROUTINE tcm_box_filter_2d_array( var, var_fil )
     5321
     5322    IMPLICIT NONE
     5323
     5324    INTEGER(iwp) ::  i    !< running index x-direction
     5325    INTEGER(iwp) ::  j    !< running index y-direction
     5326    INTEGER(iwp) ::  k    !< running index z-direction
     5327
     5328    REAL(wp)     ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)      !< variable to be filtered
     5329    REAL(wp)     ::  var_fil(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< filtered variable
     5330!
     5331!-- It is assumed that a box with a side length of 2 * dx and centered at the
     5332!-- variable's position contains one half of the four closest neigbours and one
     5333!-- forth of the diagonally closest neighbours.
     5334    DO  i = nxlg+1, nxrg-1
     5335       DO  j = nysg+1, nyng-1
     5336          DO  k = nzb, nzt+1
     5337             var_fil(k,j,i) = 0.25_wp * ( var(k,j,i) +                         &
     5338                              0.5_wp  * ( var(k,j,i+1)   + var(k,j+1,i)   +    &
     5339                                          var(k,j,i-1)   + var(k,j-1,i)     ) +&
     5340                              0.25_wp * ( var(k,j+1,i+1) + var(k,j+1,i-1) +    &
     5341                                          var(k,j-1,i+1) + var(k,j-1,i-1)   )  )
     5342          END DO
     5343       END DO
     5344    END DO
     5345
     5346 END SUBROUTINE tcm_box_filter_2d_array
    49295347
    49305348
Note: See TracChangeset for help on using the changeset viewer.