Skip to content

Commit

Permalink
added subroutines fuel_moisture and advance_moisture
Browse files Browse the repository at this point in the history
  • Loading branch information
@janmandel
janmandel committed Feb 24, 2012
1 parent c77c36f commit 7510876
Show file tree
Hide file tree
Showing 4 changed files with 232 additions and 146 deletions.
2 changes: 1 addition & 1 deletion wrfv2_fire/Registry/registry.fire
View file
Expand Up @@ -52,7 +52,7 @@ state real fcanqfx *i*j fire 1 z hr "FCANQFX"
# fuel moisture model variables
dimspec num_fmc - namelist=nfmc c fuel_moisture_classes
rconfig integer nfmc namelist,fire 1 5 - "nfmc" "number of fuel moisture classes"
state real fmc_gc ij{num_fmc} fire 1 z hr "FMC_GC" "fuel moisture contents by class" "1"
state real fmc_gc i{num_fmc}j fire 1 z hr "FMC_GC" "fuel moisture contents by class" "1"
state real rain_old ij fire 1 z hr "RAIN_OLD" "previous value of accumulated rain" "mm"
state real t2_old ij fire 1 z hr "T2_OLD" "previous value of accumulated rain" "mm"
state real q2_old ij fire 1 z hr "Q2_OLD" "previous value of accumulated rain" "mm"
Expand Down
76 changes: 0 additions & 76 deletions wrfv2_fire/phys/module_fr_sfire_atm.F
View file
Expand Up @@ -1134,81 +1134,5 @@ end subroutine coarse

end subroutine interpolate_wind2fire_height

!
!*****************************
!

subroutine interpolate_z2fire(id, & ! for debug output, <= 0 no output
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
zs, & ! atm grid arrays in
zsf) ! fire grid arrays out

implicit none
!*** purpose: interpolate height

!*** arguments
integer, intent(in)::id, &
ids,ide, jds,jde, & ! atm domain bounds
ims,ime,jms,jme, & ! atm memory bounds
ips,ipe,jps,jpe, &
its,ite,jts,jte, & ! atm tile bounds
ifds, ifde, jfds, jfde, & ! fire domain bounds
ifms, ifme, jfms, jfme, & ! fire memory bounds
ifts,ifte,jfts,jfte, & ! fire tile bounds
ir,jr ! atm/fire grid refinement ratio
real, intent(in), dimension(ims:ime, jms:jme):: zs ! terrain height at atm cell centers & ! terrain height
real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
zsf ! terrain height fire grid nodes


!*** local
real, dimension(its-2:ite+2,jts-2:jte+2):: za ! terrain height
integer:: i,j,jts1,jte1,its1,ite1,jfts1,jfte1,ifts1,ifte1,itso,jtso,iteo,jteo

! terrain height

jts1=max(jts-1,jds) ! lower loop limit by one less when at end of domain
its1=max(its-1,ids) ! ASSUMES THE HALO IS THERE if patch != domain
jte1=min(jte+1,jde)
ite1=min(ite+1,ide)
do j = jts1,jte1
do i = its1,ite1
! copy to local array
za(i,j)=zs(i,j)
enddo
enddo

call continue_at_boundary(1,1,0., & ! do x direction or y direction
its-2,ite+2,jts-2,jte+2, & ! memory dims
ids,ide,jds,jde, & ! domain dims - winds defined up to +1
ips,ipe,jps,jpe, & ! patch dims - winds defined up to +1
its1,ite1,jts1,jte1, & ! tile dims
itso,jtso,iteo,jteo, &
za) ! array

! interpolate to tile plus strip along domain boundary if at boundary
jfts1=snode(jfts,jfds,-1) ! lower loop limit by one less when at end of domain
ifts1=snode(ifts,ifds,-1)
jfte1=snode(jfte,jfde,+1)
ifte1=snode(ifte,ifde,+1)

call interpolate_2d( &
its-2,ite+2,jts-2,jte+2, & ! memory dims atm grid tile
its1-1,ite1+1,jts1-1,jte1+1, & ! where atm grid values set
ifms,ifme,jfms,jfme, & ! array dims fire grid
ifts1,ifte1,jfts1,jfte1, & ! dimensions fire grid tile
ir,jr, & ! refinement ratio
real(ids),real(jds),ifds+(ir-1)*0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
za, & ! in atm grid
zsf) ! out fire grid

end subroutine interpolate_z2fire

end module module_fr_sfire_atm
220 changes: 151 additions & 69 deletions wrfv2_fire/phys/module_fr_sfire_phys.F
View file
Expand Up @@ -12,6 +12,7 @@ module module_fr_sfire_phys
implicit none
PRIVATE


type fire_params
real,pointer,dimension(:,:):: vx,vy ! wind velocity (m/s)
real,pointer,dimension(:,:):: zsf ! terrain height (m)
Expand All @@ -35,6 +36,31 @@ module module_fr_sfire_phys
! NOTE: fcwh and fcz0 are called fwh and fz0 in read/write statements


! moisture behavior
integer, parameter::max_moisture_classes=5
integer, parameter::zm=max_moisture_classes - 3
integer:: moisture_classes=3
real, dimension(max_moisture_classes):: drying_lag,wetting_lag,saturation_moisture,saturation_rain,rain_threshold
integer, dimension(max_moisture_classes):: drying_model,wetting_model

real, dimension(7)::eq_p

data drying_lag /1., 10., 100. , zm*0./ ! time lag (h) approaching equilibrium moisture
data wetting_lag /14, 140., 1400., zm*0./ ! time lag (h) for approaching saturation in rain
data saturation_moisture /2.5, 2.5, 2.5 , zm*0./ ! saturation moisture contents (1) in rain
data saturation_rain /8.0, 8.0, 8.0 , zm*0./ ! stronger rain matters only in duration (mm/h)
data rain_threshold /0.05, 0.05, 0.05, zm*0 /! rain intensity this small is same as nothing
data drying_model /1, 1, 1, zm*0 /
data wetting_model /1, 1, 1, zm*0 /
real:: fmc_gc_default=0.09
data eq_p/ 1.035e-09, & !(3.893e-10, 1.681e-09) ! coefficients of the equilibrium fuel moisture polynomial
-2.62e-07, & !(-4.593e-07, -6.473e-08) ! fitted from the graph in Schroeder and Buck (1970)
2.507e-05, & !(2.194e-06, 4.795e-05)
-0.001107, & !(-0.002353, 0.000139)
0.02245, & !(-0.009188, 0.05409)
-0.05901, & !(-0.3721, 0.254)
3.043/ !(2.17, 3.915)

! Following table copied from module_fr_cawfe_fuel by Ned Patton with minor changes.
! Based on: Clark, T. L., J. L. Coen and D. Latham: 2004,
! "Description of a coupled atmosphere-fire model",
Expand Down Expand Up @@ -108,7 +134,8 @@ module module_fr_sfire_phys
REAL , DIMENSION( mfuelcats ), save :: windrf,weight,fgi,fci,fci_d,fct,fcbr, &
fueldepthm,fueldens,fuelmce, &
fcwh,fcz0, ffw, &
savr,st,se,adjr0,adjrw,adjrs
savr,st,se,adjr0,adjrw,adjrs, &
fmc_gw(mfuelcats,max_moisture_classes)
! =============================================================================
! Standard 13 fire behavior fuel models (for surface fires), along with some
! estimated canopy properties (for crown fire).
Expand Down Expand Up @@ -167,41 +194,95 @@ module module_fr_sfire_phys
DATA adjr0 /mfuelcats*1./
DATA adjrw /mfuelcats*1./
DATA adjrs /mfuelcats*1./

! moisture behavior
integer, parameter::max_moisture_classes=5
integer, parameter::zm=max_moisture_classes - 3
integer:: moisture_classes=3
real, dimension(max_moisture_classes):: drying_lag,wetting_lag,saturation_moisture,hysteresis,saturation_rain
integer, dimension(max_moisture_classes):: drying_model,wetting_model

real, dimension(7)::eq_p

data drying_lag /1., 10., 100. , zm*0./ ! time lag (h) approaching equilibrium moisture
data hysteresis /0., 0., 0. , zm*0./ ! equilibrium moisture increases by 1/2 of this
! when approached from above and decreases when from below
data wetting_lag /14, 140., 1400., zm*0./ ! time lag (h) for approaching saturation in rain
! 7% / hour initial rise per Fosberg and Deeming (1971)
data saturation_moisture /2.5, 2.5, 2.5 , zm*0./ ! saturation moisture contents (1) in rain
data saturation_rain /8.0, 8.0, 8.0 , zm*0./ ! stronger rain matters only in duration (mm/h)
data drying_model /1, 1, 1, zm*0 /
data wetting_model /1, 1, 1, zm*0 /
data eq_p/ 1.035e-09, & !(3.893e-10, 1.681e-09) ! coefficients of the equilibrium fuel moisture polynomial
-2.62e-07, & !(-4.593e-07, -6.473e-08) ! fitted from the graph in Schroeder and Buck (1970)
2.507e-05, & !(2.194e-06, 4.795e-05)
-0.001107, & !(-0.002353, 0.000139)
0.02245, & !(-0.009188, 0.05409)
-0.05901, & !(-0.3721, 0.254)
3.043/ !(2.17, 3.915)

integer, parameter::nz_fmc_gw = (max_moisture_classes-1)*mfuelcats
DATA fmc_gw /mfuelcats*1,nz_fmc_gw*0./ ! row k are relative weights of moisture class k for each fuel category

! =========================================================================

logical, save :: have_fuel_cats=.false.

contains

subroutine fuel_moisture( &
id, & ! for prints and maybe file names
nfmc, &
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
nfuel_cat, & ! fuel data
fmc_gc, & ! moisture contents by class on atmospheric grid
fmc_g & ! weighted fuel moisture contents on fire grid
)

!**** arguments
integer, intent(in):: &
id,nfmc, &
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

real,intent(in),dimension(ifms:ifme,jfms:jfme):: nfuel_cat ! fuel data
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:its+1,jts-1:jts+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
integer::i,j,k,n


do j=jfts,jfte
do i=ifts,ifte
fmc_g(i,j)=0. ! initialize sum over classes
enddo
enddo

do k=1,nfmc

do j=jts-1,jts+1
do i=its-1,its+1
fmc_k(i,j)=fmc_gc(i,k,j) ! copy slice to 2d array
enddo
enddo

! interpolate moisture contents in the class k to the fire mesh
call interpolate_z2fire(id, & ! for debug output, <= 0 no output
ids,ide,jds,jde, & ! atm grid dimensions
its-1,ite+1,jms-1,jme+1, &
ips,ipe,jps,jpe, &
its,ite,jts,jte, &
ifds, ifde, jfds, jfde, & ! fire grid dimensions
ifts, ifte, jfts, jfte, &
ifts,ifte, jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
fmc_k, & ! atm grid arrays in
fmc_f) ! fire grid arrays out

! add moisture contents for class k to the fuel moisture
do j=jfts,jfte
do i=ifts,ifte
n = nfuel_cat(i,j)
fmc_g(i,j)=fmc_g(i,j)+fmc_gw(k,n)*fmc_f(i,j) ! initialize sum over classes
enddo
enddo

enddo


end subroutine fuel_moisture

subroutine advance_moisture( &
initialize, & ! 0 = run, 1 = initialize previous values, 2 = initialize also the moisture
initialize, & ! 0 = run, 1 = initialize previous values, 2 = initialize also moisture
fmc_gc_init, & ! initial value of moisture
ims,ime, jms,jme, & ! memory dimensions
its,ite, jts,jte, & ! tile dimensions
Expand All @@ -211,7 +292,6 @@ subroutine advance_moisture( &
t2, q2, psfc, & ! temperature (K), vapor contents (kg/kg), pressure (Pa) at the surface
rain_old, & ! previous value of accumulated rain
t2_old, q2_old, psfc_old, & ! previous values of the atmospheric state at surface
rh_sfc, & ! relative humidity
fmc_gc & ! fuel moisture fields updated, by class, assumed set to something reasonable
)

Expand All @@ -224,52 +304,54 @@ subroutine advance_moisture( &
real, intent(in):: fmc_gc_init,timestep
real, intent(in), dimension(ims:ime,jms:jme):: t2, q2, psfc, rainc, rainnc
real, intent(inout), dimension(ims:ime,jms:jme):: t2_old, q2_old, psfc_old, rain_old
real, intent(out), dimension(ims:ime,jms:jme):: rh_sfc ! diagnostics
real, intent(inout), dimension(ims:ime,jms:jme,nfmc):: fmc_gc
real, intent(inout), dimension(ims:ime,nfmc,jms:jme):: fmc_gc

!*** local
integer:: i,j,k
real::rain_int, T, P, Q, QRS, ES, RH, tend, EMC_d, EMC_w, EMC, R
real::rain_thr=0.05
real::rain_int, T, P, Q, QRS, ES, RH, tend, EMC_d, EMC_w, EMC, R, rain_diff, fmc

!*** executable

if(initialize.ge.1)call copy2old
if(initialize.ge.2)call init_fmc_gc

do j=jts,jte
do i=its,ite
rain_int = ((rainc(i,j) + rainnc(i,j)) - rain_old(i,j))/timestep ! rain intensity
R = rain_int - rain_thr
do k=1,nfmc
if (R > 0.) then
select case(wetting_model(k))
case(1) ! saturation_moisture=2.5 wetting_lag=14h saturation_rain=8 mm/h calibrated to VanWagner&Pickett 1985 per 24 hours
tend=(saturation_moisture(k)-fmc_gc(i,j,k))/wetting_lag(k) * (1. - exp(R/saturation_rain(k)))
end select
else ! not raining
! relative humidity
T = t2(i,j)
P = psfc(i,j)
Q = q2(i,j)
ES=610.78*exp(17.269*(T-273.161)/(T-35.861));
QRS=0.622*ES/(P-0.378*ES);
RH = 100.*Q/QRS;
rh_sfc(i,j)=RH ! diagnostics
select case(drying_model(k))
case(1) ! Van Wagner formula (1972) per Viney (1991)
EMC_d=0.924*RH**0.679 + 0.000499*exp(0.1*RH) + 0.18*(21.1-T)*(1-exp(-0.115*RH)) ! equilibrium moisture for drying
EMC_w=0.618*RH**0.753 + 0.000454*exp(0.1*RH) + 0.18*(21.1-T)*(1-exp(-0.115*RH)) ! equilibrium moisture for adsorbtion
if (fmc_gc(i,j,k) > EMC_d) then
tend = (EMC_d - fmc_gc(i,j,k)) / drying_lag(k)
elseif(fmc_gc(i,j,k) < EMC_w)then
tend = (EMC_w - fmc_gc(i,j,k)) / drying_lag(k)
else
tend = 0.
endif
end select
endif
fmc_gc(i,j,k) = fmc_gc (i,j,k)+ timestep*tend
do k=1,nfmc
do i=its,ite
! old fuel moisture contents
fmc = fmc_gc(i,k,j)
! compute the rain intensity from the difference of accumulated rain
rain_diff = ((rainc(i,j) + rainnc(i,j)) - rain_old(i,j))
rain_int = rain_diff/timestep
R = rain_int - rain_threshold(k)
if (R > 0.) then
select case(wetting_model(k))
case(1) ! saturation_moisture=2.5 wetting_lag=14h saturation_rain=8 mm/h calibrated to VanWagner&Pickett 1985 per 24 hours
tend=(saturation_moisture(k)-fmc)/wetting_lag(k) * (1. - exp(R/saturation_rain(k)))
end select
else ! not raining
! average the inputs for second order accuracy
T = 0.5 * (t2_old(i,j) + t2(i,j))
P = 0.5 * (psfc_old(i,j) + psfc(i,j))
Q = 0.5 * (q2_old(i,j) + q2(i,j))
! compute the relative humidity
ES=610.78*exp(17.269*(T-273.161)/(T-35.861))
QRS=0.622*ES/(P-0.378*ES)
RH = 100.*Q/QRS
select case(drying_model(k))
case(1) ! Van Wagner (1972) per Viney (1991)
EMC_d=0.924*RH**0.679 + 0.000499*exp(0.1*RH) + 0.18*(21.1-T)*(1-exp(-0.115*RH)) ! equilibrium moisture for drying
EMC_w=0.618*RH**0.753 + 0.000454*exp(0.1*RH) + 0.18*(21.1-T)*(1-exp(-0.115*RH)) ! equilibrium moisture for adsorbtion
if (fmc > EMC_d) then
tend = (EMC_d - fmc) / drying_lag(k)
elseif(fmc < EMC_w)then
tend = (EMC_w - fmc) / drying_lag(k)
else
tend = 0.
endif
end select
endif
fmc_gc(i,k,j) = fmc + timestep*tend
enddo
enddo
enddo
Expand All @@ -293,10 +375,10 @@ subroutine copy2old
end subroutine copy2old

subroutine init_fmc_gc
do k=1,nfmc
do j=jts,jte
do j=jts,jte
do k=1,nfmc
do i=its,ite
fmc_gc(i,j,k)=0.1
fmc_gc(i,k,j)=fmc_gc_init
enddo
enddo
enddo
Expand Down

0 comments on commit 7510876

Please sign in to comment.