| 1 | !WRF:MODEL_LAYER:PHYSICS |
|---|
| 2 | |
|---|
| 3 | MODULE module_wind_fitch |
|---|
| 4 | |
|---|
| 5 | ! Calculates drag at model levels intersected by turbine blades |
|---|
| 6 | ! using wind speed dependent thrust coeff obtained from manufacturer. |
|---|
| 7 | ! Adds source of TKE from KE not extracted by turbine for power. |
|---|
| 8 | ! Output: |
|---|
| 9 | ! du, dv: horizontal velocity tendencies |
|---|
| 10 | ! qke: TKE |
|---|
| 11 | ! Input: |
|---|
| 12 | ! dz = dz between full levels (m) |
|---|
| 13 | ! phb = geopotential height |
|---|
| 14 | ! theight = hub height in m |
|---|
| 15 | ! tdiameter = turbine diameter in m |
|---|
| 16 | ! stdthrcoef = standing thrust coeff. |
|---|
| 17 | ! tpower = turbine power in MW |
|---|
| 18 | ! cospeed = cut-out speed in m/s |
|---|
| 19 | ! cispeed = cut-in speed in m/s |
|---|
| 20 | ! ewfx = x-extent of rectangular wind farm in grid cells |
|---|
| 21 | ! ewfy = y-extent of rectangular wind farm in grid cells |
|---|
| 22 | ! pwfx = x-coord of grid cell in SW corner of farm |
|---|
| 23 | ! pwfy = y-coord of grid cell in SW corner of farm |
|---|
| 24 | ! turbpercell = no. of turbines per grid cell |
|---|
| 25 | |
|---|
| 26 | USE module_wind_generic |
|---|
| 27 | USE module_driver_constants, ONLY : max_domains |
|---|
| 28 | USE module_model_constants, ONLY : g |
|---|
| 29 | |
|---|
| 30 | IMPLICIT NONE |
|---|
| 31 | |
|---|
| 32 | LOGICAL, DIMENSION(max_domains) :: inited |
|---|
| 33 | |
|---|
| 34 | PUBLIC turbine_drag |
|---|
| 35 | PRIVATE dragcof, turbine_area, inited |
|---|
| 36 | |
|---|
| 37 | CONTAINS |
|---|
| 38 | |
|---|
| 39 | SUBROUTINE turbine_drag( & |
|---|
| 40 | & id & |
|---|
| 41 | &,phb,u,v,xlat_u,xlong_u & |
|---|
| 42 | &,xlat_v,xlong_v & |
|---|
| 43 | &,dx,dz,dt,qke,du,dv & |
|---|
| 44 | &,ids,ide,jds,jde,kds,kde & |
|---|
| 45 | &,ims,ime,jms,jme,kms,kme & |
|---|
| 46 | &,its,ite,jts,jte,kts,kte & |
|---|
| 47 | &) |
|---|
| 48 | |
|---|
| 49 | INTEGER, INTENT(IN) :: id ! grid id |
|---|
| 50 | INTEGER, INTENT(IN) :: its,ite,jts,jte,kts,kte |
|---|
| 51 | INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme |
|---|
| 52 | INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde |
|---|
| 53 | REAL, INTENT(IN) :: dx,dt |
|---|
| 54 | REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: dz,u,v,phb |
|---|
| 55 | REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat_u, xlong_u |
|---|
| 56 | REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat_v, xlong_v |
|---|
| 57 | REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: du,dv,qke |
|---|
| 58 | ! Local |
|---|
| 59 | TYPE(windturbine_specs), POINTER :: p |
|---|
| 60 | INTEGER turbgridid |
|---|
| 61 | REAL hubheight,diameter,power,cutinspeed,cutoutspeed,stdthrcoef,turbpercell |
|---|
| 62 | INTEGER ewfx,ewfy,pwfx,pwfy |
|---|
| 63 | REAL blade_l_point,blade_u_point,zheightl,zheightu,z1,z2,tarea |
|---|
| 64 | REAL speed,tkecof,powcof,thrcof,wfdensity |
|---|
| 65 | INTEGER itf,jtf,ktf |
|---|
| 66 | INTEGER i,j,k,swfindx,ewfindx,swfindy,ewfindy,n,n1,n2,iturbine |
|---|
| 67 | INTEGER k_turbine_bot, k_turbine_top |
|---|
| 68 | |
|---|
| 69 | LOGICAL :: kfound |
|---|
| 70 | INTEGER :: allzero |
|---|
| 71 | |
|---|
| 72 | itf=MIN0(ite,ide-1) |
|---|
| 73 | jtf=MIN0(jte,jde-1) |
|---|
| 74 | ktf=MIN0(kte,kde-1) |
|---|
| 75 | |
|---|
| 76 | CALL nl_get_td_turbpercell(1,turbpercell) |
|---|
| 77 | CALL nl_get_td_turbgridid(1,turbgridid) |
|---|
| 78 | IF ( .NOT. inited(id) ) THEN |
|---|
| 79 | IF ( windspec .EQ. WIND_TURBINES_FROMLIST ) THEN |
|---|
| 80 | ! first check and see if xlat and xlong are all zero, if so, then use i,j directly |
|---|
| 81 | ! (just check the u variables) |
|---|
| 82 | allzero=1 |
|---|
| 83 | DO j=jts,jtf |
|---|
| 84 | DO i=its,itf |
|---|
| 85 | IF (xlat_u(i,j).NE.0. .OR. xlong_u(i,j).NE.0)allzero=0 |
|---|
| 86 | ENDDO |
|---|
| 87 | ENDDO |
|---|
| 88 | CALL wrf_dm_bcast_integer(allzero,1) |
|---|
| 89 | IF ( allzero .NE. 1 ) THEN |
|---|
| 90 | ! if there are actual lats and lons available, find i and j based on lat and lon |
|---|
| 91 | ! otherwise, it is an idealized case and the user has specified i and j in the |
|---|
| 92 | ! turbines file read in by read_windturbines_in in module_wind_generic |
|---|
| 93 | DO iturbine = 1,nwindturbines ! nwindturbines defined in module_wind_generic |
|---|
| 94 | p => windturbines(iturbine) |
|---|
| 95 | IF ( id .EQ. p%id ) THEN |
|---|
| 96 | DO j=jts,jtf |
|---|
| 97 | DO i=its,itf |
|---|
| 98 | IF (xlat_v(i,j) .LE. p%lat .AND. p%lat .LT. xlat_v(i+1,j) .AND. & |
|---|
| 99 | xlong_u(i,j).LE. p%lon .AND. p%lon .LT. xlong_u(i,j+1)) THEN |
|---|
| 100 | p%i=i |
|---|
| 101 | p%j=j |
|---|
| 102 | ENDIF |
|---|
| 103 | ENDDO |
|---|
| 104 | ENDDO |
|---|
| 105 | ENDIF |
|---|
| 106 | ENDDO |
|---|
| 107 | ENDIF |
|---|
| 108 | ELSE IF ( windspec .EQ. WIND_TURBINES_IDEAL .AND. id .EQ. turbgridid ) THEN |
|---|
| 109 | CALL nl_get_td_ewfx(1,ewfx) |
|---|
| 110 | CALL nl_get_td_ewfy(1,ewfy) |
|---|
| 111 | CALL nl_get_td_pwfx(1,pwfx) |
|---|
| 112 | CALL nl_get_td_pwfy(1,pwfy) |
|---|
| 113 | CALL nl_get_td_hubheight(1,hubheight) |
|---|
| 114 | CALL nl_get_td_diameter(1,diameter) |
|---|
| 115 | CALL nl_get_td_power(1,power) |
|---|
| 116 | CALL nl_get_td_cutinspeed(1,cutinspeed) |
|---|
| 117 | CALL nl_get_td_cutoutspeed(1,cutoutspeed) |
|---|
| 118 | CALL nl_get_td_stdthrcoef(1,stdthrcoef) |
|---|
| 119 | ! count the turbines |
|---|
| 120 | n = 0 |
|---|
| 121 | DO j = jts,jtf |
|---|
| 122 | IF ( pwfy .LE. j .AND. j .LE. (pwfy+ewfy-1) ) THEN |
|---|
| 123 | DO i = its,itf |
|---|
| 124 | IF ( pwfx .LE. i .AND. i .LE. (pwfx+ewfx-1) ) THEN |
|---|
| 125 | n = n + 1 |
|---|
| 126 | ENDIF |
|---|
| 127 | ENDDO |
|---|
| 128 | ENDIF |
|---|
| 129 | ENDDO |
|---|
| 130 | nwindturbines = n |
|---|
| 131 | ALLOCATE(windturbines(nwindturbines)) |
|---|
| 132 | ! set the turbines |
|---|
| 133 | n = 0 |
|---|
| 134 | DO j = jts,jtf |
|---|
| 135 | IF ( pwfy .LE. j .AND. j .LE. (pwfy+ewfy-1) ) THEN |
|---|
| 136 | DO i = its,itf |
|---|
| 137 | IF ( pwfx .LE. i .AND. i .LE. (pwfx+ewfx-1) ) THEN |
|---|
| 138 | n = n + 1 |
|---|
| 139 | IF ( n .GT. nwindturbines ) THEN |
|---|
| 140 | CALL wrf_error_fatal('would overrun windturbines array') |
|---|
| 141 | ENDIF |
|---|
| 142 | windturbines(n)%id = id |
|---|
| 143 | windturbines(n)%lat = 0.0 |
|---|
| 144 | windturbines(n)%lon = 0.0 |
|---|
| 145 | windturbines(n)%i = i |
|---|
| 146 | windturbines(n)%j = j |
|---|
| 147 | windturbines(n)%hubheight = hubheight |
|---|
| 148 | windturbines(n)%diameter = diameter |
|---|
| 149 | windturbines(n)%stdthrcoef = stdthrcoef |
|---|
| 150 | windturbines(n)%power = power |
|---|
| 151 | windturbines(n)%cutinspeed = cutinspeed |
|---|
| 152 | windturbines(n)%cutoutspeed = cutoutspeed |
|---|
| 153 | ENDIF |
|---|
| 154 | ENDDO |
|---|
| 155 | ENDIF |
|---|
| 156 | ENDDO |
|---|
| 157 | ENDIF |
|---|
| 158 | inited(id) = .TRUE. |
|---|
| 159 | ENDIF |
|---|
| 160 | IF ( windspec .EQ. WIND_TURBINES_FROMLIST ) THEN |
|---|
| 161 | wfdensity = 1.0/(dx*dx) ! per turbine, so numerator is 1 |
|---|
| 162 | ELSE |
|---|
| 163 | wfdensity = turbpercell/(dx*dx) |
|---|
| 164 | ENDIF |
|---|
| 165 | |
|---|
| 166 | IF (inited(id) .AND. & |
|---|
| 167 | ((windspec .EQ. WIND_TURBINES_FROMLIST) .OR. & |
|---|
| 168 | (windspec .EQ. WIND_TURBINES_IDEAL .AND. id .EQ. turbgridid ))) THEN |
|---|
| 169 | DO iturbine = 1,nwindturbines ! nwindturbines defined in module_wind_generic |
|---|
| 170 | p => windturbines(iturbine) |
|---|
| 171 | IF ( id .EQ. p%id ) THEN |
|---|
| 172 | ! vertical layers cut by turbine blades |
|---|
| 173 | k_turbine_bot=0 !bottom level |
|---|
| 174 | k_turbine_top=-1 !top level |
|---|
| 175 | i = p%i |
|---|
| 176 | j = p%j |
|---|
| 177 | |
|---|
| 178 | IF (( its .LE. i .AND. i .LE. itf ) .AND. & |
|---|
| 179 | ( jts .LE. j .AND. j .LE. jtf ) ) THEN |
|---|
| 180 | |
|---|
| 181 | blade_l_point=p%hubheight-p%diameter/2. ! height of lower blade tip above ground (m) (theight=hub height) |
|---|
| 182 | blade_u_point=p%hubheight+p%diameter/2. ! height of upper blade tip above ground (m) |
|---|
| 183 | |
|---|
| 184 | kfound = .false. |
|---|
| 185 | zheightl=0.0 |
|---|
| 186 | ! find vertical levels cut by turbine blades |
|---|
| 187 | DO k=kts,ktf |
|---|
| 188 | IF(.NOT. kfound) THEN |
|---|
| 189 | zheightu = zheightl + dz(i,k,j) ! increment height |
|---|
| 190 | |
|---|
| 191 | IF(blade_l_point .GE. zheightl .AND. blade_l_point .LE. zheightu) THEN |
|---|
| 192 | k_turbine_bot=k ! lower blade tip cuts this level |
|---|
| 193 | ENDIF |
|---|
| 194 | |
|---|
| 195 | IF(blade_u_point .GE. zheightl .AND. blade_u_point .LE. zheightu) THEN |
|---|
| 196 | k_turbine_top=k ! upper blade tip cuts this level |
|---|
| 197 | kfound = .TRUE. |
|---|
| 198 | ENDIF |
|---|
| 199 | |
|---|
| 200 | zheightl = zheightu |
|---|
| 201 | ENDIF |
|---|
| 202 | ENDDO |
|---|
| 203 | IF ( kfound ) THEN |
|---|
| 204 | DO k=k_turbine_bot,k_turbine_top ! loop over turbine blade levels |
|---|
| 205 | |
|---|
| 206 | z1=phb(i,k,j)/g-blade_l_point-phb(i,1,j)/g ! distance between k level and lower blade tip |
|---|
| 207 | z2=phb(i,k+1,j)/g-blade_l_point-phb(i,1,j)/g ! distance between k+1 level and lower blade tip |
|---|
| 208 | |
|---|
| 209 | IF(z1 .LT. 0.) z1=0.0 ! k level lower than lower blade tip |
|---|
| 210 | IF(z2 .GT. p%diameter) z2=p%diameter ! k+1 level higher than turbine upper blade tip |
|---|
| 211 | |
|---|
| 212 | ! magnitude of horizontal velocity |
|---|
| 213 | speed=sqrt(u(i,k,j)*u(i,k,j)+v(i,k,j)*v(i,k,j)) |
|---|
| 214 | |
|---|
| 215 | ! calculate TKE, power and thrust coeffs |
|---|
| 216 | CALL dragcof(tkecof,powcof,thrcof, & |
|---|
| 217 | speed,p%cutinspeed,p%cutoutspeed, & |
|---|
| 218 | p%power,p%diameter,p%stdthrcoef) |
|---|
| 219 | |
|---|
| 220 | CALL turbine_area(z1,z2,p%diameter,wfdensity,tarea) |
|---|
| 221 | |
|---|
| 222 | ! output TKE tendency |
|---|
| 223 | qke(i,k,j) = qke(i,k,j)+speed*speed*speed*tarea*tkecof*dt/dz(i,k,j) |
|---|
| 224 | ! output u tendency |
|---|
| 225 | du(i,k,j) = du(i,k,j)-.5*u(i,k,j)*thrcof*speed*tarea/dz(i,k,j) |
|---|
| 226 | ! output v tendency |
|---|
| 227 | dv(i,k,j) = dv(i,k,j)-.5*v(i,k,j)*thrcof*speed*tarea/dz(i,k,j) |
|---|
| 228 | |
|---|
| 229 | ENDDO |
|---|
| 230 | ENDIF |
|---|
| 231 | ENDIF |
|---|
| 232 | ENDIF |
|---|
| 233 | ENDDO |
|---|
| 234 | ENDIF |
|---|
| 235 | END SUBROUTINE turbine_drag |
|---|
| 236 | |
|---|
| 237 | ! calculates area of turbine between two vertical levels |
|---|
| 238 | ! Input variables : |
|---|
| 239 | ! z1 = distance between k level and lower blade tip |
|---|
| 240 | ! z2 = distance between k+1 level and lower blade tip |
|---|
| 241 | ! wfdensity = wind farm density in m^-2 |
|---|
| 242 | ! tdiameter = turbine diameter |
|---|
| 243 | ! Output variable : |
|---|
| 244 | ! tarea = area of turbine between two levels * wfdensity |
|---|
| 245 | SUBROUTINE turbine_area(z1,z2,tdiameter,wfdensity,tarea) |
|---|
| 246 | |
|---|
| 247 | REAL, INTENT(IN) ::tdiameter,wfdensity |
|---|
| 248 | REAL, INTENT(INOUT) ::z1,z2 |
|---|
| 249 | REAL, INTENT(OUT):: tarea |
|---|
| 250 | REAL r,zc1,zc2 |
|---|
| 251 | |
|---|
| 252 | r=tdiameter/2. !r = turbine radius |
|---|
| 253 | z1=r-z1 !distance of kth level from turbine center |
|---|
| 254 | z2=r-z2 !distance of k+1 th level from turbine center |
|---|
| 255 | zc1=abs(z1) |
|---|
| 256 | zc2=abs(z2) |
|---|
| 257 | !turbine area between z1 and z2 |
|---|
| 258 | IF(z1 .GT. 0. .AND. z2 .GT. 0) THEN |
|---|
| 259 | tarea=zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r)- & |
|---|
| 260 | (zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r)) |
|---|
| 261 | ELSE IF(z1 .LT. 0. .AND. z2 .LT. 0) THEN |
|---|
| 262 | tarea=zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r)- & |
|---|
| 263 | (zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r)) |
|---|
| 264 | ELSE |
|---|
| 265 | tarea=zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r)+ & |
|---|
| 266 | zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r) |
|---|
| 267 | ENDIF |
|---|
| 268 | tarea=tarea*wfdensity !turbine area * wind farm density |
|---|
| 269 | |
|---|
| 270 | END SUBROUTINE turbine_area |
|---|
| 271 | |
|---|
| 272 | ! Caculates tke, power and thrust coefficients as function of horiz wind speed |
|---|
| 273 | ! from fit to turbine power curve - needs to be changed for particular turbine used |
|---|
| 274 | |
|---|
| 275 | ! tkecof = tke coefficient |
|---|
| 276 | ! powcof = power coefficient |
|---|
| 277 | ! thrcof = thrust coefficient |
|---|
| 278 | ! cispeed = cut-in speed in m/s |
|---|
| 279 | ! cospeed = cut-out speed in m/s |
|---|
| 280 | ! tpower = turbine power in MW |
|---|
| 281 | ! speed = horiz wind speed in m/s |
|---|
| 282 | ! tdiameter = turbine diameter in m |
|---|
| 283 | ! stdthrcoef = standing thrust coefficient |
|---|
| 284 | |
|---|
| 285 | SUBROUTINE dragcof(tkecof,powcof,thrcof,speed,cispeed,cospeed, & |
|---|
| 286 | tpower,tdiameter,stdthrcoef) |
|---|
| 287 | |
|---|
| 288 | ! DISCLAIMER: The following power curve, power coefficients, and thrust |
|---|
| 289 | ! coefficients are meant for testing purposes only, and were formulated as |
|---|
| 290 | ! an approximation to a real curve. The user is strongly encouraged to |
|---|
| 291 | ! incorporate their own curves for the particular turbine of interest |
|---|
| 292 | ! to them. |
|---|
| 293 | |
|---|
| 294 | REAL, PARAMETER :: pi=22./7. |
|---|
| 295 | REAL, INTENT(IN):: speed, cispeed, cospeed, tpower,tdiameter,stdthrcoef |
|---|
| 296 | REAL, INTENT(OUT):: tkecof,powcof,thrcof |
|---|
| 297 | REAL :: power,area,mspeed,hspeed |
|---|
| 298 | |
|---|
| 299 | area=pi/4.*tdiameter*tdiameter ! area swept by turbine blades |
|---|
| 300 | |
|---|
| 301 | ! GENERIC POWER CURVE - USE AT YOUR OWN RISK |
|---|
| 302 | mspeed=0.5*(cospeed+cispeed) !average of cispeed & cospeed |
|---|
| 303 | hspeed=0.5*(cospeed-mspeed) !this regulates the transition to full power |
|---|
| 304 | power =tpower*(.5*tanh((speed - (mspeed-hspeed))/(hspeed*0.60)) + .5)*.8 |
|---|
| 305 | |
|---|
| 306 | ! GENERIC power coefficient calculation - USE AT YOUR OWN RISK |
|---|
| 307 | IF(speed .LE. cispeed .OR. speed .GE. cospeed) THEN |
|---|
| 308 | power=0. |
|---|
| 309 | powcof=0. |
|---|
| 310 | ELSE |
|---|
| 311 | powcof = power * 2.e+6 / (speed*speed*speed*area) |
|---|
| 312 | IF (speed .LT. cispeed*2.) THEN ! dampen artificial max near cispeed |
|---|
| 313 | powcof = powcof * exp(-((speed-cispeed*2.)**2./(cispeed*2.))) |
|---|
| 314 | end if |
|---|
| 315 | powcof = MIN(powcof,.55) |
|---|
| 316 | ENDIF |
|---|
| 317 | |
|---|
| 318 | ! GENERIC Thrust coefficient calculation - USE AT YOUR OWN RISK |
|---|
| 319 | IF (speed .LE. cispeed .OR. speed .GE. cospeed) THEN |
|---|
| 320 | thrcof = stdthrcoef |
|---|
| 321 | ELSE |
|---|
| 322 | !thrcof= stdthrcoef+2.3/speed**.8 |
|---|
| 323 | thrcof = powcof + powcof*0.75 |
|---|
| 324 | thrcof = MIN(thrcof,.9) |
|---|
| 325 | thrcof = MAX(thrcof,stdthrcoef) |
|---|
| 326 | ENDIF |
|---|
| 327 | |
|---|
| 328 | ! tke coefficient calculation |
|---|
| 329 | tkecof=thrcof-powcof |
|---|
| 330 | IF(tkecof .LT. 0.) tkecof=0. |
|---|
| 331 | |
|---|
| 332 | END SUBROUTINE dragcof |
|---|
| 333 | |
|---|
| 334 | SUBROUTINE init_module_wind_fitch |
|---|
| 335 | inited = .FALSE. |
|---|
| 336 | END SUBROUTINE init_module_wind_fitch |
|---|
| 337 | |
|---|
| 338 | END MODULE module_wind_fitch |
|---|
![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
