fire nodes now midpoint, model_test seems to work, new ref

wrf coupling not finished yet

other/Matlab/vis/go.m

@@ -1,2 +1,2 @@
-vc2d
+vc2d(0.5)
 exit

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

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)

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')

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
 

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

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
 

wrfv2_fire/phys/module_fr_sfire_util.F

@@ -409,17 +409,20 @@ 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)
         ie=min(ite1,io+ir)
         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)=                     &