standalone/atm.F

   1 module module_model_test
   2 use module_fr_sfire_util, only: message,crash
   3 implicit none
   4
   5 contains
   6
   7 subroutine sub_main
   8
   9 !*** purpose: standalone driver with compatible files to WRF-Fire
  10
  11 use module_fr_sfire_driver, only: sfire_driver_em, use_atm_vars
  12 use module_domain, only: domain
  13 use module_configure, only: grid_config_rec_type,read_namelist
  14 use wrf_netcdf, only : grid_info, set_info_from_file, &
  15                        create_output_file,write_vars,output_vars, &
  16                        input_vars_fire,read_vars_fire
  17 implicit none
  18
  19 !*** local
  20
  21 ! arguments to SFIRE
  22
  23 type(domain)::grid          ! all: state+inputs+outputs, compatible with wrf
  24 TYPE (grid_config_rec_type):: config_flags ! the namelist
  25 integer::  &                ! atmosphere mesh dimensions, for compatibility
  26     ids,ide, kds,kde, jds,jde,&
  27     ims,ime, kms,kme, jms,jme,&
  28     ips,ipe, kps,kpe, jps,jpe
  29 integer::  &                ! fire mesh dimensions
  30     ifds,ifde,jfds,jfde, &  ! the physical domain
  31     ifps,ifpe,jfps,jfpe, &  ! patch - assigned to one process. Here the same as domain.
  32     ifms,ifme,jfms,jfme     ! memory allocated, needs a strip around the patch
  33
  34 ! I/O interface
  35 character(len=*),parameter::inputfile='fire_input.nc'
  36 character(len=*),parameter::outputfile='fire_output.nc'
  37 type(output_vars)::output                ! output arrays
  38 type(input_vars_fire)::input                  ! input arrays_fire
  39
  40 ! other derived types
  41 type(grid_info)::info                    ! dimensions, grid controls
  42
  43 ! scalars
  44 integer:: nx,ny,nz,nfx,nfy,nfz,nsteps,itimestep,sr_x,sr_y
  45 integer::nhalo=5
  46 double precision:: dt, duration_s  ! may need more accurate time computation to get the number of timesteps right
  47 real:: time,time_step_start, dts
  48 logical::do_ouput
  49
  50 !*** executable
  51
  52 call read_namelist(config_flags)           ! read flags from namelist.input
  53 call set_info_from_file(inputfile,info)     ! get dimensions
  54
  55 ! set dimensions
  56 !nx=info%natmx
  57 !ny=info%natmy
  58 !nz=info%natmz
  59 nfx=info%nfirex
  60 nfy=info%nfirey
  61 sr_x=info%sr_x
  62 sr_y=info%sr_y
  63
  64 !write(6,2)' atmospheric domain size               ',config_flags%e_we,config_flags%e_sn,config_flags%e_vert
  65 !if(config_flags%e_we.ne.nx+1 .or. config_flags%e_sn.ne.ny+1 .or. config_flags%e_vert.ne.nz+1)then
  66 !  write(6,*)'dimensions in input file are ',nx,ny,nz,' must be 1 less than in namelist.input'
  67 !  call crash('inconsistent input files')
  68 !endif
  69
  70
  71 ! set atm domain size
  72 !ids=1
  73 !ide=nx+1
  74 !jds=1
  75 !jde=ny+1
  76 !kds=1
  77 !kde=nz+1
  78 !ims=ids-nhalo
  79 !ime=ide+nhalo
  80 !jms=jds-nhalo
  81 !jme=jde+nhalo
  82 !kms=kds
  83 !kme=kde
  84 !ips=ids
  85 !ipe=ide
  86 !jps=jds
  87 !jpe=jde
  88 !kps=kds
  89 !kpe=kde
  90 !
  91 !write(6,2)'fire_mesh refinement ratio  ',info%sr_x,info%sr_y
  92
  93 write(6,2)'fire domain size            ',nfx,nfy
  94
  95 !if(nfx.ne.sr_x*nx .or. nfy.ne.sr_y*ny)then
  96 !   call crash('fire domain size is not multiple of atmospheric domain size')
  97 !endif
  98
  99
 100 ! set fire domain size
 101 ifds=1
 102 ifde=nfx+sr_x
 103 jfds=1
 104 jfde=nfy+sr_y
 105 ifms=ifds-sr_x*nhalo
 106 ifme=ifde+sr_x*nhalo
 107 jfms=ifds-sr_y*nhalo
 108 jfme=ifde+sr_y*nhalo
 109 ifps=1
 110 ifpe=nfx
 111 jfps=1
 112 jfpe=nfy
 113
 114 !write(6,2)'atmospheric domain dimensions',ids,ide,jds,jde,kds,kde
 115 !write(6,2)'atmospheric memory dimensions',ims,ime,jms,jme,kms,kme
 116 write(6,2)'fire domain dimensions       ',ifds,ifde,jfds,jfde
 117 write(6,2)'fire memory dimensions       ',ifms,ifme,jfms,jfme
 118 2 format(a,6i6)
 119
 120 ! allocate
 121
 122 ! inputs
 123 call allocate2d(grid%uf,ifms,ifme,jfms,jfme,'uf')              ! fire winds
 124 call allocate2d(grid%vf,ifms,ifme,jfms,jfme,'vf')              ! fire winds
 125 call allocate2d(grid%zsf,ifms,ifme,jfms,jfme,'zsf')             ! terrain height
 126 call allocate2d(grid%dzdxf,ifms,ifme,jfms,jfme,'dzdxf')           ! terrain grad
 127 call allocate2d(grid%dzdyf,ifms,ifme,jfms,jfme,'dzdyf')           ! terrain grad
 128 call allocate2d(grid%fxlong,ifms,ifme,jfms,jfme,'fxlong')          !
 129 call allocate2d(grid%fxlat,ifms,ifme,jfms,jfme,'fxlat')           !
 130 call allocate2d(grid%nfuel_cat,ifms,ifme,jfms,jfme,'nfuel_cat')          !
 131
 132 ! state
 133 call allocate2d(grid%bbb,ifms,ifme,jfms,jfme,'bbb')             ! spread formula coeff
 134 call allocate2d(grid%betafl,ifms,ifme,jfms,jfme,'betafl')          ! spread formula coeff
 135 call allocate2d(grid%phiwc,ifms,ifme,jfms,jfme,'phiwc')           ! spread formula coeff
 136 call allocate2d(grid%r_0,ifms,ifme,jfms,jfme,'r_0')             ! spread formula coeff
 137 call allocate2d(grid%fgip,ifms,ifme,jfms,jfme,'fgip')            ! spread formula coeff
 138 call allocate2d(grid%ischap,ifms,ifme,jfms,jfme,'ischap')          ! spread formula coeff
 139 call allocate2d(grid%fuel_time,ifms,ifme,jfms,jfme,'fuel_time')        !
 140 call allocate2d(grid%lfn,ifms,ifme,jfms,jfme,'lfn')
 141 call allocate2d(grid%tign_g,ifms,ifme,jfms,jfme,'tign')
 142 call allocate2d(grid%fuel_frac,ifms,ifme,jfms,jfme,'fuel_frac')
 143 call allocate2d(grid%fuel_frac_burnt,ifms,ifme,jfms,jfme,'fuel_frac_burnt')
 144 call allocate2d(grid%fire_area,ifms,ifme,jfms,jfme,'fire_area')
 145 call allocate2d(grid%lfn_out,ifms,ifme,jfms,jfme,'lfn_out')
 146 call allocate2d(grid%ros,ifms,ifme,jfms,jfme,'ros')
 147
 148 ! output
 149 call allocate2d(grid%fgrnhfx,ifms,ifme,jfms,jfme,'fgrnhfx')          !
 150 call allocate2d(grid%fgrnqfx,ifms,ifme,jfms,jfme,'fgrnqfx')          !
 151 call allocate2d(grid%fcanhfx,ifms,ifme,jfms,jfme,'fcanhfx')          !
 152 call allocate2d(grid%fcanqfx,ifms,ifme,jfms,jfme,'fcanqfx')          !
 153
 154 ! atmosphere input/output compatibility arrays
 155 !call allocate2d(grid%avg_fuel_frac,ims,ime,jms,jme,'avg_fuel_frac')
 156 !call allocate2d(grid%xlong,ims,ime,jms,jme,'xlong')
 157 !call allocate2d(grid%xlat,ims,ime,jms,jme,'xlat')
 158 !call allocate2d(grid%ht,ims,ime,jms,jme,'ht')
 159 !call allocate2d(grid%z0,ims,ime,jms,jme,'z0')
 160 !call allocate2d(grid%grnhfx,ims,ime,jms,jme,'grnhfx')
 161 !call allocate2d(grid%grnqfx,ims,ime,jms,jme,'qrnqfx')
 162 !call allocate2d(grid%canhfx,ims,ime,jms,jme,'canhfx')
 163 !call allocate2d(grid%canqfx,ims,ime,jms,jme,'canqfx')
 164 !call allocate3d(grid%ph_2,ims,ime,kms,kme,jms,jme,'ph_2')
 165 !call allocate3d(grid%phb,ims,ime,kms,kme,jms,jme,'phb')
 166 !call allocate3d(grid%u_2,ims,ime,kms,kme,jms,jme,'u_2')
 167 !call allocate3d(grid%v_2,ims,ime,kms,kme,jms,jme,'v_2')
 168 !call allocate2d(grid%uah,ims,ime,jms,jme,'uah')                  !
 169 !call allocate2d(grid%vah,ims,ime,jms,jme,'vah')                  !
 170
 171 ! prepare reading input file
 172 input%nfuel_cat=>grid%nfuel_cat(1:nfx,1:nfy)
 173 input%zsf      =>grid%zsf(1:nfx,1:nfy)
 174 input%dzdxf    =>grid%dzdxf(1:nfx,1:nfy)
 175 input%dzdyf    =>grid%dzdyf(1:nfx,1:nfy)
 176 input%vf       =>grid%uf(1:nfx,1:nfy)
 177 input%uf       =>grid%vf(1:nfx,1:nfy)
 178 input%lfn      =>grid%lfn(1:nfx,1:nfy)
 179 input%tign     =>grid%tign_g(1:nfx,1:nfy)
 180 input%fuel_frac   =>grid%fuel_frac(1:nfx,1:nfy)
 181 !call allocate3d(input%ph  ,1,nx,1,ny,1,nz+1,'ph')
 182 !call allocate3d(input%phb, 1,nx,1,ny,1,nz+1,'phb')
 183 !call allocate3d(input%u,   1,nx+1,1,ny,1,nz,'u')
 184 !call allocate3d(input%v,   1,nx,1,ny+1,1,nz,'v')
 185
 186 call read_vars_fire(inputfile,info,input)
 187
 188 ! sfire uses wrf ordering of dimensions
 189 !grid%ph_2(1:nx,1:nz+1,1:ny)=reshape(input%ph, (/nx,nz+1,ny/),order=(/1,3,2/))
 190 !grid%phb (1:nx,1:nz+1,1:ny)=reshape(input%phb,(/nx,nz+1,ny/),order=(/1,3,2/))
 191 !deallocate(input%ph)
 192 !deallocate(input%phb)
 193 !grid%u_2(1:nx+1,1:nz,1:ny)=reshape(input%u,(/nx+1,nz,ny/),order=(/1,3,2/))
 194 !grid%v_2(1:nx,1:nz,1:ny+1)=reshape(input%v,(/nx,nz,ny+1/),order=(/1,3,2/))
 195 !deallocate(input%u)
 196 !deallocate(input%v)
 197 !grid%z0(1:nx,1:ny)=input%z0
 198 !deallocate(input%z0)
 199
 200 ! time control
 201 ! NOTE: dt in the netcdf input file as returned in info%dt is WRONG !!
 202 dt=config_flags%time_step
 203 if(config_flags%time_step_fract_den.ne.0)then
 204   dt=dt+dble(config_flags%time_step_fract_num)/dble(config_flags%time_step_fract_den)
 205 endif
 206 duration_s = config_flags%run_seconds           &
 207            + 60d0*(config_flags%run_minutes     &
 208            + 60d0*(config_flags%run_hours       &
 209            + 24d0*(config_flags%run_days)))
 210 nsteps = nint( duration_s / dt ) ! number of time steps
 211
 212 ! divide up for shared memory parallel execution
 213 call set_tiles(1,1,ips,ipe,jps,jpe,grid%num_tiles,grid%i_start,grid%i_end,grid%j_start,grid%j_end)
 214
 215 ! set the scalars in grid type
 216 grid%sr_x = sr_x
 217 grid%sr_y = sr_y
 218 grid%dt = dt
 219 grid%dx = info%dx
 220 grid%dy = info%dy
 221 grid%itimestep=0
 222 grid%u_frame=0.
 223 grid%v_frame=0.
 224
 225 ! start output file
 226 call create_output_file(outputfile,info)
 227
 228 ! initialize model
 229    time_step_start  =0.
 230    dts=dt
 231 call sfire_driver_em ( grid , config_flags                          &
 232             ,time_step_start,dts                                    &
 233             ,1,2,0                                                 &
 234             ,ids,ide, kds,kde, jds,jde                              &
 235             ,ims,ime, kms,kme, jms,jme                              &
 236             ,ips,ipe, kps,kpe, jps,jpe                              &
 237             ,ifds,ifde, jfds,jfde                                   &
 238             ,ifms,ifme, jfms,jfme                                   &
 239             ,ifps,ifpe, jfps,jfpe )
 240
 241 do itimestep=1,nsteps
 242 ! run  model
 243    grid%itimestep = itimestep
 244    time_step_start = itimestep*dt
 245    dts=dt
 246    call sfire_driver_em ( grid , config_flags                          &
 247             ,time_step_start,dts                                    &
 248             ,3,6,0                                                  &
 249             ,ids,ide, kds,kde, jds,jde                              &
 250             ,ims,ime, kms,kme, jms,jme                              &
 251             ,ips,ipe, kps,kpe, jps,jpe                              &
 252             ,ifds,ifde, jfds,jfde                                   &
 253             ,ifms,ifme, jfms,jfme                                   &
 254             ,ifps,ifpe, jfps,jfpe )
 255
 256     if(itimestep.le.10.or.mod(itimestep,10).eq.0)then
 257         time = dt*itimestep
 258         output%lfn=>grid%lfn(ifps:ifpe,jfps:jfpe)
 259         output%tign=>grid%tign_g(ifps:ifpe,jfps:jfpe)
 260         output%fgrnhfx=>grid%fgrnhfx(ifps:ifpe,jfps:jfpe)
 261         call write_vars(outputfile,output,time)
 262     endif
 263 enddo
 264
 265 end subroutine sub_main
 266
 267 !
 268 !******************************
 269 !
 270
 271 subroutine set_tiles(itiles,jtiles,ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 272 !*** set tiles for standalone/testing
 273 implicit none
 274 !*** arguments
 275 integer,intent(in)::itiles,jtiles,ips,ipe,jps,jpe
 276 integer,intent(out)::num_tiles
 277 integer,intent(out),dimension(itiles*jtiles)::i_start,i_end,j_start,j_end
 278 !*** local
 279 integer::i,j,istep,jstep,ij
 280 character(len=128)::msg
 281 write(msg,1)'patch',ips,':',ipe,jps,':',jpe
 282 1 format(a,5x,i6,a,2i6,a,i6)
 283 call message(msg,level=-1)
 284 if(ips.ge.ipe.or.jps.ge.jpe)call crash('bad domain bounds')
 285 num_tiles=itiles*jtiles
 286 istep=(ipe-ips+itiles)/itiles
 287 jstep=(jpe-jps+jtiles)/jtiles
 288 do i=1,itiles
 289     do j=1,jtiles
 290         ij=j+(i-1)*jtiles
 291         i_start(ij)=min(ipe,ips+(i-1)*istep)
 292         i_end(ij)  =min(ipe,ips+(i  )*istep-1)
 293         j_start(ij)=min(jpe,jps+(j-1)*jstep)
 294         j_end(ij)  =min(jpe,jps+(j  )*jstep-1)
 295     enddo
 296 enddo
 297 call check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 298 end subroutine set_tiles
 299
 300
 301 subroutine check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 302 implicit none
 303 !*** purpose: check if tiles fit
 304 !*** arguments
 305 integer,intent(in)::ips,ipe,jps,jpe,num_tiles
 306 integer,intent(in),dimension(num_tiles)::i_start,i_end,j_start,j_end
 307 !*** local
 308 character(len=128)::msg
 309 integer:: ij,ie
 310 !*** executable
 311 if(num_tiles.lt.1)call crash('check_tiles: need at least one tile')
 312 ie=0
 313 do ij=1,num_tiles
 314     if(i_start(ij).lt.ips.or.i_end(ij).gt.ipe &
 315     .or.j_start(ij).lt.jps.or.j_end(ij).gt.jpe)then
 316         write(msg,1)'patch',ips,':',ipe,jps,':',jpe
 317 1       format(a,5x,i6,a,2i6,a,i6)
 318         call message(msg,level=-1)
 319         write(msg,2)'tile',ij,i_start(ij),':',i_end(ij),j_start(ij),':',j_end(ij)
 320 2       format(a,2i6,a,2i6,a,i6)
 321         call message(msg,level=-1)
 322         call crash('bad tile bounds')
 323     endif
 324 enddo
 325 end subroutine check_tiles
 326
 327
 328 subroutine allocate2d(p,ims,ime,jms,jme,s)
 329 !*** allocate a pointer with error checking and initialization
 330 implicit none
 331 !*** arguments
 332 real, pointer, intent(out), dimension(:,:)::p
 333 integer, intent(in):: ims,ime,jms,jme
 334 character(len=*),intent(in)::s
 335 !*** local
 336 integer::err
 337 !*** executable
 338 write(6,1) ims,ime,jms,jme,trim(s)
 339 if(associated(p))call crash('already allocated')
 340 1 format('allocate2d',2(1x,i6,' :',i6),1x,a)
 341 allocate(p(ims:ime,jms:jme),stat=err)
 342 if(err.ne.0)then
 343    write(6,1)ims,ime,jms,jme,trim(s)
 344    call crash('memory allocation failed')
 345 endif
 346 p=0.
 347 end subroutine allocate2d
 348
 349 subroutine allocate3d(p,ims,ime,jms,jme,kms,kme,s)
 350 !*** allocate a pointer with error checking and initialization
 351 implicit none
 352 !*** arguments
 353 real, pointer, intent(out), dimension(:,:,:)::p
 354 integer, intent(in):: ims,ime,jms,jme,kms,kme
 355 character(len=*),intent(in)::s
 356 !*** local
 357 integer::err
 358 !*** executable
 359 write(6,1) ims,ime,jms,jme,kms,kme,trim(s)
 360 1 format('allocate3d',3(1x,i6,' :',i6),1x,a)
 361 if(associated(p))call crash('already allocated')
 362 allocate(p(ims:ime,jms:jme,kms:kme),stat=err)
 363 if(err.ne.0)then
 364    write(6,1)ims,ime,jms,jme,kms,kme,trim(s)
 365    call crash('memory allocation failed')
 366 endif
 367 p=0.
 368 end subroutine allocate3d
 369
 370
 371 end module module_model_test
 372
 373 !
 374 !******************************
 375 !
 376
 377
 378 program model_test_main
 379 use module_model_test, only: sub_main
 380 call  sub_main
 381 end program model_test_main