phys/module_wind_fitch.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2
   3 MODULE module_wind_fitch
   4 !
   5 !Represents kinetic energy extracted by wind turbines and turbulence
   6 ! (TKE) they produce at model levels within the rotor area.
   7 !
   8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   9 !! NOTICE
  10 !! The following paper should be cited whenever presenting results using this scheme
  11 !! (using either the original version or any modified versions of the scheme):
  12 !! Fitch, A. C. et al. 2012: Local and Mesoscale Impacts of Wind Farms as Parameterized in a
  13 !! Mesoscale NWP Model. Monthly Weather Review, doi:http://dx.doi.org/10.1175/MWR-D-11-00352.1
  14 !!
  15 !! Anna C. Fitch, National Center for Atmospheric Research (formerly University of Bergen)
  16 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  17 !
  18 ! History of changes:
  19 !
  20 ! WRFV3.5.1:
  21 ! WRFV3.6:   Modified by Pedro A. Jimenez to include:
  22 !             - Initialize the wind turbines in this module.
  23 !             - Introduce z_at_walls to avoid instabilities due to neglecting
  24 !                the perturbation of the geopotential height.
  25 !             - User friendly interface to introduce the technical characteritics of
  26 !                the wind turbines.
  27 !             - Only uses one set of turbine coefficients using the wind speed at hub height
  28 !             - Two standing coefficients.
  29 !             - Calculates the power produced by the wind turbines.
  30 !
  31 ! References:
  32 !
  33 ! Fitch, A. C. et al. 2012: Local and Mesoscale Impacts of Wind Farms as Parameterized in a
  34 !    Mesoscale NWP Model. Monthly Weather Review, doi:http://dx.doi.org/10.1175/MWR-D-11-00352.1
  35 ! Fitch, A. C. et al. 2013: Mesoscale Influences of Wind Farms Throughout a Diurnal Cycle.
  36 !    Monthly Weather Review, doi:http://dx.doi.org/10.1175/MWR-D-12-00185.1
  37 ! Fitch, A. C. et al. 2013: Parameterization of Wind Farms in Climate Models.
  38 !    Journal of Climate, doi:http://dx.doi.org/10.1175/JCLI-D-12-00376.1
  39 ! Jimenez, P.A., J. Navarro, A.M. Palomares and J. Dudhia:  Mesoscale modeling of offshore wind turbines
  40 !    wakes at the wind farm resolving scale: a composite-based analysis with the WRF model over Horns Rev.
  41 !    Wind Energy, (In Press.).
  42 !
  43 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  44 !
  45   USE module_driver_constants, ONLY : max_domains
  46   USE module_model_constants, ONLY :  piconst
  47 !
  48   USE module_llxy
  49   USE module_dm, ONLY : wrf_dm_min_real
  50   USE module_configure, ONLY : grid_config_rec_type
  51
  52   IMPLICIT NONE
  53
  54   INTEGER, PARAMETER :: MAXVALS  = 100
  55   INTEGER, PARAMETER :: MAXVALS2 = 100
  56 !
  57   INTEGER           :: nt
  58   INTEGER, DIMENSION(:), ALLOCATABLE :: NKIND
  59   INTEGER, DIMENSION(:,:), ALLOCATABLE :: ival,jval
  60   REAL, DIMENSION(:), ALLOCATABLE :: hubheight,diameter,stc,stc2,cutin,cutout,npower
  61 !
  62   REAL :: turbws(maxvals,maxvals2),turbtc(maxvals,maxvals2),turbpw(maxvals,maxvals2)
  63   REAL :: correction_factor
  64 !
  65 CONTAINS
  66
  67   SUBROUTINE  dragforce(                      &
  68        & id                                      &
  69        &,z_at_w,u,v                 &
  70        &,dx,dz,dt,qke                            &
  71        &,du,dv                                   &
  72        &,windfarm_opt,power                      &
  73        &,ids,ide,jds,jde,kds,kde                 &
  74        &,ims,ime,jms,jme,kms,kme                 &
  75        &,its,ite,jts,jte,kts,kte                 &
  76        &)
  77 !
  78 !
  79 !
  80   INTEGER, INTENT(IN) :: id,windfarm_opt
  81   INTEGER, INTENT(IN) :: its,ite,jts,jte,kts,kte
  82   INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme
  83   INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde
  84   REAL, INTENT(IN) :: dx,dt
  85   REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: dz,u,v,z_at_w
  86   REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: du,dv,qke
  87   REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: power
  88 !
  89 ! Local
  90 !
  91   REAL     blade_l_point,blade_u_point,zheightl,zheightu,z1,z2,tarea
  92   REAL     speed,tkecof,powcof,thrcof,wfdensity
  93   INTEGER  itf,jtf,ktf
  94   INTEGER  i,j,k,n
  95   INTEGER  k_turbine_bot, k_turbine_top
  96
  97   LOGICAL :: kfound
  98 !
  99 ! ... PAJ: more variables ...
 100 !
 101   REAL :: speedhub,speed1,speed2
 102   real :: power1,power2,area,ec
 103   INTEGER :: kbot,ktop,kt
 104
 105   itf=MIN0(ite,ide-1)
 106   jtf=MIN0(jte,jde-1)
 107   ktf=MIN0(kte,kde-1)
 108
 109     wfdensity = 1.0/(dx*dx)   !  per turbine, so numerator is 1
 110     power=0.
 111
 112     DO kt = 1,nt
 113       IF ( windfarm_opt .eq. 1 ) THEN
 114 !
 115 ! vertical layers cut by turbine blades
 116 !
 117         k_turbine_bot=0      !bottom level
 118         k_turbine_top=-1     !top level
 119         i = ival(kt,id)
 120         j = jval(kt,id)
 121 !
 122          if (i.ne.-9999.and.j.ne.-9999) then
 123         IF (( its .LE. i .AND. i .LE. itf ) .AND. &
 124             ( jts .LE. j .AND. j .LE. jtf )  ) THEN
 125 !
 126           blade_l_point=hubheight(kt)-diameter(kt)/2. ! height of lower blade tip above ground (m)
 127           blade_u_point=hubheight(kt)+diameter(kt)/2. ! height of upper blade tip above ground (m)
 128 !
 129           kfound = .false.
 130           zheightl=0.0
 131           ! find vertical levels cut by turbine blades
 132           DO k=kts,ktf
 133             IF(.NOT. kfound) THEN
 134               zheightu = zheightl + dz(i,k,j) ! increment height
 135
 136               IF(blade_l_point .GE. zheightl .AND. blade_l_point .LE. zheightu) THEN
 137                 k_turbine_bot=k ! lower blade tip cuts this level
 138               ENDIF
 139
 140               IF(blade_u_point .GE. zheightl .AND. blade_u_point .LE. zheightu) THEN
 141                 k_turbine_top=k ! upper blade tip cuts this level
 142                 kfound = .TRUE.
 143               ENDIF
 144
 145               zheightl = zheightu
 146             ENDIF
 147           ENDDO
 148           IF ( kfound ) THEN
 149 !
 150 ! ... PAJ: Changes introduced to compute only one set of turbine coefficients ...
 151 !          First computes the wind speed at the hub height.
 152 !
 153           kfound = .false.
 154           zheightl=0.
 155           ! find vertical levels (half levels) within the hub height
 156           DO k=kts,ktf
 157             IF(.NOT. kfound) THEN
 158               z2 = zheightl + 0.5*dz(i,k,j)
 159 !
 160               IF(hubheight(kt) .GE. z2 ) THEN
 161                 kbot=k
 162               ELSE
 163                 ktop=k
 164                 kfound = .TRUE.
 165               ENDIF
 166 !
 167               if (.NOT. kfound) z1=z2
 168               zheightl = z2 + 0.5*dz(i,k,j)
 169             ENDIF
 170           ENDDO
 171 !
 172           speed1=0.
 173           speed2=0.
 174           if (ktop.eq.1) then
 175            speedhub=sqrt(u(i,1,j)**2.+v(i,1,j)**2.)*hubheight(kt)/z1
 176           else
 177            speed1=sqrt(u(i,kbot,j)**2.+v(i,kbot,j)**2.)
 178            speed2=sqrt(u(i,ktop,j)**2.+v(i,ktop,j)**2.)
 179            speedhub=speed1+((speed2-speed1)/(z2-z1))*(hubheight(kt)-z1)
 180           endif
 181 !
 182 ! ... calculate TKE, power and thrust coeffs
 183 !
 184               CALL dragcof(tkecof,powcof,thrcof,               &
 185                            speedhub,cutin(kt),cutout(kt),   &
 186                            npower(kt),diameter(kt),stc(kt),stc2(kt),nkind(kt))
 187 !
 188 ! ... PAJ: Computation of power generated by the wind turbine ...
 189 !
 190           area=piconst/4.*diameter(kt)**2.          ! area swept by turbine blades
 191           power1=0.5*1.23*speedhub**3.*area*powcof
 192           power(i,j)=power1+power(i,j)
 193           power2=0.
 194 !
 195             DO k=k_turbine_bot,k_turbine_top ! loop over turbine blade levels
 196               z1=z_at_w(i,k,j)-blade_l_point-z_at_w(i,1,j)  ! distance between k level and lower blade tip
 197               z2=z_at_w(i,k+1,j)-blade_l_point-z_at_w(i,1,j) ! distance between k+1 level and lower blade tip
 198               IF(z1 .LT. 0.) z1=0.0 ! k level lower than lower blade tip
 199               IF(z2 .GT. diameter(kt)) z2=diameter(kt) ! k+1 level higher than turbine upper blade tip
 200               CALL turbine_area(z1,z2,diameter(kt),wfdensity,tarea)
 201 !
 202               speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.)
 203               power2=power2+0.5*powcof*1.23*(speed**3.)*tarea/wfdensity
 204             ENDDO
 205 !
 206 ! ... PAJ: Computes the tendencies of TKE and momentum ...
 207 !
 208             DO k=k_turbine_bot,k_turbine_top ! loop over turbine blade levels
 209               z1=z_at_w(i,k,j)-blade_l_point-z_at_w(i,1,j)  ! distance between k lev and lower blade tip
 210               z2=z_at_w(i,k+1,j)-blade_l_point-z_at_w(i,1,j) !distance between k+1 lev and lower blade tip
 211               IF(z1 .LT. 0.) z1=0.0 ! k level lower than lower blade tip
 212               IF(z2 .GT. diameter(kt)) z2=diameter(kt) ! k+1 level higher than turbine upper blade tip
 213 !
 214               CALL turbine_area(z1,z2,diameter(kt),wfdensity,tarea)
 215 !
 216               speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.)
 217 !`
 218 ! ... PAJ: normalization introduced to conserve energy ...
 219 !
 220               if (power1.eq.0.or.power2.eq.0) then
 221               ec=1.
 222               else
 223               ec=power1/power2
 224               endif
 225 !
 226               ! output TKE
 227               qke(i,k,j) = qke(i,k,j)+speed**3.*tarea*tkecof*dt/dz(i,k,j)*ec
 228               ! output u tendency
 229               du(i,k,j) = du(i,k,j)-.5*u(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec
 230               ! output v tendency
 231               dv(i,k,j) = dv(i,k,j)-.5*v(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec
 232             ENDDO
 233           ENDIF
 234         ENDIF
 235         endif
 236       ENDIF
 237     ENDDO
 238
 239   END SUBROUTINE dragforce
 240
 241 ! This subroutine calculates area of turbine between two vertical levels
 242 ! Input variables :
 243 !            z1 = distance between k level and lower blade tip
 244 !            z2 = distance between k+1 level and lower blade tip
 245 !            wfdensity = wind farm density in m^-2
 246 !     tdiameter = turbine diameter
 247 ! Output variable :
 248 !         tarea = area of turbine between two levels * wfdensity
 249   SUBROUTINE turbine_area(z1,z2,tdiameter,wfdensity,tarea)
 250
 251   REAL, INTENT(IN) ::tdiameter,wfdensity
 252   REAL, INTENT(INOUT) ::z1,z2
 253   REAL, INTENT(OUT):: tarea
 254   REAL r,zc1,zc2
 255
 256   r=tdiameter/2.              !r = turbine radius
 257   z1=r-z1                   !distance of kth level from turbine center
 258   z2=r-z2                   !distance of k+1 th level from turbine center
 259   zc1=abs(z1)
 260   zc2=abs(z2)
 261   !turbine area between z1 and z2
 262   IF(z1 .GT. 0. .AND. z2 .GT. 0.) THEN
 263      tarea=zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r)- &
 264      (zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r))
 265   ELSE IF(z1 .LT. 0. .AND. z2 .LT. 0.) THEN
 266      tarea=zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r)- &
 267      (zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r))
 268   ELSE
 269      tarea=zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r)+ &
 270      zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r)
 271   ENDIF
 272   tarea=tarea*wfdensity      !turbine area * wind farm density
 273
 274   END SUBROUTINE turbine_area
 275
 276
 277   SUBROUTINE dragcof(tkecof,powcof,thrcof,speed,cispeed,cospeed, &
 278                      tpower,tdiameter,stdthrcoef,stdthrcoef2,nkind)
 279
 280
 281   REAL, INTENT(IN):: speed, cispeed, cospeed, tpower,tdiameter,stdthrcoef,stdthrcoef2
 282   REAL, INTENT(OUT):: tkecof,powcof,thrcof
 283   REAL :: power,area,mspeed,hspeed
 284 !
 285 ! ... PAJ ...
 286 !
 287    INTEGER :: nkind,k,nu,nb
 288    LOGICAL :: vfound
 289    REAL :: fac1,fac2
 290
 291   area=piconst/4.*tdiameter**2.          ! area swept by turbine blades
 292
 293       vfound=.false.
 294       DO k=1,maxvals2
 295             IF(.NOT. vfound) THEN
 296               IF(turbws(nkind,k).GT.speed) THEN
 297                 nu=k
 298                 nb=k-1
 299                 vfound=.true.
 300               ENDIF
 301             ENDIF
 302       ENDDO
 303 !
 304   IF (speed .LE. cispeed) THEN
 305      thrcof = stdthrcoef
 306   ELSE
 307     IF (speed .GE. cospeed) THEN
 308      thrcof = stdthrcoef2
 309      ELSE
 310      thrcof = turbtc(nkind,nb)+(turbtc(nkind,nu)-turbtc(nkind,nb))/(turbws(nkind,nu)-turbws(nkind,nb))*(speed-turbws(nkind,nb))
 311     ENDIF
 312   ENDIF
 313 !
 314 ! ... power coeficient ...
 315 !
 316   IF(speed .LE. cispeed .OR. speed .GE. cospeed) THEN
 317      power=0.
 318      powcof=0.
 319   ELSE
 320       fac1=1000./(0.5*1.23*turbws(nkind,nb)**3.*area)
 321       fac2=1000./(0.5*1.23*turbws(nkind,nu)**3.*area)
 322       power = turbpw(nkind,nb)+(turbpw(nkind,nu)-turbpw(nkind,nb))/(turbws(nkind,nu)-turbws(nkind,nb)) &
 323                                *(speed-turbws(nkind,nb))
 324       powcof = turbpw(nkind,nb)*fac1+(turbpw(nkind,nu)*fac2-turbpw(nkind,nb)*fac1)/(turbws(nkind,nu)-turbws(nkind,nb)) &
 325                                      *(speed-turbws(nkind,nb))
 326   ENDIF
 327 !
 328   ! tke coefficient calculation
 329
 330   tkecof=thrcof-powcof
 331   tkecof = correction_factor * tkecof
 332   IF(tkecof .LT. 0.) tkecof=0.
 333 !
 334   END SUBROUTINE dragcof
 335 !
 336   SUBROUTINE init_module_wind_fitch(id,config_flags,xlong,xlat,windfarm_initialized,&
 337                                             ims,ime,jms,jme,its,ite,jts,jte,ids,ide,jds,jde)
 338 !
 339   IMPLICIT NONE
 340 !
 341    integer ims,ime,jms,jme,ids,ide,jds,jde
 342    integer its,ite,jts,jte
 343    REAL,     DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: xlong,xlat
 344    TYPE (grid_config_rec_type) :: config_flags
 345    TYPE (PROJ_INFO) :: ts_proj
 346    logical :: windfarm_initialized
 347 !
 348    CHARACTER*256 num,input,message_wind
 349    real lat,lon,ts_rx,ts_ry
 350    REAL :: known_lat, known_lon
 351    INTEGER i,j,nval,k,id
 352
 353    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
 354 !
 355
 356    correction_factor = config_flags%windfarm_tke_factor
 357       IF ( wrf_dm_on_monitor() ) THEN
 358 !
 359 ! ... PAJ: Opens the file with the location of the wind turbines ...
 360 !
 361         if ( config_flags%windfarm_ij .eq. 1 ) then
 362           open(70,file='windturbines-ij.txt',form='formatted',status='old')
 363         else
 364           open(70,file='windturbines.txt',form='formatted',status='old')
 365         end if
 366 !
 367 ! ... PAJ: Counts the turbines ...
 368 !
 369        nt=0
 370  10    read(70,*,end=100)
 371        nt=nt+1
 372        goto 10
 373 !
 374  100   continue
 375        rewind (70)
 376      END IF
 377 !
 378      CALL wrf_dm_bcast_integer(nt,1)
 379 !
 380 ! ... PAJ: Initializes the configuration of the wind farm(s) ...
 381 !
 382      if (.not. windfarm_initialized) then
 383        allocate (nkind(nt),ival(nt,max_domains),jval(nt,max_domains))
 384        allocate (hubheight(nt),stc(nt),stc2(nt),cutin(nt),cutout(nt),diameter(nt),npower(nt))
 385        ival=-9999
 386        jval=-9999
 387        windfarm_initialized=.true.
 388      endif
 389 !
 390      IF ( wrf_dm_on_monitor() ) THEN
 391      do k=1,nt
 392        if ( config_flags%windfarm_ij .eq. 1 ) then
 393          read(70,*) ival(k,id), jval(k,id), nkind(k)
 394          write(message_wind,*)'WINDFARM Turbine #',k,': I, J = ',ival(k,id), jval(k,id),'; Type = ',nkind(k)
 395          CALL wrf_message(message_wind)
 396
 397        else
 398
 399          read(70,*)lat,lon,nkind(k)
 400          write(message_wind,*)'WINDFARM Turbine #',k,': Lat, lon = ',lat,lon,'; Type = ',nkind(k)
 401          CALL wrf_message(message_wind)
 402
 403          CALL map_init(ts_proj)
 404
 405          known_lat = xlat(its,jts)
 406          known_lon = xlong(its,jts)
 407
 408       ! Mercator
 409       IF (config_flags%map_proj == PROJ_MERC) THEN
 410          CALL map_set(PROJ_MERC, ts_proj,               &
 411                       truelat1 = config_flags%truelat1, &
 412                       lat1     = known_lat,             &
 413                       lon1     = known_lon,             &
 414                       knowni   = REAL(its),             &
 415                       knownj   = REAL(jts),             &
 416                       dx       = config_flags%dx)
 417
 418       ! Lambert conformal
 419       ELSE IF (config_flags%map_proj == PROJ_LC) THEN
 420          CALL map_set(PROJ_LC, ts_proj,                  &
 421                       truelat1 = config_flags%truelat1,  &
 422                       truelat2 = config_flags%truelat2,  &
 423                       stdlon   = config_flags%stand_lon, &
 424                       lat1     = known_lat,              &
 425                       lon1     = known_lon,              &
 426                       knowni   = REAL(its),              &
 427                       knownj   = REAL(jts),              &
 428                       dx       = config_flags%dx)
 429 !      ! Polar stereographic
 430       ELSE IF (config_flags%map_proj == PROJ_PS) THEN
 431          CALL map_set(PROJ_PS, ts_proj,                  &
 432                       truelat1 = config_flags%truelat1,  &
 433                       stdlon   = config_flags%stand_lon, &
 434                       lat1     = known_lat,              &
 435                       lon1     = known_lon,              &
 436                       knowni   = REAL(its),              &
 437                       knownj   = REAL(jts),              &
 438                       dx       = config_flags%dx)
 439 !#if (EM_CORE == 1)
 440 !      ! Cassini (global ARW)
 441 !      ELSE IF (config_flags%map_proj == PROJ_CASSINI) THEN
 442 !         CALL map_set(PROJ_CASSINI, ts_proj,                            &
 443 !                      latinc   = grid%dy*360.0/(2.0*EARTH_RADIUS_M*PI), &
 444 !                      loninc   = grid%dx*360.0/(2.0*EARTH_RADIUS_M*PI), &
 445 !                      lat1     = known_lat,                             &
 446 !                      lon1     = known_lon,                             &
 447 !                      lat0     = config_flags%pole_lat,                 &
 448 !                      lon0     = config_flags%pole_lon,                 &
 449 !                      knowni   = 1.,                                    &
 450 !                      knownj   = 1.,                                    &
 451 !                      stdlon   = config_flags%stand_lon)
 452 !#endif
 453 !
 454 !      ! Rotated latitude-longitude
 455 !      ELSE IF (config_flags%map_proj == PROJ_ROTLL) THEN
 456 !         CALL map_set(PROJ_ROTLL, ts_proj,                      &
 457 !! I have no idea how this should work for NMM nested domains
 458 !                      ixdim    = grid%e_we-1,                   &
 459 !                      jydim    = grid%e_sn-1,                   &
 460 !                      phi      = real(grid%e_sn-2)*grid%dy/2.0, &
 461 !                      lambda   = real(grid%e_we-2)*grid%dx,     &
 462 !                      lat1     = config_flags%cen_lat,          &
 463 !                      lon1     = config_flags%cen_lon,          &
 464 !                      latinc   = grid%dy,                       &
 465 !                      loninc   = grid%dx,                       &
 466 !                      stagger  = HH)
 467 !
 468       END IF
 469 !
 470          CALL latlon_to_ij(ts_proj, lat, lon, ts_rx, ts_ry)
 471 !
 472           ival(k,id)=nint(ts_rx)
 473           jval(k,id)=nint(ts_ry)
 474           if (ival(k,id).lt.ids.and.ival(k,id).gt.ide) then
 475             ival(k,id)=-9999
 476             jval(k,id)=-9999
 477           endif
 478 !
 479 !         write(73,*) k,id,ival(k,id),jval(k,id)
 480           write(message_wind,*)'WINDFARM Turbine #',k,': Lat, lon = ',lat,lon, &
 481                                ', (i,j) = (',ival(k,id),',',jval(k,id),'); Type = ',nkind(k)
 482           CALL wrf_debug(0,message_wind)
 483 !
 484      end if
 485 !
 486      enddo
 487       close(70)
 488 !
 489 ! ... PAJ: Read the tables with the turbine's characteristics ...
 490 !
 491          turbws=0.
 492          turbtc=0.
 493          turbpw=0.
 494          DO i=1,nt
 495           write(num,*) nkind(i)
 496           num=adjustl(num)
 497           input="wind-turbine-"//trim(num)//".tbl"
 498           OPEN(file=TRIM(input),unit=19,FORM='FORMATTED',STATUS='OLD')
 499           READ (19,*,ERR=132)nval
 500           READ(19,*,ERR=132)hubheight(i),diameter(i),stc(i),npower(i)
 501             DO k=1,nval
 502               READ(19,*,ERR=132)turbws(nkind(i),k),turbtc(nkind(i),k),turbpw(nkind(i),k)
 503             ENDDO
 504           cutin(i)  = turbws(nkind(i),1)
 505           cutout(i) = turbws(nkind(i),nval)
 506           stc2(i) = turbtc(nkind(i),nval)
 507           close (19)
 508          ENDDO
 509
 510  132   continue
 511 !
 512 ! ... ...
 513 !
 514       endif
 515
 516         CALL wrf_dm_bcast_integer(ival,nt*max_domains)
 517         CALL wrf_dm_bcast_integer(jval,nt*max_domains)
 518         CALL wrf_dm_bcast_real(hubheight,nt)
 519         CALL wrf_dm_bcast_real(diameter,nt)
 520         CALL wrf_dm_bcast_real(stc,nt)
 521         CALL wrf_dm_bcast_real(npower,nt)
 522         CALL wrf_dm_bcast_real(cutin,nt)
 523         CALL wrf_dm_bcast_real(cutout,nt)
 524         CALL wrf_dm_bcast_integer(nkind,nt)
 525         CALL wrf_dm_bcast_real(turbws,maxvals*maxvals2)
 526         CALL wrf_dm_bcast_real(turbtc,maxvals*maxvals2)
 527         CALL wrf_dm_bcast_real(turbpw,maxvals*maxvals2)
 528
 529   END SUBROUTINE init_module_wind_fitch
 530
 531 END MODULE module_wind_fitch