wrfv2_fire/phys/module_wind_fitch.F

   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,j+1) .AND. &
  99                     xlong_u(i,j).LE. p%lon .AND. p%lon .LT. xlong_u(i+1,j)) 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