diff --git a/other/Matlab/vis/go.m b/other/Matlab/vis/go.m index 7644b5d4..a8237070 100755 --- a/other/Matlab/vis/go.m +++ b/other/Matlab/vis/go.m @@ -1,2 +1,2 @@ -vc2d +vc2d(0.5) exit diff --git a/other/Matlab/vis/oneframe.m b/other/Matlab/vis/oneframe.m index 86bd9056..5b0a111d 100755 --- a/other/Matlab/vis/oneframe.m +++ b/other/Matlab/vis/oneframe.m @@ -1,9 +1,11 @@ -function oneframe(k) +function oneframe(k,dt) ext=sprintf('_%5.5i.txt',k); f=read_array_m(['fgrnhfx',ext]); u=read_array_m(['uf',ext]); v=read_array_m(['vf',ext]); l=read_array_m(['lfn',ext]); sf(f,u,v,l) - title(sprintf('time %i s',k)) + if ~exist('dt','var'),dt=1,end + title(sprintf('time %i s',k*dt)) + drawnow end \ No newline at end of file diff --git a/other/Matlab/vis/vc2d.m b/other/Matlab/vis/vc2d.m index ecdca205..2c7a9bd1 100755 --- a/other/Matlab/vis/vc2d.m +++ b/other/Matlab/vis/vc2d.m @@ -1,18 +1,19 @@ -function vc2d +function vc2d(dt) j=0; t1=clock; %a=avifile('sfire.avi'); -for k=1:10:4001 +for k=[1:10,100:100:4000] j=j+1; - oneframe(k); + oneframe(k,dt); %M=getframe; %a=addframe(a,M); grid off + print('-djpeg',sprintf('frame%5.5i',j)) %M(j)=getframe(gcf); - M(j)=getframe; - if mod(j,10)==0 | j< 10, - savemovie(M) - end + %M(j)=getframe; + %if mod(j,10)==0 | j< 10, + % savemovie(M) + %end fprintf('frame %g model time %g rendering time %g s\n',j,k,etime(clock,t1)) end savemovie(M) diff --git a/wrfv2_fire/phys/model_test.m b/wrfv2_fire/phys/model_test.m index 16c562a5..6b7a8ea6 100755 --- a/wrfv2_fire/phys/model_test.m +++ b/wrfv2_fire/phys/model_test.m @@ -49,8 +49,8 @@ function vis(u,f,vx,vy,dx,dy,tNow) contour3(y,x,u,[0 0],'k') drawn=true; case '2d' - xh=[1/2:m-3/2]*dx; - yh=[1/2:n-3/2]*dy; + xh=[1/2:m-1/2]*dx; + yh=[1/2:n-1/2]*dy; h=pcolor(xh,yh,f'); % shading('interp') set(h,'edgecolor','none') diff --git a/wrfv2_fire/phys/model_test_main.F b/wrfv2_fire/phys/model_test_main.F index 5f72310c..73a28903 100644 --- a/wrfv2_fire/phys/model_test_main.F +++ b/wrfv2_fire/phys/model_test_main.F @@ -149,7 +149,7 @@ subroutine model_test( & enddo endif enddo - if(istep.le.10.or.mod(istep,50).eq.0)then + if(istep.le.10.or.mod(istep,10).eq.0)then write(1,1)1.,1.,time_start write(1,1)sm,sn,((lfn(i,j),i=ifps,ifpe),j=jfps,jfpe) write(1,1)sm,sn,((tign(i,j),i=ifps,ifpe),j=jfps,jfpe) @@ -185,7 +185,7 @@ program model_test_main ny=400 msteps=200 msteps=6 -!msteps=3 +! msteps=100 fdx=6 fdy=6 diff --git a/wrfv2_fire/phys/model_test_out.txt.ref.gz b/wrfv2_fire/phys/model_test_out.txt.ref.gz index 0d12ed5d..42f3e567 100755 Binary files a/wrfv2_fire/phys/model_test_out.txt.ref.gz and b/wrfv2_fire/phys/model_test_out.txt.ref.gz differ diff --git a/wrfv2_fire/phys/module_fr_sfire_driver.F b/wrfv2_fire/phys/module_fr_sfire_driver.F index 5da6f095..fbe9d241 100755 --- a/wrfv2_fire/phys/module_fr_sfire_driver.F +++ b/wrfv2_fire/phys/module_fr_sfire_driver.F @@ -440,17 +440,21 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no !*** local -real, dimension(its-1:ite+1, jts-1:jte+1):: ua,va, & ! atm winds, averaged over height - za ! terrain height +#define IWS its-1 +#define IWE ite+1 +#define JWS jts-1 +#define JWE jte+1 +real, dimension(IWS:IWE,JWS:JWE):: ua,va, & ! atm winds, averaged over height + za ! terrain height integer:: i,j,k,ifte1,jfte1,jts1,jte1,its1,ite1 !*** executable k=kds ! the ground - jts1=max(jts-1,jds) ! lower loop limit by one less when at end of patch or domain - its1=max(its-1,ids) ! ASSUMES THE HALO IS THERE if patch != domain - jte1=min(jte+1,jde) - ite1=min(ite+1,ide) + jts1=max(JWS,jds) ! lower loop limit by one less when at end of domain + its1=max(IWS,ids) ! ASSUMES THE HALO IS THERE if patch != domain + jte1=min(JWE,jde) + ite1=min(IWE,ide) do j = jts1,jte1 do i = its1,ite1 ! average 1st 2 layers, correct const shift @@ -460,42 +464,25 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no za(i,j)=zs(i,j) enddo enddo - ! extend beyond the domain boundary to tile+one row of cells as constant - do i=ide+1,ite+1 - do j = jts-1,jte+1 - va(i,j)=va(ide,j) - ua(i,j)=ua(ide,j) - za(i,j)=za(ide,j) - enddo - enddo - do i=its-1,ids-1 - do j = jts-1,jte+1 - va(i,j)=va(ids,j) - ua(i,j)=ua(ids,j) - za(i,j)=za(ids,j) - enddo - enddo - do j=ide+1,jte+1 - do i = its-1,ite+1 - ua(i,j)=ua(i,jde) - va(i,j)=va(i,jde) - za(i,j)=za(i,jde) - enddo - enddo - do j=jts-1,jds-1 - do i = its-1,ite+1 - ua(i,j)=ua(i,jds) - va(i,j)=va(i,jds) - za(i,j)=za(i,jds) - enddo - enddo + ! extend beyond the domain boundary + call continue_at_boundary(1,0, & ! do x direction or y direction + ims,ime,jms,jme, & ! memory dims + ids,ide,jds,jde, & ! domain dims + its,ite,jts,jte, & ! tile dims + va) ! array + + call continue_at_boundary(0,1, & ! do x direction or y direction + ims,ime,jms,jme, & ! memory dims + ids,ide,jds,jde, & ! domain dims + its,ite,jts,jte, & ! tile dims + ua) ! array if (id.gt.0) then - call write_array_m(its-1,ite+1,jts-1,jte+1,its-1,ite+1,jts-1,jte+1,ua,'ua',id) - call write_array_m(its-1,ite+1,jts-1,jte+1,its-1,ite+1,jts-1,jte+1,va,'va',id) + call write_array_m(its1,ite1,jts1,jte1,IWS,IWE,JWS,JWE,ua,'ua',id) + call write_array_m(its1,ite1,jts1,jte1,IWS,IWE,JWS,JWE,va,'va',id) endif -call print_2d_stats_vec(its-1,ite+1,jts-1,jte+1,its-1,ite+1,jts-1,jte+1,ua,va, & +call print_2d_stats_vec(its1,ite1,jts1,jte1,IWS,IWE,JWS,JWE,ua,va, & 'driver: atm wind (m/s)') ! The interpolation routine goes over each atm array tile diff --git a/wrfv2_fire/phys/module_fr_sfire_phys.F b/wrfv2_fire/phys/module_fr_sfire_phys.F index ab9b91af..41501fde 100755 --- a/wrfv2_fire/phys/module_fr_sfire_phys.F +++ b/wrfv2_fire/phys/module_fr_sfire_phys.F @@ -183,11 +183,13 @@ subroutine set_fire_params( & ifts,ifte,jfts,jfte, & ! tile dims zsfc) ! array +fdxinv=1./fdx +fdyinv=1./fdy do j=jfts,jfte do i=ifts,ifte ! central differences - dzfsdx(i,j) = (zsfc(i+1,j)-zsfc(i-1,j))/fdx - dzfsdy(i,j) = (zsfc(i,j+1)-zsfc(i,j-1))/fdy + dzfsdx(i,j) = (zsfc(i+1,j)-zsfc(i-1,j))*fdxinv + dzfsdy(i,j) = (zsfc(i,j+1)-zsfc(i,j-1))*fdyinv enddo enddo diff --git a/wrfv2_fire/phys/module_fr_sfire_util.F b/wrfv2_fire/phys/module_fr_sfire_util.F index 6125fe6a..ea0eefac 100755 --- a/wrfv2_fire/phys/module_fr_sfire_util.F +++ b/wrfv2_fire/phys/module_fr_sfire_util.F @@ -409,10 +409,10 @@ subroutine interpolate_2d_nodes2nodes( & rx=ir ry=jr -do j2=jts2,jte2-1 - jo=jp1+jr*(j2-jp2) - js=max(jts1,jo) - je=min(jte1,jo+jr) +do j2=jts2,jte2-1 ! loop over coarse mesh cells + jo=jp1+jr*(j2-jp2) ! fine grid coordinate of the coarse grid patch start + js=max(jts1,jo) ! lower bound of fine grid patch for this cell + je=min(jte1,jo+jr) ! upper bound of fine grid patch for this cell do i2=its2,ite2-1 io=ip1+ir*(i2-ip2) is=max(its1,io) @@ -420,6 +420,9 @@ subroutine interpolate_2d_nodes2nodes( & do j1=js,je ty = (j1-jo)/ry do i1=is,ie + ! in case fine grid lies on the boundary of several cells + ! the result will be written multiple times with the same value + ! up to a rounding error tx = (i1-io)/rx !print *,'coarse ',i2,j2,'to',i2+1,j2+1,' fine ',is,js,' to ',ie,je v1(i1,j1)= &