Changeset 2232 for palm/trunk/SOURCE/wind_turbine_model_mod.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90
r2153 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Adjustments to new topography concept 23 23 ! 24 24 ! Former revisions: … … 86 86 USE indices, & 87 87 ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, & 88 nzb, nz b_u_inner, nzb_v_inner, nzb_w_inner, nzt88 nzb, nzt, wall_flags_0 89 89 90 90 USE kinds … … 1358 1358 INTEGER(iwp) :: ii, jj, kk !< 1359 1359 1360 REAL(wp) :: flag !< flag to mask topography grid points 1360 1361 REAL(wp) :: sin_rot, cos_rot !< 1361 1362 REAL(wp) :: sin_yaw, cos_yaw !< … … 1971 1972 DO j = MAX( nys, j_hub(inot) - j_smear(inot) ), & 1972 1973 MIN( nyn, j_hub(inot) + j_smear(inot) ) 1973 DO k = MAX( nzb_u_inner(j,i)+1, k_hub(inot) - k_smear(inot) ), & 1974 k_hub(inot) + k_smear(inot) 1974 ! DO k = MAX( nzb_u_inner(j,i)+1, k_hub(inot) - k_smear(inot) ), & 1975 ! k_hub(inot) + k_smear(inot) 1976 DO k = nzb+1, k_hub(inot) + k_smear(inot) 1975 1977 DO ring = 1, nrings(inot) 1976 1978 DO rseg = 1, nsegs(ring,inot) 1979 ! 1980 !-- Predetermine flag to mask topography 1981 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1 ) ) 1982 1977 1983 ! 1978 1984 !-- Determine the square of the distance between the … … 1995 2001 !-- squared distances <= eps_min2: 1996 2002 IF ( dist_u_3d <= eps_min2 ) THEN 1997 thrust(k,j,i) = thrust(k,j,i) + & 1998 thrust_ring(ring,rseg) * & 1999 ( ( pol_a * dist_u_3d - pol_b ) * & 2000 dist_u_3d + 1.0_wp ) * eps_factor 2003 thrust(k,j,i) = thrust(k,j,i) + & 2004 thrust_ring(ring,rseg) * & 2005 ( ( pol_a * dist_u_3d - pol_b ) * & 2006 dist_u_3d + 1.0_wp ) * eps_factor *& 2007 flag 2001 2008 ENDIF 2002 2009 IF ( dist_v_3d <= eps_min2 ) THEN 2003 torque_y(k,j,i) = torque_y(k,j,i) + & 2004 torque_ring_y(ring,rseg) * & 2005 ( ( pol_a * dist_v_3d - pol_b ) *& 2006 dist_v_3d + 1.0_wp ) * eps_factor 2010 torque_y(k,j,i) = torque_y(k,j,i) + & 2011 torque_ring_y(ring,rseg) * & 2012 ( ( pol_a * dist_v_3d - pol_b ) * & 2013 dist_v_3d + 1.0_wp ) * eps_factor * & 2014 flag 2007 2015 ENDIF 2008 2016 IF ( dist_w_3d <= eps_min2 ) THEN 2009 torque_z(k,j,i) = torque_z(k,j,i) + & 2010 torque_ring_z(ring,rseg) * & 2011 ( ( pol_a * dist_w_3d - pol_b ) *& 2012 dist_w_3d + 1.0_wp ) * eps_factor 2017 torque_z(k,j,i) = torque_z(k,j,i) + & 2018 torque_ring_z(ring,rseg) * & 2019 ( ( pol_a * dist_w_3d - pol_b ) * & 2020 dist_w_3d + 1.0_wp ) * eps_factor * & 2021 flag 2013 2022 ENDIF 2014 2023 … … 2018 2027 ! 2019 2028 !-- Rotation of force components: 2020 rot_tend_x(k,j,i) = rot_tend_x(k,j,i) + 2029 rot_tend_x(k,j,i) = rot_tend_x(k,j,i) + ( & 2021 2030 thrust(k,j,i)*rotx(inot,1) + & 2022 2031 torque_y(k,j,i)*roty(inot,1) + & 2023 torque_z(k,j,i)*rotz(inot,1) 2032 torque_z(k,j,i)*rotz(inot,1) & 2033 ) * flag 2024 2034 2025 rot_tend_y(k,j,i) = rot_tend_y(k,j,i) + 2035 rot_tend_y(k,j,i) = rot_tend_y(k,j,i) + ( & 2026 2036 thrust(k,j,i)*rotx(inot,2) + & 2027 2037 torque_y(k,j,i)*roty(inot,2) + & 2028 torque_z(k,j,i)*rotz(inot,2) 2038 torque_z(k,j,i)*rotz(inot,2) & 2039 ) * flag 2029 2040 2030 rot_tend_z(k,j,i) = rot_tend_z(k,j,i) + 2041 rot_tend_z(k,j,i) = rot_tend_z(k,j,i) + ( & 2031 2042 thrust(k,j,i)*rotx(inot,3) + & 2032 2043 torque_y(k,j,i)*roty(inot,3) + & 2033 torque_z(k,j,i)*rotz(inot,3) 2044 torque_z(k,j,i)*rotz(inot,3) & 2045 ) * flag 2034 2046 2035 2047 ENDDO ! End of loop over k … … 2424 2436 DO i = nxlg, nxrg 2425 2437 DO j = nysg, nyng 2426 DO k = nzb _u_inner(j,i)+1, k_hub(1) + k_smear(1)2438 DO k = nzb+1, k_hub(1) + k_smear(1) 2427 2439 ! 2428 2440 !-- Calculate the thrust generated by the nacelle and the tower … … 2432 2444 SIGN( u(k,j,i)**2 , u(k,j,i) ) 2433 2445 2434 tend(k,j,i) = tend(k,j,i) - rot_tend_x(k,j,i) & 2435 - tend_nac_x - tend_tow_x 2446 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) & 2447 - tend_nac_x - tend_tow_x ) & 2448 * MERGE( 1.0_wp, 0.0_wp, & 2449 BTEST( wall_flags_0(k,j,i), 1 ) ) 2436 2450 ENDDO 2437 2451 ENDDO … … 2445 2459 DO i = nxlg, nxrg 2446 2460 DO j = nysg, nyng 2447 DO k = nzb _v_inner(j,i)+1, k_hub(1) + k_smear(1)2461 DO k = nzb+1, k_hub(1) + k_smear(1) 2448 2462 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * & 2449 2463 SIGN( v(k,j,i)**2 , v(k,j,i) ) 2450 2464 tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * & 2451 2465 SIGN( v(k,j,i)**2 , v(k,j,i) ) 2452 tend(k,j,i) = tend(k,j,i) - rot_tend_y(k,j,i) & 2453 - tend_nac_y - tend_tow_y 2466 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) & 2467 - tend_nac_y - tend_tow_y ) & 2468 * MERGE( 1.0_wp, 0.0_wp, & 2469 BTEST( wall_flags_0(k,j,i), 2 ) ) 2454 2470 ENDDO 2455 2471 ENDDO … … 2463 2479 DO i = nxlg, nxrg 2464 2480 DO j = nysg, nyng 2465 DO k = nzb_w_inner(j,i)+1, k_hub(1) + k_smear(1) 2466 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) 2481 DO k = nzb+1, k_hub(1) + k_smear(1) 2482 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) & 2483 * MERGE( 1.0_wp, 0.0_wp, & 2484 BTEST( wall_flags_0(k,j,i), 3 ) ) 2467 2485 ENDDO 2468 2486 ENDDO … … 2506 2524 IF ( simulated_time >= time_turbine_on ) THEN 2507 2525 2508 DO k = nzb _u_inner(j,i)+1, k_hub(1) + k_smear(1)2526 DO k = nzb+1, k_hub(1) + k_smear(1) 2509 2527 ! 2510 2528 !-- Calculate the thrust generated by the nacelle and the tower … … 2513 2531 tend_tow_x = 0.5_wp * tow_cd_surf(k,j,i) * & 2514 2532 SIGN( u(k,j,i)**2 , u(k,j,i) ) 2515 tend(k,j,i) = tend(k,j,i) - rot_tend_x(k,j,i) & 2516 - tend_nac_x - tend_tow_x 2533 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) & 2534 - tend_nac_x - tend_tow_x ) & 2535 * MERGE( 1.0_wp, 0.0_wp, & 2536 BTEST( wall_flags_0(k,j,i), 1 ) ) 2517 2537 ENDDO 2518 2538 ENDIF … … 2522 2542 !-- Apply the y-component of the force to the v-component of the flow: 2523 2543 IF ( simulated_time >= time_turbine_on ) THEN 2524 DO k = nzb _v_inner(j,i)+1, k_hub(1) + k_smear(1)2544 DO k = nzb+1, k_hub(1) + k_smear(1) 2525 2545 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * & 2526 2546 SIGN( v(k,j,i)**2 , v(k,j,i) ) 2527 2547 tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * & 2528 2548 SIGN( v(k,j,i)**2 , v(k,j,i) ) 2529 tend(k,j,i) = tend(k,j,i) - rot_tend_y(k,j,i) & 2530 - tend_nac_y - tend_tow_y 2549 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) & 2550 - tend_nac_y - tend_tow_y ) & 2551 * MERGE( 1.0_wp, 0.0_wp, & 2552 BTEST( wall_flags_0(k,j,i), 2 ) ) 2531 2553 ENDDO 2532 2554 ENDIF … … 2536 2558 !-- Apply the z-component of the force to the w-component of the flow: 2537 2559 IF ( simulated_time >= time_turbine_on ) THEN 2538 DO k = nzb_w_inner(j,i)+1, k_hub(1) + k_smear(1) 2539 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) 2560 DO k = nzb+1, k_hub(1) + k_smear(1) 2561 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) & 2562 * MERGE( 1.0_wp, 0.0_wp, & 2563 BTEST( wall_flags_0(k,j,i), 3 ) ) 2540 2564 ENDDO 2541 2565 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
