Wind log interpolation about wrf-sfire HOT 14 CLOSED

    loop_k: do k=kds,kdmax         ! search for layer k such that ht(k-1)<=wh<ht(k), ht(0)=z0
      
        ! interpolate height from atmospheric cell midpoints
        ht=interpolate_h(its-1,ite+1,kds,kde,jts-1,jts+1,icm,k,jcm,wicm,wjcm,z)
        ! print *,'i=',i,' j=',j,'k=',k,' ht=',ht

        if(.not. ht < wh) exit loop_k ! found layer k this point is in
        ht_last = ht
      enddo loop_k
    
      if(k .gt. kdmax) then
        goto 91  ! run out of vertical levels, this must be wrong
      endif

      kmin=min(k,kmin)
      kmax=max(k,kmax)
      ! found layer k, ht_last < wh <= ht 
      logz0 = log(z0)
      logwh= log(wh)
      loght_last = log(ht_last)
      loght = log(ht)

$ diff -wBb  interpolate_wind2fire_height interpolate_wind2can_height
1c1
< subroutine interpolate_wind2fire_height(id,       & ! to identify debugging prints and files if needed
---
> subroutine interpolate_wind2can_height(id,      & ! to identify debugging prints and files if needed
13,14c13,14
<     fz0,fwh,                                      & ! input fire arrays
<     uf,vf)                                          ! output fire arrays
---
>     wz0,can_top,                                  & ! input fire arrays
>     cuf,cvf)                                          ! output fire arrays
33,34c37,39
<     fz0,                                          & ! roughness height
<     fwh                                             ! height to read the wind from
---
>     wz0,                                          & ! WRF roughness height
>     can_top                                         ! height to read the wind from
> 
36,37c41,42
<     uf,                                           & ! atm wind at fire_wind_height, diagnostics
<     vf                                              ! 
---
>     cuf,                                           & ! atm wind at canopy top height, diagnostics
>     cvf                                              ! 
46a52,57
> 
>     !This script is nearly identical to the subroutine 'interpolate_wind2fire_height', except 
>     !that it has been repurposed to interpolate winds at the canopy top. A seperate script was
>     !necessary to prevent key variables from being overwritten since this is needed 
>     !for the massman parameterization. DVM 6/26/2019
> 
49c60
<   if(.not. have_wind_log_interpolation) call crash('interpolate_wind2fire_height: wind_log_interpolation must be set up first')
---
>     if(.not. have_wind_log_interpolation) call crash('interpolate_wind2can_height: wind_log_interpolation must be set up first')
97,98c106,107
<        uf(i,j)=0.
<        vf(i,j)=0.
---
>         cuf(i,j)=0.
>         cvf(i,j)=0.
113,115c122,123
<       z0 = fz0(i,j)                    ! roughness length 
<       wh = fwh(i,j)                    ! wind height
< 
---
>         z0 = wz0(i,j)                    ! roughness length 
>         wh = can_top(i,j)                ! wind height
168a179,188
>           !If wind height = 0, we are not dealing with a fuel type with trees
>           !so just set the canopy height winds cuf and cvf as 0.
> 
>           if(wh == 0) then 
>             cuf(i,j) = 0.
>             cvf(i,j) = 0.
>             cycle
>           end if 
>         
> 
170,171c190,191
<       uf(i,j)= uk1 + (uk - uk1) * ( logwh - loght_last) / (loght - loght_last)
<       vf(i,j)= vk1 + (vk - vk1) * ( logwh - loght_last) / (loght - loght_last)
---
>           cuf(i,j)= uk1 + (uk - uk1) * ( logwh - loght_last) / (loght - loght_last)
>           cvf(i,j)= vk1 + (vk - vk1) * ( logwh - loght_last) / (loght - loght_last)
177,179c197,198
<       uf(i,j) = 0.
<       vf(i,j) = 0.
< 
---
>           cuf(i,j) = 0.
>           cvf(i,j) = 0.
185c204
<   ! print *,'interpolate_wind2fire_height complete, id=',id
---
>     ! print *,'interpolate_wind2can_height complete, id=',id
194c213
<   91 call crash('interpolate_wind2fire_height: fire wind height too large, increase kdmax or atm height')
---
>     91 call crash('interpolate_wind2can_height: canopy wind height too large, increase kdmax or atm height')
197c216
<   write(msg,*)'fz0(',i,j,')=',fz0(i,j),'fwh(',i,j,')=',fwh(i,j)
---
>     write(msg,*)'wz0(',i,j,')=',wz0(i,j),'can_top(',i,j,')=',can_top(i,j)
200c219
<   call crash('interpolate_wind2fire_height: must have fire wind height > roughness height > 0')
---
>     call crash('interpolate_wind2can_height: must have canopy wind height > roughness height > 0')
260c279
< end subroutine interpolate_wind2fire_height
---
> end subroutine interpolate_wind2can_height
$

+           if(wh == 0) then 
+             cuf(i,j) = 0.
+             cvf(i,j) = 0.
+             cycle
+           end if 

p0=nc2struct('em_rxcadre2_orig/wrfout_d05_2012-11-11_12:00:00',{'UF','VF','U','V'},{}) 
p2=nc2struct('em_rxcadre2_devel/wrfout_d05_2012-11-11_12:00:00',{'UF','VF','U','V'},{})
t=1:139; big(p0.uf(:,:,t)-p2.uf(:,:,t))
t=1:139; big(p0.vf(:,:,t)-p2.vf(:,:,t))
!strings em_rxcadre2_orig/wrf.exe | grep commit
commit 2d744a6021ca799d812c9df1ac46c34ed46fbeea
!strings em_rxcadre2_devel/wrf.exe | grep commit
commit 85247c7f42f325e9cecd515bf23a1e57c7ea2960

long =
  -86.7300
lat =
   30.5300
names =
  1×6 cell array
    {'layer1'}    {'layer2'}    {'layer3'}    {'layer4'}    {'can_top'}    {'fwh'}
h =
    2.6581    7.9751   15.9533   26.5871   22.5000    0.6070
u =
   -3.2616   -4.2591   -5.0076   -5.5361   -3.9965   -0.5363
v =
   -0.0288   -0.0118    0.0246    0.0922   -0.0101   -0.0014

long =
  -86.7300
lat =
   30.5300
names =
  1x7 cell array
    {'layer1'}    {'layer2'}    {'layer3'}    {'layer4'}    {'can_top'}    {'uf/vf'}    {'uah/vah'}
h =
    2.6644    7.9955   16.0002   26.6716         0         0
u =
   -2.7812   -3.5387   -4.0993   -4.4913         0   -4.3654   -4.3611
v =
    0.5209    0.6878    0.8315    0.9593         0    0.9152    0.9153

names =
  1x7 cell array
    {'layer1'}    {'layer2'}    {'layer3'}    {'layer4'}    {'can_top'}    {'uf/vf'}    {'uah/vah'}
h =
    2.6670    8.0030   16.0129   26.6889   22.5000    0.6070
u =
   -2.8603   -3.4974   -3.9268   -4.2753   -4.1607   -0.5583         0
v =
    0.6379    0.7948    0.9083    1.0015    0.9734    0.1306         0