Changeset 3120 for palm/trunk/SOURCE
- Timestamp:
- Jul 11, 2018 6:30:57 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
modules.f90 (modified) (2 diffs)
-
timestep.f90 (modified) (3 diffs)
-
turbulence_closure_mod.f90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/modules.f90
r3083 r3120 25 25 ! ----------------- 26 26 ! $Id$ 27 ! +les_dynamic 28 ! 29 ! 3083 2018-06-19 14:03:12Z gronemeier 27 30 ! set dt_3d = 0.01 28 31 ! … … 1306 1309 LOGICAL :: large_scale_subsidence = .FALSE. !< namelist parameter 1307 1310 LOGICAL :: land_surface = .FALSE. !< use land surface model? 1311 LOGICAL :: les_dynamic = .FALSE. !< use dynamic subgrid model as turbulence closure for LES mode 1308 1312 LOGICAL :: les_mw = .FALSE. !< use Moeng-Wyngaard turbulence closure for LES mode 1309 1313 LOGICAL :: lsf_exception = .FALSE. !< use of lsf with buildings (temporary)? -
palm/trunk/SOURCE/timestep.f90
r3084 r3120 25 25 ! ----------------- 26 26 ! $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 27 30 ! limit increase of dt_3d only in case of RANS mode 28 31 ! … … 302 305 ! 303 306 !-- 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. 305 309 dt_diff_l = 999999.0_wp 306 310 … … 315 319 DO k = nzb+1, nzt 316 320 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 ) ) 318 323 ENDDO 319 324 ENDDO -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r3086 r3120 25 25 ! ----------------- 26 26 ! $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 27 32 ! bugfix: set rans_const_sigma(1) = 1.3 28 33 ! … … 95 100 ! -------- 96 101 ! @author Tobias Gronemeier 97 ! 102 ! @author Hauke Wurps 98 103 ! 99 104 ! Description: … … 126 131 ONLY: constant_diffusion, dt_3d, e_init, humidity, inflow_l, & 127 132 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 133 139 134 140 USE advec_ws, & … … 148 154 149 155 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, & 152 157 wall_flags_0 153 158 … … 296 301 297 302 ! 298 !-- Cal culate diffusivities303 !-- Call tcm_diffusivities_default and tcm_diffusivities_dynamic 299 304 INTERFACE tcm_diffusivities 300 305 MODULE PROCEDURE tcm_diffusivities 301 306 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 302 326 303 327 ! … … 390 414 CASE ( 'Moeng_Wyngaard' ) 391 415 les_mw = .TRUE. 416 417 CASE ( 'dynamic' ) 418 les_dynamic = .TRUE. 392 419 393 420 CASE DEFAULT … … 1945 1972 IMPLICIT NONE 1946 1973 1947 INTEGER(iwp) :: i !< grid index x-direction1948 INTEGER(iwp) :: j !< grid index y-direction1949 INTEGER(iwp) :: k !< grid index z-direction1950 INTEGER(iwp) :: m !< running index surface elements1974 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 1951 1978 1952 1979 IF ( constant_flux_layer ) THEN … … 1973 2000 !-- and km are not on the same grid. Actually, a further 1974 2001 !-- 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. 1976 2003 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 ) 1980 2007 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 ) 1984 2011 1985 2012 IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > & … … 2032 2059 !-- interpolation of km onto the u/v-grid is necessary. However, the 2033 2060 !-- 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 ) 2042 2070 2043 2071 IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > & … … 2063 2091 !-- interpolation of km onto the u/v-grid is necessary. However, the 2064 2092 !-- 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 ) 2073 2101 2074 2102 IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > & … … 4230 4258 ! 4231 4259 !-- Calculate dissipation 4232 IF ( les_ mw ) THEN4260 IF ( les_dynamic .OR. les_mw ) THEN 4233 4261 4234 4262 CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) … … 4367 4395 !-- Calculate dissipation... 4368 4396 !-- ...in case of LES 4369 IF ( les_ mw ) THEN4397 IF ( les_dynamic .OR. les_mw ) THEN 4370 4398 4371 4399 CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) … … 4678 4706 ! Description: 4679 4707 ! ------------ 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 ! ------------ 4680 4734 !> Computation of the turbulent diffusion coefficients for momentum and heat 4681 4735 !> according to Prandtl-Kolmogorov. 4682 4736 !> @todo consider non-default surfaces 4683 4737 !------------------------------------------------------------------------------! 4684 SUBROUTINE tcm_diffusivities ( var, var_reference )4738 SUBROUTINE tcm_diffusivities_default( var, var_reference ) 4685 4739 4686 4740 … … 4746 4800 ENDIF 4747 4801 4748 IF ( les_ mw ) THEN4802 IF ( les_dynamic .OR. les_mw ) THEN 4749 4803 !$OMP DO 4750 4804 DO i = nxlg, nxrg … … 4926 4980 ENDIF 4927 4981 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 4929 5347 4930 5348
Note: See TracChangeset
for help on using the changeset viewer.
