Ignore:
Timestamp:
Oct 25, 2018 1:36:23 PM (6 years ago)
Author:
maronga
Message:

added support for buildings with dynamic sgs model

File:
1 edited

Legend:

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

    r3398 r3430  
    2525! -----------------
    2626! $Id$
     27! Added support for buildings in the dynamic SGS model
     28!
     29! 3398 2018-10-22 19:30:24Z knoop
    2730! Refactored production_e and production_e_ij (removed code redundancy)
    2831!
     
    186189               simulated_time,timestep_scheme, turbulence_closure,             &
    187190               turbulent_inflow, use_upstream_for_tke, vpt_reference,          &
    188                ws_scheme_sca
     191               ws_scheme_sca, current_timestep_number
    189192
    190193    USE advec_ws,                                                              &
     
    46874690!>
    46884691!> @author Hauke Wurps
     4692!> @author Björn Maronga
    46894693!------------------------------------------------------------------------------!
    46904694 SUBROUTINE tcm_diffusivities_dynamic
     
    47144718    REAL(wp)     ::  dwdz        !< Gradient of w-component in z-direction
    47154719
     4720    REAL(wp)     ::  flag                !< topography flag
     4721   
    47164722    REAL(wp)     ::  uc(-1:1,-1:1)  !< u on grid center
    47174723    REAL(wp)     ::  vc(-1:1,-1:1)  !< v on grid center
    47184724    REAL(wp)     ::  wc(-1:1,-1:1)  !< w on grid center
    4719 !    REAL(wp)     ::  u2(-1:1,-1:1)  !< u2 on grid center
    4720 !    REAL(wp)     ::  v2(-1:1,-1:1)  !< v2 on grid center
    4721 !    REAL(wp)     ::  w2(-1:1,-1:1)  !< w2 on grid center
    4722 !    REAL(wp)     ::  uv(-1:1,-1:1)  !< u*v on grid center
    4723 !    REAL(wp)     ::  uw(-1:1,-1:1)  !< u*w on grid center
    4724 !    REAL(wp)     ::  vw(-1:1,-1:1)  !< v*w on grid center
    47254725
    47264726    REAL(wp)     ::  ut(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered u
     
    47404740    REAL(wp)     ::  sd11        !< deviatoric shear tensor
    47414741    REAL(wp)     ::  sd22        !< deviatoric shear tensor
    4742     REAL(wp)     ::  sd33        !< deviatoric shear tensor
     4742    REAL(wp)     ::  sd33        !<f deviatoric shear tensor
    47434743    REAL(wp)     ::  sd12        !< deviatoric shear tensor
    47444744    REAL(wp)     ::  sd13        !< deviatoric shear tensor
     
    47824782       DO  j = nys, nyn
    47834783          DO  k = nzb+1, nzt
     4784           
     4785             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4786             
    47844787!
    47854788!--          Compute the deviatoric shear tensor s_ij on grid centers:
     
    48854888!
    48864889!--          c* nu*^T is SGS viscosity on test filter level:
    4887              cstnust_t = -ldsd / sdt2
     4890             cstnust_t = -ldsd / ( sdt2 + 1.0E-20_wp )
    48884891!
    48894892!--          The model was only tested for an isotropic grid. The following
    48904893!--          expression was a recommendation of Stefan Heinz.
    48914894             delta = MAX( dx, dy, dzw(k) )
    4892              cst = cstnust_t / ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) )
     4895
     4896             IF ( lnn <= 0.0_wp ) THEN
     4897                cst = 0.0_wp
     4898             ELSE
     4899                cst = cstnust_t /                                              &
     4900                   ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp )
     4901             ENDIF
     4902
    48934903!
    48944904!--          Calculate border according to Mokhtarpoor and Heinz (2017)
    4895              cst_max = fac_cmax * SQRT( e(k,j,i) ) / ( delta * SQRT( 2.0_wp * sd2 ) )
    4896 
     4905             cst_max = fac_cmax * SQRT( e(k,j,i) ) /                           &
     4906                       ( delta * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )
     4907             
    48974908             IF ( ABS( cst ) > cst_max )  THEN
    48984909                cst = cst_max * cst / ABS( cst )
    48994910             ENDIF
    49004911
    4901              km(k,j,i) = cst * delta * SQRT( e(k,j,i) )
    4902 
     4912             km(k,j,i) = cst * delta * SQRT( e(k,j,i) ) * flag
     4913             
    49034914          ENDDO
    49044915       ENDDO
Note: See TracChangeset for help on using the changeset viewer.