Test file is finished, got rid of comments and added comments to core…

… file,

still left some checkings in comments

standalone/fuel_interp_test_main.F

@@ -2,24 +2,20 @@ program test
 	use module_fr_sfire_core
         implicit none
 real :: time_now
-integer :: ims,ime,jms,jme,icl,jcl,i,j,num_tests
+integer :: ims,ime,jms,jme,icl,jcl,i,j,num_tests,flag
 real,dimension(3,3)::tff,lff,lfn,tign,real_tff,real_lff
 real,dimension(3,3)::tff0,lff0,lfn0,tign0,real_tff0,real_lff0
 real,dimension(3,3)::tff1,lff1,lfn1,tign1,real_tff1,real_lff1
 real,dimension(3,3)::tff2,lff2,lfn2,tign2,real_tff2,real_lff2
 real,dimension(3,3)::tff3,lff3,lfn3,tign3,real_tff3,real_lff3
-real,dimension(3,3)::tff4,lff4,lfn4,tign4
-real,dimension(3,3)::tff5,lff5,lfn5,tign5
-real,dimension(3,3)::tff6,lff6,lfn6,tign6
-real,dimension(3,3)::tff7,lff7,lfn7,tign7
 real,dimension(3,3)::tff8,lff8,lfn8,tign8
 real,dimension(3,5)::tff9,lff9,lfn9,tign9
 real, dimension(5,5)::d
 real,dimension(2,2)::result
-real::test_err,error,test1_err,test2_err,test3_err,test4_err,test5_err
+real::test_err,error
 real::lfn00,lfn01,lfn10,lfn11,tign00,tign01,tign10,tign11,fuel_time_cell
-real:: res,Result1,Result2,Result3,Result4,Result5,result6
-real:: res1,res2,res3,res4,res5
+real:: Result1,Result2,Result3,Result4,Result5,result6
+real:: res,res1,res2,res3,res4!
 time_now=3.
 101 format(a/(3f7.3))
 ims=1
@@ -368,19 +364,33 @@ program test
 end do
 end do
 
+if (res1-result(1,1)>error.or.res2-result(1,2)>error.or.res3-result(2,1)>error.or.res4-result(2,2)>error) then
+	print*,'Test #10 is not working, one of errors '
+        print*,'For cell (1,1)',res1-result(1,1)
+        print*,'For cell (1,2)',res2-result(1,2)
+        print*,'For cell (2,1)',res3-result(2,1)
+        print*,'For cell (2,2)',res4-result(2,2)
+        print*,'is greater than 0.0001',error
+else
+num_tests=num_tests+1;
+endif
 
 print*,""
 print*,""
 !!!! Test #2
 
-write(*,*)'MAIN TEST #2'
+write(*,*)'MAIN TEST #2 (#11)'
 print*,"Check for consistency, for 3 inner cells with the same results"
 time_now=2.
 ims=1
 ime=5
 jms=1
 jme=5
-
+res1=2.9752314E-02
+res2=1.000
+res3=2.9752314E-02
+res4=1.000
+flag=1
 data lfn9/-1., -1.,-1., &
          -1., -1.,0., &
          0., 0.,0., &
@@ -396,19 +406,7 @@ program test
 do icl=2,4
         write(*,*)"icl,jcl",icl,jcl
        call tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, &
-                                    tign9,lfn9,tff,lff)
-
-d=1e20
-d(1:5:2,1:5:2)=tign9(icl-1:icl+1,1:3)
-d(2:4,2:4)=tff
-!write(*,*)'tign and tff'
-!write(*,'(5f7.3)')d
-!write(*,*)' '
-d(1:5:2,1:5:2)=lfn9(icl-1:icl+1,1:3)
-d(2:4,2:4)=lff
-!write(*,*)'lfn and lff'
-!write(*,'(5f7.3)')d
-!write(*,*)' '
+                                    tign9,lfn9,tff9,lff9)
 
 fuel_time_cell= 8.235294 ;
         write(*,*)"Calculation of fuel_frac over 4 subcells"
@@ -417,15 +415,29 @@ program test
    do j=1,2
 
 
-result6=fuel_left_cell_3(  &
-    lff(i,j),lff(i,j+1),lff(i+1,j),lff(i+1,j+1), &
-    tff(i,j),tff(i,j+1),tff(i+1,j),tff(i+1,j+1),&
+result(i,j)=fuel_left_cell_3(  &
+    lff9(i,j),lff9(i,j+1),lff9(i+1,j),lff9(i+1,j+1), &
+    tff9(i,j),tff9(i,j+1),tff9(i+1,j),tff9(i+1,j+1),&
     time_now, fuel_time_cell)
-write(*,*)'result after step i and j',result,i,j
+write(*,*)'result after step i and j',result(i,j),i,j
 
 end do
 end do
+
+if (res1-result(1,1)>error.or.res2-result(1,2)>error.or.res3-result(2,1)>error.or.res4-result(2,2)>error) then
+	print*,'Test #10 is not working, one of errors '
+        print*,'For cell (1,1)',res1-result(1,1)
+        print*,'For cell (1,2)',res2-result(1,2)
+        print*,'For cell (2,1)',res3-result(2,1)
+        print*,'For cell (2,2)',res4-result(2,2)
+        print*,'is greater than 0.0001',error
+        flag=0 
+end if
+
 end do
+if (flag.eq.1) then 
+    num_tests=num_tests+1;
+endif
 print*,"If the results for each correspondent subcell are equal in all cells"
 print*,"then the code is consistent"
 

wrfv2_fire/phys/module_fr_sfire_core.F

@@ -263,7 +263,6 @@ subroutine fuel_left( &
 real,dimension(3,3)::tff,lff
 ! help variables instead of arrays fuel_left_f and fire_area_f 
 real::lffij,lffi1j,lffij1,lffi1j1,tifij,tifi1j,tifij1,tifi1j1,tx,ty,txx,tyy
-! variables for calculation instead of lff(i,j) and tif(i,j)is lffij,tifij etc..#define IFCELLS (ite-its+1)*fuel_left_irl
 #define JFCELLS (jte-jts+1)*fuel_left_jrl
 character(len=128)::msg
 integer::m,omp_get_thread_num
@@ -306,25 +305,14 @@ subroutine fuel_left( &
 ite1=min(ite,ifde-1)
 jts1=max(jts,jfds+1)
 jte1=min(jte,jfde-1)
-!write(*,*)"fuel_left_method",fuel_left_method
 
 do icl=its1,ite1
   do jcl=jts1,jte1
     helpsum1=0
     helpsum2=0
-
-
 !   Loop over subcells in cell #(icl,jcl)
-!      write(*,*)"ifds,ifde",ifds,ifde
-!      write(*,*)"jfds,jfde",jfds,jfde
-!      write(*,*)"its,ite,jts,jte",its,ite,jts,jte
-!      write(*,*)"icl,jcl",icl,jcl
-
-
     call tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, &
                                   tign,lfn,tff,lff)
-!endif
-
     do isubcl=1,ir
       do jsubcl=1,jr 
 !we use 4 points of each subcell in fuel_left_cell_1
@@ -340,7 +328,6 @@ subroutine fuel_left( &
      lff(isubcl,jsubcl),lff(isubcl,jsubcl+1),lff(isubcl+1,jsubcl),lff(isubcl+1,jsubcl+1), &
      tff(isubcl,jsubcl),tff(isubcl,jsubcl+1),tff(isubcl+1,jsubcl),tff(isubcl+1,jsubcl+1), &
      time_now, fuel_time(icl,jcl)) 
-! dont forget to change fire_area_ff here
         else
           call crash('fuel_left: unknown fuel_left_method')
         endif
@@ -504,24 +491,26 @@ end subroutine tign_lfn_interpolation
 
 
 subroutine tign_line_interp(tign1,tign2,lfn1,lfn2,lfn_subcl,tign_subcl,time_now)
+! Description: subroutine that is calculating tign of subcell that lies on the
+! line between 2 cells, lfn_subcl is interpolated linearly look at tign_lfn_interpolation
 
 real,intent(in) :: tign1,tign2,lfn1,lfn2,lfn_subcl,time_now
 real,intent(out) :: tign_subcl
 real :: a,b,c,err
 err=0.1
 
-!      write(*,*)"lfn1,lfn2,tign1,tign2,lfn_subcl,timenow",lfn1,lfn2,tign1,tign2,lfn_subcl,time_now
 
 if((lfn1.le.0).AND.(lfn2.le.0)) then
+! both lfn1,2<0, so we do linear interpolation
    tign_subcl=0.5*(tign1+tign2)
 elseif((lfn1.gt.0).AND.(lfn2.gt.0))then
+! both lfn1,2>0, so we do linear interpolation, so it is equal to time_now
     if  ((abs(tign1-time_now).gt.err).or.(abs(tign2-time_now).gt.err)) then
-      write(*,*)"lfn1,lfn2,tign1,tign2,timenow",lfn1,lfn2,tign1,tign2,time_now
       call crash('line_interp: when lfn1(2)>0 their tign1(2) should = time_now')
    else
       tign_subcl=time_now;
    endif
-elseif(lfn_subcl.gt.0) then
+elseif(lfn_subcl.gt.0) then !So the subcell is not on fire
    if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err)) then
       call crash('tign_line_interp one of tign1,2 should be equal time_now')
    else
@@ -555,7 +544,9 @@ end subroutine tign_line_interp
 
 subroutine tign_four_pnts_interp(tign1,tign2,tign3,tign4,lfn1,lfn2,  &
 lfn3,lfn4,lfn_subcl,tign_subcl,time_now)
-
+! Description: subroutine that is calculating tign of subcell that lies in the
+! middle of the square with 4 cells in the endpoints, lfn_subcl is interpolated linearly
+! look at tign_lfn_interpolation
 real,intent(in) :: tign1,tign2,tign3,tign4
 real,intent(in) :: lfn1,lfn2,lfn3,lfn4,lfn_subcl,time_now
 real,intent(out) :: tign_subcl
@@ -564,14 +555,16 @@ subroutine tign_four_pnts_interp(tign1,tign2,tign3,tign4,lfn1,lfn2,  &
 err=0.0001
 
 if((lfn1.le.0).AND.(lfn2.le.0).AND.(lfn3.le.0).AND.(lfn4.le.0)) then
+! all 4 cells:i:1,4 lfn(i)<0, so we do linear interpolation
    tign_subcl=0.25*(tign1+tign2+tign3+tign4)
 elseif((lfn1.gt.0).AND.(lfn2.gt.0).AND.(lfn3.gt.0).AND.(lfn4.gt.0))then
+! all 4 cells:i:1,4 lfn(i)>0, not on fire so tign_subcl=time_now
    !if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err).OR.(abs(tign3-time_now).gt.err).OR.(abs(tign4-time_now).gt.err)) then
    !   call crash('tign_four_pnts_interp: when lfn1(2,3,4)>0 their tign1(2,3,4) should = time_now')
    !else
       tign_subcl=time_now;
    !endif
-elseif(lfn_subcl.gt.0) then
+elseif(lfn_subcl.gt.0) then  ! so the subcell is not on fire
   ! if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err).OR.(abs(tign3-time_now).gt.err).OR.(abs(tign4-time_now).gt.err)) then
   !    call crash('tign_line_interp one of tign1(2,3,4) should be equal time_now')
   ! else
@@ -583,6 +576,7 @@ subroutine tign_four_pnts_interp(tign1,tign2,tign3,tign4,lfn1,lfn2,  &
 ! tign_subcl~=c*lfn_subcl+time_now;
 a=0; 
 b=0;
+!Calculating weights
    if (lfn1.le.0) then
 !       if (tign1.gt.time_now)
 !           call crash('tign_four_pnts_interp tign1 should be less then time_now');
@@ -936,14 +930,9 @@ real function fuel_left_cell_3(  &
 external dggglm
 !executable statements
 
-! a very crude approximation - replace by a better code
-! call check_mesh_2dim(ids,ide+1,jds,jde+1,ims,ime,jms,jme)
 
 dx=1
 dy=1
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!comment - changed tign-time_now to time_now-tign
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 t00=time_now-tign00
 if(lfn00>=0. .or. t00<0.)t00=0.
 t01=time_now-tign01
@@ -954,14 +943,10 @@ real function fuel_left_cell_3(  &
 if(lfn11>=0. .or. t11<0.)t11=0.
 !*** case0 Do nothing
 if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then
-!print*,"Case 0"
     out = 1.0 !  fuel_left, no burning
-!*** case4 all four coners are burning
+!*** case4 all four corners are burning
 else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then
-!!!!!!!!!!!
-
-!print*,"Case 4"
-! All burning
+! All burning, use Finite Differences
 ! T=u(1)*x+u(2)*y+u(3)
 ! t(0,0)=tign(1,1)
 ! t(0,fd(2))=tign(1,2)
@@ -977,41 +962,28 @@ real function fuel_left_cell_3(  &
 ! u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy)
 
  tign_middle=(t00+t01+t10+t11)/4
-! write(*,*)"tign_middle",tign_middle
-! write(*,*)"t00,t01,t10,t11",t00,t01,t10,t11
 
  ! since mesh_size is 1 we replace fd(1) and fd(2) by 1
  dt_dx=(t10-t00+t11-t01)/2
  dt_dy=(t01-t00+t11-t10)/2
 
-! write(*,*)"dt_dx,dt_dy",dt_dx,dt_dy
-! write(*,*)"dx,dy",dx,dy
-
  u(1)=dt_dx
  u(2)=dt_dy
  u(3)=tign_middle-(dt_dx+dt_dy)/2
 
-!Write(*,*)"coefficients of the plane u=",u(1),u(2),u(3)
-
-    ! integrate
-    u(1)=-u(1)/fuel_time_cell
-    u(2)=-u(2)/fuel_time_cell
-    u(3)=-u(3)/fuel_time_cell
-!write(*,*)"u/fuel_time_cell",u(1),u(2),u(3)
-
-!write(*,*)"intexp(u(1)),intexp(u(2))",intexp(u(1)*dx),intexp(u(2)*dy)
-    !fuel_burn(i,j)=1-exp(u(3))*intexp(u(1)*dx)*intexp(u(2)*dy)
+ ! integrate
+ u(1)=-u(1)/fuel_time_cell
+ u(2)=-u(2)/fuel_time_cell
+ u(3)=-u(3)/fuel_time_cell
     s1=u(1)
     s2=u(2)            
     out=1-exp(u(3))*intexp(s1)*intexp(s2)
-!    print *,intexp(s1),intexp(s2),out
     if ( out<0 .or. out>1.0 ) then
         print *,'case4, out should be between 0 and 1'
     end if
-!*** case 1,2,3
+!*** case 1,2,3- other cases
 
 else
-!print*,"Case 123"
     ! set xx, yy for the coner points
     ! move these values out of i and j loop to speed up
     xx(1) = -0.5
@@ -1065,33 +1037,26 @@ real function fuel_left_cell_3(  &
     xylist(npoint+1,2)=xylist(1,2)
 
 
-
-    !print *,'after LS in case3'pproximation of the plane for lfn using finite
-    !differences
+!approximation of the plane for lfn using finite differences, look at "case 4" explanations
 ! approximate L(x,y)=u(1)*x+u(2)*y+u(3)
  lfn_middle=(lfn00+lfn01+lfn10+lfn11)/4
-! print *,"lfn_middle",lfn_middle
  dt_dx=(lfn10-lfn00+lfn11-lfn01)/2
  dt_dy=(lfn01-lfn00+lfn11-lfn10)/2
-!print *,"dt_dx,dt_dy",dt_dx,dt_dy
  u(1)=dt_dx
  u(2)=dt_dy
  u(3)=lfn_middle-(dt_dx+dt_dy)/2
 ! finding the coefficient c, reminder we work over one subcell only
 ! T(x,y)=c*L(x,y)+time_now
-!write(*,*)"vector u before c",u(1),u(2),u(3)
     a=0
     b=0
-
+!Calculating weights
     if (lfn00 <= 0) then
                         a=a+lfn00*lfn00
                         if (t00 < 0) then
                         call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
                         else
                         b=b+t00*lfn00
                         end if
-                     !   write(*,*)"t00,time_now,lfn00",t00,time_now,lfn00
-                     !   write(*,*)"a,b",a,b
     end if
 
 
@@ -1102,7 +1067,6 @@ real function fuel_left_cell_3(  &
                         else
                         b=b+t01*lfn01
                         end if
-                     !   write(*,*)"a,b",a,b
     end if
 
 
@@ -1113,7 +1077,6 @@ real function fuel_left_cell_3(  &
                         else
                         b=b+t10*lfn10
                         end if
-                      !  write(*,*)"a,b",a,b
   end if
 
     if (lfn11<=0) then
@@ -1123,7 +1086,6 @@ real function fuel_left_cell_3(  &
                         else
                         b=b+t11*lfn11
                         end if
-                      !  write(*,*)"a,b",a,b
     end if
 
 
@@ -1136,12 +1098,6 @@ real function fuel_left_cell_3(  &
                         u(2)=u(2)*c
                         u(3)=u(3)*c
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-    !print *,'vecX from LS',vecX
-    !print *,'tign inputed',tign00,tign10,tign11,tign01
-
-!write(*,*)"vector u",u(1),u(2),u(3)
     ! rotate to gradient on x only
     nr = sqrt(u(1)**2+u(2)**2)
     if(.not.nr.gt.eps)then