From 75108761b97f96bb64ed2285329167a5ec0ab082 Mon Sep 17 00:00:00 2001 From: Jan Mandel Date: Fri, 24 Feb 2012 06:56:34 -0700 Subject: [PATCH] added subroutines fuel_moisture and advance_moisture --- wrfv2_fire/Registry/registry.fire | 2 +- wrfv2_fire/phys/module_fr_sfire_atm.F | 76 --------- wrfv2_fire/phys/module_fr_sfire_phys.F | 220 +++++++++++++++++-------- wrfv2_fire/phys/module_fr_sfire_util.F | 80 +++++++++ 4 files changed, 232 insertions(+), 146 deletions(-) diff --git a/wrfv2_fire/Registry/registry.fire b/wrfv2_fire/Registry/registry.fire index c93199ec96..be77adfda0 100644 --- a/wrfv2_fire/Registry/registry.fire +++ b/wrfv2_fire/Registry/registry.fire @@ -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" diff --git a/wrfv2_fire/phys/module_fr_sfire_atm.F b/wrfv2_fire/phys/module_fr_sfire_atm.F index 9c8a22d09c..737077b468 100644 --- a/wrfv2_fire/phys/module_fr_sfire_atm.F +++ b/wrfv2_fire/phys/module_fr_sfire_atm.F @@ -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 diff --git a/wrfv2_fire/phys/module_fr_sfire_phys.F b/wrfv2_fire/phys/module_fr_sfire_phys.F index 40bd96430d..942022be4a 100644 --- a/wrfv2_fire/phys/module_fr_sfire_phys.F +++ b/wrfv2_fire/phys/module_fr_sfire_phys.F @@ -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) @@ -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", @@ -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). @@ -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 @@ -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 ) @@ -224,13 +304,11 @@ 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 @@ -238,38 +316,42 @@ subroutine advance_moisture( & 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 @@ -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 diff --git a/wrfv2_fire/phys/module_fr_sfire_util.F b/wrfv2_fire/phys/module_fr_sfire_util.F index 16d1a009e8..19c9a955b0 100644 --- a/wrfv2_fire/phys/module_fr_sfire_util.F +++ b/wrfv2_fire/phys/module_fr_sfire_util.F @@ -85,6 +85,86 @@ module module_fr_sfire_util contains +! +!***************************** +! + +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 + +! +!*********************************************************************** +! subroutine fire_ignition_convert (config_flags,ignition) USE module_configure