Changeset 3065 for palm/trunk/SOURCE/wind_turbine_model_mod.f90
- Timestamp:
- Jun 12, 2018 7:03:02 AM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90
r3049 r3065 26 26 ! ----------------- 27 27 ! $Id$ 28 ! dz was replaced by dz(1), error message concerning grid stretching has been 29 ! introduced 30 ! 31 ! 3049 2018-05-29 13:52:36Z Giersch 28 32 ! Error messages revised 29 33 ! … … 132 136 !> automatically). 133 137 !> 138 !> @todo Replace dz(1) appropriatly to account for grid stretching 134 139 !> @todo Revise code according to PALM Coding Standard 135 140 !> @todo Implement ADM and ALM turbine models … … 740 745 delta_r_factor = segment_width 741 746 delta_t_factor = segment_length 742 delta_r_init = delta_r_factor * MIN( dx, dy, dz )747 delta_r_init = delta_r_factor * MIN( dx, dy, dz(1)) 743 748 744 749 DO inot = 1, nturbines … … 937 942 938 943 944 USE control_parameters, & 945 ONLY: dz_stretch_level_start 946 939 947 IMPLICIT NONE 940 948 … … 992 1000 pi**( 1.0_wp / 6.0_wp ) * eps_kernel 993 1001 ! 1002 !-- Stretching (non-uniform grid spacing) is not considered in the wind 1003 !-- turbine model. Therefore, vertical stretching has to be applied above 1004 !-- the area where the wtm is active. ABS (...) is required because the 1005 !-- default value of dz_stretch_level_start is -9999999.9_wp (negative). 1006 IF ( ABS( dz_stretch_level_start(1) ) <= MAXVAL(rcz) + MAXVAL(rr) + & 1007 eps_min) THEN 1008 WRITE( message_string, * ) 'The lowest level where vertical ', & 1009 'stretching is applied &have to be ', & 1010 'greater than ', MAXVAL(rcz) + & 1011 MAXVAL(rr) + eps_min 1012 CALL message( 'init_grid', 'PA0484', 1, 2, 0, 6, 0 ) 1013 ENDIF 1014 ! 994 1015 !-- Square of eps_min: 995 1016 eps_min2 = eps_min**2 … … 1022 1043 i_hub(inot) = INT( rcx(inot) / dx ) 1023 1044 j_hub(inot) = INT( ( rcy(inot) + 0.5_wp * dy ) / dy ) 1024 k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz ) / dz)1045 k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz(1) ) / dz(1) ) 1025 1046 1026 1047 ! … … 1031 1052 i_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dx ) 1032 1053 j_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dy ) 1033 k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz )1054 k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz(1) ) 1034 1055 1035 1056 ENDDO … … 1121 1142 !-- surface. All points between z=0 and z=dz/s would already be contained 1122 1143 !-- in grid box 1. 1123 index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz ) + 11124 index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz ) + 11144 index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz(1) ) + 1 1145 index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz(1) ) + 1 1125 1146 1126 1147 ! … … 1148 1169 IF ( j == tower_s ) THEN 1149 1170 tow_cd_surf(k,j,i) = ( rcz(inot) - & 1150 ( k_hub(inot) * dz - 0.5_wp * dz ) ) *& ! extension in z-direction1171 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction 1151 1172 ( ( tower_s + 1.0_wp + 0.5_wp ) * dy - & 1152 1173 ( rcy(inot) - 0.5_wp * dtow(inot) ) ) * & ! extension in y-direction … … 1154 1175 ELSEIF ( j == tower_n ) THEN 1155 1176 tow_cd_surf(k,j,i) = ( rcz(inot) - & 1156 ( k_hub(inot) * dz - 0.5_wp * dz ) ) *& ! extension in z-direction1177 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction 1157 1178 ( ( rcy(inot) + 0.5_wp * dtow(inot) ) - & 1158 1179 ( tower_n + 0.5_wp ) * dy ) * & ! extension in y-direction … … 1163 1184 ELSE 1164 1185 tow_cd_surf(k,j,i) = ( rcz(inot) - & 1165 ( k_hub(inot) * dz - 0.5_wp * dz ) ) *&1186 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& 1166 1187 dy * turb_cd_tower(inot) 1167 1188 ENDIF … … 1169 1190 !-- tower lies completely within one grid box: 1170 1191 ELSE 1171 tow_cd_surf(k,j,i) = ( rcz(inot) -&1172 ( k_hub(inot) * dz - 0.5_wp * dz ) ) *&1192 tow_cd_surf(k,j,i) = ( rcz(inot) - ( k_hub(inot) * & 1193 dz(1) - 0.5_wp * dz(1) ) ) * & 1173 1194 dtow(inot) * turb_cd_tower(inot) 1174 1195 ENDIF … … 1182 1203 !-- leftmost and rightmost grid box: 1183 1204 IF ( j == tower_s ) THEN 1184 tow_cd_surf(k,j,i) = dz * (&1205 tow_cd_surf(k,j,i) = dz(1) * ( & 1185 1206 ( tower_s + 1 + 0.5_wp ) * dy - & 1186 1207 ( rcy(inot) - 0.5_wp * dtow(inot) ) & 1187 1208 ) * turb_cd_tower(inot) 1188 1209 ELSEIF ( j == tower_n ) THEN 1189 tow_cd_surf(k,j,i) = dz * (&1210 tow_cd_surf(k,j,i) = dz(1) * ( & 1190 1211 ( rcy(inot) + 0.5_wp * dtow(inot) ) - & 1191 1212 ( tower_n + 0.5_wp ) * dy & … … 1195 1216 !-- (where tow_cd_surf = grid box area): 1196 1217 ELSE 1197 tow_cd_surf(k,j,i) = dz * dy * turb_cd_tower(inot) 1218 tow_cd_surf(k,j,i) = dz(1) * dy * & 1219 turb_cd_tower(inot) 1198 1220 ENDIF 1199 1221 ! 1200 1222 !-- tower lies completely within one grid box: 1201 1223 ELSE 1202 tow_cd_surf(k,j,i) = dz * dtow(inot) *&1224 tow_cd_surf(k,j,i) = dz(1) * dtow(inot) * & 1203 1225 turb_cd_tower(inot) 1204 1226 ENDIF ! end if larger than grid box … … 1302 1324 ENDDO ! end of loop over turbines 1303 1325 1304 tow_cd_surf = tow_cd_surf / ( dx * dy * dz )! Normalize tower drag1305 nac_cd_surf = nac_cd_surf / ( dx * dy * dz ) ! Normalize nacelle drag1326 tow_cd_surf = tow_cd_surf / ( dx * dy * dz(1) ) ! Normalize tower drag 1327 nac_cd_surf = nac_cd_surf / ( dx * dy * dz(1) ) ! Normalize nacelle drag 1306 1328 1307 1329 CALL wtm_read_blade_tables … … 1746 1768 ii = rbx(ring,rseg) * ddx 1747 1769 jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy 1748 kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz1770 kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 1749 1771 ! 1750 1772 !-- Interpolate only if all required information is available on … … 1776 1798 1777 1799 u_int_1_l(inot,ring,rseg) = u_int_l + & 1778 ( rbz(ring,rseg) - zu(kk) ) / dz *&1800 ( rbz(ring,rseg) - zu(kk) ) / dz(1) * & 1779 1801 ( u_int_u - u_int_l ) 1780 1802 … … 1791 1813 ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx 1792 1814 jj = rby(ring,rseg) * ddy 1793 kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz1815 kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 1794 1816 ! 1795 1817 !-- Interpolate only if all required information is available on … … 1821 1843 1822 1844 v_int_1_l(inot,ring,rseg) = v_int_l + & 1823 ( rbz(ring,rseg) - zu(kk) ) / dz *&1845 ( rbz(ring,rseg) - zu(kk) ) / dz(1) * & 1824 1846 ( v_int_u - v_int_l ) 1825 1847 … … 1836 1858 ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx 1837 1859 jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy 1838 kk = rbz(ring,rseg) / dz 1860 kk = rbz(ring,rseg) / dz(1) 1839 1861 ! 1840 1862 !-- Interpolate only if all required information is available on … … 1866 1888 1867 1889 w_int_1_l(inot,ring,rseg) = w_int_l + & 1868 ( rbz(ring,rseg) - zw(kk) ) / dz *&1890 ( rbz(ring,rseg) - zw(kk) ) / dz(1) * & 1869 1891 ( w_int_u - w_int_l ) 1870 1892 ELSE … … 2273 2295 dist_u_3d = ( i * dx - rbx(ring,rseg) )**2 + & 2274 2296 ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & 2275 ( k * dz - 0.5_wp * dz- rbz(ring,rseg) )**22297 ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2 2276 2298 dist_v_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & 2277 2299 ( j * dy - rby(ring,rseg) )**2 + & 2278 ( k * dz - 0.5_wp * dz- rbz(ring,rseg) )**22300 ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2 2279 2301 dist_w_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & 2280 2302 ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & 2281 ( k * dz - rbz(ring,rseg) )**22303 ( k * dz(1) - rbz(ring,rseg) )**2 2282 2304 2283 2305 !
Note: See TracChangeset
for help on using the changeset viewer.
