Ignore:
Timestamp:
Oct 18, 2017 8:03:45 AM (6 years ago)
Author:
Giersch
Message:

Changed programming style and bugfixes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/wind_turbine_model_mod.f90

    r2410 r2553  
    2626! -----------------
    2727! $Id$
     28! Bugfix of vertical loop in wtm_tendencies to account for different turbine
     29! heights, bugfix of the interpolation of the u-component concerning the
     30! vertical index and further small adjustments of the programming style
     31!
     32! 2410 2017-09-06 08:16:30Z Giersch
    2833! Revise error message PA0462
    2934!
     
    9398
    9499    USE arrays_3d,                                                             &
    95         ONLY:  dd2zu, tend, u, v, w, zu, zw
    96 
    97     USE constants
     100        ONLY:  tend, u, v, w, zu, zw
     101
     102    USE constants,                                                             &
     103        ONLY:  pi
    98104
    99105    USE control_parameters,                                                    &
     
    140146
    141147    REAL(wp), DIMENSION(1:100) ::  dtow             = 0.0_wp  !< tower diameter [m]
    142     REAL(wp), DIMENSION(1:100) ::  omega_rot        = 0.0_wp  !< inital or constant rotor speed [rad/s]
     148    REAL(wp), DIMENSION(1:100) ::  omega_rot        = 0.9_wp  !< inital or constant rotor speed [rad/s]
    143149    REAL(wp), DIMENSION(1:100) ::  phi_yaw          = 0.0_wp  !< yaw angle [degree] ( clockwise, 0 = facing west )
    144150    REAL(wp), DIMENSION(1:100) ::  pitch_add        = 0.0_wp  !< constant pitch angle
     
    240246     
    241247    REAL(wp) ::  cur_r                       !<
    242     REAL(wp) ::  delta_t                     !<  tangential segment length
    243248    REAL(wp) ::  phi_rotor                   !<
    244249    REAL(wp) ::  pre_factor                  !< 
     
    383388    INTERFACE wtm_forces
    384389       MODULE PROCEDURE wtm_forces
     390    END INTERFACE wtm_forces
     391
     392    INTERFACE wtm_yawcontrol
    385393       MODULE PROCEDURE wtm_yawcontrol
    386     END INTERFACE wtm_forces
     394    END INTERFACE wtm_yawcontrol
    387395   
    388396    INTERFACE wtm_rotate_rotor
     
    402410   
    403411    PUBLIC wtm_check_parameters, wtm_forces, wtm_init, wtm_init_arrays,        &
    404            wtm_parin, wtm_tendencies, wtm_tendencies_ij, wind_turbine
     412           wtm_parin, wtm_tendencies, wind_turbine
    405413
    406414 CONTAINS
     
    524532       delta_t_factor = segment_length
    525533       delta_r_init   = delta_r_factor * MIN( dx, dy, dz)
    526        delta_t        = delta_t_factor * MIN( dx, dy, dz)
    527534
    528535       DO inot = 1, nturbines
     
    14841491                   ii =   rbx(ring,rseg) * ddx
    14851492                   jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy
    1486                    kk = ( rbz(ring,rseg) - 0.5_wp * dz ) / dz
     1493                   kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz
    14871494!
    14881495!--                Interpolate only if all required information is available on
     
    17101717!
    17111718!--                In the following the 3D-velocity field is projected its
    1712 !--                components perpedicular and parallel to the rotor area
     1719!--                components perpendicular and parallel to the rotor area
    17131720!--                The calculation of forces will be done in the rotor-
    17141721!--                coordinates y' and z.
     
    22302237    SUBROUTINE wtm_yawcontrol( inot )
    22312238   
    2232        USE constants
    22332239       USE kinds
    22342240               
     
    24612467             DO  i = nxlg, nxrg
    24622468                DO  j = nysg, nyng
    2463                    DO  k = nzb+1, k_hub(1) + k_smear(1)
     2469                   DO  k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear)
    24642470!
    24652471!--                   Calculate the thrust generated by the nacelle and the tower
    24662472                      tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) *               &
    2467                                          SIGN( u(k,j,i)**2 , u(k,j,i) )     
     2473                                      SIGN( u(k,j,i)**2 , u(k,j,i) )     
    24682474                      tend_tow_x   = 0.5_wp * tow_cd_surf(k,j,i) *             &
    24692475                                         SIGN( u(k,j,i)**2 , u(k,j,i) )
    2470                                                    
     2476                                               
    24712477                      tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i)        &
    24722478                                  - tend_nac_x - tend_tow_x )                  &
     
    24842490             DO  i = nxlg, nxrg
    24852491                DO  j = nysg, nyng
    2486                    DO  k = nzb+1, k_hub(1) + k_smear(1)
     2492                   DO  k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear)
    24872493                      tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) *               &
    24882494                                         SIGN( v(k,j,i)**2 , v(k,j,i) )     
     
    25042510             DO  i = nxlg, nxrg
    25052511                DO  j = nysg, nyng
    2506                    DO  k = nzb+1,  k_hub(1) + k_smear(1)
     2512                   DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
    25072513                      tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i)            &
    25082514                                      * MERGE( 1.0_wp, 0.0_wp,                 &
     
    25482554!--       Apply the x-component of the force to the u-component of the flow:
    25492555          IF ( simulated_time >= time_turbine_on )  THEN
    2550 
    2551              DO  k = nzb+1,  k_hub(1) + k_smear(1)
     2556             DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
    25522557!
    25532558!--             Calculate the thrust generated by the nacelle and the tower
     
    25592564                            - tend_nac_x - tend_tow_x )                        &
    25602565                                      * MERGE( 1.0_wp, 0.0_wp,                 &
    2561                                                BTEST( wall_flags_0(k,j,i), 1 ) )
     2566                                            BTEST( wall_flags_0(k,j,i), 1 ) )
    25622567             ENDDO
    25632568          ENDIF
     
    25672572!--       Apply the y-component of the force to the v-component of the flow:
    25682573          IF ( simulated_time >= time_turbine_on )  THEN
    2569              DO  k = nzb+1,  k_hub(1) + k_smear(1)
     2574             DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
    25702575                tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) *                     &
    25712576                                   SIGN( v(k,j,i)**2 , v(k,j,i) )     
     
    25762581                                      * MERGE( 1.0_wp, 0.0_wp,                 &
    25772582                                               BTEST( wall_flags_0(k,j,i), 2 ) )
    2578              ENDDO
     2583                ENDDO
    25792584          ENDIF
    25802585
     
    25832588!--       Apply the z-component of the force to the w-component of the flow:
    25842589          IF ( simulated_time >= time_turbine_on )  THEN
    2585              DO  k = nzb+1,  k_hub(1) + k_smear(1)
     2590             DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
    25862591                tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i)                  &
    25872592                                      * MERGE( 1.0_wp, 0.0_wp,                 &
Note: See TracChangeset for help on using the changeset viewer.