debugging NDWI test problem

wrfv2_fire/phys/module_fr_sfire_phys.F

@@ -275,7 +275,8 @@ subroutine fuel_moisture(                &
                                        nwdi_f      ! inerpolation of nwdi to the fire grid
 integer::i,j,k,n
 integer::ibs,ibe,jbs,jbe
-real::f1,w1,w2,f2
+real::f1,w1,w2,f2,fa,fc
+
 character(len=128)::msg
 
 call check_mesh_2dim(ifts,ifte,jfts,jfte,ifds,ifde,jfds,jfde) ! check if fire tile fits into domain
@@ -332,16 +333,20 @@ subroutine fuel_moisture(                &
         do i=ifts,ifte
             n = nfuel_cat(i,j)
             if(n > 0)then
-                f1=fmc_f(i,j)
-                if(k .eq. kfmc_ndwi)then  ! if live moisture, assimilate ndwi
+                if(k .ne. kfmc_ndwi)then 
+                   fmc_g(i,j)=fmc_g(i,j)+fmc_gw(n,k)*fmc_f(i,j)      ! add to sum over classes
+                else                                                    ! if live moisture, assimilate ndwi
+                   f1=fmc_f(i,j)
                    w1 = fmc_gl_stdev(n)                            
                    w1 = 1./(w1*w1)                                      ! weight of forecast
                    w2 = fmc_gl_ndwi_stdev(n)
                    w2 = 1./(w2*w2)                                      ! weight of update
                    f2 = fmc_gl_ndwi_0(n) + fmc_gl_ndwi_rate(n) * fndwi(i,j)   ! from regression
-                   f1 = (w1*f1 + w2*f2) / (w1 + w2)                     ! updated value
+                   fa = (w1*f1 + w2*f2) / (w1 + w2)                     ! updated value
+                   fc = fmc_gw(n,k)*fa                                  ! times proportion of live fuel
+                   fmc_g(i,j)=fmc_g(i,j)+fc      ! add to sum over classes
+                   ! write(*,*)'NDWI:',i,j,f1,f2,w1,w2,f1,fa,fmc_gw(n,k),fc,fmc_g(i,j)
                 endif
-                fmc_g(i,j)=fmc_g(i,j)+fmc_gw(n,k)*f1      ! add to sum over classes
             endif
         enddo
     enddo

wrfv2_fire/test/em_fire/ndwi/namelist.fire

@@ -61,8 +61,16 @@ no_fuel_cat = 14                   ! extra category for no fuel
           1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
   adjrs = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
           1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
-  fmc_gw01 = 12*0.5 ,
-  fmc_gw05 = 12*0.5 ,
+  fmc_gw01 = 12*0.2 ,
+  fmc_gw02 = 12*0.3 ,
+  fmc_gw03 = 12*0.0 ,
+  fmc_gw04 = 12*0.0 ,
+  fmc_gw05 = 12*0.5 ,      ! fuel is consists of 0.5 live fuel
+  fmc_gl_stdev = 12*0.1,   ! live fuel moisture prior has error +- 0.1
+  ! observation:  FMC = 0.1 + 0.6 * NDWI +- 0.2
+  fmc_gl_ndwi_0 = 12*0.1,
+  fmc_gl_ndwi_rate = 10*0.6,
+  fmc_gl_ndwi_stdev = 0.2,
  /
 
 &moisture

wrfv2_fire/test/em_fire/ndwi/readme.txt

@@ -0,0 +1,8 @@
+Test example for assimilation of NDWI
+
+1. run ideal.exe
+2. in matlab, cd up directory (em_fire), startup, cd ndwi
+3. still in matlab, set_wrfinput, note the value of fmc_g at the end
+3. look at output - the values of fmc_g should be the same
+4. in matlab, confirm fmc_g=ncread('wrfout_d01_0001-01-01_00:00:00','FMG_C') 
+for fmc_g(:,:,2) (fmc_g is set only after the initial wrfout is made) 

wrfv2_fire/test/em_fire/ndwi/set_wrfinput.m

@@ -0,0 +1,31 @@
+disp('Set wrfinput for NDWI assimilation testing')
+format compact
+
+f='wrfinput_d01'
+
+ndwi = 0.5
+live_fmc_prior=0.1
+
+w=ncread(f,'NDWI');
+w=ndwi * ones(size(w));
+ncreplace(f,'NDWI',w)
+
+% set FMC_GC(:,:,5)=0.1
+u=ncread(f,'FMC_GC');
+u(:,:,5)=live_fmc_prior;
+ncreplace(f,'FMC_GC',u);
+
+% in namelist.fire we have set
+fmc_gw05 = 0.5 ,      % fuel is consists of 0.5 live fuel
+fmc_gl_stdev = 0.1,   % live fuel moisture prior has error +- 0.1
+% ! observation:  FMC = 0.1 + 0.6 * NDWI +- 0.2
+fmc_gl_ndwi_0 = 0.1,
+fmc_gl_ndwi_rate = 0.6,
+fmc_gl_ndwi_stdev = 0.2,
+
+% so we should see
+live_fmc_obs = fmc_gl_ndwi_0 + fmc_gl_ndwi_rate * ndwi
+w_prior = 1/(fmc_gl_stdev^2)
+w_obs   = 1/(fmc_gl_ndwi_stdev^2)
+fmc_live = (live_fmc_prior*w_prior + live_fmc_obs*w_obs)/(w_prior+w_obs)
+fmc_g = fmc_live * fmc_gw05    % contribution to fmc_g on fire mesh