- Timestamp:
- Mar 4, 2019 1:06:12 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/turbulence_closure_mod.f90 (modified) (71 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/turbulence_closure_mod.f90
r3775 r3776 25 25 ! ----------------- 26 26 ! $Id$ 27 ! write out preprocessor directives; remove tailing whitespaces 28 ! 29 ! 3775 2019-03-04 12:40:20Z gronemeier 27 30 ! removed unused variables 28 ! 31 ! 29 32 ! 3724 2019-02-06 16:28:23Z kanani 30 33 ! Correct double-used log_point_s units 31 ! 34 ! 32 35 ! 3719 2019-02-06 13:10:18Z kanani 33 ! Changed log_point to log_point_s, otherwise this overlaps with 36 ! Changed log_point to log_point_s, otherwise this overlaps with 34 37 ! 'all progn.equations' cpu measurement. 35 ! 38 ! 36 39 ! 3684 2019-01-20 20:20:58Z knoop 37 40 ! Remove unused variable simulated_time 38 ! 41 ! 39 42 ! 3636 2018-12-19 13:48:34Z raasch 40 43 ! nopointer option removed 41 ! 44 ! 42 45 ! 3634 2018-12-18 12:31:28Z knoop 43 46 ! OpenACC port for SPEC 44 ! 47 ! 45 48 ! 3550 2018-11-21 16:01:01Z gronemeier 46 49 ! - calculate diss production same in vector and cache optimization 47 50 ! - move boundary condition for e and diss to boundary_conds 48 ! 51 ! 49 52 ! 3545 2018-11-21 11:19:41Z gronemeier 50 53 ! - Set rans_mode according to value of turbulence_closure 51 54 ! - removed debug output 52 ! 55 ! 53 56 ! 3430 2018-10-25 13:36:23Z maronga 54 57 ! Added support for buildings in the dynamic SGS model 55 ! 58 ! 56 59 ! 3398 2018-10-22 19:30:24Z knoop 57 60 ! Refactored production_e and production_e_ij (removed code redundancy) 58 ! 61 ! 59 62 ! 3386 2018-10-19 16:28:22Z gronemeier 60 63 ! Renamed tcm_prognostic to tcm_prognostic_equations 61 ! 64 ! 62 65 ! 3385 2018-10-19 14:52:29Z knoop 63 66 ! Restructured loops and ifs in production_e to ease vectorization on GPUs 64 ! 67 ! 65 68 ! 3300 2018-10-02 14:16:54Z gronemeier 66 69 ! - removed global array wall_flags_0_global, hence reduced accuracy of l_wall … … 68 71 ! - removed maxloc call as this produced different results for different 69 72 ! compiler options 70 ! 73 ! 71 74 ! 3294 2018-10-01 02:37:10Z raasch 72 75 ! changes concerning modularization of ocean option 73 ! 76 ! 74 77 ! 3274 2018-09-24 15:42:55Z knoop 75 78 ! Modularization of all bulk cloud physics code components 76 ! 79 ! 77 80 ! 3245 2018-09-13 14:08:16Z knoop 78 81 ! unused variables removed, shortest_distance has wp now 79 ! 82 ! 80 83 ! 3183 2018-07-27 14:25:55Z suehring 81 84 ! Rename variables and remove unused variable from USE statement 82 ! 85 ! 83 86 ! 3182 2018-07-27 13:36:03Z suehring 84 87 ! Use MOST for km only in RANS mode 85 ! 88 ! 86 89 ! 3130 2018-07-16 11:08:55Z gronemeier 87 90 ! - move boundary condition of km and kh to tcm_diffusivities 88 91 ! - calculate km at boundaries according to MOST 89 92 ! - move phi_m to surface_layer_fluxes_mod 90 ! 93 ! 91 94 ! 3129 2018-07-16 07:45:13Z gronemeier 92 95 ! - move limitation of diss to boundary_conds … … 96 99 ! - resort output after case select -> reduced code duplication 97 100 ! - when using rans_tke_e and 1d-model, do not use e1d, km1d and diss1d 98 ! 101 ! 99 102 ! 3121 2018-07-11 18:46:49Z gronemeier 100 103 ! - created the function phi_m 101 104 ! - implemented km = u* * kappa * zp / phi_m in production_e_init for all 102 105 ! surfaces 103 ! 106 ! 104 107 ! 3120 2018-07-11 18:30:57Z gronemeier 105 108 ! - changed tcm_diffusivities to tcm_diffusivities_default 106 109 ! - created subroutine tcm_diffusivities that calls tcm_diffusivities_default 107 110 ! and tcm_diffusivities_dynamic 108 ! 111 ! 109 112 ! 3086 2018-06-25 09:08:04Z gronemeier 110 113 ! bugfix: set rans_const_sigma(1) = 1.3 111 ! 114 ! 112 115 ! 3083 2018-06-19 14:03:12Z gronemeier 113 116 ! - set limits of diss at the end of prognostic equations … … 130 133 ! - calculate km based on l_wall 131 134 ! - initialize diss if 1D model is not used 132 ! 135 ! 133 136 ! 3045 2018-05-28 07:55:41Z Giersch 134 137 ! Error message revised 135 ! 138 ! 136 139 ! 3014 2018-05-09 08:42:38Z maronga 137 140 ! Bugfix: nzb_do and nzt_do were not used for 3d data output 138 ! 141 ! 139 142 ! 3004 2018-04-27 12:33:25Z Giersch 140 143 ! Further allocation checks implemented 141 ! 144 ! 142 145 ! 2938 2018-03-27 15:52:42Z suehring 143 146 ! Further todo's 144 ! 147 ! 145 148 ! 2936 2018-03-27 14:49:27Z gronemeier 146 149 ! - defined l_grid only within this module … … 153 156 ! - Moved init_mixing_length from init_grid.f90 and 154 157 ! renamed it to tcm_init_mixing_length 155 ! 158 ! 156 159 ! 2764 2018-01-22 09:25:36Z gronemeier 157 160 ! Bugfix: remove duplicate SAVE statements 158 ! 161 ! 159 162 ! 2746 2018-01-15 12:06:04Z suehring 160 163 ! Move flag plant canopy to modules 161 ! 164 ! 162 165 ! 2718 2018-01-02 08:49:38Z maronga 163 166 ! Corrected "Former revisions" section 164 ! 167 ! 165 168 ! 2701 2017-12-15 15:40:50Z suehring 166 169 ! Changes from last commit documented 167 ! 170 ! 168 171 ! 2698 2017-12-14 18:46:24Z suehring 169 172 ! Bugfix in get_topography_top_index … … 172 175 ! Initial revision 173 176 ! 174 ! 175 ! 177 ! 178 ! 176 179 ! 177 180 ! Authors: … … 191 194 !------------------------------------------------------------------------------! 192 195 MODULE turbulence_closure_mod 193 196 194 197 195 198 USE arrays_3d, & … … 264 267 REAL(wp), DIMENSION(:), ALLOCATABLE :: l_black !< mixing length according to Blackadar 265 268 266 REAL(wp), DIMENSION(:), ALLOCATABLE :: l_grid !< geometric mean of grid sizes dx, dy, dz 269 REAL(wp), DIMENSION(:), ALLOCATABLE :: l_grid !< geometric mean of grid sizes dx, dy, dz 267 270 268 271 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: l_wall !< near-wall mixing length … … 309 312 310 313 ! 311 !-- Initialization actions 314 !-- Initialization actions 312 315 INTERFACE tcm_init 313 316 MODULE PROCEDURE tcm_init … … 476 479 c_0 = 0.1_wp !according to Lilly (1967) and Deardorff (1980) 477 480 478 dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead 481 dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead 479 482 !of K_e which is used in RANS mode 480 483 … … 489 492 !------------------------------------------------------------------------------! 490 493 SUBROUTINE tcm_check_data_output( var, unit ) 491 494 492 495 IMPLICIT NONE 493 496 … … 519 522 !------------------------------------------------------------------------------! 520 523 SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) 521 524 522 525 IMPLICIT NONE 523 526 … … 567 570 !------------------------------------------------------------------------------! 568 571 SUBROUTINE tcm_3d_data_averaging( mode, variable ) 569 572 570 573 571 574 USE averaging, & … … 616 619 617 620 CASE ( 'diss' ) 618 IF ( ALLOCATED( diss_av ) ) THEN 621 IF ( ALLOCATED( diss_av ) ) THEN 619 622 DO i = nxlg, nxrg 620 623 DO j = nysg, nyng … … 662 665 DO j = nysg, nyng 663 666 DO k = nzb, nzt+1 664 diss_av(k,j,i) = diss_av(k,j,i) & 667 diss_av(k,j,i) = diss_av(k,j,i) & 665 668 / REAL( average_count_3d, KIND=wp ) 666 669 ENDDO … … 674 677 DO j = nysg, nyng 675 678 DO k = nzb, nzt+1 676 kh_av(k,j,i) = kh_av(k,j,i) & 679 kh_av(k,j,i) = kh_av(k,j,i) & 677 680 / REAL( average_count_3d, KIND=wp ) 678 681 ENDDO … … 686 689 DO j = nysg, nyng 687 690 DO k = nzb, nzt+1 688 km_av(k,j,i) = km_av(k,j,i) & 691 km_av(k,j,i) = km_av(k,j,i) & 689 692 / REAL( average_count_3d, KIND=wp ) 690 693 ENDDO … … 707 710 SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, & 708 711 nzb_do, nzt_do ) 709 712 710 713 USE averaging, & 711 714 ONLY: diss_av, kh_av, km_av … … 734 737 735 738 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable 736 739 737 740 found = .TRUE. 738 741 resorted = .FALSE. … … 791 794 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 792 795 REAL( fill_value, KIND = wp ), & 793 BTEST( wall_flags_0(k,j,i), flag_nr ) ) 796 BTEST( wall_flags_0(k,j,i), flag_nr ) ) 794 797 ENDDO 795 798 ENDDO 796 799 ENDDO 797 800 ENDIF 798 801 799 802 END SUBROUTINE tcm_data_output_2d 800 803 801 804 802 805 !------------------------------------------------------------------------------! 803 806 ! Description: … … 806 809 !------------------------------------------------------------------------------! 807 810 SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 808 811 809 812 810 813 USE averaging, & … … 873 876 to_be_resorted => km_av 874 877 ENDIF 875 878 876 879 CASE DEFAULT 877 880 found = .FALSE. … … 923 926 924 927 ! 925 !-- Allocate arrays required for dissipation. 928 !-- Allocate arrays required for dissipation. 926 929 !-- Please note, if it is a nested run, arrays need to be allocated even if 927 !-- they do not necessarily need to be transferred, which is attributed to 928 !-- the design of the model coupler which allocates memory for each variable. 930 !-- they do not necessarily need to be transferred, which is attributed to 931 !-- the design of the model coupler which allocates memory for each variable. 929 932 IF ( rans_mode .OR. use_sgs_for_particles .OR. wang_kernel .OR. & 930 933 collision_turbulence .OR. nested_run ) THEN … … 1025 1028 ELSE 1026 1029 IF ( .NOT. ocean_mode ) THEN 1027 kh = 0.01_wp ! there must exist an initial diffusion, because 1030 kh = 0.01_wp ! there must exist an initial diffusion, because 1028 1031 km = 0.01_wp ! otherwise no TKE would be produced by the 1029 1032 ! production terms, as long as not yet … … 1069 1072 ELSE 1070 1073 IF ( .NOT. ocean_mode ) THEN 1071 kh = 0.01_wp ! there must exist an initial diffusion, because 1074 kh = 0.01_wp ! there must exist an initial diffusion, because 1072 1075 km = 0.01_wp ! otherwise no TKE would be produced by the 1073 1076 ! production terms, as long as not yet … … 1112 1115 ! 1113 1116 !-- In case of complex terrain and cyclic fill method as initialization, 1114 !-- shift initial data in the vertical direction for each point in the 1117 !-- shift initial data in the vertical direction for each point in the 1115 1118 !-- x-y-plane depending on local surface height 1116 1119 IF ( complex_terrain .AND. & … … 1142 1145 mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e 1143 1146 ! 1144 !-- In case of complex terrain, determine vertical displacement at inflow 1147 !-- In case of complex terrain, determine vertical displacement at inflow 1145 1148 !-- boundary and adjust mean inflow profiles 1146 1149 IF ( complex_terrain ) THEN … … 1344 1347 SQRT( 0.25_wp * dy**2 + & 1345 1348 ( zw(k) - zu(k) )**2 ) ) 1346 ! 1349 ! 1347 1350 !-- xz-distance (vertical edges, down) 1348 1351 IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 ) .OR. & … … 1384 1387 SQRT( 0.25_wp * ( dx**2 + dy**2 ) & 1385 1388 + ( zw(k) - zu(k) )**2 ) ) 1386 1389 1387 1390 ENDIF 1388 1391 ENDDO … … 1426 1429 ! 1427 1430 !-- Limit mixing length to either nearest wall or Blackadar mixing length. 1428 !-- For that, analyze each grid point (i/j/k) ("analysed grid point") and 1431 !-- For that, analyze each grid point (i/j/k) ("analysed grid point") and 1429 1432 !-- search within its vicinity for the shortest distance to a wall by cal- 1430 1433 !-- culating the distance between the analysed grid point and the "viewed … … 1459 1462 ! 1460 1463 !-- When analysed grid point lies above maximum topography, set search 1461 !-- radius to 0 if the distance between the analysed grid point and max 1464 !-- radius to 0 if the distance between the analysed grid point and max 1462 1465 !-- topography height is larger than the maximum search radius 1463 1466 IF ( zu(k-rad_k_b) > zu(k_max_topo) ) rad_k_b = 0 … … 1628 1631 ! Description: 1629 1632 ! ------------ 1630 !> Calculate the shortest distance between position (i/j/k)=(0/0/0) and 1633 !> Calculate the shortest distance between position (i/j/k)=(0/0/0) and 1631 1634 !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array' 1632 1635 !> closest to the origin (0/0) of 'array'. … … 1694 1697 ENDDO 1695 1698 ENDIF 1696 1699 1697 1700 END FUNCTION 1698 1701 … … 1700 1703 ! Description: 1701 1704 ! ------------ 1702 !> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point 1705 !> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point 1703 1706 !> (kp,jp,ip) containing the first bit of wall_flags_0 into the array 1704 1707 !> 'vicinity'. Only copy first bit as this indicates the presence of topography. … … 1753 1756 1754 1757 USE surface_mod, & 1755 ONLY : surf_def_h, surf_lsm_h, surf_usm_h 1758 ONLY : surf_def_h, surf_lsm_h, surf_usm_h 1756 1759 1757 1760 IMPLICIT NONE … … 1761 1764 INTEGER(iwp) :: k !< grid index z-direction 1762 1765 INTEGER(iwp) :: m !< running index surface elements 1763 1766 1764 1767 REAL(wp) :: km_sfc !< diffusion coefficient, used to compute virtual velocities 1765 1768 … … 1780 1783 DO m = 1, surf_def_h(0)%ns 1781 1784 1782 i = surf_def_h(0)%i(m) 1785 i = surf_def_h(0)%i(m) 1783 1786 j = surf_def_h(0)%j(m) 1784 1787 k = surf_def_h(0)%k(m) 1785 1788 ! 1786 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1787 !-- and km are not on the same grid. Actually, a further 1789 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1790 !-- and km are not on the same grid. Actually, a further 1788 1791 !-- interpolation of km onto the u/v-grid is necessary. However, the 1789 1792 !-- effect of this error is negligible. … … 1798 1801 drho_air_zw(k-1) * & 1799 1802 ( zu(k+1) - zu(k-1) ) / & 1800 ( km_sfc + 1.0E-20_wp ) 1803 ( km_sfc + 1.0E-20_wp ) 1801 1804 1802 1805 IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > & … … 1816 1819 DO m = 1, surf_def_h(1)%ns 1817 1820 1818 i = surf_def_h(1)%i(m) 1821 i = surf_def_h(1)%i(m) 1819 1822 j = surf_def_h(1)%j(m) 1820 1823 k = surf_def_h(1)%k(m) 1821 1824 ! 1822 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1823 !-- and km are not on the same grid. Actually, a further 1825 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1826 !-- and km are not on the same grid. Actually, a further 1824 1827 !-- interpolation of km onto the u/v-grid is necessary. However, the 1825 1828 !-- effect of this error is negligible. … … 1827 1830 drho_air_zw(k-1) * & 1828 1831 ( zu(k+1) - zu(k-1) ) / & 1829 ( km(k,j,i) + 1.0E-20_wp ) 1832 ( km(k,j,i) + 1.0E-20_wp ) 1830 1833 surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * & 1831 1834 drho_air_zw(k-1) * & 1832 1835 ( zu(k+1) - zu(k-1) ) / & 1833 ( km(k,j,i) + 1.0E-20_wp ) 1836 ( km(k,j,i) + 1.0E-20_wp ) 1834 1837 1835 1838 IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > & … … 1853 1856 k = surf_lsm_h%k(m) 1854 1857 ! 1855 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1856 !-- and km are not on the same grid. Actually, a further 1858 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1859 !-- and km are not on the same grid. Actually, a further 1857 1860 !-- interpolation of km onto the u/v-grid is necessary. However, the 1858 1861 !-- effect of this error is negligible. … … 1863 1866 drho_air_zw(k-1) * & 1864 1867 ( zu(k+1) - zu(k-1) ) / & 1865 ( km_sfc + 1.0E-20_wp ) 1868 ( km_sfc + 1.0E-20_wp ) 1866 1869 surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & 1867 1870 drho_air_zw(k-1) * & … … 1889 1892 k = surf_usm_h%k(m) 1890 1893 ! 1891 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1892 !-- and km are not on the same grid. Actually, a further 1894 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 1895 !-- and km are not on the same grid. Actually, a further 1893 1896 !-- interpolation of km onto the u/v-grid is necessary. However, the 1894 1897 !-- effect of this error is negligible. … … 2365 2368 REAL(wp) :: theta !< virtual potential temperature 2366 2369 REAL(wp) :: temp !< theta * Exner-function 2367 REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation 2370 REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation 2368 2371 REAL(wp) :: usvs !< momentum flux u"v" 2369 2372 REAL(wp) :: vsus !< momentum flux v"u" … … 3669 3672 SUBROUTINE diffusion_e( var, var_reference ) 3670 3673 3671 #if def _OPENACC3674 #if defined( _OPENACC ) 3672 3675 USE arrays_3d, & 3673 3676 ONLY: ddzu, dd2zu, ddzw, drho_air, rho_air_zw … … 3703 3706 REAL(wp) :: ll !< adjusted l 3704 3707 REAL(wp) :: var_reference !< reference temperature 3705 #if def _OPENACC3708 #if defined( _OPENACC ) 3706 3709 ! 3707 3710 !-- From mixing_length_les: … … 3742 3745 IF ( les_dynamic .OR. les_mw ) THEN 3743 3746 3744 #if def _OPENACC3747 #if defined( _OPENACC ) 3745 3748 ! 3746 3749 !-- Cannot call subroutine mixing_length_les because after adding all required … … 3777 3780 ELSEIF ( rans_tke_l ) THEN 3778 3781 3779 #if def _OPENACC3782 #if defined( _OPENACC ) 3780 3783 ! 3781 3784 !-- Same reason as above... … … 3857 3860 !-- Upward facing surfaces 3858 3861 DO m = 1, bc_h(0)%ns 3859 i = bc_h(0)%i(m) 3862 i = bc_h(0)%i(m) 3860 3863 j = bc_h(0)%j(m) 3861 3864 k = bc_h(0)%k(m) … … 3865 3868 !-- Downward facing surfaces 3866 3869 DO m = 1, bc_h(1)%ns 3867 i = bc_h(1)%i(m) 3870 i = bc_h(1)%i(m) 3868 3871 j = bc_h(1)%j(m) 3869 3872 k = bc_h(1)%k(m) … … 3889 3892 USE grid_variables, & 3890 3893 ONLY: ddx2, ddy2 3891 3894 3892 3895 USE bulk_cloud_model_mod, & 3893 3896 ONLY: collision_turbulence … … 3984 3987 !-- For each surface type determine start and end index (in case of elevated 3985 3988 !-- topography several up/downward facing surfaces may exist. 3986 surf_s = bc_h(0)%start_index(j,i) 3987 surf_e = bc_h(0)%end_index(j,i) 3989 surf_s = bc_h(0)%start_index(j,i) 3990 surf_e = bc_h(0)%end_index(j,i) 3988 3991 DO m = surf_s, surf_e 3989 3992 k = bc_h(0)%k(m) … … 3992 3995 ! 3993 3996 !-- Downward facing surfaces 3994 surf_s = bc_h(1)%start_index(j,i) 3995 surf_e = bc_h(1)%end_index(j,i) 3997 surf_s = bc_h(1)%start_index(j,i) 3998 surf_e = bc_h(1)%end_index(j,i) 3996 3999 DO m = surf_s, surf_e 3997 4000 k = bc_h(1)%k(m) … … 4298 4301 kh(k,j,i) = 1.35_wp * km(k,j,i) 4299 4302 ENDDO 4300 ! 4303 ! 4301 4304 !-- Natural surfaces 4302 4305 !$OMP PARALLEL DO PRIVATE(i,j,k,m) … … 4309 4312 kh(k,j,i) = 1.35_wp * km(k,j,i) 4310 4313 ENDDO 4311 ! 4314 ! 4312 4315 !-- Urban surfaces 4313 4316 !$OMP PARALLEL DO PRIVATE(i,j,k,m) … … 4320 4323 kh(k,j,i) = 1.35_wp * km(k,j,i) 4321 4324 ENDDO 4322 4323 ! 4325 4326 ! 4324 4327 !-- North-, south-, west and eastward facing surfaces 4325 4328 !-- Do not consider stratification at these surfaces. 4326 4329 DO n = 0, 3 4327 ! 4330 ! 4328 4331 !-- Default surfaces 4329 4332 !$OMP PARALLEL DO PRIVATE(i,j,k,m) … … 4335 4338 kh(k,j,i) = 1.35_wp * km(k,j,i) 4336 4339 ENDDO 4337 ! 4340 ! 4338 4341 !-- Natural surfaces 4339 4342 !$OMP PARALLEL DO PRIVATE(i,j,k,m) … … 4345 4348 kh(k,j,i) = 1.35_wp * km(k,j,i) 4346 4349 ENDDO 4347 ! 4350 ! 4348 4351 !-- Urban surfaces 4349 4352 !$OMP PARALLEL DO PRIVATE(i,j,k,m) … … 4356 4359 ENDDO 4357 4360 ENDDO 4358 4361 4359 4362 CALL exchange_horiz( km, nbgp ) 4360 4363 CALL exchange_horiz( kh, nbgp ) … … 4424 4427 ! Description: 4425 4428 ! ------------ 4426 !> Computation of the turbulent diffusion coefficients for momentum and heat 4429 !> Computation of the turbulent diffusion coefficients for momentum and heat 4427 4430 !> according to Prandtl-Kolmogorov. 4428 4431 !------------------------------------------------------------------------------! 4429 4432 SUBROUTINE tcm_diffusivities_default( var, var_reference ) 4430 4431 #if def _OPENACC4433 4434 #if defined( _OPENACC ) 4432 4435 USE arrays_3d, & 4433 4436 ONLY: dd2zu … … 4454 4457 REAL(wp) :: ll !< adjusted mixing length 4455 4458 REAL(wp) :: var_reference !< reference temperature 4456 #if def _OPENACC4459 #if defined( _OPENACC ) 4457 4460 ! 4458 4461 !-- From mixing_length_les: … … 4493 4496 ! 4494 4497 !-- Determine the mixing length for LES closure 4495 #if def _OPENACC4498 #if defined( _OPENACC ) 4496 4499 ! 4497 4500 !-- Cannot call subroutine mixing_length_les because after adding all required … … 4602 4605 ! Description: 4603 4606 ! ------------ 4604 !> Calculates the eddy viscosity dynamically using the linear dynamic model 4605 !> according to 4606 !> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via 4607 !> stochastic analysis." 4607 !> Calculates the eddy viscosity dynamically using the linear dynamic model 4608 !> according to 4609 !> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via 4610 !> stochastic analysis." 4608 4611 !> Monte Carlo Methods and Applications 14.4 (2008): 311-329. 4609 4612 !> 4610 !> Furthermore dynamic bounds are used to limit the absolute value of c* as 4613 !> Furthermore dynamic bounds are used to limit the absolute value of c* as 4611 4614 !> described in 4612 !> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation: 4615 !> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation: 4613 4616 !> Stability via realizability." Physics of Fluids 29.10 (2017): 105104. 4614 4617 !> … … 4617 4620 !------------------------------------------------------------------------------! 4618 4621 SUBROUTINE tcm_diffusivities_dynamic 4619 4622 4620 4623 USE arrays_3d, & 4621 4624 ONLY: ddzw, dzw, dd2zu, w, ug, vg 4622 4625 4623 4626 USE grid_variables, & 4624 4627 ONLY : ddx, ddy, dx, dy … … 4643 4646 4644 4647 REAL(wp) :: flag !< topography flag 4645 4648 4646 4649 REAL(wp) :: uc(-1:1,-1:1) !< u on grid center 4647 4650 REAL(wp) :: vc(-1:1,-1:1) !< v on grid center … … 4706 4709 DO j = nys, nyn 4707 4710 DO k = nzb+1, nzt 4708 4711 4709 4712 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 4710 4713 4711 4714 ! 4712 4715 !-- Compute the deviatoric shear tensor s_ij on grid centers: … … 4829 4832 cst_max = fac_cmax * SQRT( e(k,j,i) ) / & 4830 4833 ( delta * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp ) 4831 4834 4832 4835 IF ( ABS( cst ) > cst_max ) THEN 4833 4836 cst = cst_max * cst / ABS( cst ) … … 4835 4838 4836 4839 km(k,j,i) = cst * delta * SQRT( e(k,j,i) ) * flag 4837 4840 4838 4841 ENDDO 4839 4842 ENDDO … … 4850 4853 !------------------------------------------------------------------------------! 4851 4854 SUBROUTINE tcm_box_filter_2d_single( var, var_fil ) 4852 4855 4853 4856 IMPLICIT NONE 4854 4857 … … 4856 4859 REAL(wp) :: var_fil !< filtered variable 4857 4860 ! 4858 !-- It is assumed that a box with a side length of 2 * dx and centered at the 4861 !-- It is assumed that a box with a side length of 2 * dx and centered at the 4859 4862 !-- variable*s position contains one half of the four closest neigbours and one 4860 4863 !-- forth of the diagonally closest neighbours. … … 4884 4887 REAL(wp) :: var_fil(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< filtered variable 4885 4888 ! 4886 !-- It is assumed that a box with a side length of 2 * dx and centered at the 4889 !-- It is assumed that a box with a side length of 2 * dx and centered at the 4887 4890 !-- variable's position contains one half of the four closest neigbours and one 4888 4891 !-- forth of the diagonally closest neighbours.
Note: See TracChangeset
for help on using the changeset viewer.
