Changeset 3130 for palm/trunk/SOURCE/turbulence_closure_mod.f90
 Timestamp:
 Jul 16, 2018 11:08:55 AM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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 20180716 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 noncyclic 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, upwardfacing 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/vgrid 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(k1) * & 1943 1941 ( zu(k+1)  zu(k1) ) / & 1944 ( km _sfc+ 1.0E20_wp )1942 ( km(k,j,i) + 1.0E20_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(k1) * & 1947 1945 ( zu(k+1)  zu(k1) ) / & 1948 ( km _sfc+ 1.0E20_wp )1946 ( km(k,j,i) + 1.0E20_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/vgrid is necessary. However, the 1969 ! effect of this error is negligible. 1968 1970 surf_def_h(1)%u_0(m) = u(k1,j,i)  surf_def_h(1)%usws(m) * & 1969 1971 drho_air_zw(k1) * & … … 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/vgrid 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(k1) * & 2005 ( zu(k+1)  zu(k1) ) /&2006 ( km _sfc+ 1.0E20_wp )2003 ( zu(k+1)  zu(k1) ) / & 2004 ( km(k,j,i) + 1.0E20_wp ) 2007 2005 surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & 2008 2006 drho_air_zw(k1) * & 2009 2007 ( zu(k+1)  zu(k1) ) / & 2010 ( km _sfc+ 1.0E20_wp )2008 ( km(k,j,i) + 1.0E20_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/vgrid 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(k1) * & 2039 2034 ( zu(k+1)  zu(k1) ) / & 2040 ( km _sfc+ 1.0E20_wp )2035 ( km(k,j,i) + 1.0E20_wp ) 2041 2036 surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * & 2042 2037 drho_air_zw(k1) * & 2043 2038 ( zu(k+1)  zu(k1) ) / & 2044 ( km _sfc+ 1.0E20_wp )2039 ( km(k,j,i) + 1.0E20_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 PrandtlKolmogorov.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 TKEl 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.0E30_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.0E30_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 boundarycondition 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 wallboundary 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 Prandtllayer where the boundary 4806 ! values of the diffusivities are not needed. 4807 IF ( .NOT. rans_tke_e ) THEN 4808 ! 4809 ! Upwardfacing 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(k1,j,i) = km(k,j,i) 4816 kh(k1,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 PrandtlKolmogorov. 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 TKEl 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.0E30_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.0E30_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 boundarycondition 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 ucomponent in xdirection … … 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 wallboundary 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 constantflux layer where the5189 ! boundary values of the diffusivities are not needed.5190 5191 !5192 ! Upwardfacing5193 !$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(k1,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 ! noncyclic lateral boundaries5220 IF ( outflow_l ) THEN5221 km(:,:,nxl1) = km(:,:,nxl)5222 ENDIF5223 IF ( outflow_r ) THEN5224 km(:,:,nxr+1) = km(:,:,nxr)5225 ENDIF5226 IF ( outflow_s ) THEN5227 km(:,nys1,:) = 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.