diff --git a/LICENSE.txt b/LICENSE.txt old mode 100644 new mode 100755 index e0c9342282..c524bbbba6 --- a/LICENSE.txt +++ b/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. diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 4d133c9bfa..a6266bc4e8 100644 --- a/Registry/Registry.EM_COMMON +++ b/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 diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index 9a56bb2414..39ea81d6be 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/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 & diff --git a/phys/module_pbl_driver.F b/phys/module_pbl_driver.F index 24c2174556..78907f4482 100644 --- a/phys/module_pbl_driver.F +++ b/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 & diff --git a/phys/module_wind_fitch.F b/phys/module_wind_fitch.F index 3adaa8274c..743fc9fa92 100644 --- a/phys/module_wind_fitch.F +++ b/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,8 +100,13 @@ 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) @@ -108,6 +114,7 @@ SUBROUTINE dragforce( & 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 ... !