Fitch bug fix: axial induction factor

LICENSE.txt

@@ -1,3 +1,21 @@
-WRF was developed at the National Center for Atmospheric Research (NCAR) which is operated by the University Corporation for Atmospheric Research (UCAR). NCAR and UCAR make no proprietary claims, either statutory or otherwise, to this version and release of WRF and consider WRF to be in the public domain for use by any person or entity for any purpose without any fee or charge. UCAR requests that any WRF user include this notice on any partial or full copies of WRF. WRF is provided on an "AS IS" basis and any warranties, either express or implied, including but not limited to implied warranties of non-infringement, originality, merchantability and fitness for a particular purpose, are disclaimed. In no event shall UCAR be liable for any damages, whatsoever, whether direct, indirect, consequential or special, that arise out of or in connection with the access, use or performance of WRF, including infringement actions.
+MIT License
 
-WRF® is a registered trademark of the University Corporation for Atmospheric Research (UCAR).
+Copyright (c) 2022 FraunhoferIWES
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

Registry/Registry.EM_COMMON

@@ -411,7 +411,7 @@ state   real    bldtacttime    -        -         -         -     r         "bld
 state   real    cudtacttime    -        -         -         -     r         "cudtacttime"              "CUDTACTTIME"         "CPS   ACTIVATION TIME in s"
 state   real    ltngacttime    -        -         -         -     r         "ltngacttime"              "LTNGACTTIME"         "LTNG  ACTIVATION TIME in s"
 state   real    power          ij       misc      1         -     irh       "Power"                   "Power production"         "W"
-
+state   real    thrust         ij       misc      1         -     irh       "Thrust"                   "Thrust force"         "N"
 
 # State for derived time quantities.  
 state   integer itimestep      -        -          -         -     rh         "itimestep"             ""         ""      
@@ -3389,7 +3389,7 @@ package   wrfhydro       wrf_hydro==1                -            state:SOLDRAIN
 
 #WRF Windfarm
 package   no_windfarm    windfarm_opt==0             -            -
-package   fitchscheme    windfarm_opt==1             -            state:power
+package   fitchscheme    windfarm_opt==1             -            state:power,thrust
 # Yulong add for WLM
 package   mavscheme      windfarm_opt==2             -            state:power
 

dyn_em/module_first_rk_step_part1.F

@@ -1113,7 +1113,7 @@ SUBROUTINE first_rk_step_part1 (   grid , config_flags              &
       CALL pbl_driver(                                                              &
      &         AKHS=grid%akhs          ,AKMS=grid%akms                              &
      &        ,BL_PBL_PHYSICS=config_flags%bl_pbl_physics                           &
-     &        ,WINDFARM_OPT=config_flags%windfarm_opt,power=grid%power              &
+     &        ,WINDFARM_OPT=config_flags%windfarm_opt,power=grid%power,thrust=grid%thrust              &
      &        ,windfarm_wake_model=config_flags%windfarm_wake_model          & ! Yulong add for WLM
      &        ,windfarm_overlap_method=config_flags%windfarm_overlap_method  & ! Yulong add for WLM
      &        ,BLDT=grid%bldt, CURR_SECS=curr_secs, ADAPT_STEP_FLAG=adapt_step_flag &

phys/module_pbl_driver.F

@@ -27,7 +27,7 @@ SUBROUTINE pbl_driver(                                          &
                  ,stepbl,warm_rain                                 &
                  ,kpbl,mixht,ct,lh,snow,xice                       &
                  ,znu, znw, mut, p_top                             &
-                 ,ctopo,ctopo2,windfarm_opt,power                  &
+                 ,ctopo,ctopo2,windfarm_opt,power,thrust           &
                  ,windfarm_wake_model, windfarm_overlap_method     &
                  ,ysu_topdown_pblmix                               &
                  ,shinhong_tke_diag                                &
@@ -542,6 +542,7 @@ SUBROUTINE pbl_driver(                                          &
                                                             VOCE, &
                                                               T2, &
                                                            POWER, &
+                                                           THRUST, &
                                                             WSPD
    REAL,       DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(INOUT)    ::                     HPBL_HOLD
@@ -2057,6 +2058,7 @@ SUBROUTINE pbl_driver(                                          &
                      &,QKE=qke                                           &
                      &,DU=rublten,DV=rvblten                             &
                      &,WINDFARM_OPT=windfarm_opt,POWER=power             &
+                     &,THRUST=thrust                                     &
                      &,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde   &
                      &,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme   &
                      &,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte   &

phys/module_wind_fitch.F

@@ -69,7 +69,7 @@ SUBROUTINE  dragforce(                      &
        &,z_at_w,u,v                 &
        &,dx,dz,dt,qke                            &
        &,du,dv                                   &
-       &,windfarm_opt,power                      &
+       &,windfarm_opt,power,thrust               &
        &,ids,ide,jds,jde,kds,kde                 &
        &,ims,ime,jms,jme,kms,kme                 &
        &,its,ite,jts,jte,kts,kte                 &
@@ -84,12 +84,13 @@ SUBROUTINE  dragforce(                      &
   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, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: power,thrust
 !
 ! Local
 !
   REAL     blade_l_point,blade_u_point,zheightl,zheightu,z1,z2,tarea
   REAL     speed,tkecof,powcof,thrcof,wfdensity
+  REAL     a_ind ! BAMS
   INTEGER  itf,jtf,ktf
   INTEGER  i,j,k,n
   INTEGER  k_turbine_bot, k_turbine_top
@@ -99,15 +100,21 @@ SUBROUTINE  dragforce(                      &
 ! ... PAJ: more variables ...
 !
   REAL :: speedhub,speed1,speed2
-  real :: power1,power2,area,ec
+  real :: power1,power2,area,ec,thrust1
   INTEGER :: kbot,ktop,kt
+  INTEGER :: ii ! BAMS: running index
+  INTEGER, DIMENSION(:), ALLOCATABLE :: idx ! BAMS: temporary array of indices
+  INTEGER :: nturb ! BAMS: number of turbines in grid cell
+
+  CHARACTER*256 num,input,message_wind
 
   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.
+    thrust=0.
 
     DO kt = 1,nt  
       IF ( windfarm_opt .eq. 1 ) THEN
@@ -118,6 +125,11 @@ SUBROUTINE  dragforce(                      &
         k_turbine_top=-1     !top level
         i = ival(kt,id)
         j = jval(kt,id)
+
+        ! BAMS: Number of turbines in current grid cell
+        idx = PACK([(ii, ii=1,nt)], ival(:,id) == i) ! indices of values equal to i in ival
+        idx = PACK([(ii, ii=1,SIZE(jval(idx,id)))], jval(idx,id) == j) ! indices of values equal j in jval(idx)
+        nturb = SIZE(idx) ! number of turbines in grid cell (i,j)
 !
          if (i.ne.-9999.and.j.ne.-9999) then
         IF (( its .LE. i .AND. i .LE. itf ) .AND. &
@@ -169,27 +181,58 @@ SUBROUTINE  dragforce(                      &
             ENDIF
           ENDDO
 !
+! ... BAMS: first iteration with induction first estimate
           speed1=0.
           speed2=0.
+          a_ind=0.01
           if (ktop.eq.1) then
            speedhub=sqrt(u(i,1,j)**2.+v(i,1,j)**2.)*hubheight(kt)/z1
+           speedhub=speedhub/((1-a_ind)**nturb) ! BAMS
           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)
+           speedhub=speedhub/((1-a_ind)**nturb) ! BAMS
           endif
 !
 ! ... 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))
+          CALL dragcof(tkecof,powcof,thrcof,               &
+                        speedhub,cutin(kt),cutout(kt),   &
+                        npower(kt),diameter(kt),stc(kt),stc2(kt),nkind(kt))
+!
+! ... BAMS: second iteration with induction
+          area = piconst/4.*diameter(kt)**2.       ! area swept by turbine blades; BAMS
+          a_ind = 0.5 * (1 - sqrt(1 - thrcof))
+          a_ind = a_ind * area / ( diameter(kt) * dx * &
+                  MIN(ABS(1/COS(ATAN( ((v(i,ktop,j) + v(i,kbot,j))/2) / ((u(i,ktop,j) + u(i,kbot,j))/2) ))), &
+                      ABS(1/SIN(ATAN( ((v(i,ktop,j) + v(i,kbot,j))/2) / ((u(i,ktop,j) + u(i,kbot,j))/2) )))) )     ! Corrected induction factor BAMS
+
+!          write(message_wind,*)'Turbine #',kt,': ',thrcof,0.5 * (1 - sqrt(1 - thrcof)),a_ind
+!          CALL wrf_message(message_wind)
+
+          speed1=0.
+          speed2=0.
+          if (ktop.eq.1) then
+          speedhub=sqrt(u(i,1,j)**2.+v(i,1,j)**2.)*hubheight(kt)/z1
+          speedhub=speedhub/((1-a_ind)**nturb) ! BAMS
+          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)
+          speedhub=speedhub/((1-a_ind)**nturb) ! BAMS
+          endif
+
+         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)
+          thrust1=0.5*1.23*speedhub**2.*area*thrcof
+          thrust(i,j)=thrust1+thrust(i,j)
           power2=0.
 !
             DO k=k_turbine_bot,k_turbine_top ! loop over turbine blade levels
@@ -199,7 +242,7 @@ SUBROUTINE  dragforce(                      &
               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)
 !
-              speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.)
+              speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.) / ((1-a_ind)**nturb) ! BAMS
               power2=power2+0.5*powcof*1.23*(speed**3.)*tarea/wfdensity
             ENDDO
 !
@@ -213,7 +256,7 @@ SUBROUTINE  dragforce(                      &
 !
               CALL turbine_area(z1,z2,diameter(kt),wfdensity,tarea)
 !
-              speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.)
+              speed=sqrt(u(i,k,j)**2.+v(i,k,j)**2.)  / ((1-a_ind)**nturb) ! BAMS
 !`
 ! ... PAJ: normalization introduced to conserve energy ...
 !
@@ -226,9 +269,9 @@ SUBROUTINE  dragforce(                      &
               ! output TKE
               qke(i,k,j) = qke(i,k,j)+speed**3.*tarea*tkecof*dt/dz(i,k,j)*ec
               ! output u tendency
-              du(i,k,j) = du(i,k,j)-.5*u(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec
+              du(i,k,j) = du(i,k,j)-.5*(u(i,k,j)/((1-a_ind)**nturb))*thrcof*speed*tarea/dz(i,k,j)*ec ! BAMS
               ! output v tendency
-              dv(i,k,j) = dv(i,k,j)-.5*v(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec
+              dv(i,k,j) = dv(i,k,j)-.5*(v(i,k,j)/((1-a_ind)**nturb))*thrcof*speed*tarea/dz(i,k,j)*ec ! BAMS
             ENDDO
           ENDIF
         ENDIF
@@ -310,6 +353,7 @@ SUBROUTINE dragcof(tkecof,powcof,thrcof,speed,cispeed,cospeed, &
      thrcof = turbtc(nkind,nb)+(turbtc(nkind,nu)-turbtc(nkind,nb))/(turbws(nkind,nu)-turbws(nkind,nb))*(speed-turbws(nkind,nb))
     ENDIF
   ENDIF
+  thrcof = MAX(MIN(thrcof,0.9999),0.0001) ! BAMS
 !
 ! ... power coeficient ...
 !