Skip to content

Commit

Permalink
adding update of live fuel moisture by ndwi
Browse files Browse the repository at this point in the history
  • Loading branch information
Jan Mandel committed Jul 6, 2015
1 parent c862db8 commit 254a102
Show file tree
Hide file tree
Showing 4 changed files with 72 additions and 22 deletions.
6 changes: 4 additions & 2 deletions wrfv2_fire/Registry/registry.fire
Expand Up @@ -53,7 +53,7 @@ state real fgrnhfx *i*j fire 1 z hr "FGRNHFX"
state real fgrnqfx *i*j fire 1 z hr "FGRNQFX" "moisture flux from ground fire" "W/m^2"
state real fcanhfx *i*j fire 1 z hr "FCANHFX" "heat flux from crown fire" "W/m^2"
state real fcanqfx *i*j fire 1 z hr "FCANQFX" "moisture flux from crown fire" "W/m^2"
state real fndwi *i*j fire 1 z hr "NDWI" "Normalized Difference Water Index on fire grid" "1"
state real fndwi *i*j fire 1 z hr "FNDWI" "Normalized Difference Water Index on fire grid" "1"

# fuel moisture model section
dimspec num_fmc - namelist=nfmc z fuel_moisture_classes
Expand All @@ -74,11 +74,13 @@ rconfig logical fmoist_run namelist,fire max_domains .false. h
rconfig logical fmoist_interp namelist,fire max_domains .false. hr "interpolate moisture from the model or the input to fuels on the fire grid"
rconfig logical fmoist_only namelist,fire max_domains .false. hr "only run moisture model, skip fire"
rconfig integer fmoist_freq namelist,fire max_domains 0 hr "fmoist_freq" "frequency to run moisture model 0: use fmoist_dt, k>0: every k timesteps" "1"
rconfig integer kfmc_ndwi namelist,fire 1 0 hr "KFMC_NDWI number of moisture class to update from NDWI, or zero"
rconfig integer kfmc_ndwi namelist,fire 1 0 hr "KFMC_NDWI" number of moisture class to update from NDWI, or zero"
rconfig integer fire_ndwi_from_atm namelist,fire 1 1 hr "FIRE_NDWI_FROM_ATM" "number of moisture class to update from NDWI, or zero"
rconfig real fmoist_dt namelist,fire max_domains 600 hr "fmoist_dt " "moisture model time step" "s"
rconfig real fmep_decay_tlag namelist,fire 1 0.01 hr "fmep_decay_tlag" "time constant of assimilated adjustments of equilibria decay" "1"
halo HALO_FIRE_MFG dyn_em 24:fmc_g
halo HALO_FIRE_MAG dyn_em 8:fmc_gc
halo HALO_FIRE_NDWI dyn_em 8:ndwi

# diagnostics
# for the actual modeled fire
Expand Down
45 changes: 33 additions & 12 deletions wrfv2_fire/phys/module_fr_sfire_driver.F
Expand Up @@ -121,7 +121,8 @@ subroutine sfire_driver_em ( grid , config_flags &
USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
USE module_comm_dm , ONLY : halo_fire_fuel_sub, halo_fire_tign_sub, halo_fire_wind_f_sub, &
halo_fire_wind_a_sub, halo_fire_ph_sub, halo_fire_zsf_sub, halo_fire_longlat_sub, &
halo_fire_phb_sub, halo_fire_z0_sub, halo_fire_lfn_sub, HALO_FIRE_MAG_sub, HALO_FIRE_MFG_sub
halo_fire_phb_sub, halo_fire_z0_sub, halo_fire_lfn_sub, HALO_FIRE_MAG_sub, HALO_FIRE_MFG_sub, &
halo_fire_ndwi_sub
#endif

use module_fr_sfire_atm, only: read_emissions_table, add_fire_emissions
Expand Down Expand Up @@ -361,6 +362,7 @@ subroutine sfire_driver_em ( grid , config_flags &
! base geopotential and roughness
#include "HALO_FIRE_PHB.inc"
#include "HALO_FIRE_Z0.inc"
#include "HALO_FIRE_NDWI.inc"

elseif(fire_ifun.eq.2)then
! halo exchange on zsf width 2
Expand Down Expand Up @@ -450,7 +452,9 @@ subroutine sfire_driver_em ( grid , config_flags &
grid%fmep, & ! fuel moisture extended model parameters
grid%fmc_equi, & ! fuel moisture fields updated, by class, equilibrium diagnostic
grid%fmc_lag, & ! fuel moisture fields updated, by class, tendency diagnostic
fp%fmc_g ) ! write-only alias. need to exit before using fp again
fp%fmc_g, & ! write-only alias. need to exit before using fp again
grid%ndwi, &
grid%fndwi)


#ifdef DM_PARALLEL
Expand Down Expand Up @@ -551,7 +555,9 @@ subroutine sfire_driver_phys (ifun, &
fmep, & ! fuel moisture extended model parameters
fmc_equi, & ! fuel moisture fields updated, by class equilibrium diagnostic
fmc_lag, & ! fuel moisture fields updated, by class tendency diagnostic
fmc_g) ! fuel moisture, alias of fp%fmc_g
fmc_g, & ! fuel moisture, alias of fp%fmc_g
ndwi, &
fndwi) ! ndwi on fire grid


implicit none
Expand Down Expand Up @@ -600,11 +606,11 @@ subroutine sfire_driver_phys (ifun, &
uah, & ! atm wind at fire_wind_height, diagnostics
vah ! atm wind at fire_wind_height, diagnostics

real, dimension(ims:ime, jms:jme), intent(inout)::xlong, xlat ! inout because of extension at bdry
real, dimension(ims:ime, jms:jme), intent(inout)::xlong, xlat, ndwi ! inout because of extension at bdry

real, intent(inout), dimension(ifms:ifme,jfms:jfme):: & ! fuel data; can be also set inside (cell based, fire grid)
fz0,fwh, &
nfuel_cat
nfuel_cat,fndwi

real, intent(inout), dimension(ifms:ifme, jfms:jfme):: &
lfn,tign,fuel_frac, & ! state: level function, ign time, fuel left
Expand Down Expand Up @@ -713,6 +719,7 @@ subroutine sfire_driver_phys (ifun, &
if(itimestep.le.fire_print_file.or.mod(itimestep,fire_print_file).eq.0)pid=itimestep ! print 1-fire_print_file then every fire_print_file-th
endif


if(ifun.eq.1)then
call init_fuel_cats(fmoist_run .or. fmoist_interp) ! properties of fuel categories and moisture classes from namelist.fire
endif
Expand Down Expand Up @@ -806,13 +813,28 @@ subroutine sfire_driver_phys (ifun, &

if(fire_run)then

if(ifun.eq.1)then ! set terrain
if(ifun.eq.2)then ! interpolate

if(restart)then

call message('restart - topo initialization skipped')
call message('restart - interpolation skipped')

else
if(kfmc_ndwi > 0 .and. fire_ndwi_from_atm.eq.1)then
call print_2d_stats(ips,ipe,jps,jpe,ims,ime,jms,jme,fndwi,'driver:ndwi')
call interpolate_z2fire(id,0, & ! for debug output, <= 0 no output, extend strip
ids,ide, jds,jde, & ! atm grid dimensions
ims,ime, jms,jme, &
ips,ipe,jps,jpe, &
its,ite,jts,jte, &
ifds, ifde, jfds, jfde, & ! fire grid dimensions
ifms, ifme, jfms, jfme, &
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
ndwi, & ! atm grid arrays in
fndwi) ! fire grid arrays out
call print_2d_stats(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fndwi,'driver:fndwi')
endif
!
! call print_2d_stats(ips,ipe,jps,jpe,ims,ime,jms,jme,zs,'driver:zs')
!
Expand Down Expand Up @@ -856,16 +878,13 @@ subroutine sfire_driver_phys (ifun, &
call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,xlat,'xlat',id)
call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,xlong,'xlong',id)
#endif
endif

endif

elseif(ifun.eq.2)then

! after the loop where zsf created exited and all synced
call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%zsf,'driver_phys:zsf')

! cannot initialize moisture model because T2 Q2 PSFC are not set yet
endif
endif

elseif(ifun.eq.3)then ! interpolate winds to the fire grid

Expand Down Expand Up @@ -948,6 +967,7 @@ subroutine sfire_driver_phys (ifun, &
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
nfuel_cat, & ! fuel data
fndwi, & ! satellite sensing interpolated on fire grid
fmc_gc, & ! moisture contents by class on atmospheric grid
fmc_g & ! weighted fuel moisture contents on fire grid
)
Expand Down Expand Up @@ -1429,6 +1449,7 @@ subroutine set_flags(config_flags)
fire_use_windrf = config_flags%fire_use_windrf
fire_fmc_read = config_flags%fire_fmc_read
kfmc_ndwi = config_flags%kfmc_ndwi
fire_ndwi_from_atm = config_flags%fire_ndwi_from_atm

end subroutine set_flags

Expand Down
40 changes: 33 additions & 7 deletions wrfv2_fire/phys/module_fr_sfire_phys.F
Expand Up @@ -161,7 +161,8 @@ module module_fr_sfire_phys
fueldepthm,fueldens,fuelmce, &
fcwh,fcz0, ffw, &
savr,st,se,adjr0,adjrw,adjrs, &
fmc_gw(mfuelcats,max_moisture_classes)
fmc_gl_stdev,fmc_gl_ndwi_0,fmc_gl_ndwi_rate,fmc_gl_ndwi_stdev
REAL , DIMENSION( mfuelcats , max_moisture_classes), save :: fmc_gw
! =============================================================================
! Standard 13 fire behavior fuel models (for surface fires), along with some
! estimated canopy properties (for crown fire).
Expand Down Expand Up @@ -217,6 +218,10 @@ module module_fr_sfire_phys
!DATA fcz0 /0.3, 0.3, 0.75, 18., 0.6, 0.75, 0.75, &
! 0.06, 0.06, 0.3, 0.3, 0.69, 0.9, 0.9, zf*0 / ! fuel roughness height, added jm 2/23/11
DATA ffw /nf* 0.9, zf*0/
DATA fmc_gl_ndwi_0 /nf*0.1, zf*0./
DATA fmc_gl_ndwi_rate /nf*0.6, zf*0./
DATA fmc_gl_ndwi_stdev /nf*0.2, zf*0./
DATA fmc_gl_stdev /nf*0.2, zf*0./
DATA adjr0 /mfuelcats*1./
DATA adjrw /mfuelcats*1./
DATA adjrs /mfuelcats*1./
Expand All @@ -239,6 +244,7 @@ subroutine fuel_moisture( &
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
nfuel_cat, & ! fuel data
fndwi, & ! satellite sensing on fire grid
fmc_gc, & ! moisture contents by class on atmospheric grid
fmc_g & ! weighted fuel moisture contents on fire grid
)
Expand All @@ -258,15 +264,18 @@ subroutine fuel_moisture( &
ir,jr ! atm/fire grid ratio


real,intent(in),dimension(ifms:ifme,jfms:jfme):: nfuel_cat ! fuel data
real,intent(in),dimension(ifms:ifme,jfms:jfme):: nfuel_cat, & ! fuel data
fndwi ! satellite sensing interpolated to fire grid
real,intent(inout),dimension(ims:ime,nfmc,jms:jme):: fmc_gc
real,intent(out),dimension(ifms:ifme,jfms:jfme):: fmc_g ! fuel data

!**** local
real, dimension(its-1:ite+1,jts-1:jte+1):: fmc_k ! copy of fmc_gc(:,k,:)
real, dimension(ifts:ifte,jfts:jfte):: fmc_f ! interpolation of fmc_gc(:,k,:) to the fire grid
real, dimension(ifts:ifte,jfts:jfte):: fmc_f, & ! interpolation of fmc_gc(:,k,:) to the fire grid
nwdi_f ! inerpolation of nwdi to the fire grid
integer::i,j,k,n
integer::ibs,ibe,jbs,jbe
real::f1,w1,w2,f2
character(len=128)::msg

call check_mesh_2dim(ifts,ifte,jfts,jfte,ifds,ifde,jfds,jfde) ! check if fire tile fits into domain
Expand Down Expand Up @@ -312,20 +321,35 @@ subroutine fuel_moisture( &

call print_2d_stats(ifts,ifte,jfts,jfte,ifts,ifte,jfts,jfte,fmc_f,'fuel_moisture: fmc_f')

if(k .eq. kfmc_ndwi)then ! if live moisture, assimilate ndwi
call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fndwi,'fuel_moisture: fndwi')
write(msg,'(a,i4)')'Assimilating NDWI in fuel moisture class ',k
call message(msg)
endif

! add moisture contents for class k to the fuel moisture
do j=jfts,jfte
do i=ifts,ifte
n = nfuel_cat(i,j)
if(n > 0)then
fmc_g(i,j)=fmc_g(i,j)+fmc_gw(n,k)*fmc_f(i,j) ! add to sum over classes
f1=fmc_f(i,j)
if(k .eq. kfmc_ndwi)then ! if live moisture, assimilate ndwi
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
endif
fmc_g(i,j)=fmc_g(i,j)+fmc_gw(n,k)*f1 ! add to sum over classes
endif
enddo
enddo

! call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fmc_g,'fuel_moisture: fmc_g')

enddo

! call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fmc_g,'fuel_moisture: fmc_g')


end subroutine fuel_moisture

Expand Down Expand Up @@ -654,7 +678,9 @@ subroutine init_fuel_cats(init_fuel_moisture)
! read
namelist /fuel_scalars/ cmbcnst,hfgl,fuelmc_g,fuelmc_c,nfuelcats,no_fuel_cat,fire_wind_height,ibeh
namelist /fuel_categories/ fuel_name,windrf,fgi,fueldepthm,savr, &
fuelmce,fueldens,st,se,weight,fci_d,fct,ichap,fwh,fz0,ffw,adjr0,adjrw,adjrs,fmc_gw01,fmc_gw02,fmc_gw03,fmc_gw04,fmc_gw05
fuelmce,fueldens,st,se,weight,fci_d,fct,ichap,fwh,fz0,ffw, &
fmc_gl_ndwi_0,fmc_gl_ndwi_rate,fmc_gl_ndwi_stdev, fmc_gl_stdev, &
adjr0,adjrw,adjrs,fmc_gw01,fmc_gw02,fmc_gw03,fmc_gw04,fmc_gw05

namelist /moisture/ moisture_classes,drying_lag,wetting_lag,saturation_moisture,saturation_rain,rain_threshold, &
drying_model,wetting_model, moisture_class_name,fmc_gc_initialization
Expand Down
3 changes: 2 additions & 1 deletion wrfv2_fire/phys/module_fr_sfire_util.F
Expand Up @@ -42,7 +42,8 @@ module module_fr_sfire_util
fire_fmc_read=1, & ! fuel moisture: 0 from wrfinput, 1 from namelist.fire, 2 read from file in ideal
fire_hfx_given=0, & ! "0=no, run normally, 1=from wrfinput, 2=from file input_hfx in ideal, 4=by parameters" ""
fire_hfx_num_lines=1, & ! number of heatflux parameter sets defining the heaflux lines (must be 1)
kfmc_ndwi=0 ! if set, number of fuel moisture class to update from ndwi
fire_ndwi_from_atm=1, & ! interpolate ndwi from atmosphere resolution
kfmc_ndwi=0 ! if>0 , number of the fuel moisture class to update from ndwi


real, save:: &
Expand Down

0 comments on commit 254a102

Please sign in to comment.