diff --git a/wrfv2_fire/phys/Makefile b/wrfv2_fire/phys/Makefile index 453fcff5..3b5ec439 100644 --- a/wrfv2_fire/phys/Makefile +++ b/wrfv2_fire/phys/Makefile @@ -61,11 +61,10 @@ MODULES = \ module_fr_cawfe_fuel.o \ module_fire_debug_output.o \ # module_fr_spread_model.o \ - module_fr_spread_core.o \ - module_fr_propagate.o \ - module_fr_fuelused.o \ - module_fr_spread_util.o \ - module_fr_wrfstubs.o + module_fr_sfire_core.o \ + module_fr_sfire_prop.o \ + module_fr_sfire_fuel.o \ + module_fr_sfire_util.o OBJS = @@ -271,20 +270,13 @@ module_fr_spread_model.o: \ ../share/module_model_constants.o \ ../frame/module_wrf_error.o \ module_fr_cawfe_fuel.o \ - module_fr_wrfstubs.o \ - module_fr_spread_core.o \ - module_fr_spread_util.o - -module_fr_spread_util.o: \ - module_fr_wrfstubs.o \ - -module_fr_wrfstubs.o: \ - ../frame/module_wrf_error.o \ + module_fr_sfire_core.o \ + module_fr_sfire_util.o module_fr_spread_core.o: \ - module_fr_wrfstubs.o \ - module_fr_fuelused.o \ - module_fr_propagate.o + module_fr_sfire_util.o \ + module_fr_sfire_fuel.o \ + module_fr_sfire_prop.o module_fire_debug_output.o: \ ../frame/module_domain.o \ diff --git a/wrfv2_fire/phys/module_fr_sfire_core.F b/wrfv2_fire/phys/module_fr_sfire_core.F new file mode 100644 index 00000000..83f8ea03 --- /dev/null +++ b/wrfv2_fire/phys/module_fr_sfire_core.F @@ -0,0 +1,151 @@ +module module_fr_spread_core +! +!*** Jan Mandel August 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com +! +contains + +subroutine fire_spread_core( ids,ide,jds,jde,ims,ime,jms,jme,& +time_start,time_end,fire_dx,fire_dy,fuel_time, & +vx,vy,r, & +phi,tign,frac_lost,frac_end) + +!*** purpose +! +! Dynamical core of the fire spread model, insulated from the physics. +! References the fire grid ONLY. Keep separated from the coupling +! with atm. model and data input. +! + +use module_fr_sfire_fuel +use module_fr_sfire_prop + +implicit none + +!*** arguments + +integer, intent(in) :: ids,ide,jds,jde,ims,ime,jms,jme +real, intent(inout), dimension(ims:ime,jms:jme)::phi,tign ! state +real, intent(in), dimension(ims:ime,jms:jme)::vx,vy,r +real, intent(in), dimension(ims:ime-1,jms:jme-1)::fuel_time +real, intent(out), dimension(ims:ime-1,jms:jme-1)::frac_lost,frac_end +real, intent(in):: time_start,time_end,fire_dx,fire_dy + +! argument intent description (unit) lives at +! +! ids,ide,jds,jde in mesh domain dimensions (1) +! ims,ime,jms,jme in mesh aray dimensions (1) +! time_start in the starting time (s) +! time_end in the ending time (s) +! fire_dx,fire_dy in fire mesh spacings (m) +! fuel_time in time fuel burns down to 1/e (s) cells +! vx,vy in the wind velocity (m/s) nodes +! r in the spread rate w/o wind (m/s) nodes +! phi inout level function (state) (1) nodes +! tign inout ignition time (state) (s) nodes +! frac_lost out the fuel fraction decrease (1) cells +! frac_end out the fuel fraction at the end (1) cells + +! ***NOTE: vx,vy,r will be replaced by a call to a speed function in future.*** +! The speed function will be called for the spread rate at points of fireline. +! A version of such level set code is in progress. + +!*** description + +! This is a dynamical core of the fire spread model, insulated from the physics. +! The physics should be done in the pre- and postprocessing, and in +! the speed function for the fireline propagation (future). +! +! The state of the model is the level function phi, which determines the fire +! area, and the ignition time tign, both interpolated from values at nodes. +! The fire area is the level set where phi <= 0. The fireline is where phi=0. +! The array tign outside of the fire area is not set or referenced. +! The state should be preserved between the calls, and it can be modified by +! data assimilation. All other quantities are derived from the state in each call. +! +! The level function evolves the fireline with the speed in the normal direction given by +! the spread rate r and the normal component of the wind. The level set method +! takes care of of various special cases automagically, such as ignition of a cell +! surrounded by cells that all completely burning, and merging of approaching firelines. +! +! The fuel fraction is estimated from the the ignition times assuming +! exponential decrease since ignition with decrease of fule fraction to 1/e in +! fuel_tim. The ignition times at nodes are interpolated linearly from the +! evolving level function at the start and at the end, i.e. by assuming that the +! value of the level function at a point varies linearly with time. +! +! ***NOTE: If a narrow band scheme is used to advance the level function in time +! then the firelines at time_start and time_end must fit within the band +! at either time, and level function values away from the band should be set +! to some large positive and negative constants to assure that the ignition times +! in the area between the firelines are set reasonably. This is also important +! when the level function is used for data assimilation. *** +! +! It is the responsibility of the caller to: +! +! before the call +! - polulate all cells with the proper fuel_time coefficient +! - interpolate and correct atmospheric winds +! +! after the call +! - compute the fluxes from the fuel fraction burned +! and sum up the the fluxes over atmopheric grid cells +! +!*** local +real, dimension(ims:ime-1,jms:jme-1)::frac_start +real, dimension(ims:ime,jms:jme)::phi_start +real::dummy +integer::i,j +real, dimension(ids:ide,jds:jde)::phi_d,vx_d,vy_d,r_d + + +!*** executable + +! compute the fuel fraction at time_start +call fuel_left(ids,ide,jds,jde,ims,ime,jms,jme, & + phi,tign,fuel_time,time_start,frac_start) + +! save the level function for interpolation between the old and the new one later +do i=ids,ide + do j=jds,jde + phi_start(i,j)=phi(i,j) + enddo +enddo + +! advance the level function from time_start to time_end +! prop_ls will be replaced anyway, do not bother adding overlaps +phi_d=phi(ids:ide,jds:jde) +vx_d=vx(ids:ide,jds:jde) +vy_d=vy(ids:ide,jds:jde) +r_d=r(ids:ide,jds:jde) +call prop_ls(ids,ide,jds,jde,phi_d,dummy,time_start,time_end,vx_d,vy_d,r_d,fire_dx,fire_dy) +phi(ids:ide,jds:jde)=phi_d + +! set ignition time at the nodes the fireline moved over +do j=jds,jde + do i=ids,ide + if(phi_start(i,j)>0 .and. .not. phi(i,j)>0)then + ! node was not burning at start but it is burning at end + ! interpolate from the level functions at start and at end + ! phi_start(i,j) is the level function value at time_start + ! phi(i,j) is the level function value at time_end + ! 0 is the level function value at tign(i,j) + tign(i,j)=time_start + (time_end-time_start) & + * phi_start(i,j) / (phi_start(i,j) - phi(i,j)) + endif + enddo +enddo + +! compute the fuel fraction at time_end +call fuel_left(ids,ide,jds,jde,ims,ime,jms,jme, & + phi,tign,fuel_time,time_end,frac_end) + +! compute the fuel fraction lost +do j=jds,jde-1 + do i=ids,ide-1 + frac_lost(i,j)=frac_end(i,j)-frac_start(i,j) + enddo +enddo + +end subroutine fire_spread_core + +end module module_fr_spread_core diff --git a/wrfv2_fire/phys/module_fr_sfire_fuel.F b/wrfv2_fire/phys/module_fr_sfire_fuel.F new file mode 100644 index 00000000..bc980333 --- /dev/null +++ b/wrfv2_fire/phys/module_fr_sfire_fuel.F @@ -0,0 +1,126 @@ +module module_fr_sfire_fuel + +private +public:: fuel_left + +contains + +subroutine fuel_left(ids,ide,jds,jde,ims,ime,jms,jme,phi,tign,w,tnow,fuel_frac) +implicit none + +!*** purpose: determine fraction of fuel remaining + +!*** Jan Mandel August 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com + +!*** arguments + +integer, intent(in) :: ids,ide,jds,jde,ims,ime,jms,jme +real, intent(in), dimension(ims:ime,jms:jme)::phi,tign +real, intent(in), dimension(ims:ime-1,jms:jme-1)::w +real, intent(in):: tnow +real, intent(out), dimension(ims:ime-1,jms:jme-1)::fuel_frac + +! ids,ide,jds,jde in mesh domain dimensions (1) +! ims,ime,jms,jme in mesh aray dimensions (1) +! phi in level function, at nodes +! tign in ignition time, at nodes +! w in time constant of fuel, per cell +! tnow in time now +! fuel_frac out fraction of fuel remaining + +!*** Description +! The area burning is given by the condition P(x,y) <= 0, where the function P is +! interpolated from the values of phi at mesh nodes, +! P(dx*(i-1),dy*(j-1))=phi(i,j). +! +! The time since ignition in location (x,y) is the function T, interpolated in +! each mesh cell from the values T(dx*(i-1),dy*(j-1))=tign(i,j) at the nodes +! where phi(i,j)<=0, and T(x,y)=tnow on all points on the grid lines where P(x,y) = 0. +! The values of tign(i,j) where phi(i,j)>0 are ignored. +! +! The subroutine computes for each mesh cell [dx*(i-1),dx*i] by [dy*(j-1),dy*j] +! an approximation of the average of exp(-T(x,y)/w(i,j)) over the burning area +! in the cell, that is an approximation of the integral +! +! /\ +! 1 | T(x,y)-tnow +! fuel_frac(i,j) = ----- | exp( - ------------ ) dxdy +! dx*dy | w(i,j) +! \/ +! dx*(i-1)=0), then fuel_frac(i,j)=1. +! Because of symmetries, the result should not depend on the mesh spacing dx dy +! so dx and dy are not in the argument list. +! +! Example: +! +! phi<0 phi>0 +! (i,j+1)-----O--(i+1,j+1) O = points on the fireline, T=tnow +! | \ | A = the burning area for computing +! | \| fuel_frac(i,j) +! | A O +! | | +! | | +! (i,j)---------(i+1,j) +! phi<0 phi<0 +! +! Approximations allowed: +! The fireline can be approximated by straight line(s). +! When all cell is burning, approximation by 1 point Gaussian quadrature is OK. +! +! Requirements: +! The output should be a continuous function of the arrays phi and tign whenever +! phi(i,j)=0 implies tign(i,j)=tnow. +! The output should be invariant to the symmetries of the input in each cell. +! Arbitrary combinations of the signs of phi(i,j) should work. +! The result should be at least 1st order accurate, i.e. exact when +! exp(T) is replaced by a linear function equal to one at the fireline +! +! IMPORTANT: follow WRF coding conventions +! http://www.mmm.ucar.edu/wrf/WG2/WRF_conventions.html + +!*** local + +integer::i,j +real,dimension(ims:ime,jms:jme)::t,ap +real:: ta,pf,aps,ps,a + +! a very crude approximation - replace by a better code +! just an idea - not debugged yet + +do j=jds,jde ! note the order of indices for fast execution + do i=ids,ide + if (phi(i,j)>0) then + t(i,j)=1.0 ! add missing values - should be > tnow what the hell + else + t(i,j)=exp(tign(i,j)-tnow) + endif + ap(i,j)=abs(phi(i,j)) + enddo +enddo + +do j=jds,jde-1 + do i=ids,ide-1 ! it is OK to introduce extra scalars, just as fast + + ! average unscaled fuel fraction since ignition + ta=0.25d0*(t(i+1,j+1)+t(i+1,j)+t(i,j+1)+t(i,j)) + ps=phi(i+1,j+1)+phi(i+1,j)+phi(i,j+1)+phi(i,j) + aps=ap(i+1,j+1)+ ap(i+1,j)+ ap(i,j+1)+ ap(i,j) + + ! a=0 if all phi>0, 1 if all <0, transition except when all phi=0 + if (aps>0.0 .or. aps<0.0) then ! never compare =0.0, slow + a=(-ps/aps+1.0)*0.5d0 + else ! for fast code, the else clause should happen less often + a=0.5d0 ! just to have something, if all phi=0 it is junk anyway + endif + + fuel_frac(i,j)=a*ta**(1.0/w(i,j))+(1.0-a) + enddo +enddo + +end subroutine fuel_left + +end module module_fr_sfire_fuel diff --git a/wrfv2_fire/phys/module_fr_sfire_level.F b/wrfv2_fire/phys/module_fr_sfire_level.F new file mode 100644 index 00000000..63e0c89f --- /dev/null +++ b/wrfv2_fire/phys/module_fr_sfire_level.F @@ -0,0 +1,6 @@ +module module_fr_level +! Minjeong Kim, Min.Kim@cudenver.edu + +! the narrow band level set code comes here + +end module module_fr_level diff --git a/wrfv2_fire/phys/module_fr_sfire_model.F b/wrfv2_fire/phys/module_fr_sfire_model.F new file mode 100644 index 00000000..cee513d7 --- /dev/null +++ b/wrfv2_fire/phys/module_fr_sfire_model.F @@ -0,0 +1,94 @@ +module module_fr_spread_model +! +!*** Jan Mandel September 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com +! +contains + +subroutine fire_spread_model & + (ids,ide, kds,kde, jds,jde, & ! incoming + ims,ime, kms,kme, jms,jme, & + its,ite, kts,kte, jts,jte, & + ifds,ifde, kfds,kfde, jfds,jfde, & + ifms,ifme, kfms,kfme, jfms,jfme, & + itimestep,dt,dx,dy, & + grid_id,cen_lat,cen_lon,lat_ll,lon_ll, & + moad_cen_lat,moad_cen_lon, & + moad_lat_ll,moad_lon_ll,moad_dx,moad_dy, & + moad_s_we,moad_e_we,moad_s_sn,moad_e_sn, & + nfrx,nfry, & + tlat_stf,tlon_stf,t_ignite,ishape,ibeh, & + z1can,alfg,alfc,ifuelread,nfuel_cat0, & + z,z_at_w,dz8w,zs,u,v,mu,rho, & + nfuel_cat,nfl,nfl_t,nfl_c,ncod, & ! in and out + in1,in2,ixb,iyb,icn, & + fg,fc,r_0,bbb,betafl,phiwc,area,area2, & + zf,zsf,tign_g,tign_c,tign_crt, & + xfg,yfg,xcd,ycd,xcn,ycn,sprdx,sprdy, & + rthfrten,rqvfrten, & ! outgoing + grnhfx,grnqfx,canhfx,canqfx) ! temp? + + +use module_fr_spread_core +use module_fr_spread_util +implicit none + +!*** arguments + +INTEGER, INTENT(in) :: ids,ide, kds,kde, jds,jde ! atmosphere domain indices +INTEGER, INTENT(in) :: ims,ime, kms,kme, jms,jme ! atmosphere memory indices +INTEGER, INTENT(in) :: its,ite, kts,kte, jts,jte ! atmosphere tile indices + +INTEGER, INTENT(in) :: ifds,ifde, jfds,jfde, kfds,kfde ! fire domain indices +INTEGER, INTENT(in) :: ifms,ifme, jfms,jfme, kfms,kfme ! fire memory indices + +INTEGER, INTENT(in) :: itimestep ! current time step (cumulative) +REAL, INTENT(in) :: dt ! time step +REAL, INTENT(in) :: dx,dy ! dx,dy on innermost atm mesh + +! incoming atmos. winds (m/s at arakawa-c grid locations) +REAL, INTENT(in), DIMENSION( ims:ime, kms:kme, jms:jme ) :: u,v +! height above sea level of w points (m) +REAL, INTENT(in), DIMENSION( ims:ime, kms:kme, jms:jme ) :: z,z_at_w + + +! local +integer::if_st,if_en,jf_st,jf_en + +write(*,7001) & + ids,ide, kds,kde, jds,jde, & ! incoming + ims,ime, kms,kme, jms,jme, & + its,ite, kts,kte, jts,jte, & + ifds,ifde, kfds,kfde, jfds,jfde, & + ifms,ifme, kfms,kfme, jfms,jfme +7001 format('FR: indices to cawfe:'/'ids ',6i8/'ims ',6i8/'idfs ',6i8/'imfs ',6i8) + + CALL fire_winds(u,v, & ! in + ids,ide, kds,kde, jds,jde, & + ims,ime, kms,kme, jms,jme, & + its,ite, kts,kte, jts,jte, & + sfcu,sfcv) ! out + +! ----- get the time from model start (assumes non-variable dt) + + time = FLOAT(itimestep) * dt + +! ----- set indicies over which fire grid exists +! +! these indicies are needed so that we are properly handling +! fire calculations within tiles that butt up against domain +! boundaries where the halo information is not available. +! therefore the fire exists on one less atmospheric grid point +! than the innermost domain is dimensioned. + + if_st = MAX( (its-1)*nfrx+1, ids*nfrx+1 ) + if_en = MIN( (ite)*nfrx , (ide-1)*nfrx ) + jf_st = MAX( (jts-1)*nfry+1, jds*nfry+1 ) + jf_en = MIN( (jte)*nfry , (jde-1)*nfry ) + +! ----- begin initialization for all runs, including restart + + + +end subroutine fire_spread_model + +end module module_fr_spread_model diff --git a/wrfv2_fire/phys/module_fr_sfire_prop.F b/wrfv2_fire/phys/module_fr_sfire_prop.F new file mode 100644 index 00000000..7673a6ce --- /dev/null +++ b/wrfv2_fire/phys/module_fr_sfire_prop.F @@ -0,0 +1,217 @@ +module module_fr_sfire_prop + +contains + +subroutine prop_ls(m1,m,n1,n,phi,t,ts,te,vx,vy,r,dx,dy) + implicit none + + !*** purpose: advance level function in time + + ! Jan Mandel and Minjeong Kim, August 2007 + + !*** description + ! + ! Propagation of closed curve by a level function method. The level function + ! phi is defined by its values at the nodes of a rectangular grid. + ! The area where phi < 0 is inside the curve. The curve is + ! described implicitly by phi=0. Points where the curve intersects gridlines + ! can be found by linear interpolation from nodes. + ! + ! The level function is advanced from time ts to time te (up to rounding). + ! At exit, the actual time t is returned, which should essentially equal ts. + ! + ! The level function should be initialized to (an approximation of) the signed + ! distance from the curve. If the initial curve is a circle, the initial level + ! function is simply the distance from the center minus the radius. + ! + ! The curve moves outside with speed r in normal direction and by + ! advection by velocity field (vx,vy), given by values at grid nodes. + ! Thus, the speed of the curve evolution in the normal direction is + ! r + (vx,vy) dot normal vector. + ! + ! Currently r is a scalar but it could be easily given by a grid array. All that is + ! needed is to declare it so. + ! + ! Based on S. Osher and R. Fedkiw, Level set methods and dynamic implicit surfaces, + ! Springer, 2003, Sec. 6.4, as implemented in toolboxLS for Matlab by + ! I. Mitchell, A toolbox of Level Set Methods (Version 1.1), TR-2007-11, + ! Dept. Computer Science, University of British Columbia, 2007 + ! http://www.cs.ubc.ca/\~mitchell/ToolboxLS + ! + ! This code can reproduce (up to rounding) spinStarDemo from + ! toolboxLS with accuracy='low'. See prop_ls_readme.txt how. + ! + ! Method: Upwinding based on 1st order one-sided differences in space, + ! Euler method in time. Godunov method for the normal motion combined with + ! simple upwinding for advection. The timestep is set automatically based on + ! CFL condition. For a straight segment in a constant field and locally linear + ! level function, the method reduces to the exact normal motion. The advantage of + ! the level set method is that it treats automatically special cases such as + ! the curve approaching itself and merging components of the area inside the curve. + ! + ! IMPORTANT NOTE: This code passed all test except for a circular curve + ! moving in a circular field with center outside of the curve at the distance + ! of several diameters. The curve should travel on a circular path and grow; + ! instead it slowly shrinks. This seems to be an artefact of the low order + ! of the method. In toolboxLS, one has to set accuracy='high' to make it grow. + ! But this involves a high-order Runge-Kutta scheme in time and a quite complicated + ! high-order ENO method in space, and is much more expensive. + + !*** arguments + + ! m1,m,n1,n grid bounds (usually n1=m1=1) + ! phi the level function + ! t actual end time + ! ts star time + ! te end time + ! vx,vy velocity for advection + ! r normal speed, must be positive. + ! dx,dy grid spacing + + integer,intent(in)::m1,m,n1,n + real,dimension(m1:m,n1:n),intent(inout)::phi + real,intent(out)::t + real,dimension(m1:m,n1:n),intent(in)::vx,vy,r + real,intent(in)::dx,dy,ts,te + + !*** local arrays, same dimension as phi + + real,dimension(m1:m,n1:n)::tend,grad, & + diff2x,diff2y, & + diffLx,diffLy,diffRx,diffRy, & + flowLx,flowRx,flowLy,flowRy, & + diffCx,diffCy, & + tend_n,tend_a, & + tbound_np + + !*** local variables + + integer::istep,mstep + real::dt,tol,tbound,tbound_n,tbound_a,zero,one + parameter(zero=0.0,one=1.0) + + ! f90 intrinsic function + + intrinsic max,min,maxval,sqrt,lbound,ubound + + !*** executable + + !print *,'prop_ls:m1,m,n1,n,ts,te,r,dx,dy=', & + ! m1,m,n1,n,ts,te,r,dx,dy + + t=ts + tol=6.*10**(-14) ! 300*2.22*10^-16 + istep = 0 + mstep = 1000 + + do while ( t < (te-tol) .and. istep < mstep ) + + istep=istep+1 + + ! one sided differences + ! we waste a little and store them separately to make the code more + ! readable, and to allow for higher order scheme in future if needed + ! allow for general array bounds n1:n,m1:m - may be useful in parallel + + diffLx(m1+1:m,:) = (phi(m1+1:m,:)-phi(m1:m-1,:))/dx + diffRx(m1:m-1,:) = diffLx(m1+1:m,:) + diffLx(m1 ,:) = diffLx(m1+1,:) + diffRx(m ,:) = diffLx(m,:) + diffLy(:,n1+1:n) = (phi(:,n1+1:n)-phi(:,n1:n-1))/dy + diffRy(:,n1:n-1) = diffLy(:,n1+1:n) + diffLy(: ,n1) = diffLy(:,n1+1) + diffRy(: ,n) = diffLy(:,n) + + ! 2 times central differences + + diffCx=diffLx + diffRx + diffCy=diffLy + diffRy + + ! Godunov scheme: choose the upwind direction, L or R or none + ! always test on > or < never = , much faster because of IEEE + + ! flowLx = (diffLx >=0 & diffCx >=0) + + where (.not. diffLx<0 .and. .not. diffCx<0) + flowLx=1 + elsewhere + flowLx=0 + end where + + ! flowRx=(diffRx<=0 & diffCx<0) + + where (.not. diffRx>0 .and. diffCx<0) + flowRx=1 + elsewhere + flowRx=0 + end where + + ! flowLy = (diffLy >=0 & diffCy >=0) + + where (.not. diffLy<0 .and. .not. diffCy<0) + flowLy=1 + elsewhere + flowLy=0 + end where + + ! flowRy=(diffRy<=0 & diffCy<0) + + where (.not. diffRy>0 .and. diffCy<0) + flowRy=1 + elsewhere + flowRy=0 + end where + + ! Godunov scheme: choose the proper one sided difference or zero + + diff2x=diffLx*flowLx + diffRx*flowRx + diff2y=diffLy*flowLy + diffRy*flowRy + + ! magnitude of the gradient + ! do not use **2 who knows if the compiler is smart enough + + grad=sqrt(diff2x*diff2x + diff2y*diff2y) + + ! time step bound - CFL condition + + ! contribution of the normal term + + where (grad > 0) tbound_np = r*(abs(diff2x)/dx+abs(diff2y)/dy)/grad + tbound_n = maxval(tbound_np, grad>0) ! do not use undefined entries + + ! contribution of the advection term + + tbound_a = maxval(abs(vx))/dx+maxval(abs(vy))/dy + + ! the final CFL bound and the timestep + + tbound = 1/(tbound_n+tbound_a+tol) + dt = min(te-t, 0.5*tbound) + + ! the rhs of the diff eq, a.k.a. phi dot, a.k.a. the phi "tendency" + ! 0. because max requires both arguments same type (real) + + tend_n = -r*grad ! normal term + tend_a = -(diffLx*max(vx,zero)+diffRx*min(vx,zero)+ & ! advection term + diffLy*max(vy,zero)+diffRy*min(vy,zero)) ! advection term + tend=tend_n + tend_a + + ! trailing edge correction - do not allow fireline to go backwards + ! THIS IS DIFFERENT FROM THE TOOLBOX!! + where (phi<=0.0) tend=min(tend,-0.5*one) + + !print *,'prop_ls:step,time,dt=',istep,t,dt + + ! advance one step + + phi = phi + dt*tend + t = t + dt + + end do !while + + !if (istep >= mstep) print *,'max iteration reached' + !print *,'prop_ls:step,ts,te=',istep,ts,te + +end subroutine prop_ls + +end module module_fr_sfire_prop diff --git a/wrfv2_fire/phys/module_fr_sfire_util.F b/wrfv2_fire/phys/module_fr_sfire_util.F new file mode 100755 index 00000000..dd5e5f0a --- /dev/null +++ b/wrfv2_fire/phys/module_fr_sfire_util.F @@ -0,0 +1,324 @@ +module module_fr_sfire_util + +contains + +! +!**************** +! +subroutine crash(s) +use module_wrf_error +implicit none +character(len=*), intent(in)::s +character(len=128)msg +msg='SFIRE:'//s +call wrf_error_fatal(msg) +end subroutine crash + +! +!**************** +! + +subroutine message(s) +use module_wrf_error +implicit none +character(len=*), intent(in)::s +character(len=128)msg +msg='SFIRE:'//s +call wrf_message(msg) +end subroutine message + +subroutine sum_2dmesh_cells( & + ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2,& ! input + ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1) ! output +implicit none + +!*** purpose +! sum cell values in mesh2 to cell values of coarser mesh1 + +!*** arguments +! the dimensions are in cells, not nodes! + +integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1 +real, intent(out)::v1(ims1:ime1,jms1:jme1) +integer, intent(in)::ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2 +real, intent(in)::v2(ims2:ime2,jms2:jme2) + +!*** local +integer:: i1,i2,j1,j2,ir,jr,isz1,isz2,jsz1,jsz2,ioff,joff +real t +character(len=128)msg + +!*** executable + +!check mesh dimensions and domain dimensions +call check_mesh_dim(ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1) +call check_mesh_dim(ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2) + +! compute mesh sizes +isz1 = ide1-ids1+1 +jsz1 = jde1-jds1+1 +isz2 = ide2-ids2+1 +jsz2 = jde2-jds2+1 + +! check mesh sizes +if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)then + call message('all mesh sizes must be positive') + goto 9 +endif +if(mod(isz2,isz1).ne.0.or.mod(jsz2,jsz1).ne.0)then + call message('input mesh size must be multiple of output mesh size') + goto 9 +endif + +! compute mesh ratios +ir=isz2/isz1 +jr=jsz2/jsz1 + +! v1 = v1 + sum(v2) +do j1=jds1,jde1 + do i1=ids1,ide1 + t=0 + do joff=0,jr-1 + do ioff=0,ir-1 + i2=ioff+ids2+ir*(i1-ids1) + j2=joff+jds2+jr*(j1-jds1) + t=t+v2(i2,j2) + enddo + enddo + v1(i1,j1)=t + enddo +enddo + +return +9 continue +write(msg,91)ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2 +call message(msg) +write(msg,91)ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1 +call message(msg) +write(msg,92)'input mesh size:',isz2,jsz2 +call message(msg) +91 format('dimensions: ',8i8) +write(msg,92)'output mesh size:',isz1,jsz1 +call message(msg) +92 format(a,2i8) +call crash('module_fr_spread_util:sum_mesh_cells: bad mesh sizes') +end subroutine sum_2dmesh_cells + +! +!**************** +! + +subroutine interpolate_2dmesh_nodes( & + ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2,v2, & ! in + ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,v1 ) ! out +implicit none + +!*** purpose +! interpolate nodal values in mesh2 to nodal values in mesh1 +! input mesh 2 is coarse output mesh 1 is fine + +!*** arguments + +integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1 +real, intent(out)::v1(ims1:ime1,jms1:jme1) +integer, intent(in)::ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2 +real, intent(in)::v2(ims2:ime2,jms2:jme2) + +!*** local +integer:: i1,i2,j1,j2,ir,jr,isz1,isz2,jsz1,jsz2,ioff,joff +real:: tx,ty,rx,ry +character(len=128)msg + +!*** executable + +!check mesh dimensions and domain dimensions +call check_mesh_dim(ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1) +call check_mesh_dim(ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2) + +! compute mesh sizes +isz1 = ide1-ids1 +jsz1 = jde1-jds1 +isz2 = ide2-ids2 +jsz2 = jde2-jds2 + +! check mesh sizes +if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)goto 9 +if(mod(isz1,isz2).ne.0.or.mod(jsz1,jsz2).ne.0)goto 9 + +! compute mesh ratios +ir=isz1/isz2 +jr=jsz1/jsz2 +rx=ir +ry=jr + +! v1 = v1 + interpolate(v2) +! this loop goes over coarse mesh lines twice but avoid complications +do j2=jds2,jde2-1 + do i2=ids2,ide2-1 + do ioff=0,ir + do joff=0,jr + ! compute fine mesh coordinate + i1=ioff+ids1+ir*(i2-ids2) + j1=joff+jds1+jr*(j2-jds2) + ! weights + tx = ioff/rx + ty = joff/ry + ! interpolation + v1(i1,j1)= & + (1-tx)*(1-ty)*v2(i2,j2) & + + tx*(1-ty) *v2(i2,j2+1) & + + (1-tx)*ty *v2(i2+1,j2) & + + tx*ty *v2(i2+1,j2+1) + !print *,'coarse ',i2,j2,' fine ',i1,j1, 'weights ',tx,ty, & + ! 'in ',v2(i2,j2),v2(i2,j2+1),v2(i2+1,j2),v2(i2+1,j2+1),' out ',v1(i1,j1) + enddo + enddo + enddo +enddo + +return +9 continue +write(msg,91)ids2,ide2,jds2,jde2,ims2,ime2,jms2,jme2 +call message(msg) +write(msg,91)ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1 +call message(msg) +write(msg,92)'input mesh size:',isz2,jsz2 +call message(msg) +91 format('dimensions: ',8i8) +write(msg,92)'output mesh size:',isz1,jsz1 +call message(msg) +92 format(a,2i8) +call crash("module_fr_spread_util:interpolate_mesh_nodes: bad mesh sizes") +end subroutine interpolate_2dmesh_nodes + +! +!**************** +! + +subroutine check_mesh_dim(ids,ide,jds,jde,ims,ime,jms,jme) +integer, intent(in)::ids,ide,jds,jde,ims,ime,jms,jme +if(idsime.or.jdsjme)call crash('domain dimensions too large') +end subroutine check_mesh_dim + +! +!**************** +! + +real function interp(ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,ip1,jp1,x1,y1,dx1,dy1,v1,x,y) +implicit none +!*** purpose +! general interpolation in a mesh aligned with coordinates + +!*** arguments +! the mesh is given by the position (x1,y1) of the node (ip1,ip2) and the spacing (dx1,dy1) +! the coordinates to interpolate to are x y +! returns the interpolation from v1 at mesh nodes + +integer, intent(in)::ids1,ide1,jds1,jde1,ims1,ime1,jms1,jme1,ip1,jp1 +real, intent(in)::x1,y1,dx1,dy1,v1(ims1:ime1,jms1:jme1),x,y + +!*** calls +intrinsic floor + +!*** local +integer i1,j1,i2,j2 +real tx,ty,xx,yy + +! executable + +! indices of the lower left corner of the cell in mesh 1 that contains (x,y) +i1 = ip1+floor((x - x1)/dx1) +i1=max(min(i1,ide1-1),ids1) +j1 = jp1+floor((y - y1)/dy1) +j1=max(min(j1,jde1-1),jds1) + +! position of node (i1,j1) +xx = x1 + dx1*(i1-ip1) +yy = y1 + dy1*(j1-jp1) + +! the leftover +tx = x - xx +ty = y - yy +if(tx.gt.1.0)call crash('bad tx') +if(ty.gt.1.0)call crash('bad ty') + +! interpolate the values +interp = & + (1-tx)*(1-ty)*v1(i1,j1) & + + tx*(1-ty) *v1(i1+1,j1) & + + (1-tx)*ty *v1(i1,j1+1) & + + tx*ty *v1(i1+1,j1+1) + +print *,'x,y=',x,y,'i1,i2=',i1,j1,'tx,ty=',tx,ty,' interp=',interp +end function interp + +! +!**************** +! + +SUBROUTINE fire_winds(u,v, & ! incoming + ids,ide, kds,kde, jds,jde, & + ims,ime, kms,kme, jms,jme, & + its,ite, kts,kte, jts,jte, & + u_i,v_i) ! outgoing + +! This subroutine is taken from module_fr_cawfe by Ned Patton + +! --- this routine takes u and v from the arakawa c-grid and interpolates +! them horizontally and upward to the w-level (i.e. to the grid cube corners) +! as desired by the fire code. note that the final values are two +! dimensional arrays that are six grid points tall valid at the w-levels +! and that the exterior single grid point on all four edges of the domain +! are not filled. +! +! v(1,2) u(1,2) u(2,2) +! ----*---- v(1,2) *--------* v(2,2) +! | | | | +! u(1,1) * * u(2,1) ===> | | +! | | | | +! ^ y ----*---- u(1,1) *--------* u(2,1) +! | v(1,1) v(1,1) v(2,1) +! | +! *----> x and shifted up to w-level +! + IMPLICIT NONE + + INTEGER , INTENT(in) :: ids,ide, kds,kde, jds,jde, & + ims,ime, kms,kme, jms,jme, & + its,ite, kts,kte, jts,jte + + REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: u,v + + REAL, INTENT(out), DIMENSION( ims:ime,jms:jme,6 ) :: u_i,v_i + + INTEGER :: i,j,k + INTEGER :: i_st,i_en + INTEGER :: j_st,j_en + + ! --- set indicies + + i_st = MAX(its,ids+1) + i_en = MIN(ite,ide-1) + j_st = MAX(jts,jds+1) + j_en = MIN(jte,jde-1) + + ! --- get velocities + + DO k = 1,6 + DO j = j_st,j_en + DO i = i_st,i_en + u_i(i,j,k) = .25*( u(i,k,j) + u(i,k,j+1) + u(i,k+1,j) + u(i,k+1,j+1) ) + v_i(i,j,k) = .25*( v(i-1,k,j) + v(i,k,j) + v(i-1,k+1,j) + v(i,k+1,j) ) + END DO + END DO + END DO + + RETURN + +END SUBROUTINE fire_winds + +! +!**************** +! + +end module module_fr_sfire_util