wrfv2_fire/phys/model_netcdf_main.F

   1 ! the main program is at the end because of fortran limitations
   2 #ifndef NF_TYPE
   3 #define NF_TYPE NF90_FLOAT
   4 #endif
   5
   6 module jb_utils
   7 contains
   8
   9 subroutine check(status)
  10 use netcdf
  11   integer, intent ( in) :: status
  12
  13   if(status /= nf90_noerr) then
  14     print *, nf90_strerror(status)
  15     stop "Stopped"
  16   end if
  17 end subroutine check
  18
  19 end module jb_utils
  20
  21
  22 module module_model_test
  23 use module_fr_sfire_util
  24 use module_fr_sfire_model
  25 use module_fr_sfire_phys
  26 use jb_utils
  27 use netcdf
  28
  29 contains
  30
  31 subroutine main_sub
  32 implicit none
  33
  34 !*** declarations
  35
  36 integer:: nx,ny, msteps  ! problem dimension, in cells, number of steps
  37 real:: fdx,fdy, dt,wind,alpha       ! fire mesh spacing (m), time step (s), wind azimuth
  38 integer::rfx,rfy ! refinement, convenience only
  39
  40 nx=400
  41 ny=400
  42 msteps=200
  43 msteps=6
  44 ! msteps=100
  45 fdx=6
  46 fdy=6
  47
  48 rfx=1
  49 rfy=1
  50
  51 nx=nx*rfx
  52 ny=ny*rfy
  53 fdx=fdx/rfx
  54 fdy=fdy/rfy
  55
  56 dt=0.5
  57 wind=10
  58 alpha=0.0
  59
  60 print *,'nx=',nx,' ny=',ny,' msteps=',msteps,' fdx=',fdx,' fdy=',fdy,' dt=',dt,' wind=',wind,' alpha=',alpha
  61 call model_test(1,nx,1,ny, &
  62 !   -1,nx+3,-2,ny+2, &
  63    0,nx+1,0,ny+1, &
  64    1,nx,1,ny, &
  65    fdx,fdy,wind,alpha,dt,msteps)
  66
  67 end subroutine main_sub
  68
  69 !
  70 !****************
  71 !
  72
  73
  74 subroutine model_test(  &
  75     ifds,ifde,jfds,jfde, &
  76     ifms,ifme,jfms,jfme, &
  77     ifps,ifpe,jfps,jfpe, &
  78     fdx,fdy,wind,alpha,dt,msteps)
  79 implicit none
  80
  81 !*** arguments
  82 integer, intent(in):: &
  83     ifds,ifde,jfds,jfde, &
  84     ifps,ifpe,jfps,jfpe, &
  85     ifms,ifme,jfms,jfme, msteps
  86 real, intent(in)::fdx,fdy,dt,wind,alpha
  87
  88 !*** local
  89 real, dimension(ifms:ifme,jfms:jfme):: zsf,     &
  90                  lfn,tign,fuel_frac,                    &
  91                  grnhfx,grnqfx,lfn_out
  92 integer::  num_ignitions,i,j,ifuelread,istep,nfuel_cat0,ifun,ifun_start
  93 logical::need_lfn_update
  94 real:: t0,time_start,sm,sn
  95 integer, dimension(ifms:ifme,jfms:jfme)::nfuel_cat,ischap
  96 real, dimension(ifms:ifme,jfms:jfme)::fuel_time,vx,vy,dzfsdx,dzfsdy,bbb,betafl,phiwc,r_0,fgip
  97 integer, parameter::max_tiles=10
  98 integer::num_tiles,ij,ifts,ifte,jfts,jfte
  99 integer, dimension(max_tiles)::if_start,if_end,jf_start,jf_end
 100 integer, parameter :: max_ignitions=10
 101 real, dimension(max_ignitions) :: ignition_start_x,ignition_start_y, &
 102     ignition_end_x,ignition_end_y,ignition_radius,ignition_time
 103 character (len=NF90_MAX_NAME), parameter :: iofile = "fire_io.nc"
 104 integer :: ncid,xid,yid,tid,lfnid,vfxid,vfyid,tignid,grnhftid,itmp,tdim,tvalsid
 105 integer :: xsid,ysid  ! staggered dims
 106 integer :: ostat,stepnum,grnqftid,fuel_fracid,zsfid,fueltimeid
 107 real :: tvals
 108 integer,parameter :: ndims = 3
 109 integer,dimension(ndims):: dimids,count,start,atmp,counts,starts
 110 integer::sstep
 111
 112 !*** executable
 113
 114 call set_tiles(2,2,ifps,ifpe,jfps,jfpe,num_tiles,if_start,if_end,jf_start,jf_end)
 115
 116 print *,'mesh size in cells: ',ifps,ifpe,jfps,jfpe
 117 print *,'array allocation:   ',ifms,ifme,jfms,jfme
 118
 119 t0=0                                        ! starting time, arbitrary
 120 ostat=nf90_open(iofile,NF90_WRITE,ncid)
 121 print*,ostat,nf90_noerr
 122 if(ostat.ne.NF90_NOERR)then
 123
 124
 125 ! populate the arrays somehow
 126 do j=jfps,jfpe
 127     do i=ifps,ifpe
 128         zsf(i,j)=1000   ! flat ground
 129         vx(i,j)=wind*cos(alpha)    ! constant wind
 130         vy(i,j)=wind*sin(alpha)
 131     enddo
 132 enddo
 133 call check( nf90_create(iofile,NF90_NOCLOBBER,ncid) )
 134 call check( nf90_def_dim(ncid,"x",ifpe-ifps+1,xid) )
 135 call check( nf90_def_dim(ncid,"y",jfpe-jfps+1,yid) )
 136 call check( nf90_def_dim(ncid,"x_stag",ifpe-ifps+2,xsid) )
 137 call check( nf90_def_dim(ncid,"y_stag",jfpe-jfps+2,ysid) )
 138 call check( nf90_def_dim(ncid,"Time",NF90_UNLIMITED,tid) )
 139 call check( nf90_put_att(ncid,NF90_GLOBAL,"DX",fdx) )
 140 call check( nf90_put_att(ncid,NF90_GLOBAL,"DY",fdy) )
 141 stepnum=0
 142 call check( nf90_put_att(ncid,NF90_GLOBAL,"stepnum",stepnum) )
 143
 144 dimids=(/xsid,ysid,tid/)
 145 call check( nf90_def_var(ncid,"lfn",NF_TYPE,dimids,lfnid) )
 146 call check( nf90_def_var(ncid,"vx",NF_TYPE,dimids,vfxid) )
 147 call check( nf90_def_var(ncid,"vy",NF_TYPE,dimids,vfyid) )
 148 call check( nf90_def_var(ncid,"Times",NF_TYPE,(/tid/),tvalsid) )
 149
 150 dimids=(/xid,yid,tid/)
 151 call check( nf90_def_var(ncid,"tign",NF_TYPE,dimids,tignid) )
 152 call check( nf90_def_var(ncid,"grnhfx",NF_TYPE,dimids,grnhftid) )
 153 call check( nf90_def_var(ncid,"grnqfx",NF_TYPE,dimids,grnqftid) )
 154 call check( nf90_def_var(ncid,"fuelfrac",NF_TYPE,dimids,fuel_fracid) )
 155 call check( nf90_def_var(ncid,"fueltime",NF_TYPE,dimids,fueltimeid) )
 156 call check( nf90_enddef(ncid) )
 157 print*,'askjdfksdjf'
 158 tvals=0.0
 159 tdim=0
 160
 161 else
 162
 163 call check( nf90_inq_dimid(ncid,"x_stag",xsid) )
 164 call check( nf90_inquire_dimension(ncid,xsid,len=itmp) )
 165 if(itmp.ne.ifpe-ifps+2)call check(NF90_EVARSIZE)
 166 call check( nf90_inq_dimid(ncid,"y_stag",ysid) )
 167 call check( nf90_inquire_dimension(ncid,ysid,len=itmp) )
 168 if(itmp.ne.jfpe-jfps+2)call check(NF90_EVARSIZE)
 169
 170 call check( nf90_inq_dimid(ncid,"x",xid) )
 171 call check( nf90_inquire_dimension(ncid,xid,len=itmp) )
 172 if(itmp.ne.ifpe-ifps+1)call check(NF90_EVARSIZE)
 173 call check( nf90_inq_dimid(ncid,"y",yid) )
 174 call check( nf90_inquire_dimension(ncid,yid,len=itmp) )
 175 if(itmp.ne.jfpe-jfps+1)call check(NF90_EVARSIZE)
 176 call check( nf90_inq_dimid(ncid,"Time",tid) )
 177 call check( nf90_inquire(ncid,unlimitedDimID=itmp) )
 178 if(itmp.ne.tid)call check(NF90_EVARSIZE)
 179
 180 call check( nf90_inq_varid(ncid,"lfn",lfnid) )
 181 call check( nf90_inq_varid(ncid,"vx",vfxid) )
 182 call check( nf90_inq_varid(ncid,"vy",vfyid) )
 183 call check( nf90_inq_varid(ncid,"tign",tignid) )
 184 call check( nf90_inq_varid(ncid,"grnhfx",grnhftid) )
 185 call check( nf90_inq_varid(ncid,"grnqfx",grnqftid) )
 186 call check( nf90_inq_varid(ncid,"fuelfrac",fuel_fracid) )
 187 call check( nf90_inq_varid(ncid,"fueltime",fueltimeid) )
 188 call check( nf90_inquire_dimension(ncid,tid,len=tdim) )
 189 call check( nf90_inq_varid(ncid,"Times",tvalsid) )
 190 call check( nf90_get_var(ncid,tvalsid,atmp,start=(/tdim/),count=(/1/)) )
 191 tvals=atmp(1)
 192 dimids=(/xid,yid,tid/)
 193 call check( nf90_get_att(ncid,nf90_global,"stepnum",stepnum) )
 194
 195 count=(/ifpe-ifps+1,jfpe-jfps+1,1/)
 196 counts=(/ifpe-ifps+2,jfpe-jfps+2,1/)
 197 start=(/1,1,tdim/)
 198 starts=(/1,1,tdim/)
 199
 200 #define _SD (ifps:ifpe+1,jfps:jfpe+1)
 201 #define _D  (ifps:ifpe,jfps:jfpe)
 202 call check( nf90_get_var(ncid,lfnid,lfn _SD,start=starts,count=counts) )
 203 call check( nf90_get_var(ncid,tignid,tign _D,start=starts,count=count) )
 204 call check( nf90_get_var(ncid,vfxid,vx _SD,start=starts,count=counts) )
 205 call check( nf90_get_var(ncid,vfyid,vy _SD,start=starts,count=counts) )
 206 call check( nf90_get_var(ncid,grnhftid,grnhfx _D,start=start,count=count) )
 207 call check( nf90_get_var(ncid,grnqftid,grnqfx _D,start=start,count=count) )
 208 call check( nf90_get_var(ncid,fuel_fracid,fuel_frac _D,start=start,count=count) )
 209 call check( nf90_get_var(ncid,fueltimeid,fuel_time _D,start=start,count=count) )
 210 endif
 211
 212
 213 count=(/ifpe-ifps+1,jfpe-jfps+1,1/)
 214 counts=(/ifpe-ifps+2,jfpe-jfps+2,1/)
 215
 216 ! fuel data
 217 ifuelread=0
 218 nfuel_cat0=3
 219
 220 ! for matlab
 221 open(1,file='model_test_out.txt',form='formatted')
 222 1   format(e25.12e3)
 223 sm=ifpe-ifps+1
 224 sn=jfpe-jfps+1
 225 write(1,1)1.,1.,fdx
 226 write(1,1)1.,1.,fdy
 227
 228 time_start=t0
 229 num_ignitions=2
 230 ignition_start_x(1)=0.5*fdx*(ifde-ifds)
 231 ignition_start_y(1)=0.5*fdy*(jfde-jfds)
 232 ignition_end_x(1)=0.5*fdx*(ifde-ifds)*0.9999999
 233 ignition_end_y(1)=0.5*fdy*(jfde-jfds)*1.0000001
 234 ignition_radius(1) = 0.5*max(5.0,6*max(fdx,fdy))
 235 ignition_time(1)=1
 236 ignition_start_x(2)=1000
 237 ignition_start_y(2)=500
 238 ignition_end_x(2)=1500
 239 ignition_end_y(2)=1500
 240 ! at least 6 by 6 cells but no less than 5 m
 241 ignition_radius(2) = 0.5*max(5.0,6*max(fdx,fdy))
 242 ignition_time(2)=2
 243
 244 t0=tvals
 245 time_start=t0
 246 stepnum=stepnum+1
 247 sstep=stepnum
 248
 249 do istep=sstep,msteps
 250     ifun_start=1
 251     if(istep.ne.1)ifun_start=3
 252     do ifun=ifun_start,6
 253 !OMP    PARALLEL DO PRIVATE(ij,ifts,ifte,jfts,jfte)
 254         do ij=1,num_tiles
 255             ifts= if_start(ij)
 256             ifte= if_end(ij)
 257             jfts= jf_start(ij)
 258             jfte= jf_end(ij)
 259
 260             call   sfire_model (10*istep+ij,ifun,         &
 261                 need_lfn_update, num_ignitions,              &
 262                 ifuelread,nfuel_cat0,                   &
 263                 ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
 264                 ifms,ifme,jfms,jfme,                    & ! fire memory dims - how declared
 265                 ifts,ifte,jfts,jfte,                    & ! fire tile dims  - this thread
 266                 time_start,dt,                          & ! time and increment
 267                 fdx,fdy,                                & ! fire mesh spacing
 268                 2,                                      & ! boundary_guard
 269                 ignition_start_x,ignition_start_y,      &
 270                 ignition_end_x,ignition_end_y,          &
 271                 ignition_radius,                        &
 272                 ignition_time,                          &
 273                 zsf,                                    & ! terrain height (for gradient)
 274                 vx,vy,                                  & ! input: wind
 275                 lfn,lfn_out,tign,fuel_frac,                     & ! state: level function, ign time, fuel left
 276                 grnhfx,grnqfx,                          & ! output: heat fluxes
 277                 nfuel_cat,                              & ! fuel data per point
 278                 fuel_time,                              & ! save derived internal data
 279                 dzfsdx,dzfsdy,bbb,betafl,phiwc,r_0,fgip,ischap &
 280             )
 281 tdim=tdim+1
 282 start=(/1,1,tdim/)
 283 atmp(1)=tvals+dt
 284 call check( &
 285 nf90_put_var(ncid,tvalsid,atmp,start=(/start(3)/),count=(/count(3)/)) &
 286            )
 287 call check( nf90_put_var(ncid,lfnid,lfn _SD,start=start,count=counts) )
 288 call check( nf90_put_var(ncid,tignid,tign _D,start=start,count=count) )
 289 call check( nf90_put_var(ncid,vfxid,vx _SD,start=start,count=counts) )
 290 call check( nf90_put_var(ncid,vfyid,vy _SD,start=start,count=counts) )
 291 call check( nf90_put_var(ncid,grnhftid,grnhfx _D,start=start,count=count) )
 292 call check( nf90_put_var(ncid,grnqftid,grnqfx _D,start=start,count=count) )
 293 call check( nf90_put_var(ncid,fuel_fracid,fuel_frac _D,start=start,count=count) )
 294 call check( nf90_put_var(ncid,tvalsid,atmp,start=(/tdim/),count=(/1/)) )
 295 call check( nf90_put_var(ncid,fueltimeid,fuel_time _D,start=start,count=count) )
 296 call check( nf90_put_att(ncid,nf90_global,"stepnum",istep) )
 297 call check( nf90_sync(ncid) )
 298 time_start=time_start+dt
 299
 300 !close(1)
 301 call check( nf90_close(ncid) )
 302
 303         enddo
 304
 305     enddo !OMP PARALLEL
 306
 307     if(istep.le.10.or.mod(istep,10).eq.0)then
 308         write(1,1)1.,1.,time_start
 309         write(1,1)sm,sn,((lfn(i,j),i=ifps,ifpe),j=jfps,jfpe)
 310         write(1,1)sm,sn,((tign(i,j),i=ifps,ifpe),j=jfps,jfpe)
 311         write(1,1)sm,sn,((vx(i,j),i=ifps,ifpe),j=jfps,jfpe)
 312         write(1,1)sm,sn,((vy(i,j),i=ifps,ifpe),j=jfps,jfpe)
 313         write(1,1)sm,sn,((grnhfx(i,j),i=ifps,ifpe),j=jfps,jfpe)
 314     endif
 315     print *,'test_main: step ',istep,' of ',msteps,' time ',time_start
 316     time_start=time_start+dt
 317 enddo
 318
 319 close(1)
 320
 321 end subroutine model_test
 322
 323 !
 324 !******************************
 325 !
 326
 327 subroutine set_tiles(itiles,jtiles,ids,ide,jds,jde,num_tiles,i_start,i_end,j_start,j_end)
 328 !*** set tiles for standalone/testing
 329 implicit none
 330 !*** arguments
 331 integer,intent(in)::itiles,jtiles,ids,ide,jds,jde
 332 integer,intent(out)::num_tiles
 333 integer,intent(out),dimension(itiles*jtiles)::i_start,i_end,j_start,j_end
 334 !*** local
 335 integer::i,j,istep,jstep,ij
 336 num_tiles=itiles*jtiles
 337 istep=(ide-ids+itiles)/itiles
 338 jstep=(jde-jds+jtiles)/jtiles
 339 do i=1,itiles
 340     do j=1,jtiles
 341         ij=j+(i-1)*jtiles
 342         i_start(ij)=min(ide,ids+(i-1)*istep)
 343         i_end(ij)  =min(ide,ids+(i  )*istep-1)
 344         j_start(ij)=min(jde,jds+(j-1)*jstep)
 345         j_end(ij)  =min(jde,jds+(j  )*jstep-1)
 346     enddo
 347 enddo
 348 call check_tiles(ids,ide,jds,jde,num_tiles,i_start,i_end,j_start,j_end)
 349 end subroutine set_tiles
 350
 351
 352 subroutine check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 353 implicit none
 354 !*** purpose: check if tiles fit
 355 !*** arguments
 356 integer,intent(in)::ips,ipe,jps,jpe,num_tiles
 357 integer,intent(in),dimension(num_tiles)::i_start,i_end,j_start,j_end
 358 !*** local
 359 character(len=128)::msg
 360 integer:: ij,ie
 361 !*** executable
 362 if(num_tiles.lt.1)call crash('check_tiles: need at least one tile')
 363 ie=0
 364 if (num_tiles.eq.1) then
 365     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
 366 else
 367     do ij=1,num_tiles
 368         if(i_start(ij).lt.ips.or.i_end(ij).gt.ipe &
 369         .or.j_start(ij).lt.jps.or.j_end(ij).gt.jpe)ie=ij
 370     enddo
 371 endif
 372 if(ie.ne.0)then
 373     write(msg,*)'bad tile ',ie
 374     call message(msg)
 375     write(msg,*)'patch dimensions:',ips,ipe,jps,jpe
 376     call message(msg)
 377     do ij=1,num_tiles
 378         write(msg,*)'tile',ij,i_start(ij),i_end(ij),j_start(ij),j_end(ij)
 379         call message(msg)
 380     enddo
 381     call crash('bad tile bounds')
 382 endif
 383 end subroutine check_tiles
 384
 385
 386 end module module_model_test
 387
 388 !
 389 !******************************
 390 !
 391
 392
 393 program model_test_main
 394 use module_model_test
 395 call  main_sub
 396 end program model_test_main