Changeset 3430 for palm/trunk/SOURCE/turbulence_closure_mod.f90

Timestamp:

Oct 25, 2018 1:36:23 PM (6 years ago)

Author:

maronga

Message:

added support for buildings with dynamic sgs model

File:

• palm/trunk/SOURCE/turbulence_closure_mod.f90 (modified) (7 diffs)

palm/trunk/SOURCE/turbulence_closure_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Added support for buildings in the dynamic SGS model
+! 
+! 3398 2018-10-22 19:30:24Z knoop
 ! Refactored production_e and production_e_ij (removed code redundancy)
 ! 
@@ -186 +189 @@
                simulated_time,timestep_scheme, turbulence_closure,             &
                turbulent_inflow, use_upstream_for_tke, vpt_reference,          &
-               ws_scheme_sca
+               ws_scheme_sca, current_timestep_number
 
     USE advec_ws,                                                              &
@@ -4687 +4690 @@
 !>
 !> @author Hauke Wurps
+!> @author Björn Maronga
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_diffusivities_dynamic
@@ -4714 +4718 @@
     REAL(wp)     ::  dwdz        !< Gradient of w-component in z-direction
 
+    REAL(wp)     ::  flag                !< topography flag
+    
     REAL(wp)     ::  uc(-1:1,-1:1)  !< u on grid center
     REAL(wp)     ::  vc(-1:1,-1:1)  !< v on grid center
     REAL(wp)     ::  wc(-1:1,-1:1)  !< w on grid center
-!    REAL(wp)     ::  u2(-1:1,-1:1)  !< u2 on grid center
-!    REAL(wp)     ::  v2(-1:1,-1:1)  !< v2 on grid center
-!    REAL(wp)     ::  w2(-1:1,-1:1)  !< w2 on grid center
-!    REAL(wp)     ::  uv(-1:1,-1:1)  !< u*v on grid center
-!    REAL(wp)     ::  uw(-1:1,-1:1)  !< u*w on grid center
-!    REAL(wp)     ::  vw(-1:1,-1:1)  !< v*w on grid center
 
     REAL(wp)     ::  ut(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered u
@@ -4740 +4740 @@
     REAL(wp)     ::  sd11        !< deviatoric shear tensor
     REAL(wp)     ::  sd22        !< deviatoric shear tensor
-    REAL(wp)     ::  sd33        !< deviatoric shear tensor
+    REAL(wp)     ::  sd33        !<f deviatoric shear tensor
     REAL(wp)     ::  sd12        !< deviatoric shear tensor
     REAL(wp)     ::  sd13        !< deviatoric shear tensor
@@ -4782 +4782 @@
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
+            
+             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+             
 !
 !--          Compute the deviatoric shear tensor s_ij on grid centers:
@@ -4885 +4888 @@
 !
 !--          c* nu*^T is SGS viscosity on test filter level:
-             cstnust_t = -ldsd / sdt2
+             cstnust_t = -ldsd / ( sdt2 + 1.0E-20_wp )
 !
 !--          The model was only tested for an isotropic grid. The following
 !--          expression was a recommendation of Stefan Heinz.
              delta = MAX( dx, dy, dzw(k) )
-             cst = cstnust_t / ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) )
+
+             IF ( lnn <= 0.0_wp ) THEN
+                cst = 0.0_wp
+             ELSE
+                cst = cstnust_t /                                              &
+                   ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp )
+             ENDIF
+
 !
 !--          Calculate border according to Mokhtarpoor and Heinz (2017)
-             cst_max = fac_cmax * SQRT( e(k,j,i) ) / ( delta * SQRT( 2.0_wp * sd2 ) )
-
+             cst_max = fac_cmax * SQRT( e(k,j,i) ) /                           &
+                       ( delta * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )
+             
              IF ( ABS( cst ) > cst_max )  THEN
                 cst = cst_max * cst / ABS( cst )
              ENDIF
 
-             km(k,j,i) = cst * delta * SQRT( e(k,j,i) )
-
+             km(k,j,i) = cst * delta * SQRT( e(k,j,i) ) * flag
+             
           ENDDO
        ENDDO