standalone/old/model_test_main.F

   1 ! the main program is at the end because of fortran limitations
   2
   3 module module_model_test
   4 use module_fr_sfire_util
   5 use module_fr_sfire_model
   6 use module_fr_sfire_phys
   7
   8
   9 contains
  10
  11 subroutine main_sub
  12 implicit none
  13
  14 !*** declarations
  15
  16 integer:: nx,ny, msteps  ! problem dimension, in cells, number of steps
  17 real:: fdx,fdy, dt,wind,alpha       ! fire mesh spacing (m), time step (s), wind azimuth
  18 integer::rfx,rfy ! refinement, convenience only
  19
  20 nx=400
  21 ny=400
  22 msteps=200
  23 msteps=6
  24 msteps=100
  25 fdx=6
  26 fdy=6
  27
  28 rfx=1
  29 rfy=1
  30
  31 nx=nx*rfx
  32 ny=ny*rfy
  33 fdx=fdx/rfx
  34 fdy=fdy/rfy
  35
  36 dt=0.5
  37 wind=10
  38 alpha=0.0
  39
  40 print *,'nx=',nx,' ny=',ny,' msteps=',msteps,' fdx=',fdx,' fdy=',fdy,' dt=',dt,' wind=',wind,' alpha=',alpha
  41 call model_test(1,nx,1,ny, &
  42 !   -1,nx+3,-2,ny+2, &
  43    0,nx+1,0,ny+1, &
  44    1,nx,1,ny, &
  45    fdx,fdy,wind,alpha,dt,msteps)
  46
  47 end subroutine main_sub
  48
  49 !
  50 !****************
  51 !
  52
  53
  54 subroutine model_test(  &
  55     ifds,ifde,jfds,jfde, &
  56     ifms,ifme,jfms,jfme, &
  57     ifps,ifpe,jfps,jfpe, &
  58     fdx,fdy,wind,alpha,dt,msteps)
  59 implicit none
  60
  61 !*** arguments
  62 integer, intent(in):: &
  63     ifds,ifde,jfds,jfde, &
  64     ifps,ifpe,jfps,jfpe, &
  65     ifms,ifme,jfms,jfme, msteps
  66 real, intent(in)::fdx,fdy,dt,wind,alpha
  67
  68 !*** local
  69 real, dimension(ifms:ifme,jfms:jfme):: zsf,     &
  70                  lfn,tign,fuel_frac,fire_area,                    &
  71                  grnhfx,grnqfx,lfn_out
  72 integer::  num_ignitions,i,j,ifuelread,istep,nfuel_cat0,ifun,ifun_start
  73 logical::need_lfn_update
  74 real:: t0,time_start,sm,sn
  75 integer, dimension(ifms:ifme,jfms:jfme)::nfuel_cat,ischap
  76 real, dimension(ifms:ifme,jfms:jfme)::fuel_time,vx,vy,dzfsdx,dzfsdy,bbb,betafl,phiwc,r_0,fgip, &
  77        coord_xf,coord_yf
  78 real:: unit_xf,unit_yf
  79 integer, parameter::max_tiles=10
  80 integer::num_tiles,ij,ifts,ifte,jfts,jfte
  81 integer, dimension(max_tiles)::if_start,if_end,jf_start,jf_end
  82 integer, parameter :: max_ignitions=10
  83 real, dimension(max_ignitions) :: ignition_start_x,ignition_start_y, &
  84     ignition_end_x,ignition_end_y,ignition_radius,ignition_time
  85
  86 !*** executable
  87
  88 call set_tiles(2,2,ifps,ifpe,jfps,jfpe,num_tiles,if_start,if_end,jf_start,jf_end)
  89
  90 print *,'mesh size in cells: ',ifps,ifpe,jfps,jfpe
  91 print *,'array allocation:   ',ifms,ifme,jfms,jfme
  92
  93 t0=0                                        ! starting time, arbitrary
  94
  95 ! populate the arrays somehow
  96 do j=jfps,jfpe
  97     do i=ifps,ifpe
  98         zsf(i,j)=1000   ! flat ground
  99         vx(i,j)=wind*cos(alpha)    ! constant wind
 100         vy(i,j)=wind*sin(alpha)
 101     enddo
 102 enddo
 103
 104 ! fuel data
 105 ifuelread=0
 106 nfuel_cat0=3
 107
 108 ! for matlab
 109 open(1,file='model_test_out.txt',form='formatted')
 110 1   format(e25.12e3)
 111 sm=ifpe-ifps+1
 112 sn=jfpe-jfps+1
 113 write(1,1)1.,1.,fdx
 114 write(1,1)1.,1.,fdy
 115
 116 time_start=t0
 117 num_ignitions=2
 118 ignition_start_x(1)=0.5*fdx*(ifde-ifds)
 119 ignition_start_y(1)=0.5*fdy*(jfde-jfds)
 120 ignition_end_x(1)=0.5*fdx*(ifde-ifds)*0.9999999
 121 ignition_end_y(1)=0.5*fdy*(jfde-jfds)*1.0000001
 122 ignition_radius(1) = 0.5*max(5.0,6*max(fdx,fdy))
 123 ignition_time(1)=1
 124 ignition_start_x(2)=1000
 125 ignition_start_y(2)=500
 126 ignition_end_x(2)=1500
 127 ignition_end_y(2)=1500
 128 ! at least 6 by 6 cells but no less than 5 m
 129 ignition_radius(2) = 0.5*max(5.0,6*max(fdx,fdy))
 130 ignition_time(2)=2
 131
 132
 133 unit_xf=1
 134 unit_yf=1
 135
 136 do istep=1,msteps
 137     ifun_start=1
 138     if(istep.ne.1)ifun_start=3
 139     do ifun=ifun_start,6
 140 !OMP    PARALLEL DO PRIVATE(ij,ifts,ifte,jfts,jfte)
 141         do ij=1,num_tiles
 142             ifts= if_start(ij)
 143             ifte= if_end(ij)
 144             jfts= jf_start(ij)
 145             jfte= jf_end(ij)
 146
 147
 148             call set_ideal_coord( fdx,fdy, &
 149                 ifds,ifde,jfds,jfde,  &
 150                 ifms,ifme,jfms,jfme,  &
 151                 ifts,ifte,jfts,jfte,  &
 152                 coord_xf,coord_yf     &
 153             )
 154
 155             call   sfire_model (10*istep+ij,ifun,         &
 156                 need_lfn_update, num_ignitions,              &
 157                 ifuelread,nfuel_cat0,                   &
 158                 ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
 159                 ifms,ifme,jfms,jfme,                    & ! fire memory dims - how declared
 160                 ifds,ifde,jfds,jfde,                    & ! patch=domain
 161                 ifts,ifte,jfts,jfte,                    & ! fire tile dims  - this thread
 162                 time_start,dt,                          & ! time and increment
 163                 fdx,fdy,                                & ! fire mesh spacing
 164                 ignition_start_x,ignition_start_y,      &
 165                 ignition_end_x,ignition_end_y,          &
 166                 ignition_radius,                        &
 167                 ignition_time,                          &
 168                 coord_xf,coord_yf,unit_xf,unit_yf,      & ! fire mesh coordinates
 169                 zsf,                                    & ! terrain height (for gradient)
 170                 vx,vy,                                  & ! input: wind
 171                 lfn,lfn_out,tign,fuel_frac,fire_area,   & ! state: level function, ign time, fuel left
 172                 grnhfx,grnqfx,                          & ! output: heat fluxes
 173                 nfuel_cat,                              & ! fuel data per point
 174                 fuel_time,                              & ! save derived internal data
 175                 bbb,betafl,phiwc,r_0,fgip,ischap &
 176             )
 177         enddo
 178
 179     enddo !OMP PARALLEL
 180
 181     if(istep.le.10.or.mod(istep,10).eq.0)then
 182         write(1,1)1.,1.,time_start
 183         write(1,1)sm,sn,((lfn(i,j),i=ifps,ifpe),j=jfps,jfpe)
 184         write(1,1)sm,sn,((tign(i,j),i=ifps,ifpe),j=jfps,jfpe)
 185         write(1,1)sm,sn,((vx(i,j),i=ifps,ifpe),j=jfps,jfpe)
 186         write(1,1)sm,sn,((vy(i,j),i=ifps,ifpe),j=jfps,jfpe)
 187         write(1,1)sm,sn,((grnhfx(i,j),i=ifps,ifpe),j=jfps,jfpe)
 188     endif
 189     print *,'test_main: step ',istep,' of ',msteps,' time ',time_start
 190     time_start=time_start+dt
 191 enddo
 192
 193 close(1)
 194
 195 end subroutine model_test
 196
 197 !
 198 !******************************
 199 !
 200
 201 subroutine set_tiles(itiles,jtiles,ids,ide,jds,jde,num_tiles,i_start,i_end,j_start,j_end)
 202 !*** set tiles for standalone/testing
 203 implicit none
 204 !*** arguments
 205 integer,intent(in)::itiles,jtiles,ids,ide,jds,jde
 206 integer,intent(out)::num_tiles
 207 integer,intent(out),dimension(itiles*jtiles)::i_start,i_end,j_start,j_end
 208 !*** local
 209 integer::i,j,istep,jstep,ij
 210 num_tiles=itiles*jtiles
 211 istep=(ide-ids+itiles)/itiles
 212 jstep=(jde-jds+jtiles)/jtiles
 213 do i=1,itiles
 214     do j=1,jtiles
 215         ij=j+(i-1)*jtiles
 216         i_start(ij)=min(ide,ids+(i-1)*istep)
 217         i_end(ij)  =min(ide,ids+(i  )*istep-1)
 218         j_start(ij)=min(jde,jds+(j-1)*jstep)
 219         j_end(ij)  =min(jde,jds+(j  )*jstep-1)
 220     enddo
 221 enddo
 222 call check_tiles(ids,ide,jds,jde,num_tiles,i_start,i_end,j_start,j_end)
 223 end subroutine set_tiles
 224
 225
 226 subroutine check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 227 implicit none
 228 !*** purpose: check if tiles fit
 229 !*** arguments
 230 integer,intent(in)::ips,ipe,jps,jpe,num_tiles
 231 integer,intent(in),dimension(num_tiles)::i_start,i_end,j_start,j_end
 232 !*** local
 233 character(len=128)::msg
 234 integer:: ij,ie
 235 !*** executable
 236 if(num_tiles.lt.1)call crash('check_tiles: need at least one tile')
 237 ie=0
 238 if (num_tiles.eq.1) then
 239     if(i_start(1).ne.ips.or.i_end(1).ne.ipe.or.j_start(1).ne.jps.or.j_end(1).ne.jpe)ie=1
 240 else
 241     do ij=1,num_tiles
 242         if(i_start(ij).lt.ips.or.i_end(ij).gt.ipe &
 243         .or.j_start(ij).lt.jps.or.j_end(ij).gt.jpe)ie=ij
 244     enddo
 245 endif
 246 if(ie.ne.0)then
 247     write(msg,*)'bad tile ',ie
 248     call message(msg)
 249     write(msg,*)'patch dimensions:',ips,ipe,jps,jpe
 250     call message(msg)
 251     do ij=1,num_tiles
 252         write(msg,*)'tile',ij,i_start(ij),i_end(ij),j_start(ij),j_end(ij)
 253         call message(msg)
 254     enddo
 255     call crash('bad tile bounds')
 256 endif
 257 end subroutine check_tiles
 258
 259
 260 end module module_model_test
 261
 262 !
 263 !******************************
 264 !
 265
 266
 267 program model_test_main
 268 use module_model_test
 269 call  main_sub
 270 end program model_test_main