!WRF:MODEL_LAYER:PHYSICS

MODULE module_wind_fitch
!
!Represents kinetic energy extracted by wind turbines and turbulence
! (TKE) they produce at model levels within the rotor area.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! NOTICE
!! The following paper should be cited whenever presenting results using this scheme
!! (using either the original version or any modified versions of the scheme):
!! Fitch, A. C. et al. 2012: Local and Mesoscale Impacts of Wind Farms as Parameterized in a
!! Mesoscale NWP Model. Monthly Weather Review, doi:http://dx.doi.org/10.1175/MWR-D-11-00352.1
!!
!! Anna C. Fitch, National Center for Atmospheric Research (formerly University of Bergen)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! History of changes:
!
! WRFV3.5.1: 
! WRFV3.6:   Modified by Pedro A. Jimenez to include:
!             - Initialize the wind turbines in this module.
!             - Introduce z_at_walls to avoid instabilities due to neglecting
!                the perturbation of the geopotential height.
!             - User friendly interface to introduce the technical characteritics of
!                the wind turbines.
!             - Only uses one set of turbine coefficients using the wind speed at hub height
!             - Two standing coefficients.
!             - Calculates the power produced by the wind turbines.
!
! References:
!
! Fitch, A. C. et al. 2012: Local and Mesoscale Impacts of Wind Farms as Parameterized in a
!    Mesoscale NWP Model. Monthly Weather Review, doi:http://dx.doi.org/10.1175/MWR-D-11-00352.1
! Fitch, A. C. et al. 2013: Mesoscale Influences of Wind Farms Throughout a Diurnal Cycle.
!    Monthly Weather Review, doi:http://dx.doi.org/10.1175/MWR-D-12-00185.1
! Fitch, A. C. et al. 2013: Parameterization of Wind Farms in Climate Models.
!    Journal of Climate, doi:http://dx.doi.org/10.1175/JCLI-D-12-00376.1
! Jimenez, P.A., J. Navarro, A.M. Palomares and J. Dudhia:  Mesoscale modeling of offshore wind turbines
!    wakes at the wind farm resolving scale: a composite-based analysis with the WRF model over Horns Rev. 
!    Wind Energy, (In Press.).
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
  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
!
! Local
!
  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
!
! ... PAJ: more variables ...
!
  REAL :: speedhub,speed1,speed2
  real :: power1,power2,area,ec
  INTEGER :: kbot,ktop,kt
  REAL, DIMENSION(:), ALLOCATABLE :: kspeed
  integer  tt !网格间距=tt*风机间距（tt=times）
  integer sn !风机前遮挡风机的个数 （sn=superposition number）
  integer yy !风机编号除以倍数关系的余数
  integer ll !循环变量
  real aa,bb2,bb3 !尾流效应修正系数
  real diameterr !风机直径
  real rr
  real pp !尾流膨胀系数
  real dd !风机间距
  real xx,nn
  real uabs,vabs,vel_z
  real theta

  itf=MIN0(ite,ide-1)
  jtf=MIN0(jte,jde-1)
  ktf=MIN0(kte,kde-1)

    wfdensity = 1.0/(dx*dx)   !  per turbine, so numerator is 1
    power=0.
    allocate(kspeed(nt))    !supplement code
    DO kt = 1,nt 
!supplement code
    open(file='parainput.txt',unit=615,form='formatted',status='old')
    read(615,*)tt,diameterr,pp,dd
    close(615)
    yy=mod(kt,tt)
    if(yy==0)then
       sn=tt-1
    else
       sn=yy-1
    end if

    if(sn==0)then 
    aa = 1
    else 
    aa = 0
    rr = diameterr/2
    xx = (real(kt)-1.0)/real(tt)
    nn = int(xx)*tt + 1
    Do LL = sn , 1 , -1
    CALL dragcof(tkecof,powcof,thrcof,  & 
                 kspeed(nn),cutin(kt),cutout(kt),  &
                 npower(kt),diameter(kt),stc(kt),stc2(kt),nkind(kt))
        aa = aa + ((1-sqrt(1-thrcof))*(rr/(rr+LL*dd*pp))**2)**2
        nn = nn + 1
    end Do  
    aa=1-sqrt(aa)
    end if

!supplement code
      IF ( windfarm_opt .eq. 1 ) THEN
!
! vertical layers cut by turbine blades
!
        k_turbine_bot=0      !bottom level
        k_turbine_top=-1     !top level
        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. ! height of lower blade tip above ground (m)
          blade_u_point=hubheight(kt)+diameter(kt)/2. ! height of upper blade tip above ground (m)
!
          kfound = .false.
          zheightl=0.0
          ! find vertical levels cut by turbine blades
          DO k=kts,ktf
            IF(.NOT. kfound) THEN
              zheightu = zheightl + dz(i,k,j) ! increment height

              IF(blade_l_point .GE. zheightl .AND. blade_l_point .LE. zheightu) THEN
                k_turbine_bot=k ! lower blade tip cuts this level
              ENDIF

              IF(blade_u_point .GE. zheightl .AND. blade_u_point .LE. zheightu) THEN
                k_turbine_top=k ! upper blade tip cuts this level
                kfound = .TRUE.
              ENDIF

              zheightl = zheightu
            ENDIF
          ENDDO
          IF ( kfound ) THEN
!
! ... PAJ: Changes introduced to compute only one set of turbine coefficients ...
!          First computes the wind speed at the hub height.
!
          kfound = .false.
          zheightl=0.
          ! find vertical levels (half levels) within the hub height
          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
!velocity adjust        
	     uabs=abs(u(i,1,j))
	     vabs=abs(v(i,1,j))
	     if (uabs>vabs) then
		   theta = atan(vabs/uabs)*360.0/2.0/3.14159
	     else
		   theta = atan(uabs/vabs)*360.0/2.0/3.14159
	     endif
	     vel_z = sqrt(uabs**2+vabs**2)
	     bb3 = 0.4992-0.03641*theta + 0.5409*vel_z + 0.01042*theta**2-0.003057*theta*vel_z-0.2048*vel_z**2-0.0007402*theta**3-0.0006517*theta**2*vel_z&
	     + 0.00225*theta*vel_z**2+0.03835*vel_z**3+2.053e-05*theta**4+5.856e-05*theta**3*vel_z-4.497e-05*theta**2*vel_z**2-0.0002557*theta*vel_z**3-0.003962*vel_z**4& 
	     - 1.107e-07*theta**5-2.032e-06*theta**4*vel_z-5.569e-07*theta**3*vel_z**2 + 6.971e-06*theta**2*vel_z**3 + 9.526e-06*theta*vel_z**4 + 0.00023*vel_z**5&
	     - 4.582e-09*theta**6 + 3.564e-08*theta**5*vel_z + 4.596e-09*theta**4*vel_z**2 + 1.594e-08*theta**3*vel_z**3-3.679e-07*theta**2*vel_z**4 + 1.925e-08*theta*vel_z**5& 
	     - 7.023e-06*vel_z**6+5.931e-11*theta**7-2.426e-10*theta**6*vel_z-3.01e-11*theta**5*vel_z**2-6.279e-12*theta**4*vel_z**3-2.726e-10*theta**3*vel_z**4& 
	     + 6.696e-09*theta**2*vel_z**5-5.215e-09*theta*vel_z**6 + 8.779e-08*vel_z**7
	     
           speedhub=sqrt(u(i,1,j)**2.+v(i,1,j)**2.)*hubheight(kt)/z1*aa*bb3
          else
 !velocity adjust          
	     uabs=abs(u(i,kbot,j))
	     vabs=abs(v(i,kbot,j))
	     if (uabs>vabs) then
		   theta = atan(vabs/uabs)*360.0/2.0/3.14159
	     else
	 	   theta = atan(uabs/vabs)*360.0/2.0/3.14159
	     endif
	     vel_z = sqrt(uabs**2+vabs**2)
	     bb3 = 0.4992-0.03641*theta + 0.5409*vel_z + 0.01042*theta**2-0.003057*theta*vel_z-0.2048*vel_z**2-0.0007402*theta**3-0.0006517*theta**2*vel_z&
	     + 0.00225*theta*vel_z**2+0.03835*vel_z**3+2.053e-05*theta**4+5.856e-05*theta**3*vel_z-4.497e-05*theta**2*vel_z**2-0.0002557*theta*vel_z**3-0.003962*vel_z**4& 
	     - 1.107e-07*theta**5-2.032e-06*theta**4*vel_z-5.569e-07*theta**3*vel_z**2 + 6.971e-06*theta**2*vel_z**3 + 9.526e-06*theta*vel_z**4 + 0.00023*vel_z**5&
	     - 4.582e-09*theta**6 + 3.564e-08*theta**5*vel_z + 4.596e-09*theta**4*vel_z**2 + 1.594e-08*theta**3*vel_z**3-3.679e-07*theta**2*vel_z**4 + 1.925e-08*theta*vel_z**5& 
	     - 7.023e-06*vel_z**6+5.931e-11*theta**7-2.426e-10*theta**6*vel_z-3.01e-11*theta**5*vel_z**2-6.279e-12*theta**4*vel_z**3-2.726e-10*theta**3*vel_z**4& 
	     + 6.696e-09*theta**2*vel_z**5-5.215e-09*theta*vel_z**6 + 8.779e-08*vel_z**7
                 
           speed1=sqrt(u(i,kbot,j)**2.+v(i,kbot,j)**2.)*aa*bb3
!velocity adjust           
	     uabs=abs(u(i,ktop,j))
	     vabs=abs(v(i,ktop,j))
	     if (uabs>vabs) then
		   theta = atan(vabs/uabs)*360.0/2.0/3.14159
	     else
	 	   theta = atan(uabs/vabs)*360.0/2.0/3.14159
	     endif
	     vel_z = sqrt(uabs**2+vabs**2)
	     bb3 = 0.4992-0.03641*theta + 0.5409*vel_z + 0.01042*theta**2-0.003057*theta*vel_z-0.2048*vel_z**2-0.0007402*theta**3-0.0006517*theta**2*vel_z&
	     + 0.00225*theta*vel_z**2+0.03835*vel_z**3+2.053e-05*theta**4+5.856e-05*theta**3*vel_z-4.497e-05*theta**2*vel_z**2-0.0002557*theta*vel_z**3-0.003962*vel_z**4& 
	     - 1.107e-07*theta**5-2.032e-06*theta**4*vel_z-5.569e-07*theta**3*vel_z**2 + 6.971e-06*theta**2*vel_z**3 + 9.526e-06*theta*vel_z**4 + 0.00023*vel_z**5&
	     - 4.582e-09*theta**6 + 3.564e-08*theta**5*vel_z + 4.596e-09*theta**4*vel_z**2 + 1.594e-08*theta**3*vel_z**3-3.679e-07*theta**2*vel_z**4 + 1.925e-08*theta*vel_z**5& 
	     - 7.023e-06*vel_z**6+5.931e-11*theta**7-2.426e-10*theta**6*vel_z-3.01e-11*theta**5*vel_z**2-6.279e-12*theta**4*vel_z**3-2.726e-10*theta**3*vel_z**4& 
	     + 6.696e-09*theta**2*vel_z**5-5.215e-09*theta*vel_z**6 + 8.779e-08*vel_z**7
                 
           speed2=sqrt(u(i,ktop,j)**2.+v(i,ktop,j)**2.)*aa*bb3
           speedhub=speed1+((speed2-speed1)/(z2-z1))*(hubheight(kt)-z1)
          endif
!supplement code
          kspeed(kt)=speedhub
!
! ... calculate TKE, power and thrust coeffs
!
              CALL dragcof(tkecof,powcof,thrcof,               &
                           speedhub,cutin(kt),cutout(kt),   &
                           npower(kt),diameter(kt),stc(kt),stc2(kt),nkind(kt))
!
! ... PAJ: Computation of power generated by the wind turbine ...
!
          area=piconst/4.*diameter(kt)**2.          ! area swept by turbine blades
          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 ! loop over turbine blade levels
              z1=z_at_w(i,k,j)-blade_l_point-z_at_w(i,1,j)  ! distance between k level and lower blade tip
              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
              IF(z1 .LT. 0.) z1=0.0 ! k level lower than lower blade tip
              IF(z2 .GT. diameter(kt)) z2=diameter(kt) ! k+1 level higher than turbine upper blade tip
              CALL turbine_area(z1,z2,diameter(kt),wfdensity,tarea)
!
!velocity adjust           
	     uabs=abs(u(i,k,j))
	     vabs=abs(v(i,k,j))
	     if (uabs>vabs) then
		   theta = atan(vabs/uabs)*360.0/2.0/3.14159
	     else  
	 	   theta = atan(uabs/vabs)*360.0/2.0/3.14159
	     endif
	     vel_z = sqrt(uabs**2+vabs**2)
	     bb3 = 0.4992-0.03641*theta + 0.5409*vel_z + 0.01042*theta**2-0.003057*theta*vel_z-0.2048*vel_z**2-0.0007402*theta**3-0.0006517*theta**2*vel_z&
	     + 0.00225*theta*vel_z**2+0.03835*vel_z**3+2.053e-05*theta**4+5.856e-05*theta**3*vel_z-4.497e-05*theta**2*vel_z**2-0.0002557*theta*vel_z**3-0.003962*vel_z**4& 
	     - 1.107e-07*theta**5-2.032e-06*theta**4*vel_z-5.569e-07*theta**3*vel_z**2 + 6.971e-06*theta**2*vel_z**3 + 9.526e-06*theta*vel_z**4 + 0.00023*vel_z**5&
	     - 4.582e-09*theta**6 + 3.564e-08*theta**5*vel_z + 4.596e-09*theta**4*vel_z**2 + 1.594e-08*theta**3*vel_z**3-3.679e-07*theta**2*vel_z**4 + 1.925e-08*theta*vel_z**5& 
	     - 7.023e-06*vel_z**6+5.931e-11*theta**7-2.426e-10*theta**6*vel_z-3.01e-11*theta**5*vel_z**2-6.279e-12*theta**4*vel_z**3-2.726e-10*theta**3*vel_z**4& 
	     + 6.696e-09*theta**2*vel_z**5-5.215e-09*theta*vel_z**6 + 8.779e-08*vel_z**7

              speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.)*aa*bb3
              power2=power2+0.5*powcof*1.23*(speed**3.)*tarea/wfdensity
            ENDDO
!
! ... PAJ: Computes the tendencies of TKE and momentum ...
!
            DO k=k_turbine_bot,k_turbine_top ! loop over turbine blade levels
              z1=z_at_w(i,k,j)-blade_l_point-z_at_w(i,1,j)  ! distance between k lev and lower blade tip
              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
              IF(z1 .LT. 0.) z1=0.0 ! k level lower than lower blade tip
              IF(z2 .GT. diameter(kt)) z2=diameter(kt) ! k+1 level higher than turbine upper blade tip
!
              CALL turbine_area(z1,z2,diameter(kt),wfdensity,tarea)
!
!velocity adjust           
             uabs=abs(u(i,k,j))
	     vabs=abs(v(i,k,j))
	     if (uabs>vabs) then
		   theta = atan(vabs/uabs)*360.0/2.0/3.14159
	     else  
	 	   theta = atan(uabs/vabs)*360.0/2.0/3.14159
	     endif
	     vel_z = sqrt(uabs**2+vabs**2)
	     bb2 = 0.6455-0.01021*theta+0.3905*vel_z+0.009437*theta**2-0.01169*theta*vel_z-0.1554*vel_z**2-0.0007517*theta**3-0.0002167*theta**2*vel_z&
	     + 0.003552*theta*vel_z**2+0.03078*vel_z**3+2.136e-05*theta**4+5.26e-05*theta**3*vel_z-0.0001078*theta**2*vel_z**2-0.0003365*theta*vel_z**3-0.00336*vel_z**4&
	     - 9.719e-08*theta**5-2.091e-06*theta**4*vel_z + 2.055e-07*theta**3*vel_z**2 + 1.17e-05*theta**2*vel_z**3+9.421e-06*theta*vel_z**4+0.000205*vel_z**5&
	     - 5.516e-09*theta**6+3.911e-08*theta**5*vel_z-4.403e-09*theta**4*vel_z**2-3.984e-09*theta**3*vel_z**3-5.582e-07*theta**2*vel_z**4 + 2.103e-07*theta*vel_z**5&
	     - 6.544e-06*vel_z**6 + 6.924e-11*theta**7-2.76e-10*theta**6*vel_z + 2.421e-11*theta**5*vel_z**2 + 4.886e-11*theta**4*vel_z**3 + 3.825e-11*theta**3*vel_z**4& 
	     + 9.626e-09*theta**2*vel_z**5-1.021e-08*theta*vel_z**6+8.501e-08*vel_z**7
                 
	     bb3 = 0.4992-0.03641*theta + 0.5409*vel_z + 0.01042*theta**2-0.003057*theta*vel_z-0.2048*vel_z**2-0.0007402*theta**3-0.0006517*theta**2*vel_z&
	     + 0.00225*theta*vel_z**2+0.03835*vel_z**3+2.053e-05*theta**4+5.856e-05*theta**3*vel_z-4.497e-05*theta**2*vel_z**2-0.0002557*theta*vel_z**3-0.003962*vel_z**4& 
	     - 1.107e-07*theta**5-2.032e-06*theta**4*vel_z-5.569e-07*theta**3*vel_z**2 + 6.971e-06*theta**2*vel_z**3 + 9.526e-06*theta*vel_z**4 + 0.00023*vel_z**5&
	     - 4.582e-09*theta**6 + 3.564e-08*theta**5*vel_z + 4.596e-09*theta**4*vel_z**2 + 1.594e-08*theta**3*vel_z**3-3.679e-07*theta**2*vel_z**4 + 1.925e-08*theta*vel_z**5& 
	     - 7.023e-06*vel_z**6+5.931e-11*theta**7-2.426e-10*theta**6*vel_z-3.01e-11*theta**5*vel_z**2-6.279e-12*theta**4*vel_z**3-2.726e-10*theta**3*vel_z**4& 
	     + 6.696e-09*theta**2*vel_z**5-5.215e-09*theta*vel_z**6 + 8.779e-08*vel_z**7

              speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.)
!`
! ... PAJ: normalization introduced to conserve energy ...
!
              if (power1.eq.0.or.power2.eq.0) then
              ec=1.
              else
              ec=power1/power2
              endif
!
              ! output TKE
              qke(i,k,j) = qke(i,k,j)+(aa*bb3*speed)**3.*tarea*tkecof*dt/dz(i,k,j)*ec
              ! output u tendency
              du(i,k,j) = du(i,k,j)-.5*aa*bb2*u(i,k,j)*thrcof*aa*bb2*speed*tarea/dz(i,k,j)*ec
              ! output v tendency
              dv(i,k,j) = dv(i,k,j)-.5*aa*bb2*v(i,k,j)*thrcof*aa*bb2*speed*tarea/dz(i,k,j)*ec
            ENDDO
          ENDIF
        ENDIF
        endif
      ENDIF
    ENDDO

  END SUBROUTINE dragforce

! This subroutine calculates area of turbine between two vertical levels
! Input variables : 
!            z1 = distance between k level and lower blade tip
!            z2 = distance between k+1 level and lower blade tip
!            wfdensity = wind farm density in m^-2
!     tdiameter = turbine diameter
! Output variable :
!         tarea = area of turbine between two levels * wfdensity
  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.              !r = turbine radius
  z1=r-z1                   !distance of kth level from turbine center 
  z2=r-z2                   !distance of k+1 th level from turbine center
  zc1=abs(z1)
  zc2=abs(z2)
  !turbine area between z1 and 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      !turbine area * wind farm density 

  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
!
! ... PAJ ...
!
   INTEGER :: nkind,k,nu,nb
   LOGICAL :: vfound
   REAL :: fac1,fac2

  area=piconst/4.*tdiameter**2.          ! area swept by turbine blades

      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
!
! ... power coeficient ...
!
  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
!
  ! tke coefficient calculation 

  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
!
! ... PAJ: Opens the file with the location of the wind turbines ...
!
        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
!
! ... PAJ: Counts the turbines ...
!
       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)
!
! ... PAJ: Initializes the configuration of the wind farm(s) ...
!
     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)

      ! Mercator
      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)

      ! Lambert conformal
      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)
!      ! Polar stereographic
      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)
!#if (EM_CORE == 1)
!      ! Cassini (global ARW)
!      ELSE IF (config_flags%map_proj == PROJ_CASSINI) THEN
!         CALL map_set(PROJ_CASSINI, ts_proj,                            &
!                      latinc   = grid%dy*360.0/(2.0*EARTH_RADIUS_M*PI), &
!                      loninc   = grid%dx*360.0/(2.0*EARTH_RADIUS_M*PI), &
!                      lat1     = known_lat,                             &
!                      lon1     = known_lon,                             &
!                      lat0     = config_flags%pole_lat,                 &
!                      lon0     = config_flags%pole_lon,                 &
!                      knowni   = 1.,                                    &
!                      knownj   = 1.,                                    &
!                      stdlon   = config_flags%stand_lon)
!#endif
!
!      ! Rotated latitude-longitude
!      ELSE IF (config_flags%map_proj == PROJ_ROTLL) THEN
!         CALL map_set(PROJ_ROTLL, ts_proj,                      &
!! I have no idea how this should work for NMM nested domains
!                      ixdim    = grid%e_we-1,                   &
!                      jydim    = grid%e_sn-1,                   &
!                      phi      = real(grid%e_sn-2)*grid%dy/2.0, &
!                      lambda   = real(grid%e_we-2)*grid%dx,     &
!                      lat1     = config_flags%cen_lat,          &
!                      lon1     = config_flags%cen_lon,          &
!                      latinc   = grid%dy,                       &
!                      loninc   = grid%dx,                       &
!                      stagger  = HH)
!
      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(73,*) k,id,ival(k,id),jval(k,id)
          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)
!
! ... PAJ: Read the tables with the turbine's characteristics ...
!
         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
