Changeset 3776
 Timestamp:
 Mar 4, 2019 1:06:12 PM (5 years ago)
 File:

 1 edited
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 20190304 12:40:20Z gronemeier 27 30 ! removed unused variables 28 ! 31 ! 29 32 ! 3724 20190206 16:28:23Z kanani 30 33 ! Correct doubleused log_point_s units 31 ! 34 ! 32 35 ! 3719 20190206 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 20190120 20:20:58Z knoop 37 40 ! Remove unused variable simulated_time 38 ! 41 ! 39 42 ! 3636 20181219 13:48:34Z raasch 40 43 ! nopointer option removed 41 ! 44 ! 42 45 ! 3634 20181218 12:31:28Z knoop 43 46 ! OpenACC port for SPEC 44 ! 47 ! 45 48 ! 3550 20181121 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 20181121 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 20181025 13:36:23Z maronga 54 57 ! Added support for buildings in the dynamic SGS model 55 ! 58 ! 56 59 ! 3398 20181022 19:30:24Z knoop 57 60 ! Refactored production_e and production_e_ij (removed code redundancy) 58 ! 61 ! 59 62 ! 3386 20181019 16:28:22Z gronemeier 60 63 ! Renamed tcm_prognostic to tcm_prognostic_equations 61 ! 64 ! 62 65 ! 3385 20181019 14:52:29Z knoop 63 66 ! Restructured loops and ifs in production_e to ease vectorization on GPUs 64 ! 67 ! 65 68 ! 3300 20181002 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 20181001 02:37:10Z raasch 72 75 ! changes concerning modularization of ocean option 73 ! 76 ! 74 77 ! 3274 20180924 15:42:55Z knoop 75 78 ! Modularization of all bulk cloud physics code components 76 ! 79 ! 77 80 ! 3245 20180913 14:08:16Z knoop 78 81 ! unused variables removed, shortest_distance has wp now 79 ! 82 ! 80 83 ! 3183 20180727 14:25:55Z suehring 81 84 ! Rename variables and remove unused variable from USE statement 82 ! 85 ! 83 86 ! 3182 20180727 13:36:03Z suehring 84 87 ! Use MOST for km only in RANS mode 85 ! 88 ! 86 89 ! 3130 20180716 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 20180716 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 1dmodel, do not use e1d, km1d and diss1d 98 ! 101 ! 99 102 ! 3121 20180711 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 20180711 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 20180625 09:08:04Z gronemeier 110 113 ! bugfix: set rans_const_sigma(1) = 1.3 111 ! 114 ! 112 115 ! 3083 20180619 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 20180528 07:55:41Z Giersch 134 137 ! Error message revised 135 ! 138 ! 136 139 ! 3014 20180509 08:42:38Z maronga 137 140 ! Bugfix: nzb_do and nzt_do were not used for 3d data output 138 ! 141 ! 139 142 ! 3004 20180427 12:33:25Z Giersch 140 143 ! Further allocation checks implemented 141 ! 144 ! 142 145 ! 2938 20180327 15:52:42Z suehring 143 146 ! Further todo's 144 ! 147 ! 145 148 ! 2936 20180327 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 20180122 09:25:36Z gronemeier 157 160 ! Bugfix: remove duplicate SAVE statements 158 ! 161 ! 159 162 ! 2746 20180115 12:06:04Z suehring 160 163 ! Move flag plant canopy to modules 161 ! 164 ! 162 165 ! 2718 20180102 08:49:38Z maronga 163 166 ! Corrected "Former revisions" section 164 ! 167 ! 165 168 ! 2701 20171215 15:40:50Z suehring 166 169 ! Changes from last commit documented 167 ! 170 ! 168 171 ! 2698 20171214 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 !< nearwall 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 ! xyplane 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 ! xzdistance (vertical edges, down) 1348 1351 IF ( .NOT. BTEST( wall_flags_0(k1,j,i1), 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(krad_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 zdirection 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/vgrid is necessary. However, the 1789 1792 ! effect of this error is negligible. … … 1798 1801 drho_air_zw(k1) * & 1799 1802 ( zu(k+1)  zu(k1) ) / & 1800 ( km_sfc + 1.0E20_wp ) 1803 ( km_sfc + 1.0E20_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/vgrid is necessary. However, the 1825 1828 ! effect of this error is negligible. … … 1827 1830 drho_air_zw(k1) * & 1828 1831 ( zu(k+1)  zu(k1) ) / & 1829 ( km(k,j,i) + 1.0E20_wp ) 1832 ( km(k,j,i) + 1.0E20_wp ) 1830 1833 surf_def_h(1)%v_0(m) = v(k1,j,i)  surf_def_h(1)%vsws(m) * & 1831 1834 drho_air_zw(k1) * & 1832 1835 ( zu(k+1)  zu(k1) ) / & 1833 ( km(k,j,i) + 1.0E20_wp ) 1836 ( km(k,j,i) + 1.0E20_wp ) 1834 1837 1835 1838 IF ( ABS( surf_def_h(1)%u_0(m)  u(k1,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/vgrid is necessary. However, the 1858 1861 ! effect of this error is negligible. … … 1863 1866 drho_air_zw(k1) * & 1864 1867 ( zu(k+1)  zu(k1) ) / & 1865 ( km_sfc + 1.0E20_wp ) 1868 ( km_sfc + 1.0E20_wp ) 1866 1869 surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & 1867 1870 drho_air_zw(k1) * & … … 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/vgrid 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 * Exnerfunction 2367 REAL(wp) :: sign_dir !< sign of walltke flux, depending on wall orientation 2370 REAL(wp) :: sign_dir !< sign of walltke 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 PrandtlKolmogorov. 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 subgridscale 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 subgridscale stress models via 4610 !> stochastic analysis." 4608 4611 !> Monte Carlo Methods and Applications 14.4 (2008): 311329. 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.0E20_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.