- Timestamp:
- Mar 27, 2020 12:56:41 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
modules.f90 (modified) (3 diffs)
-
turbulence_closure_mod.f90 (modified) (45 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/modules.f90
r4473 r4476 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Renamed variables for subgrids-scale model steering 28 ! 29 ! 4473 2020-03-25 21:04:07Z gronemeier 27 30 ! moved wall_adjustment_factor to turbulence_closure_mod 28 31 ! … … 548 551 CHARACTER (LEN=20) :: reference_state = 'initial_profile' !< namelist parameter 549 552 CHARACTER (LEN=20) :: timestep_scheme = 'runge-kutta-3' !< namelist parameter 550 CHARACTER (LEN=20) :: turbulence_closure = ' Moeng_Wyngaard'!< namelist parameter553 CHARACTER (LEN=20) :: turbulence_closure = '1.5-order' !< namelist parameter 551 554 CHARACTER (LEN=40) :: topography = 'flat' !< namelist parameter 552 555 CHARACTER (LEN=64) :: host = '????' !< configuration identifier as given by palmrun option -c, ENVPAR namelist parameter provided by palmrun … … 750 753 LOGICAL :: large_scale_subsidence = .FALSE. !< namelist parameter 751 754 LOGICAL :: land_surface = .FALSE. !< use land surface model? 755 LOGICAL :: les_dai = .FALSE. !< use Dai et al. turbulence closure (modified 1.5-order closure) for LES mode. Shall replace the default 1.5-order closure 752 756 LOGICAL :: les_dynamic = .FALSE. !< use dynamic subgrid model as turbulence closure for LES mode 753 LOGICAL :: les_ mw = .FALSE. !< use Moeng-Wyngaardturbulence closure for LES mode757 LOGICAL :: les_default = .FALSE. !< use 1.5-order default turbulence closure for LES mode 754 758 LOGICAL :: lsf_exception = .FALSE. !< use of lsf with buildings (temporary)? 755 759 LOGICAL :: lsf_surf = .TRUE. !< use surface forcing (large scale forcing)? -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r4473 r4476 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - added new LES closure after Dai et al. (2020), which provides much better grid 28 ! convergence in stable boundary layer runs. The implementation is experimental 29 ! at the moment and should be used with special care. The new SGS closure can be 30 ! switched on via turbulence_closure = '1.5-order-dai' 31 ! - variable ml_wall_adjusted renamed to delta as it represents a grid size and 32 ! not a mixing length (see Equations 14 and 18 in Maronga et al. 2015, GMD) 33 ! - nameing of turbulence closures revised: 34 ! 'Moeng_Wyngaard' to '1.5-order' 35 ! 'TKE-l' to 'tke-l' 36 ! 'TKE-e' to 'tke-e' 37 ! - LOGICAL steering variable renamed: 38 ! les_mw to les_default 39 ! 40 ! 4473 2020-03-25 21:04:07Z gronemeier 27 41 ! - rename l-grid to gridsize-geometric-mean 28 42 ! l-wall to ml-wall-adjusted … … 140 154 initializing_actions, intermediate_timestep_count, & 141 155 intermediate_timestep_count_max, km_constant, & 142 les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, & 156 les_dai, les_dynamic, les_default, & 157 ocean_mode, plant_canopy, prandtl_number, & 143 158 pt_reference, rans_mode, rans_tke_e, rans_tke_l, & 144 159 timestep_scheme, turbulence_closure, & … … 210 225 211 226 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ml_blackadar !< mixing length according to Blackadar 212 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ml_wall_adjusted !< mixing length adjusted according to near-by walls 227 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: delta !< grid size, possibly limited by wall adjustment factor 228 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: distance_to_wall !< distance to the surface/wall 213 229 214 230 ! … … 584 600 SELECT CASE ( TRIM( turbulence_closure ) ) 585 601 602 CASE ( '1.5-order' ) 603 les_default = .TRUE. 604 605 CASE ( '1.5-order-dai' ) 606 les_dai = .TRUE. 607 586 608 CASE ( 'dynamic' ) 587 609 les_dynamic = .TRUE. 588 610 589 CASE ( 'Moeng_Wyngaard' ) 590 les_mw = .TRUE. 591 592 CASE ( 'TKE-l' ) 611 CASE ( 'tke-l' ) 593 612 rans_tke_l = .TRUE. 594 613 rans_mode = .TRUE. 595 614 596 CASE ( ' TKE-e' )615 CASE ( 'tke-e' ) 597 616 rans_tke_e = .TRUE. 598 617 rans_mode = .TRUE. … … 1157 1176 DO j = nysg, nyng 1158 1177 DO k = nzb+1, nzt 1159 km(k,j,i) = c_0 * SQRT( e_init ) * MIN( ml_wall_adjusted(k,j,i), &1178 km(k,j,i) = c_0 * SQRT( e_init ) * MIN( delta(k,j,i), & 1160 1179 ml_blackadar(k) ) 1161 1180 ENDDO … … 1202 1221 DO j = nysg, nyng 1203 1222 DO k = nzb+1, nzt 1204 km(k,j,i) = c_0 * SQRT( e_init ) * ml_wall_adjusted(k,j,i)1223 km(k,j,i) = c_0 * SQRT( e_init ) * delta(k,j,i) 1205 1224 ENDDO 1206 1225 ENDDO … … 1423 1442 REAL(wp) :: distance_edge_yz_down !< distance of grid box center to its boundary along yz diagonal (down) 1424 1443 REAL(wp) :: distance_edge_yz_up !< distance of grid box center to its boundary along yz diagonal (up) 1425 REAL(wp) :: distance_edge_xz_down !< distance of grid box center to its boundary along xz diagonal (down)1444 REAL(wp) :: distance_edge_xz_down !< distance of grid box center to its boundary along xz diagonal 1426 1445 REAL(wp) :: distance_edge_xz_up !< distance of grid box center to its boundary along xz diagonal (up) 1427 1446 REAL(wp) :: distance_edge_xy !< distance of grid box center to its boundary along xy diagonal … … 1432 1451 REAL(wp), DIMENSION(nzb:nzt+1) :: gridsize_geometric_mean !< geometric mean of grid sizes dx, dy, dz 1433 1452 1434 ! ALLOCATE( gridsize_geometric_mean(nzb:nzt+1) ) 1435 ALLOCATE( ml_wall_adjusted(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1453 1454 ALLOCATE( delta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1455 1456 1457 IF ( les_dai ) THEN 1458 ALLOCATE ( distance_to_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1459 ENDIF 1436 1460 ! 1437 1461 !-- Initialize the mixing length in case of an LES-simulation … … 1454 1478 1455 1479 ! 1456 !-- Initialize near-wall mixing length ml_wall_adjusted1480 !-- Initialize the grid size variable (delta) and the distance to the wall, if needed 1457 1481 DO k = nzb, nzt+1 1458 ml_wall_adjusted(k,:,:) = gridsize_geometric_mean(k) 1482 delta(k,:,:) = gridsize_geometric_mean(k) 1483 1484 ! 1485 !-- If Dai et al. (2020) closure is used, the distance to the wall must be stored 1486 IF ( les_dai ) THEN 1487 distance_to_wall(k,:,:) = zu(k) 1488 ENDIF 1459 1489 ENDDO 1460 1490 … … 1569 1599 distance_corners_up = 9999999.9_wp 1570 1600 ENDIF 1571 ! 1572 !-- Adjust mixing length by wall-adjustment factor 1573 ml_wall_adjusted(k,j,i) = wall_adjustment_factor * MIN( & 1601 1602 ! 1603 !-- Calculate the minimum distance from the wall and store it 1604 !-- temporarily in the array delta 1605 delta(k,j,i) = MIN( & 1574 1606 distance_up, distance_down, distance_ns, distance_lr, & 1575 1607 distance_edge_yz_down, distance_edge_yz_up, & … … 1578 1610 distance_corners_down, distance_corners_up ) 1579 1611 1612 ! 1613 !-- If Dai et al. (2020) closure is used, the distance to the wall 1614 !-- must be permanently stored 1615 IF ( les_dai .AND. delta(k,j,i) /= 9999999.9_wp ) THEN 1616 distance_to_wall(k,j,i) = delta(k,j,i) 1617 ENDIF 1618 1619 ! 1620 !-- Close to the surface, the maximum mixing length is limited to 1.8 z 1621 delta(k,j,i) = wall_adjustment_factor * delta(k,j,i) 1622 1580 1623 ENDIF !if grid point belongs to atmosphere 1581 1624 1582 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), & 1583 gridsize_geometric_mean(k) ) 1625 1626 1627 ! 1628 !-- The grid size (delta) is defined as the the minimum of the distance to 1629 !-- the nearest wall * 1.8 and the geometric mean grid size. 1630 delta(k,j,i) = MIN( delta(k,j,i), gridsize_geometric_mean(k) ) 1631 1584 1632 1585 1633 ENDDO !k loop … … 1622 1670 ENDDO 1623 1671 1624 ml_wall_adjusted(nzb,:,:) = ml_blackadar(nzb)1625 ml_wall_adjusted(nzt+1,:,:) = ml_blackadar(nzt+1)1672 delta(nzb,:,:) = ml_blackadar(nzb) 1673 delta(nzt+1,:,:) = ml_blackadar(nzt+1) 1626 1674 ! 1627 1675 !-- Limit mixing length to either nearest wall or Blackadar mixing length. … … 1634 1682 radius = ml_blackadar(k) ! radius within walls are searched 1635 1683 ! 1636 !-- Set ml_wall_adjustedto its default maximum value (ml_blackadar)1637 ml_wall_adjusted(k,:,:) = radius1684 !-- Set delta to its default maximum value (ml_blackadar) 1685 delta(k,:,:) = radius 1638 1686 1639 1687 ! … … 1730 1778 !-- Search for walls within octant (+++) 1731 1779 vic_yz = vicinity(0:rad_k+1,0:rad_j,ii) 1732 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &1780 delta(k,j,i) = MIN( delta(k,j,i), & 1733 1781 shortest_distance( vic_yz, .TRUE., ii ) ) 1734 1782 ! … … 1738 1786 !-- shortest_distance"). 1739 1787 vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii) 1740 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &1788 delta(k,j,i) = MIN( delta(k,j,i), & 1741 1789 shortest_distance( vic_yz, .TRUE., ii ) ) 1742 1790 … … 1745 1793 !-- Search for walls within octant (+--) 1746 1794 vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii) 1747 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &1795 delta(k,j,i) = MIN( delta(k,j,i), & 1748 1796 shortest_distance( vic_yz, .FALSE., ii ) ) 1749 1797 ! 1750 1798 !-- Search for walls within octant (++-) 1751 1799 vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii) 1752 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &1800 delta(k,j,i) = MIN( delta(k,j,i), & 1753 1801 shortest_distance( vic_yz, .FALSE., ii ) ) 1754 1802 ! 1755 1803 !-- Reduce search along x by already found distance 1756 dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )1804 dist_dx = CEILING( delta(k,j,i) / dx ) 1757 1805 1758 1806 ENDDO … … 1767 1815 !-- Search for walls within octant (-++) 1768 1816 vic_yz = vicinity(0:rad_k+1,0:rad_j,ii) 1769 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &1817 delta(k,j,i) = MIN( delta(k,j,i), & 1770 1818 shortest_distance( vic_yz, .TRUE., -ii ) ) 1771 1819 ! … … 1775 1823 !-- shortest_distance"). 1776 1824 vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii) 1777 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &1825 delta(k,j,i) = MIN( delta(k,j,i), & 1778 1826 shortest_distance( vic_yz, .TRUE., -ii ) ) 1779 1827 … … 1782 1830 !-- Search for walls within octant (---) 1783 1831 vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii) 1784 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &1832 delta(k,j,i) = MIN( delta(k,j,i), & 1785 1833 shortest_distance( vic_yz, .FALSE., -ii ) ) 1786 1834 ! 1787 1835 !-- Search for walls within octant (-+-) 1788 1836 vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii) 1789 ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &1837 delta(k,j,i) = MIN( delta(k,j,i), & 1790 1838 shortest_distance( vic_yz, .FALSE., -ii ) ) 1791 1839 ! 1792 1840 !-- Reduce search along x by already found distance 1793 dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )1841 dist_dx = CEILING( delta(k,j,i) / dx ) 1794 1842 1795 1843 ENDDO … … 1799 1847 ELSE !Check if (i,j,k) belongs to atmosphere 1800 1848 1801 ml_wall_adjusted(k,j,i) = ml_blackadar(k)1849 delta(k,j,i) = ml_blackadar(k) 1802 1850 1803 1851 ENDIF … … 1818 1866 1819 1867 ! 1820 !-- Set lateral boundary conditions for ml_wall_adjusted 1821 CALL exchange_horiz( ml_wall_adjusted, nbgp ) 1822 1823 !$ACC ENTER DATA COPYIN(ml_wall_adjusted(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) 1868 !-- Set lateral boundary conditions for delta 1869 CALL exchange_horiz( delta, nbgp ) 1870 !$ACC ENTER DATA COPYIN(delta(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) 1871 1872 IF ( les_dai ) THEN 1873 CALL exchange_horiz( distance_to_wall, nbgp ) 1874 !$ACC ENTER DATA COPYIN(distance_to_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) 1875 ENDIF 1876 1877 1824 1878 1825 1879 CONTAINS … … 4032 4086 ! 4033 4087 !-- Calculate the dissipation 4034 IF ( les_dynamic .OR. les_mw) THEN4088 IF( les_default .OR. les_dynamic ) THEN 4035 4089 4036 4090 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 4037 4091 !$ACC PRIVATE(ml, ml_stratification, dvar_dz) & 4038 4092 !$ACC PRESENT(diss, e, var, wall_flags_total_0) & 4039 !$ACC PRESENT(dd2zu, ml_wall_adjusted)4093 !$ACC PRESENT(dd2zu, delta) 4040 4094 DO i = nxl, nxr 4041 4095 DO j = nys, nyn … … 4053 4107 ENDIF 4054 4108 ELSE 4055 ml_stratification(k) = ml_wall_adjusted(k,j,i)4109 ml_stratification(k) = delta(k,j,i) 4056 4110 ENDIF 4057 4111 4058 4112 ENDDO 4059 4113 4114 ! 4115 !-- ATTENTION: Don't merge the following loop with the previous one, because this would prohibit proper vectorization by 4116 !-- the Intel18 compiler. This behaviour might change for future compiler versions. 4060 4117 !$ACC LOOP PRIVATE(k) 4061 4118 !DIR$ IVDEP 4062 4119 DO k = nzb+1, nzt 4063 4120 4064 ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )4065 4066 diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / ml_wall_adjusted(k,j,i) ) &4121 ml = MIN( delta(k,j,i), ml_stratification(k) ) 4122 4123 diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) ) & 4067 4124 * e(k,j,i) * SQRT( e(k,j,i) ) / ml & 4068 4125 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) … … 4073 4130 ENDDO 4074 4131 4132 ELSEIF ( les_dai ) THEN 4133 4134 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 4135 !$ACC PRIVATE(ml, ml_stratification, dvar_dz) & 4136 !$ACC PRESENT(diss, e, var, wall_flags_total_0) & 4137 !$ACC PRESENT(dd2zu, delta) 4138 DO i = nxl, nxr 4139 DO j = nys, nyn 4140 !$ACC LOOP PRIVATE(k) 4141 DO k = nzb+1, nzt 4142 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 4143 IF ( dvar_dz > 0.0_wp ) THEN 4144 ! 4145 !-- The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is 4146 !-- the stratification term. 1E-5 is added as l is zero at the beginning of 4147 !-- the simulation. 4148 IF ( use_single_reference_value ) THEN 4149 ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) & 4150 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 4151 ELSE 4152 ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) & 4153 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 4154 ENDIF 4155 ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) & 4156 + 1.0_wp / ml_stratification(k) ) 4157 ELSE 4158 ml_stratification(k) = delta(k,j,i) 4159 ENDIF 4160 4161 ENDDO 4162 4163 ! 4164 !-- ATTENTION: Don't merge the following loop with the previous one, because this would prohibit proper vectorization by 4165 !-- the Intel18 compiler. This behaviour might change for future compiler versions. 4166 !$ACC LOOP PRIVATE(k) 4167 !DIR$ IVDEP 4168 DO k = nzb+1, nzt 4169 4170 ml = ml_stratification(k) 4171 diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) ) & 4172 * e(k,j,i) * SQRT( e(k,j,i) ) / ml & 4173 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4174 4175 ENDDO 4176 4177 ENDDO 4178 ENDDO 4075 4179 ELSEIF ( rans_tke_l ) THEN 4076 4180 … … 4078 4182 !$ACC PRIVATE(ml_stratification, duv2_dz2, rif, dvar_dz) & 4079 4183 !$ACC PRESENT(diss, e, u, v, var, wall_flags_total_0) & 4080 !$ACC PRESENT(dd2zu, ml_blackadar, ml_wall_adjusted)4184 !$ACC PRESENT(dd2zu, ml_blackadar, delta) 4081 4185 DO i = nxl, nxr 4082 4186 DO j = nys, nyn … … 4084 4188 !-- Calculate Richardson-flux number 4085 4189 IF ( use_single_reference_value ) THEN 4190 4086 4191 !$ACC LOOP PRIVATE(k) 4087 4192 DO k = nzb+1, nzt … … 4096 4201 ENDDO 4097 4202 ELSE 4203 ! 4204 !-- ATTENTION: Don't merge the following loops with the previous one, because this would prohibit proper vectorization 4205 !-- by the Intel18 compiler. This behaviour might change for future compiler versions. 4098 4206 !$ACC LOOP PRIVATE(k) 4099 4207 DO k = nzb+1, nzt … … 4124 4232 4125 4233 diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) & 4126 / MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) )&4234 / MIN( ml_stratification(k), delta(k,j,i) ) & 4127 4235 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4128 4236 … … 4243 4351 !-- Calculate the mixing length and dissipation... 4244 4352 !-- ...in case of LES 4245 IF ( les_d ynamic .OR. les_mw) THEN4353 IF ( les_default .OR. les_dynamic ) THEN 4246 4354 4247 4355 DO k = nzb+1, nzt … … 4256 4364 ENDIF 4257 4365 ELSE 4258 ml_stratification(k) = ml_wall_adjusted(k,j,i)4366 ml_stratification(k) = delta(k,j,i) 4259 4367 ENDIF 4260 4368 ENDDO … … 4262 4370 !DIR$ IVDEP 4263 4371 DO k = nzb+1, nzt 4264 ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )4265 4266 diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / ml_wall_adjusted(k,j,i) )&4372 ml = MIN( delta(k,j,i), ml_stratification(k) ) 4373 4374 diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) ) & 4267 4375 * e(k,j,i) * SQRT( e(k,j,i) ) / ml & 4268 4376 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4269 4377 ENDDO 4378 4379 ELSEIF ( les_dai ) THEN 4380 4381 DO k = nzb+1, nzt 4382 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 4383 IF ( dvar_dz > 0.0_wp ) THEN 4384 ! 4385 !-- The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is 4386 !-- the stratification term. 1E-5 is added as l is zero at the beginning of 4387 !-- the simulation. 4388 IF ( use_single_reference_value ) THEN 4389 ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) & 4390 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 4391 ELSE 4392 ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) & 4393 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 4394 ENDIF 4395 4396 ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) & 4397 + 1.0_wp / ml_stratification(k) ) 4398 4399 ELSE 4400 ml_stratification(k) = delta(k,j,i) 4401 ENDIF 4402 ENDDO 4403 4404 !DIR$ IVDEP 4405 DO k = nzb+1, nzt 4406 ml = ml_stratification(k) 4407 4408 diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) ) & 4409 * e(k,j,i) * SQRT( e(k,j,i) ) / ml & 4410 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4411 ENDDO 4412 4413 4414 4270 4415 4271 4416 ! … … 4310 4455 DO k = nzb+1, nzt 4311 4456 diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) & 4312 / MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) )&4457 / MIN( ml_stratification(k), delta(k,j,i) ) & 4313 4458 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4314 4459 ENDDO … … 4730 4875 !$ tn = omp_get_thread_num() 4731 4876 4732 IF ( les_d ynamic .OR. les_mw) THEN4877 IF ( les_default .OR. les_dynamic ) THEN 4733 4878 !$OMP DO 4734 4879 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 4735 4880 !$ACC PRIVATE(dvar_dz, ml, ml_stratification, ml_local_profile) & 4736 !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, ml_wall_adjusted) &4881 !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, delta) & 4737 4882 !$ACC PRESENT(kh, km, sums_l_l, rmask) 4738 4883 DO i = nxlg, nxrg … … 4745 4890 !-- due to errors when using OpenACC directives. The execution 4746 4891 !-- crashes reliably if a subroutine is called at this point (the 4747 !-- reason gfor this behaviour is unknown, however).4892 !-- reason for this behaviour is unknown, however). 4748 4893 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 4749 4894 IF ( dvar_dz > 0.0_wp ) THEN … … 4756 4901 ENDIF 4757 4902 ELSE 4758 ml_stratification(k) = ml_wall_adjusted(k,j,i)4903 ml_stratification(k) = delta(k,j,i) 4759 4904 ENDIF 4760 4905 … … 4765 4910 DO k = nzb+1, nzt 4766 4911 4767 ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )&4912 ml = MIN( delta(k,j,i), ml_stratification(k) ) & 4768 4913 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4769 4914 ml_local_profile(k) = ml … … 4771 4916 !-- Compute diffusion coefficients for momentum and heat 4772 4917 km(k,j,i) = c_0 * ml * SQRT( e(k,j,i) ) 4773 kh(k,j,i) = ( 1.0_wp + 2.0_wp * ml / ml_wall_adjusted(k,j,i) ) * km(k,j,i)4918 kh(k,j,i) = ( 1.0_wp + 2.0_wp * ml / delta(k,j,i) ) * km(k,j,i) 4774 4919 4775 4920 ENDDO … … 4786 4931 ENDDO 4787 4932 4933 ELSEIF ( les_dai ) THEN 4934 !$OMP DO 4935 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 4936 !$ACC PRIVATE(dvar_dz, ml, ml_stratification, ml_local_profile) & 4937 !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, delta) & 4938 !$ACC PRESENT(kh, km, sums_l_l, rmask) 4939 DO i = nxlg, nxrg 4940 DO j = nysg, nyng 4941 !$ACC LOOP PRIVATE(k) 4942 DO k = nzb+1, nzt 4943 4944 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 4945 IF ( dvar_dz > 0.0_wp ) THEN 4946 ! 4947 !-- The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is 4948 !-- the stratification term. 1E-5 is added as l is zero at the beginning of 4949 !-- the simulation. 4950 IF ( use_single_reference_value ) THEN 4951 ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) & 4952 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 4953 ELSE 4954 ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) & 4955 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 4956 ENDIF 4957 4958 ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) & 4959 + 1.0_wp / ml_stratification(k) ) 4960 ELSE 4961 ml_stratification(k) = delta(k,j,i) 4962 ENDIF 4963 4964 ENDDO 4965 4966 !$ACC LOOP PRIVATE(k) 4967 !DIR$ IVDEP 4968 DO k = nzb+1, nzt 4969 4970 ml_local_profile(k) = ml_stratification(k) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4971 ml = ml_local_profile(k) 4972 ! 4973 !-- Compute diffusion coefficients for momentum and heat 4974 km(k,j,i) = c_0 * ml * SQRT( e(k,j,i) ) 4975 4976 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 4977 IF ( dvar_dz > 0.0_wp ) THEN 4978 kh(k,j,i) = km(k,j,i) 4979 ELSE 4980 kh(k,j,i) = 3.0_wp * km(k,j,i) 4981 ENDIF 4982 4983 ENDDO 4984 ! 4985 !-- Summation for averaged profile (cf. flow_statistics) 4986 !$ACC LOOP PRIVATE(sr, k) 4987 DO sr = 0, statistic_regions 4988 DO k = nzb+1, nzt 4989 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr) 4990 ENDDO 4991 ENDDO 4992 4993 ENDDO 4994 ENDDO 4995 4788 4996 ELSEIF ( rans_tke_l ) THEN 4789 4997 … … 4791 4999 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 4792 5000 !$ACC PRIVATE(dvar_dz, duv2_dz2, ml_stratification, ml_local_profile, rif) & 4793 !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, ml_wall_adjusted, ml_blackadar) &5001 !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, delta, ml_blackadar) & 4794 5002 !$ACC PRESENT(kh, km, sums_l_l, rmask) 4795 5003 DO i = nxlg, nxrg … … 4839 5047 !DIR$ IVDEP 4840 5048 DO k = nzb+1, nzt 4841 ml_local_profile(k) = MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) ) &5049 ml_local_profile(k) = MIN( ml_stratification(k), delta(k,j,i) ) & 4842 5050 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4843 5051 km(k,j,i) = c_0 * ml_local_profile(k) * SQRT( e(k,j,i) ) … … 4993 5201 REAL(wp) :: ldsd !< sum: ld_ij*sd_ij 4994 5202 4995 REAL(wp) :: delta !< grid size5203 REAL(wp) :: delta_max !< maximum of the grid spacings 4996 5204 REAL(wp) :: cst !< c* 4997 5205 REAL(wp) :: cstnust_t !< product c*nu* … … 5119 5327 !-- The model was only tested for an isotropic grid. The following 5120 5328 !-- expression was a recommendation of Stefan Heinz. 5121 delta = MAX( dx, dy, dzw(k) )5329 delta_max = MAX( dx, dy, dzw(k) ) 5122 5330 5123 5331 IF ( lnn <= 0.0_wp ) THEN … … 5125 5333 ELSE 5126 5334 cst = cstnust_t / & 5127 ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp )5335 ( 4.0_wp * delta_max * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp ) 5128 5336 ENDIF 5129 5337 … … 5131 5339 !-- Calculate border according to Mokhtarpoor and Heinz (2017) 5132 5340 cst_max = fac_cmax * SQRT( e(k,j,i) ) / & 5133 ( delta * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )5341 ( delta_max * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp ) 5134 5342 5135 5343 IF ( ABS( cst ) > cst_max ) THEN … … 5137 5345 ENDIF 5138 5346 5139 km(k,j,i) = cst * delta * SQRT( e(k,j,i) ) * flag5347 km(k,j,i) = cst * delta_max * SQRT( e(k,j,i) ) * flag 5140 5348 5141 5349 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
