Ignore:
Timestamp:
Jun 12, 2018 7:03:02 AM (3 years ago)
Author:
Giersch
Message:

New vertical stretching procedure has been introduced

File:
1 edited

Legend:

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

    r3049 r3065  
    2626! -----------------
    2727! $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
    2832! Error messages revised
    2933!
     
    132136!> automatically).
    133137!>
     138!> @todo Replace dz(1) appropriatly to account for grid stretching
    134139!> @todo Revise code according to PALM Coding Standard
    135140!> @todo Implement ADM and ALM turbine models
     
    740745       delta_r_factor = segment_width
    741746       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))
    743748
    744749       DO inot = 1, nturbines
     
    937942
    938943   
     944       USE control_parameters,                                                 &
     945           ONLY:  dz_stretch_level_start
     946   
    939947       IMPLICIT NONE
    940948
     
    9921000                 pi**( 1.0_wp / 6.0_wp ) * eps_kernel
    9931001!
     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!
    9941015!--    Square of eps_min:
    9951016       eps_min2 = eps_min**2
     
    10221043          i_hub(inot) = INT(   rcx(inot)                 / dx )
    10231044          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) )
    10251046
    10261047!
     
    10311052          i_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dx )
    10321053          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) )
    10341055       
    10351056       ENDDO
     
    11211142!--       surface. All points between z=0 and z=dz/s would already be contained
    11221143!--       in grid box 1.
    1123           index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz ) + 1
    1124           index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz ) + 1
     1144          index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz(1) ) + 1
     1145          index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz(1) ) + 1
    11251146
    11261147!
     
    11481169                            IF ( j == tower_s )  THEN
    11491170                               tow_cd_surf(k,j,i) = ( rcz(inot) -              &
    1150                                     ( k_hub(inot) * dz - 0.5_wp * dz ) )  *    & ! extension in z-direction
     1171                                    ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction
    11511172                                  ( ( tower_s + 1.0_wp + 0.5_wp ) * dy    -    &
    11521173                                    ( rcy(inot) - 0.5_wp * dtow(inot) ) ) *    & ! extension in y-direction
     
    11541175                            ELSEIF ( j == tower_n )  THEN
    11551176                               tow_cd_surf(k,j,i) = ( rcz(inot)            -   &
    1156                                     ( k_hub(inot) * dz - 0.5_wp * dz ) )  *    & ! extension in z-direction
     1177                                    ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction
    11571178                                  ( ( rcy(inot) + 0.5_wp * dtow(inot) )   -    &
    11581179                                    ( tower_n + 0.5_wp ) * dy )           *    & ! extension in y-direction
     
    11631184                            ELSE
    11641185                               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) ) )*&
    11661187                                    dy * turb_cd_tower(inot)
    11671188                            ENDIF
     
    11691190!--                      tower lies completely within one grid box:
    11701191                         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) ) ) *            &
    11731194                                       dtow(inot) * turb_cd_tower(inot)
    11741195                         ENDIF
     
    11821203!--                         leftmost and rightmost grid box:
    11831204                            IF ( j == tower_s )  THEN                         
    1184                                tow_cd_surf(k,j,i) = dz * (                     &
     1205                               tow_cd_surf(k,j,i) = dz(1) * (                  &
    11851206                                      ( tower_s + 1 + 0.5_wp ) * dy         -  &
    11861207                                      ( rcy(inot) - 0.5_wp * dtow(inot) )      &
    11871208                                                        ) * turb_cd_tower(inot)
    11881209                            ELSEIF ( j == tower_n )  THEN
    1189                                tow_cd_surf(k,j,i) = dz * (                     &
     1210                               tow_cd_surf(k,j,i) = dz(1) * (                  &
    11901211                                      ( rcy(inot) + 0.5_wp * dtow(inot) )   -  &
    11911212                                      ( tower_n + 0.5_wp ) * dy                &
     
    11951216!--                         (where tow_cd_surf = grid box area):
    11961217                            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)
    11981220                            ENDIF
    11991221!
    12001222!--                         tower lies completely within one grid box:
    12011223                         ELSE
    1202                             tow_cd_surf(k,j,i) = dz * dtow(inot) *             &
     1224                            tow_cd_surf(k,j,i) = dz(1) * dtow(inot) *          &
    12031225                                                turb_cd_tower(inot)
    12041226                         ENDIF ! end if larger than grid box
     
    13021324       ENDDO   ! end of loop over turbines
    13031325
    1304        tow_cd_surf   = tow_cd_surf   / ( dx * dy * dz )      ! Normalize tower drag
    1305        nac_cd_surf = nac_cd_surf / ( dx * dy * dz )      ! Normalize nacelle drag
     1326       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
    13061328
    13071329       CALL wtm_read_blade_tables
     
    17461768                   ii =   rbx(ring,rseg) * ddx
    17471769                   jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy
    1748                    kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz
     1770                   kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1)
    17491771!
    17501772!--                Interpolate only if all required information is available on
     
    17761798
    17771799                         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) *     &
    17791801                                     ( u_int_u - u_int_l )
    17801802
     
    17911813                   ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx
    17921814                   jj =   rby(ring,rseg)                 * ddy
    1793                    kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz
     1815                   kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1)
    17941816!
    17951817!--                Interpolate only if all required information is available on
     
    18211843
    18221844                         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) *     &
    18241846                                     ( v_int_u - v_int_l )
    18251847
     
    18361858                   ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx
    18371859                   jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy
    1838                    kk =   rbz(ring,rseg)                 / dz
     1860                   kk =   rbz(ring,rseg)                 / dz(1)
    18391861!
    18401862!--                Interpolate only if all required information is available on
     
    18661888
    18671889                         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) *     &
    18691891                                     ( w_int_u - w_int_l )
    18701892                      ELSE
     
    22732295                            dist_u_3d = ( i * dx               - rbx(ring,rseg) )**2 + &
    22742296                                        ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + &
    2275                                         ( k * dz - 0.5_wp * dz - rbz(ring,rseg) )**2
     2297                                        ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2
    22762298                            dist_v_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + &
    22772299                                        ( j * dy               - rby(ring,rseg) )**2 + &
    2278                                         ( k * dz - 0.5_wp * dz - rbz(ring,rseg) )**2
     2300                                        ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2
    22792301                            dist_w_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + &
    22802302                                        ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + &
    2281                                         ( k * dz               - rbz(ring,rseg) )**2
     2303                                        ( k * dz(1)               - rbz(ring,rseg) )**2
    22822304
    22832305!
Note: See TracChangeset for help on using the changeset viewer.