MODULE module_wind_fitch USE module_driver_constants, ONLY : max_domains USE module_model_constants, ONLY : piconst USE module_llxy USE module_dm, ONLY : wrf_dm_min_real USE module_configure, ONLY : grid_config_rec_type IMPLICIT NONE INTEGER, PARAMETER :: MAXVALS = 100 INTEGER, PARAMETER :: MAXVALS2 = 100 INTEGER :: nt INTEGER, DIMENSION(:), ALLOCATABLE :: NKIND INTEGER, DIMENSION(:,:), ALLOCATABLE :: ival,jval REAL, DIMENSION(:), ALLOCATABLE :: hubheight,diameter,stc,stc2,cutin,cutout,npower REAL :: turbws(maxvals,maxvals2),turbtc(maxvals,maxvals2),turbpw(maxvals,maxvals2) CONTAINS SUBROUTINE dragforce( & & id & &,z_at_w,u,v & &,dx,dz,dt,qke & &,du,dv & &,windfarm_opt,power & &,ids,ide,jds,jde,kds,kde & &,ims,ime,jms,jme,kms,kme & &,its,ite,jts,jte,kts,kte & &) INTEGER, INTENT(IN) :: id,windfarm_opt INTEGER, INTENT(IN) :: its,ite,jts,jte,kts,kte INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde REAL, INTENT(IN) :: dx,dt REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: dz,u,v,z_at_w REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: du,dv,qke REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: power REAL blade_l_point,blade_u_point,zheightl,zheightu,z1,z2,tarea REAL speed,tkecof,powcof,thrcof,wfdensity INTEGER itf,jtf,ktf INTEGER i,j,k,n INTEGER k_turbine_bot, k_turbine_top LOGICAL :: kfound REAL :: speedhub,speed1,speed2 real :: power1,power2,area,ec INTEGER :: kbot,ktop,kt itf=MIN0(ite,ide-1) jtf=MIN0(jte,jde-1) ktf=MIN0(kte,kde-1) wfdensity = 1.0/(dx*dx) power=0. DO kt = 1,nt IF ( windfarm_opt .eq. 1 ) THEN k_turbine_bot=0 k_turbine_top=-1 i = ival(kt,id) j = jval(kt,id) if (i.ne.-9999.and.j.ne.-9999) then IF (( its .LE. i .AND. i .LE. itf ) .AND. & ( jts .LE. j .AND. j .LE. jtf ) ) THEN blade_l_point=hubheight(kt)-diameter(kt)/2. blade_u_point=hubheight(kt)+diameter(kt)/2. kfound = .false. zheightl=0.0 DO k=kts,ktf IF(.NOT. kfound) THEN zheightu = zheightl + dz(i,k,j) IF(blade_l_point .GE. zheightl .AND. blade_l_point .LE. zheightu) THEN k_turbine_bot=k ENDIF IF(blade_u_point .GE. zheightl .AND. blade_u_point .LE. zheightu) THEN k_turbine_top=k kfound = .TRUE. ENDIF zheightl = zheightu ENDIF ENDDO IF ( kfound ) THEN kfound = .false. zheightl=0. DO k=kts,ktf IF(.NOT. kfound) THEN z2 = zheightl + 0.5*dz(i,k,j) IF(hubheight(kt) .GE. z2 ) THEN kbot=k ELSE ktop=k kfound = .TRUE. ENDIF if (.NOT. kfound) z1=z2 zheightl = z2 + 0.5*dz(i,k,j) ENDIF ENDDO speed1=0. speed2=0. if (ktop.eq.1) then speedhub=sqrt(u(i,1,j)**2.+v(i,1,j)**2.)*hubheight(kt)/z1 else speed1=sqrt(u(i,kbot,j)**2.+v(i,kbot,j)**2.) speed2=sqrt(u(i,ktop,j)**2.+v(i,ktop,j)**2.) speedhub=speed1+((speed2-speed1)/(z2-z1))*(hubheight(kt)-z1) endif CALL dragcof(tkecof,powcof,thrcof, & speedhub,cutin(kt),cutout(kt), & npower(kt),diameter(kt),stc(kt),stc2(kt),nkind(kt)) area=piconst/4.*diameter(kt)**2. power1=0.5*1.23*speedhub**3.*area*powcof power(i,j)=power1+power(i,j) power2=0. DO k=k_turbine_bot,k_turbine_top z1=z_at_w(i,k,j)-blade_l_point-z_at_w(i,1,j) z2=z_at_w(i,k+1,j)-blade_l_point-z_at_w(i,1,j) IF(z1 .LT. 0.) z1=0.0 IF(z2 .GT. diameter(kt)) z2=diameter(kt) CALL turbine_area(z1,z2,diameter(kt),wfdensity,tarea) speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.) power2=power2+0.5*powcof*1.23*(speed**3.)*tarea/wfdensity ENDDO DO k=k_turbine_bot,k_turbine_top z1=z_at_w(i,k,j)-blade_l_point-z_at_w(i,1,j) z2=z_at_w(i,k+1,j)-blade_l_point-z_at_w(i,1,j) IF(z1 .LT. 0.) z1=0.0 IF(z2 .GT. diameter(kt)) z2=diameter(kt) CALL turbine_area(z1,z2,diameter(kt),wfdensity,tarea) speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.) if (power1.eq.0.or.power2.eq.0) then ec=1. else ec=power1/power2 endif qke(i,k,j) = qke(i,k,j)+speed**3.*tarea*tkecof*dt/dz(i,k,j)*ec du(i,k,j) = du(i,k,j)-.5*u(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec dv(i,k,j) = dv(i,k,j)-.5*v(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec ENDDO ENDIF ENDIF endif ENDIF ENDDO END SUBROUTINE dragforce SUBROUTINE turbine_area(z1,z2,tdiameter,wfdensity,tarea) REAL, INTENT(IN) ::tdiameter,wfdensity REAL, INTENT(INOUT) ::z1,z2 REAL, INTENT(OUT):: tarea REAL r,zc1,zc2 r=tdiameter/2. z1=r-z1 z2=r-z2 zc1=abs(z1) zc2=abs(z2) IF(z1 .GT. 0. .AND. z2 .GT. 0.) THEN tarea=zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r)- & (zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r)) ELSE IF(z1 .LT. 0. .AND. z2 .LT. 0.) THEN tarea=zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r)- & (zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r)) ELSE tarea=zc2*sqrt(r*r-zc2*zc2)+r*r*asin(zc2/r)+ & zc1*sqrt(r*r-zc1*zc1)+r*r*asin(zc1/r) ENDIF tarea=tarea*wfdensity END SUBROUTINE turbine_area SUBROUTINE dragcof(tkecof,powcof,thrcof,speed,cispeed,cospeed, & tpower,tdiameter,stdthrcoef,stdthrcoef2,nkind) REAL, INTENT(IN):: speed, cispeed, cospeed, tpower,tdiameter,stdthrcoef,stdthrcoef2 REAL, INTENT(OUT):: tkecof,powcof,thrcof REAL :: power,area,mspeed,hspeed INTEGER :: nkind,k,nu,nb LOGICAL :: vfound REAL :: fac1,fac2 area=piconst/4.*tdiameter**2. vfound=.false. DO k=1,maxvals2 IF(.NOT. vfound) THEN IF(turbws(nkind,k).GT.speed) THEN nu=k nb=k-1 vfound=.true. ENDIF ENDIF ENDDO IF (speed .LE. cispeed) THEN thrcof = stdthrcoef ELSE IF (speed .GE. cospeed) THEN thrcof = stdthrcoef2 ELSE thrcof = turbtc(nkind,nb)+(turbtc(nkind,nu)-turbtc(nkind,nb))/(turbws(nkind,nu)-turbws(nkind,nb))*(speed-turbws(nkind,nb)) ENDIF ENDIF IF(speed .LE. cispeed .OR. speed .GE. cospeed) THEN power=0. powcof=0. ELSE fac1=1000./(0.5*1.23*turbws(nkind,nb)**3.*area) fac2=1000./(0.5*1.23*turbws(nkind,nu)**3.*area) power = turbpw(nkind,nb)+(turbpw(nkind,nu)-turbpw(nkind,nb))/(turbws(nkind,nu)-turbws(nkind,nb)) & *(speed-turbws(nkind,nb)) powcof = turbpw(nkind,nb)*fac1+(turbpw(nkind,nu)*fac2-turbpw(nkind,nb)*fac1)/(turbws(nkind,nu)-turbws(nkind,nb)) & *(speed-turbws(nkind,nb)) ENDIF tkecof=thrcof-powcof IF(tkecof .LT. 0.) tkecof=0. END SUBROUTINE dragcof SUBROUTINE init_module_wind_fitch(id,config_flags,xlong,xlat,windfarm_initialized,& ims,ime,jms,jme,its,ite,jts,jte,ids,ide,jds,jde) IMPLICIT NONE integer ims,ime,jms,jme,ids,ide,jds,jde integer its,ite,jts,jte REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: xlong,xlat TYPE (grid_config_rec_type) :: config_flags TYPE (PROJ_INFO) :: ts_proj logical :: windfarm_initialized CHARACTER*256 num,input,message_wind real lat,lon,ts_rx,ts_ry REAL :: known_lat, known_lon INTEGER i,j,nval,k,id LOGICAL, EXTERNAL :: wrf_dm_on_monitor IF ( wrf_dm_on_monitor() ) THEN if ( config_flags%windfarm_ij .eq. 1 ) then open(70,file='windturbines-ij.txt',form='formatted',status='old') else open(70,file='windturbines.txt',form='formatted',status='old') end if nt=0 10 read(70,*,end=100) nt=nt+1 goto 10 100 continue rewind (70) END IF CALL wrf_dm_bcast_integer(nt,1) if (.not. windfarm_initialized) then allocate (nkind(nt),ival(nt,max_domains),jval(nt,max_domains)) allocate (hubheight(nt),stc(nt),stc2(nt),cutin(nt),cutout(nt),diameter(nt),npower(nt)) ival=-9999 jval=-9999 windfarm_initialized=.true. endif IF ( wrf_dm_on_monitor() ) THEN do k=1,nt if ( config_flags%windfarm_ij .eq. 1 ) then read(70,*) ival(k,id), jval(k,id), nkind(k) write(message_wind,*)'WINDFARM Turbine #',k,': I, J = ',ival(k,id), jval(k,id),'; Type = ',nkind(k) CALL wrf_message(message_wind) else read(70,*)lat,lon,nkind(k) write(message_wind,*)'WINDFARM Turbine #',k,': Lat, lon = ',lat,lon,'; Type = ',nkind(k) CALL wrf_message(message_wind) CALL map_init(ts_proj) known_lat = xlat(its,jts) known_lon = xlong(its,jts) IF (config_flags%map_proj == PROJ_MERC) THEN CALL map_set(PROJ_MERC, ts_proj, & truelat1 = config_flags%truelat1, & lat1 = known_lat, & lon1 = known_lon, & knowni = REAL(its), & knownj = REAL(jts), & dx = config_flags%dx) ELSE IF (config_flags%map_proj == PROJ_LC) THEN CALL map_set(PROJ_LC, ts_proj, & truelat1 = config_flags%truelat1, & truelat2 = config_flags%truelat2, & stdlon = config_flags%stand_lon, & lat1 = known_lat, & lon1 = known_lon, & knowni = REAL(its), & knownj = REAL(jts), & dx = config_flags%dx) ELSE IF (config_flags%map_proj == PROJ_PS) THEN CALL map_set(PROJ_PS, ts_proj, & truelat1 = config_flags%truelat1, & stdlon = config_flags%stand_lon, & lat1 = known_lat, & lon1 = known_lon, & knowni = REAL(its), & knownj = REAL(jts), & dx = config_flags%dx) END IF CALL latlon_to_ij(ts_proj, lat, lon, ts_rx, ts_ry) ival(k,id)=nint(ts_rx) jval(k,id)=nint(ts_ry) if (ival(k,id).lt.ids.and.ival(k,id).gt.ide) then ival(k,id)=-9999 jval(k,id)=-9999 endif write(message_wind,*)'WINDFARM Turbine #',k,': Lat, lon = ',lat,lon, & ', (i,j) = (',ival(k,id),',',jval(k,id),'); Type = ',nkind(k) CALL wrf_debug(0,message_wind) end if enddo close(70) turbws=0. turbtc=0. turbpw=0. DO i=1,nt write(num,*) nkind(i) num=adjustl(num) input="wind-turbine-"//trim(num)//".tbl" OPEN(file=TRIM(input),unit=19,FORM='FORMATTED',STATUS='OLD') READ (19,*,ERR=132)nval READ(19,*,ERR=132)hubheight(i),diameter(i),stc(i),npower(i) DO k=1,nval READ(19,*,ERR=132)turbws(nkind(i),k),turbtc(nkind(i),k),turbpw(nkind(i),k) ENDDO cutin(i) = turbws(nkind(i),1) cutout(i) = turbws(nkind(i),nval) stc2(i) = turbtc(nkind(i),nval) close (19) ENDDO 132 continue endif CALL wrf_dm_bcast_integer(ival,nt*max_domains) CALL wrf_dm_bcast_integer(jval,nt*max_domains) CALL wrf_dm_bcast_real(hubheight,nt) CALL wrf_dm_bcast_real(diameter,nt) CALL wrf_dm_bcast_real(stc,nt) CALL wrf_dm_bcast_real(npower,nt) CALL wrf_dm_bcast_real(cutin,nt) CALL wrf_dm_bcast_real(cutout,nt) CALL wrf_dm_bcast_integer(nkind,nt) CALL wrf_dm_bcast_real(turbws,maxvals*maxvals2) CALL wrf_dm_bcast_real(turbtc,maxvals*maxvals2) CALL wrf_dm_bcast_real(turbpw,maxvals*maxvals2) END SUBROUTINE init_module_wind_fitch END MODULE module_wind_fitch