- Timestamp:
- Jul 16, 2018 11:08:55 AM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
Makefile (modified) (2 diffs)
-
surface_layer_fluxes_mod.f90 (modified) (3 diffs)
-
turbulence_closure_mod.f90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r3129 r3130 25 25 # ----------------- 26 26 # $Id$ 27 # add surface_layer_fluxes_mod to turbulence_closure_mod 28 # 29 # 3129 2018-07-16 07:45:13Z gronemeier 27 30 # add turbulence_closure_mod to parin 28 31 # … … 1531 1534 plant_canopy_model_mod.o \ 1532 1535 pmc_interface_mod.o \ 1536 surface_layer_fluxes_mod.o \ 1533 1537 surface_mod.o \ 1534 1538 user_actions.o -
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r3045 r3130 26 26 ! ----------------- 27 27 ! $Id$ 28 ! move phi_m from turbulence_closure_mod 29 ! 30 ! 3045 2018-05-28 07:55:41Z Giersch 28 31 ! Error message revised 29 32 ! … … 287 290 PRIVATE 288 291 289 PUBLIC init_surface_layer_fluxes, surface_layer_fluxes292 PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes 290 293 291 294 INTERFACE init_surface_layer_fluxes 292 295 MODULE PROCEDURE init_surface_layer_fluxes 293 296 END INTERFACE init_surface_layer_fluxes 297 298 INTERFACE phi_m 299 MODULE PROCEDURE phi_m 300 END INTERFACE phi_m 294 301 295 302 INTERFACE surface_layer_fluxes … … 2293 2300 END FUNCTION psi_h 2294 2301 2302 2303 !------------------------------------------------------------------------------! 2304 ! Description: 2305 ! ------------ 2306 !> Calculates stability function for momentum 2307 !> 2308 !> @author Hauke Wurps 2309 !------------------------------------------------------------------------------! 2310 FUNCTION phi_m( zeta ) 2311 2312 IMPLICIT NONE 2313 2314 REAL(wp) :: phi_m !< Value of the function 2315 REAL(wp) :: zeta !< Stability parameter z/L 2316 2317 REAL(wp), PARAMETER :: a = 16.0_wp !< constant 2318 REAL(wp), PARAMETER :: c = 5.0_wp !< constant 2319 2320 IF ( zeta < 0.0_wp ) THEN 2321 phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) ) 2322 ELSE 2323 phi_m = 1.0_wp + c * zeta 2324 ENDIF 2325 2326 END FUNCTION phi_m 2327 2295 2328 END MODULE surface_layer_fluxes_mod -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r3129 r3130 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - move boundary condition of km and kh to tcm_diffusivities 28 ! - calculate km at boundaries according to MOST 29 ! - move phi_m to surface_layer_fluxes_mod 30 ! 31 ! 3129 2018-07-16 07:45:13Z gronemeier 27 32 ! - move limitation of diss to boundary_conds 28 33 ! - move boundary conditions for e and diss to boundary_conds … … 1910 1915 INTEGER(iwp) :: m !< running index surface elements 1911 1916 1912 REAL(wp) :: km_sfc !< km according to MOST1913 1914 1917 IF ( constant_flux_layer ) THEN 1915 1918 ! … … 1922 1925 !-- surface winds). 1923 1926 !-- not available in case of non-cyclic boundary conditions. 1924 !-- WARNING: the exact analytical solution would require the determination1925 !-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.1926 1927 !-- Default surfaces, upward-facing 1927 1928 !$OMP PARALLEL DO PRIVATE(i,j,k,m) … … 1932 1933 k = surf_def_h(0)%k(m) 1933 1934 ! 1934 !-- Note, calculation eof u_0 and v_0 is not fully accurate, as u/v1935 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1935 1936 !-- and km are not on the same grid. Actually, a further 1936 1937 !-- interpolation of km onto the u/v-grid is necessary. However, the 1937 1938 !-- effect of this error is negligible. 1938 km_sfc = surf_def_h(0)%us(m) * kappa * surf_def_h(0)%z_mo(m) / &1939 phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )1940 1941 1939 surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * & 1942 1940 drho_air_zw(k-1) * & 1943 1941 ( zu(k+1) - zu(k-1) ) / & 1944 ( km _sfc+ 1.0E-20_wp )1942 ( km(k,j,i) + 1.0E-20_wp ) 1945 1943 surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * & 1946 1944 drho_air_zw(k-1) * & 1947 1945 ( zu(k+1) - zu(k-1) ) / & 1948 ( km _sfc+ 1.0E-20_wp )1946 ( km(k,j,i) + 1.0E-20_wp ) 1949 1947 1950 1948 IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > & … … 1965 1963 j = surf_def_h(1)%j(m) 1966 1964 k = surf_def_h(1)%k(m) 1967 1965 ! 1966 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1967 !-- and km are not on the same grid. Actually, a further 1968 !-- interpolation of km onto the u/v-grid is necessary. However, the 1969 !-- effect of this error is negligible. 1968 1970 surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * & 1969 1971 drho_air_zw(k-1) * & … … 1989 1991 DO m = 1, surf_lsm_h%ns 1990 1992 1991 i = surf_lsm_h%i(m) 1993 i = surf_lsm_h%i(m) 1992 1994 j = surf_lsm_h%j(m) 1993 1995 k = surf_lsm_h%k(m) 1994 1996 ! 1995 !-- Note, calculation eof u_0 and v_0 is not fully accurate, as u/v1997 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1996 1998 !-- and km are not on the same grid. Actually, a further 1997 1999 !-- interpolation of km onto the u/v-grid is necessary. However, the 1998 !-- effect of this error is negligible. 1999 !-- effect of this error is negligible. 2000 km_sfc = surf_lsm_h%us(m) * kappa * surf_lsm_h%z_mo(m) / & 2001 phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) ) 2002 2000 !-- effect of this error is negligible. 2003 2001 surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * & 2004 2002 drho_air_zw(k-1) * & 2005 ( zu(k+1) - zu(k-1) ) /&2006 ( km _sfc+ 1.0E-20_wp )2003 ( zu(k+1) - zu(k-1) ) / & 2004 ( km(k,j,i) + 1.0E-20_wp ) 2007 2005 surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & 2008 2006 drho_air_zw(k-1) * & 2009 2007 ( zu(k+1) - zu(k-1) ) / & 2010 ( km _sfc+ 1.0E-20_wp )2008 ( km(k,j,i) + 1.0E-20_wp ) 2011 2009 2012 2010 IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > & … … 2024 2022 DO m = 1, surf_usm_h%ns 2025 2023 2026 i = surf_usm_h%i(m) 2024 i = surf_usm_h%i(m) 2027 2025 j = surf_usm_h%j(m) 2028 2026 k = surf_usm_h%k(m) 2029 2027 ! 2030 !-- Note, calculation eof u_0 and v_0 is not fully accurate, as u/v2028 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 2031 2029 !-- and km are not on the same grid. Actually, a further 2032 2030 !-- interpolation of km onto the u/v-grid is necessary. However, the 2033 !-- effect of this error is negligible. 2034 km_sfc = surf_usm_h%us(m) * kappa * surf_usm_h%z_mo(m) / & 2035 phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) ) 2036 2031 !-- effect of this error is negligible. 2037 2032 surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * & 2038 2033 drho_air_zw(k-1) * & 2039 2034 ( zu(k+1) - zu(k-1) ) / & 2040 ( km _sfc+ 1.0E-20_wp )2035 ( km(k,j,i) + 1.0E-20_wp ) 2041 2036 surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * & 2042 2037 drho_air_zw(k-1) * & 2043 2038 ( zu(k+1) - zu(k-1) ) / & 2044 ( km _sfc+ 1.0E-20_wp )2039 ( km(k,j,i) + 1.0E-20_wp ) 2045 2040 2046 2041 IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > & … … 2057 2052 2058 2053 END SUBROUTINE production_e_init 2059 2060 2061 !------------------------------------------------------------------------------!2062 ! Description:2063 ! ------------2064 !> Calculates stability function for momentum2065 !>2066 !> @author Hauke Wurps2067 !------------------------------------------------------------------------------!2068 FUNCTION phi_m( zeta )2069 2070 IMPLICIT NONE2071 2072 REAL(wp) :: phi_m !< Value of the function2073 REAL(wp) :: zeta !< Stability parameter z/L2074 2075 REAL(wp), PARAMETER :: a = 16.0_wp !< constant2076 REAL(wp), PARAMETER :: b = -0.25_wp !< constant2077 REAL(wp), PARAMETER :: c = 5.0_wp !< constant2078 2079 IF ( zeta < 0.0_wp ) THEN2080 phi_m = ( 1.0_wp - a * zeta )**b2081 ELSE2082 phi_m = 1.0_wp + c * zeta2083 ENDIF2084 2085 END FUNCTION phi_m2086 2054 2087 2055 … … 4626 4594 SUBROUTINE tcm_diffusivities( var, var_reference ) 4627 4595 4596 USE control_parameters, & 4597 ONLY: e_min, outflow_l, outflow_n, outflow_r, outflow_s 4598 4599 USE surface_layer_fluxes_mod, & 4600 ONLY: phi_m 4601 4602 USE surface_mod, & 4603 ONLY : bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, & 4604 surf_usm_h, surf_usm_v 4605 4606 INTEGER(iwp) :: i !< loop index 4607 INTEGER(iwp) :: j !< loop index 4608 INTEGER(iwp) :: k !< loop index 4609 INTEGER(iwp) :: m !< loop index 4610 INTEGER(iwp) :: n !< loop index 4611 4628 4612 REAL(wp) :: var_reference !< reference temperature 4629 4613 … … 4633 4617 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature 4634 4618 #endif 4635 !4636 !-- Call default diffusivities routine. This is always used to calculate kh.4637 CALL tcm_diffusivities_default( var, var_reference )4638 !4639 !-- Call dynamic subgrid model to calculate km.4640 IF ( les_dynamic ) THEN4641 CALL tcm_diffusivities_dynamic4642 ENDIF4643 4644 END SUBROUTINE tcm_diffusivities4645 4646 4647 !------------------------------------------------------------------------------!4648 ! Description:4649 ! ------------4650 !> Computation of the turbulent diffusion coefficients for momentum and heat4651 !> according to Prandtl-Kolmogorov.4652 !------------------------------------------------------------------------------!4653 SUBROUTINE tcm_diffusivities_default( var, var_reference )4654 4655 4656 USE control_parameters, &4657 ONLY: e_min, outflow_l, outflow_n, outflow_r, outflow_s4658 4659 USE grid_variables, &4660 ONLY: dx, dy4661 4662 USE statistics, &4663 ONLY : rmask, sums_l_l4664 4665 USE surface_mod, &4666 ONLY : bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, &4667 surf_usm_h, surf_usm_v4668 4669 IMPLICIT NONE4670 4671 INTEGER(iwp) :: i !< loop index4672 INTEGER(iwp) :: j !< loop index4673 INTEGER(iwp) :: k !< loop index4674 INTEGER(iwp) :: m !< loop index4675 INTEGER(iwp) :: n !< loop index4676 INTEGER(iwp) :: omp_get_thread_num !< opemmp function to get thread number4677 INTEGER(iwp) :: sr !< statistic region4678 INTEGER(iwp) :: tn !< thread number4679 4680 REAL(wp) :: flag !< topography flag4681 REAL(wp) :: l !< mixing length4682 REAL(wp) :: ll !< adjusted mixing length4683 REAL(wp) :: var_reference !< reference temperature4684 4685 #if defined( __nopointer )4686 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature4687 #else4688 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature4689 #endif4690 4691 !4692 !-- Default thread number in case of one thread4693 tn = 04694 4695 !4696 !-- Initialization for calculation of the mixing length profile4697 sums_l_l = 0.0_wp4698 4699 !4700 !-- Compute the turbulent diffusion coefficient for momentum4701 !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)4702 !$ tn = omp_get_thread_num()4703 4619 4704 4620 ! 4705 4621 !-- Introduce an optional minimum tke 4706 4622 IF ( e_min > 0.0_wp ) THEN 4707 !$OMP DO4623 !$OMP PARALLEL DO PRIVATE(i,j,k) 4708 4624 DO i = nxlg, nxrg 4709 4625 DO j = nysg, nyng … … 4716 4632 ENDIF 4717 4633 4718 IF ( les_dynamic .OR. les_mw ) THEN 4719 !$OMP DO 4720 DO i = nxlg, nxrg 4721 DO j = nysg, nyng 4722 DO k = nzb+1, nzt 4723 4724 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 4725 4726 ! 4727 !-- Determine the mixing length for LES closure 4728 CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) 4729 ! 4730 !-- Compute diffusion coefficients for momentum and heat 4731 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag 4732 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag 4733 ! 4734 !-- Summation for averaged profile (cf. flow_statistics) 4735 DO sr = 0, statistic_regions 4736 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) & 4737 * flag 4738 ENDDO 4739 4740 ENDDO 4741 ENDDO 4742 ENDDO 4743 4744 ELSEIF ( rans_tke_l ) THEN 4745 4746 !$OMP DO 4747 DO i = nxlg, nxrg 4748 DO j = nysg, nyng 4749 DO k = nzb+1, nzt 4750 4751 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 4752 ! 4753 !-- Mixing length for RANS mode with TKE-l closure 4754 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference ) 4755 ! 4756 !-- Compute diffusion coefficients for momentum and heat 4757 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag 4758 kh(k,j,i) = km(k,j,i) / prandtl_number * flag 4759 ! 4760 !-- Summation for averaged profile (cf. flow_statistics) 4761 DO sr = 0, statistic_regions 4762 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) & 4763 * flag 4764 ENDDO 4765 4766 ENDDO 4767 ENDDO 4768 ENDDO 4769 4770 ELSEIF ( rans_tke_e ) THEN 4771 4772 !$OMP DO 4773 DO i = nxlg, nxrg 4774 DO j = nysg, nyng 4775 DO k = nzb+1, nzt 4776 4777 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 4778 ! 4779 !-- Compute diffusion coefficients for momentum and heat 4780 km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag 4781 kh(k,j,i) = km(k,j,i) / prandtl_number * flag 4782 ! 4783 !-- Summation for averaged profile of mixing length (cf. flow_statistics) 4784 DO sr = 0, statistic_regions 4785 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + & 4786 c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) / & 4787 ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag 4788 ENDDO 4789 4790 ENDDO 4791 ENDDO 4792 ENDDO 4793 4634 ! 4635 !-- Call default diffusivities routine. This is always used to calculate kh. 4636 CALL tcm_diffusivities_default( var, var_reference ) 4637 ! 4638 !-- Call dynamic subgrid model to calculate km. 4639 IF ( les_dynamic ) THEN 4640 CALL tcm_diffusivities_dynamic 4794 4641 ENDIF 4795 4642 4796 sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn) ! quasi boundary-condition for 4797 ! data output 4798 !$OMP END PARALLEL 4799 4800 ! 4801 !-- Set vertical boundary values (Neumann conditions both at upward- and 4802 !-- downward facing walls. To set wall-boundary values, the surface data type 4803 !-- is applied. 4804 !-- Horizontal boundary conditions at vertical walls are not set because 4805 !-- so far vertical surfaces require usage of a Prandtl-layer where the boundary 4806 !-- values of the diffusivities are not needed. 4807 IF ( .NOT. rans_tke_e ) THEN 4808 ! 4809 !-- Upward-facing 4810 !$OMP PARALLEL DO PRIVATE( i, j, k ) 4811 DO m = 1, bc_h(0)%ns 4812 i = bc_h(0)%i(m) 4813 j = bc_h(0)%j(m) 4814 k = bc_h(0)%k(m) 4815 km(k-1,j,i) = km(k,j,i) 4816 kh(k-1,j,i) = kh(k,j,i) 4817 ENDDO 4818 ! 4819 !-- Downward facing surfaces 4820 !$OMP PARALLEL DO PRIVATE( i, j, k ) 4821 DO m = 1, bc_h(1)%ns 4822 i = bc_h(1)%i(m) 4823 j = bc_h(1)%j(m) 4824 k = bc_h(1)%k(m) 4825 km(k+1,j,i) = km(k,j,i) 4826 kh(k+1,j,i) = kh(k,j,i) 4827 ENDDO 4828 ELSE 4829 ! 4830 !-- Up- and downward facing surfaces 4643 ! 4644 !-- Set boundary values (Neumann conditions at downward facing surfaces, 4645 !-- according to MOST at upward facing surfaces and vertical surfaces). 4646 ! 4647 !-- Downward facing surfaces 4648 !$OMP PARALLEL DO PRIVATE(i,j,k) 4649 DO m = 1, bc_h(1)%ns 4650 i = bc_h(1)%i(m) 4651 j = bc_h(1)%j(m) 4652 k = bc_h(1)%k(m) 4653 km(k+1,j,i) = km(k,j,i) 4654 kh(k+1,j,i) = kh(k,j,i) 4655 ENDDO 4656 4657 !-- Upward facing surfaces 4658 !-- Default surfaces 4659 n = 0 4660 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4661 DO m = 1, surf_def_h(n)%ns 4662 i = surf_def_h(n)%i(m) 4663 j = surf_def_h(n)%j(m) 4664 k = surf_def_h(n)%k(m) 4665 km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m) / & 4666 phi_m( surf_def_h(n)%z_mo(m) / surf_def_h(n)%ol(m) ) 4667 kh(k,j,i) = 1.35_wp * km(k,j,i) 4668 ENDDO 4669 ! 4670 !-- Natural surfaces 4671 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4672 DO m = 1, surf_lsm_h%ns 4673 i = surf_lsm_h%i(m) 4674 j = surf_lsm_h%j(m) 4675 k = surf_lsm_h%k(m) 4676 km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / & 4677 phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) ) 4678 kh(k,j,i) = 1.35_wp * km(k,j,i) 4679 ENDDO 4680 ! 4681 !-- Urban surfaces 4682 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4683 DO m = 1, surf_usm_h%ns 4684 i = surf_usm_h%i(m) 4685 j = surf_usm_h%j(m) 4686 k = surf_usm_h%k(m) 4687 km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / & 4688 phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) ) 4689 kh(k,j,i) = 1.35_wp * km(k,j,i) 4690 ENDDO 4691 4692 ! 4693 !-- North-, south-, west and eastward facing surfaces 4694 DO n = 0, 3 4695 ! 4831 4696 !-- Default surfaces 4832 DO n = 0, 1 4833 DO m = 1, surf_def_h(n)%ns 4834 i = surf_def_h(n)%i(m) 4835 j = surf_def_h(n)%j(m) 4836 k = surf_def_h(n)%k(m) 4837 km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m) 4838 kh(k,j,i) = 1.35_wp * km(k,j,i) 4839 ENDDO 4840 ENDDO 4841 ! 4842 !-- Upward facing surfaces 4843 !-- Natural surfaces 4844 DO m = 1, surf_lsm_h%ns 4845 i = surf_lsm_h%i(m) 4846 j = surf_lsm_h%j(m) 4847 k = surf_lsm_h%k(m) 4848 km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) 4697 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4698 DO m = 1, surf_def_v(n)%ns 4699 i = surf_def_v(n)%i(m) 4700 j = surf_def_v(n)%j(m) 4701 k = surf_def_v(n)%k(m) 4702 km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m) / & 4703 phi_m( surf_def_v(n)%z_mo(m) / surf_def_v(n)%ol(m) ) 4849 4704 kh(k,j,i) = 1.35_wp * km(k,j,i) 4850 4705 ENDDO 4851 4706 ! 4852 !-- Urban surfaces 4853 DO m = 1, surf_usm_h%ns 4854 i = surf_usm_h%i(m) 4855 j = surf_usm_h%j(m) 4856 k = surf_usm_h%k(m) 4857 km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) 4707 !-- Natural surfaces 4708 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4709 DO m = 1, surf_lsm_v(n)%ns 4710 i = surf_lsm_v(n)%i(m) 4711 j = surf_lsm_v(n)%j(m) 4712 k = surf_lsm_v(n)%k(m) 4713 km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m) / & 4714 phi_m( surf_lsm_v(n)%z_mo(m) / surf_lsm_v(n)%ol(m) ) 4858 4715 kh(k,j,i) = 1.35_wp * km(k,j,i) 4859 4716 ENDDO 4860 4861 ! 4862 !-- North-, south-, west and eastward facing surfaces 4863 DO n = 0, 3 4864 ! 4865 !-- Default surfaces 4866 DO m = 1, surf_def_v(n)%ns 4867 i = surf_def_v(n)%i(m) 4868 j = surf_def_v(n)%j(m) 4869 k = surf_def_v(n)%k(m) 4870 km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m) 4871 kh(k,j,i) = 1.35_wp * km(k,j,i) 4872 ENDDO 4873 ! 4874 !-- Natural surfaces 4875 DO m = 1, surf_lsm_v(n)%ns 4876 i = surf_lsm_v(n)%i(m) 4877 j = surf_lsm_v(n)%j(m) 4878 k = surf_lsm_v(n)%k(m) 4879 km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m) 4880 kh(k,j,i) = 1.35_wp * km(k,j,i) 4881 ENDDO 4882 ! 4883 !-- Urban surfaces 4884 DO m = 1, surf_usm_v(n)%ns 4885 i = surf_usm_v(n)%i(m) 4886 j = surf_usm_v(n)%j(m) 4887 k = surf_usm_v(n)%k(m) 4888 km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m) 4889 kh(k,j,i) = 1.35_wp * km(k,j,i) 4890 ENDDO 4717 ! 4718 !-- Urban surfaces 4719 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 4720 DO m = 1, surf_usm_v(n)%ns 4721 i = surf_usm_v(n)%i(m) 4722 j = surf_usm_v(n)%j(m) 4723 k = surf_usm_v(n)%k(m) 4724 km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m) / & 4725 phi_m( surf_usm_v(n)%z_mo(m) / surf_usm_v(n)%ol(m) ) 4726 kh(k,j,i) = 1.35_wp * km(k,j,i) 4891 4727 ENDDO 4892 4893 CALL exchange_horiz( km, nbgp ) 4894 CALL exchange_horiz( kh, nbgp ) 4895 4896 ENDIF 4728 ENDDO 4729 4730 CALL exchange_horiz( km, nbgp ) 4731 CALL exchange_horiz( kh, nbgp ) 4897 4732 4898 4733 ! … … 4926 4761 ENDIF 4927 4762 4763 END SUBROUTINE tcm_diffusivities 4764 4765 4766 !------------------------------------------------------------------------------! 4767 ! Description: 4768 ! ------------ 4769 !> Computation of the turbulent diffusion coefficients for momentum and heat 4770 !> according to Prandtl-Kolmogorov. 4771 !------------------------------------------------------------------------------! 4772 SUBROUTINE tcm_diffusivities_default( var, var_reference ) 4773 4774 USE statistics, & 4775 ONLY : rmask, sums_l_l 4776 4777 IMPLICIT NONE 4778 4779 INTEGER(iwp) :: i !< loop index 4780 INTEGER(iwp) :: j !< loop index 4781 INTEGER(iwp) :: k !< loop index 4782 INTEGER(iwp) :: omp_get_thread_num !< opemmp function to get thread number 4783 INTEGER(iwp) :: sr !< statistic region 4784 INTEGER(iwp) :: tn !< thread number 4785 4786 REAL(wp) :: flag !< topography flag 4787 REAL(wp) :: l !< mixing length 4788 REAL(wp) :: ll !< adjusted mixing length 4789 REAL(wp) :: var_reference !< reference temperature 4790 4791 #if defined( __nopointer ) 4792 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature 4793 #else 4794 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature 4795 #endif 4796 4797 ! 4798 !-- Default thread number in case of one thread 4799 tn = 0 4800 4801 ! 4802 !-- Initialization for calculation of the mixing length profile 4803 sums_l_l = 0.0_wp 4804 4805 ! 4806 !-- Compute the turbulent diffusion coefficient for momentum 4807 !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag) 4808 !$ tn = omp_get_thread_num() 4809 4810 IF ( les_dynamic .OR. les_mw ) THEN 4811 !$OMP DO 4812 DO i = nxlg, nxrg 4813 DO j = nysg, nyng 4814 DO k = nzb+1, nzt 4815 4816 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 4817 4818 ! 4819 !-- Determine the mixing length for LES closure 4820 CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) 4821 ! 4822 !-- Compute diffusion coefficients for momentum and heat 4823 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag 4824 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag 4825 ! 4826 !-- Summation for averaged profile (cf. flow_statistics) 4827 DO sr = 0, statistic_regions 4828 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) & 4829 * flag 4830 ENDDO 4831 4832 ENDDO 4833 ENDDO 4834 ENDDO 4835 4836 ELSEIF ( rans_tke_l ) THEN 4837 4838 !$OMP DO 4839 DO i = nxlg, nxrg 4840 DO j = nysg, nyng 4841 DO k = nzb+1, nzt 4842 4843 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 4844 ! 4845 !-- Mixing length for RANS mode with TKE-l closure 4846 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference ) 4847 ! 4848 !-- Compute diffusion coefficients for momentum and heat 4849 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag 4850 kh(k,j,i) = km(k,j,i) / prandtl_number * flag 4851 ! 4852 !-- Summation for averaged profile (cf. flow_statistics) 4853 DO sr = 0, statistic_regions 4854 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) & 4855 * flag 4856 ENDDO 4857 4858 ENDDO 4859 ENDDO 4860 ENDDO 4861 4862 ELSEIF ( rans_tke_e ) THEN 4863 4864 !$OMP DO 4865 DO i = nxlg, nxrg 4866 DO j = nysg, nyng 4867 DO k = nzb+1, nzt 4868 4869 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 4870 ! 4871 !-- Compute diffusion coefficients for momentum and heat 4872 km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag 4873 kh(k,j,i) = km(k,j,i) / prandtl_number * flag 4874 ! 4875 !-- Summation for averaged profile of mixing length (cf. flow_statistics) 4876 DO sr = 0, statistic_regions 4877 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + & 4878 c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) / & 4879 ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag 4880 ENDDO 4881 4882 ENDDO 4883 ENDDO 4884 ENDDO 4885 4886 ENDIF 4887 4888 sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn) ! quasi boundary-condition for 4889 ! data output 4890 !$OMP END PARALLEL 4891 4928 4892 END SUBROUTINE tcm_diffusivities_default 4929 4893 … … 4949 4913 USE arrays_3d, & 4950 4914 ONLY: ddzw, dzw, dd2zu, w, ug, vg 4951 4952 USE control_parameters, &4953 ONLY: e_min, outflow_l, outflow_n, outflow_r, outflow_s4954 4915 4955 4916 USE grid_variables, & 4956 4917 ONLY : ddx, ddy, dx, dy 4957 4958 USE surface_mod, &4959 ONLY : bc_h4960 4918 4961 4919 IMPLICIT NONE … … 4966 4924 INTEGER(iwp) :: l !< running index 4967 4925 INTEGER(iwp) :: m !< running index 4968 INTEGER(iwp) :: n !< running index4969 4926 4970 4927 REAL(wp) :: dudx !< Gradient of u-component in x-direction … … 5037 4994 REAL(wp), PARAMETER :: fac_cmax = 23.0_wp/(24.0_wp*sqrt(3.0_wp)) !< constant 5038 4995 5039 !5040 !-- Introduce an optional minimum tke5041 IF ( e_min > 0.0_wp ) THEN5042 !$OMP DO5043 DO i = nxlg, nxrg5044 DO j = nysg, nyng5045 DO k = nzb+1, nzt5046 e(k,j,i) = MAX( e(k,j,i), e_min ) * &5047 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )5048 ENDDO5049 ENDDO5050 ENDDO5051 ENDIF5052 4996 ! 5053 4997 !-- velocities on grid centers: … … 5181 5125 ENDDO 5182 5126 ENDDO 5183 !5184 !-- Set vertical boundary values (Neumann conditions both at upward- and5185 !-- downward facing walls. To set wall-boundary values, the surface data type5186 !-- is applied.5187 !-- Horizontal boundary conditions at vertical walls are not set because5188 !-- so far vertical surfaces require usage of a constant-flux layer where the5189 !-- boundary values of the diffusivities are not needed.5190 5191 !5192 !-- Upward-facing5193 !$OMP PARALLEL DO PRIVATE( i, j, k )5194 DO m = 1, bc_h(0)%ns5195 i = bc_h(0)%i(m)5196 j = bc_h(0)%j(m)5197 k = bc_h(0)%k(m)5198 km(k-1,j,i) = km(k,j,i)5199 ENDDO5200 !5201 !-- Downward facing surfaces5202 !$OMP PARALLEL DO PRIVATE( i, j, k )5203 DO m = 1, bc_h(1)%ns5204 i = bc_h(1)%i(m)5205 j = bc_h(1)%j(m)5206 k = bc_h(1)%k(m)5207 km(k+1,j,i) = km(k,j,i)5208 ENDDO5209 !5210 !-- Model top5211 !$OMP PARALLEL DO5212 DO i = nxlg, nxrg5213 DO j = nysg, nyng5214 km(nzt+1,j,i) = km(nzt,j,i)5215 ENDDO5216 ENDDO5217 !5218 !-- Set Neumann boundary conditions at the outflow boundaries in case of5219 !-- non-cyclic lateral boundaries5220 IF ( outflow_l ) THEN5221 km(:,:,nxl-1) = km(:,:,nxl)5222 ENDIF5223 IF ( outflow_r ) THEN5224 km(:,:,nxr+1) = km(:,:,nxr)5225 ENDIF5226 IF ( outflow_s ) THEN5227 km(:,nys-1,:) = km(:,nys,:)5228 ENDIF5229 IF ( outflow_n ) THEN5230 km(:,nyn+1,:) = km(:,nyn,:)5231 ENDIF5232 5127 5233 5128 END SUBROUTINE tcm_diffusivities_dynamic
Note: See TracChangeset
for help on using the changeset viewer.
