femwind/fortran/femwind_wrfout.f90

   1 program femwind_wrfout
   2
   3 use module_femwind
   4 use module_utils
   5 use module_common
   6 use module_netcdf
   7 use module_wrfout
   8
   9 implicit none
  10
  11 integer::write_debug = 0,call_femwind = 0
  12
  13 integer ::                          &
  14     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  15     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  16     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  17     ifts, ifte, kfts, kfte, jfts, jfte                          ! fire tile bounds
  18
  19 ! declarations of A, msize included from module femwind
  20
  21 ! variables read from wrfout
  22 real, pointer:: u0_fmw(:,:,:), v0_fmw(:,:,:), w0_fmw(:,:,:), zsf(:,:), ht_fmw(:)
  23
  24 ! variables allocated and computed here
  25 real, pointer:: u0(:,:,:), v0(:,:,:), w0(:,:,:)
  26 real, pointer:: u(:,:,:), v(:,:,:), w(:,:,:), uf(:,:), vf(:,:)
  27
  28 integer:: i, j, k, n(3), nx, ny, nz
  29 integer::ncid,frame,sr(2),frame0_fmw,mframe=100,dims(3)
  30 real:: rate,A(3,3),dx,dy,zx,zy
  31 character(len=256)::filename, msg
  32
  33 !*** executable
  34
  35 filename = "wrf.nc" ! file to read from and write to
  36 frame = 1           ! frame in the file
  37
  38 !************************************
  39 ! read static data
  40 !************************************
  41
  42 call ncopen(filename,nf90_nowrite,ncid)
  43
  44 ! horizontal height at cell centers
  45 call netcdf_read_array_wrf(ncid,"ZSF",frame=frame,a2d=zsf)
  46
  47 ! height of the vertical layers above the terrain
  48 call netcdf_read_array_wrf(ncid,"HT_FMW",frame=frame,a1d=ht_fmw)
  49
  50 !  fire subgrid refinement ratio, 3d grid dimensions cell centered
  51 call get_wrf_dims(ncid,sr,dims)
  52 ! horizontal mesh spacing
  53 dx = netcdf_read_att(ncid,"DX")/sr(1)
  54 dy = netcdf_read_att(ncid,"DY")/sr(2)
  55
  56 call ncclose(ncid)
  57
  58 !************************************
  59 ! process static data
  60 !************************************
  61
  62 !construct mg data
  63
  64 params%A = reshape((/1.0,0.0,0.0, &
  65                      0.0,1.0,0.0, &
  66                      0.0,0.0,1.0 /), &
  67                    (/3,3/))
  68
  69 ! finest level 1
  70 mg(1)%nx = dims(1) + 1  ! dimensions are vertex centered
  71 mg(1)%ny = dims(2) + 1  ! dimensions are vertex centered
  72 mg(1)%nz = dims(3) + 1  ! dimensions are vertex centered
  73
  74 call get_mg_dims(mg(1), &                 !  set bounds compatible with WRF
  75     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
  76     ifms, ifme, kfms,kfme, jfms, jfme,            &
  77     ifps, ifpe, kfps,kfpe, jfps, jfpe,           & ! fire patch bounds
  78     ifts, ifte, kfts,kfte, jfts, jfte)
  79
  80 ! allocations persistent over the whole run
  81
  82 allocate(mg(1)%X(ifms:ifme,kfms:kfme,jfms:jfme))
  83 allocate(mg(1)%Y(ifms:ifme,kfms:kfme,jfms:jfme))
  84 allocate(mg(1)%Z(ifms:ifme,kfms:kfme,jfms:jfme))
  85 allocate(u0(ifms:ifme,kfms:kfme,jfms:jfme)) ! vector components called u v W
  86 allocate(v0(ifms:ifme,kfms:kfme,jfms:jfme))
  87 allocate(w0(ifms:ifme,kfms:kfme,jfms:jfme))
  88 allocate(u(ifms:ifme,kfms:kfme,jfms:jfme)) ! vector components called u v W
  89 allocate(v(ifms:ifme,kfms:kfme,jfms:jfme))
  90 allocate(w(ifms:ifme,kfms:kfme,jfms:jfme))
  91 allocate(uf(ifts:ifte,jfts:jfte))
  92 allocate(vf(ifts:ifte,jfts:jfte))
  93   write(msg,*)"shape(uf)=",shape(uf)
  94   call message(msg)
  95   write(msg,*)"shape(vf)=",shape(vf)
  96   call message(msg)
  97
  98
  99 ! X Y Z are corner based, upper bound larger by one
 100 ! first interpolate/extrapolate zsf to corners
 101 do j=jfts,jfte+1
 102   do i=ifts,ifte+1
 103     if(i.eq.ifts)then
 104       zx = zsf(i,j) - 0.5*(zsf(i+1,j)-zsf(i,j))  ! extrapolate down
 105     elseif(i.eq.ifte+1)then
 106       zx = zsf(i,j) - 0.5*(zsf(i-1,j)-zsf(i,j))  ! extrapolate up
 107     else
 108       zx = 0.5*(zsf(i,j) + zsf(i+1,j))           ! interpolate
 109     endif
 110     if(j.eq.jfts)then
 111       zy = zsf(i,j) - 0.5*(zsf(i,j+1)-zsf(i,j))  ! extrapolate down
 112     elseif(j.eq.jfte+1)then
 113       zy = zsf(i,j) - 0.5*(zsf(i,j-1)-zsf(i,j))  ! extrapolate up
 114     else
 115       zy = 0.5*(zsf(i,j) + zsf(i,j+1))           ! interpolate
 116     endif
 117     mg(1)%Z(i,kfts,j) = 0.5*(zx + zy)
 118   enddo
 119 enddo
 120
 121 ! (X,Y) is the same uniform grid on all levels
 122 do j=jfts,jfte+1
 123   do k=kfts,kfte+1
 124     do i=ifts,ifte+1
 125         mg(1)%X(i,k,j) = (i-1)*dx
 126         mg(1)%Y(i,k,j) = (j-1)*dy
 127     enddo
 128   enddo
 129 enddo
 130
 131 ! add top of levels above terrain to Z(:,kfts,:)
 132 print *,'debug ht_fmw=',ht_fmw
 133 do j=jfts,jfte+1
 134   do k=kfts+1,kfte+1
 135     do i=ifts,ifte+1
 136         mg(1)%Z(i,k,j) = mg(1)%Z(i,kfts,j) + ht_fmw(k - kfts)
 137     enddo
 138   enddo
 139 enddo
 140
 141 ! mesh spacing
 142 mg(1)%dx = dx
 143 mg(1)%dy = dy
 144 allocate(mg(1)%dz(mg(1)%nz-1))
 145 do k=kfds,kfte
 146     mg(1)%dz(k)=mg(1)%Z(1,k+1,1)-mg(1)%Z(1,k,1)
 147 enddo
 148
 149 deallocate(ht_fmw,zsf)  ! no longer needed
 150
 151 if(call_femwind.gt.0)then
 152 call message('calling femwind_setup')
 153 call femwind_setup(mg)
 154 call message('femwind_setup done')
 155 endif
 156
 157 if(write_debug.gt.0)then
 158   ! write coordinates for debug/display, even if the caller knows
 159   call write_average_to_center(mg(1)%X,'X_c')
 160   call write_average_to_center(mg(1)%Y,'Y_c')
 161   call write_average_to_center(mg(1)%Z,'Z_c')
 162 endif
 163
 164 !************************************
 165 ! loop over time steps
 166 !************************************
 167
 168 do frame0_fmw=1,mframe
 169
 170   ! read initial velocity field, loop until delivered
 171   call read_initial_wind(filename,u0_fmw,v0_fmw,w0_fmw,frame0_fmw,frame=1)
 172
 173   ! copy the input data to tile sized bounds
 174   ! initial wind is at cell centers, indexing was already switched to ikj in reading
 175   u0(ifts:ifte,kfts:kfte,jfts:jfte) = u0_fmw
 176   v0(ifts:ifte,kfts:kfte,jfts:jfte) = v0_fmw
 177   w0(ifts:ifte,kfts:kfte,jfts:jfte) = w0_fmw
 178
 179   ! save memory while solver is running
 180   deallocate(u0_fmw,v0_fmw,w0_fmw)
 181
 182   ! compute/update mass consistent flow
 183
 184 if(call_femwind.gt.0)then
 185   write(*,'(a)')'calling femwind_solve'
 186   call femwind_solve(  mg,&
 187     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 188     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 189     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 190     ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 191     u0, v0, w0,                                  & ! input arrays
 192     u, v, w,                                                  & ! output arrays
 193     rate)
 194   write(*,'(a)')'femwind_solve returned OK'
 195 else
 196   ! interpolate to fire height
 197   call message('TESTING ONLY: copying lowest level of input to uf vf')
 198   write(msg,*)"shape(uf)=",shape(uf)
 199   call message(msg)
 200   write(msg,*)"shape(vf)=",shape(vf)
 201   call message(msg)
 202   uf = u0(ifts:ifte,1,jfts:jfte)
 203   vf = v0(ifts:ifte,1,jfts:jfte)
 204   call write_fire_wind(filename,uf,vf,frame0_fmw,frame=1)
 205 endif
 206
 207   if(write_debug.gt.0)then
 208     ! write output as is in 3D but with tight dimensions
 209     call write_array(u(ifts:ifte,kfts:kfte,jfts:jfte),'u')
 210     call write_array(v(ifts:ifte,kfts:kfte,jfts:jfte),'v')
 211     call write_array(w(ifts:ifte,kfts:kfte,jfts:jfte),'w')
 212     call write_scalar(rate,'rate')
 213   endif
 214
 215 enddo
 216
 217 contains
 218
 219 subroutine write_average_to_center(F,name)
 220 ! average and return with tight bounds
 221 implicit none
 222 !*** arguments
 223 real, intent(in)::F(ifms:ifme, kfms:kfme, jfms:jfme)
 224 character(len=*)::name
 225 !*** local
 226 real, allocatable::Fm(:,:,:)
 227 integer::ii,kk,jj
 228 real::s
 229 !*** executable
 230 allocate(Fm(ifts:ifte,kfts:kfte,jfts:jfte))
 231 do j=jfts,jfte
 232   do k=kfts,kfte
 233     do i=ifts,ifte
 234       s = 0.
 235       do ii=0,1
 236         do kk=0,1
 237           do jj=0,1
 238             s = s + F(i+ii,k+kk,j+jj)
 239           enddo
 240         enddo
 241       enddo
 242       Fm(i,k,j) = s/8
 243     enddo
 244   enddo
 245 enddo
 246 call write_array(Fm,name)
 247 end subroutine write_average_to_center
 248
 249 end program femwind_wrfout