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
  12 integer ::                          &
  13     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  14     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  15     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  16     ifts, ifte, kfts, kfte, jfts, jfte                          ! fire tile bounds
  17
  18 ! declarations of A, msize included from module femwind
  19
  20 ! variables read from wrfout
  21 real, pointer:: u0_fmw(:,:,:), v0_fmw(:,:,:), w0_fmw(:,:,:), zsfe(:,:), ht_fmw(:), zsf_fmw(:,:), fwh_fmw(:,:)
  22
  23 ! variables allocated and computed here
  24 real, pointer:: u0(:,:,:), v0(:,:,:), w0(:,:,:), hth(:)
  25 real, pointer:: u(:,:,:), v(:,:,:), w(:,:,:), uf(:,:), vf(:,:), wh(:,:)
  26 integer, pointer:: kh(:,:)
  27
  28 integer:: i, j, k, n(3), nx, ny, nz, nfx, nfy, nfz, iu=10
  29 integer::ncid,sr(2),frame0_fmw,dims(3)
  30 real:: rate,A(3,3),dx,dy,zx,zy
  31 character(len=256)::filename, msg
  32
  33 ! params
  34 integer ::     &
  35   frame = 1,             & ! frame to use in the communication file
  36   frame_terminate = -99, &
  37   write_uvw_fmw = 1,     &
  38   mframe = 10,           & ! send code to terminate and exit after this frame
  39   write_debug = 0,       &
  40   call_femwind = 1,      &
  41   interpolate_uf_vf = 1
  42
  43
  44 namelist /femwind_nl/    &
  45   frame ,                &
  46   frame_terminate ,      &
  47   write_uvw_fmw ,        &
  48   mframe ,               & ! send code to terminate and exit after this frame
  49   write_debug ,          &
  50   call_femwind,          &
  51   interpolate_uf_vf
  52
  53 !*** executable
  54
  55 filename = "wrf.nc" ! file to read from and write to
  56
  57 ! params from namelist
  58 open(iu,file="namelist.femwind")
  59 read(iu, femwind_nl)
  60 close(iu)
  61 write(*,femwind_nl)
  62
  63 !************************************
  64 ! read static data
  65 !************************************
  66
  67 call ncopen(filename,nf90_nowrite,ncid)
  68
  69 !  fire subgrid refinement ratio, 3d grid dimensions cell centered
  70 call get_wrf_dims(ncid,sr,dims)
  71 ! horizontal mesh spacing
  72 dx = netcdf_read_att(ncid,"DX")/sr(1)
  73 dy = netcdf_read_att(ncid,"DY")/sr(2)
  74
  75 ! height of the vertical layers above the terrain
  76 call netcdf_read_array_wrf(ncid,"HT_FMW",frame=frame,a1d=ht_fmw)
  77
  78 ! fire mesh size in cells
  79 nfx = dims(1)*sr(1)
  80 nfy = dims(2)*sr(2)
  81 nfz = size(ht_fmw,1)
  82
  83 ! terrain height at cell centers
  84 call netcdf_read_array_wrf(ncid,"ZSF",frame=frame,a2d=zsf_fmw)
  85
  86 ! fire wind height to interpolate to at cell centers
  87 call netcdf_read_array_wrf(ncid,"FWH",frame=frame,a2d=fwh_fmw)
  88
  89 call ncclose(ncid)
  90
  91 !************************************
  92 ! process static data
  93 !************************************
  94
  95 !construct mg data
  96
  97 params%A = reshape((/1.0,0.0,0.0, &
  98                      0.0,1.0,0.0, &
  99                      0.0,0.0,1.0 /), &
 100                    (/3,3/))
 101
 102 ! finest level 1
 103 mg(1)%nx = nfx + 1  ! dimensions are vertex centered
 104 mg(1)%ny = nfy + 1  ! dimensions are vertex centered
 105 mg(1)%nz = nfz + 1  ! dimensions are vertex centered
 106
 107 call get_mg_dims(mg(1), &                 !  set bounds compatible with WRF
 108     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
 109     ifms, ifme, kfms,kfme, jfms, jfme,            &
 110     ifps, ifpe, kfps,kfpe, jfps, jfpe,           & ! fire patch bounds
 111     ifts, ifte, kfts,kfte, jfts, jfte)
 112
 113 ! allocations persistent over the whole run
 114
 115 allocate(mg(1)%X(ifms:ifme,kfms:kfme,jfms:jfme))
 116 allocate(mg(1)%Y(ifms:ifme,kfms:kfme,jfms:jfme))
 117 allocate(mg(1)%Z(ifms:ifme,kfms:kfme,jfms:jfme))
 118 allocate(u0(ifms:ifme,kfms:kfme,jfms:jfme)) ! vector components called u v W
 119 allocate(v0(ifms:ifme,kfms:kfme,jfms:jfme))
 120 allocate(w0(ifms:ifme,kfms:kfme,jfms:jfme))
 121 allocate(u(ifms:ifme,kfms:kfme,jfms:jfme)) ! vector components called u v W
 122 allocate(v(ifms:ifme,kfms:kfme,jfms:jfme))
 123 allocate(w(ifms:ifme,kfms:kfme,jfms:jfme))
 124 allocate(uf(ifts:ifte,jfts:jfte))
 125 allocate(vf(ifts:ifte,jfts:jfte))
 126 allocate(kh(ifts:ifte,jfts:jfte))
 127 allocate(wh(ifts:ifte,jfts:jfte))
 128 allocate(zsfe(ifts-1:ifte+1,jfts-1:jfte+1))
 129
 130   write(msg,*)"shape(uf)=",shape(uf)
 131   call message(msg)
 132   write(msg,*)"shape(vf)=",shape(vf)
 133   call message(msg)
 134
 135 ! X Y Z are corner based, upper bound larger by one
 136
 137 ! copy interior
 138 zsfe(ifts:ifte,jfts:jfte)=zsf_fmw
 139 ! extrapolate on sides
 140 zsfe(ifts-1,jfts:jfte)=1.5*zsfe(ifts,jfts:jfte)-0.5*zsfe(ifts+1,jfts:jfte)
 141 zsfe(ifte+1,jfts:jfte)=1.5*zsfe(ifte,jfts:jfte)-0.5*zsfe(ifte-1,jfts:jfte)
 142 zsfe(ifts:ifte,jfts-1)=1.5*zsfe(ifts:ifte,jfts)-0.5*zsfe(ifts:ifte,jfts+1)
 143 zsfe(ifts:ifte,jfte+1)=1.5*zsfe(ifts:ifte,jfte)-0.5*zsfe(ifts:ifte,jfte-1)
 144 ! extrapolate at domain corners
 145 zsfe(ifts-1,jfts-1)=1.5*zsfe(ifts,jfts)-0.5*zsfe(ifts+1,jfts+1)
 146 zsfe(ifte+1,jfts-1)=1.5*zsfe(ifte,jfts)-0.5*zsfe(ifte-1,jfts+1)
 147 zsfe(ifts-1,jfte+1)=1.5*zsfe(ifts,jfte)-0.5*zsfe(ifts+1,jfte-1)
 148 zsfe(ifte+1,jfte+1)=1.5*zsfe(ifte,jfte)-0.5*zsfe(ifte-1,jfte-1)
 149 ! interpolate to cell corners
 150 mg(1)%Z(ifts:ifte+1,kfts,jfts:jfte+1) = 0.25 * (    &
 151     zsfe(ifts:ifte+1,jfts:jfte+1) +         &
 152     zsfe(ifts-1:ifte,jfts:jfte+1) +         &
 153     zsfe(ifts:ifte+1,jfts-1:jfte) +         &
 154     zsfe(ifts-1:ifte,jfts-1:jfte) )
 155
 156 ! (X,Y) is the same uniform grid on all levels
 157 do j=jfts,jfte+1
 158   do k=kfts,kfte+1
 159     do i=ifts,ifte+1
 160         mg(1)%X(i,k,j) = (i-1)*dx
 161         mg(1)%Y(i,k,j) = (j-1)*dy
 162     enddo
 163   enddo
 164 enddo
 165
 166 ! add top of levels above terrain to Z(:,kfts,:)
 167 print *,'debug ht_fmw=',ht_fmw
 168 do j=jfts,jfte+1
 169   do k=kfts+1,kfte+1
 170     do i=ifts,ifte+1
 171         mg(1)%Z(i,k,j) = mg(1)%Z(i,kfts,j) + ht_fmw(k - kfts)
 172     enddo
 173   enddo
 174 enddo
 175
 176 ! mesh spacing
 177 mg(1)%dx = dx
 178 mg(1)%dy = dy
 179 allocate(mg(1)%dz(mg(1)%nz-1))
 180 do k=kfds,kfte
 181     mg(1)%dz(k)=mg(1)%Z(1,k+1,1)-mg(1)%Z(1,k,1)
 182 enddo
 183
 184 ! heights of wind nodes in cell centers
 185 allocate(hth(0:kfte))
 186 hth(0)=0
 187 hth(1)=ht_fmw(1)*0.5
 188 do k=2,kfte
 189   hth(k)=0.5*(ht_fmw(k-1)+ht_fmw(k))
 190 enddo
 191
 192 ! indices to interpolate to fwh_fmw
 193 do j=jfts,jfte
 194   do i=ifts,ifte
 195     do k=1,kfte
 196       if (hth(k) .ge. fwh_fmw(i,j))then
 197         kh(i,j)=k
 198         wh(i,j) = (hth(k) - fwh_fmw(i,j))/(hth(k) - hth(k-1)) ! weight 0 if at the top
 199         ! write(msg,*)'i,j=',i,j,'k=',k,'hth=',hth(k),hth(k-1),'fwh=',fwh_fmw(i,j),'wh=',wh(i,j)
 200         ! call message(msg)
 201         goto 1
 202       endif
 203     enddo
 204     call crash('did not find interpolation layer')
 205 1   continue
 206   enddo
 207 enddo
 208
 209 deallocate(hth,ht_fmw,zsfe)  ! no longer needed
 210
 211 if(call_femwind.gt.0)then
 212   call message('calling femwind_setup')
 213   call femwind_setup(mg)
 214   call message('femwind_setup done')
 215 endif
 216
 217 if(write_debug.gt.0)then
 218   ! write coordinates for debug/display, even if the caller knows
 219   call write_average_to_center(mg(1)%X,'X_c')
 220   call write_average_to_center(mg(1)%Y,'Y_c')
 221   call write_average_to_center(mg(1)%Z,'Z_c')
 222 endif
 223
 224 !************************************
 225 ! loop over time steps
 226 !************************************
 227
 228 do frame0_fmw=1,mframe
 229
 230   ! read initial velocity field, loop until delivered
 231   call read_initial_wind(filename,u0_fmw,v0_fmw,w0_fmw,frame0_fmw,frame=frame)
 232
 233   ! copy the input data to tile sized bounds
 234   ! initial wind is at cell centers, indexing was already switched to ikj in reading
 235   u0(ifts:ifte,kfts:kfte,jfts:jfte) = u0_fmw
 236   v0(ifts:ifte,kfts:kfte,jfts:jfte) = v0_fmw
 237   w0(ifts:ifte,kfts:kfte,jfts:jfte) = w0_fmw
 238
 239   ! save memory while solver is running
 240   deallocate(u0_fmw,v0_fmw,w0_fmw)
 241
 242   ! compute/update mass consistent flow
 243
 244   if(call_femwind.gt.0)then
 245     write(*,'(a)')'calling femwind_solve'
 246     call femwind_solve(  mg,&
 247       ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 248       ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 249       ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 250       ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 251       u0, v0, w0,                                  & ! input arrays
 252       u, v, w,                                                  & ! output arrays
 253       rate)
 254     write(*,'(a)')'femwind_solve done'
 255   else
 256     u=u0
 257     v=v0
 258     w=w0
 259     write(*,'(a)')'femwind_solve skipped'
 260   endif
 261
 262   if(interpolate_uf_vf.gt.0)then
 263     ! interpolate u,v to fire height
 264     do j=jfts,jfte
 265       do i=ifts,ifte
 266          k=kh(i,j)
 267          if(k>1)then
 268            uf(i,j) = wh(i,j) * u(i,k,j) + (1.0 - wh(i,j)) * u(i,k-1,j)
 269            vf(i,j) = wh(i,j) * v(i,k,j) + (1.0 - wh(i,j)) * v(i,k-1,j)
 270          else
 271            uf(i,j) = u(i,k,j)
 272            vf(i,j) = v(i,k,j)
 273          endif
 274       enddo
 275     enddo
 276   else
 277     call message('TESTING ONLY: copying lowest level of input to uf vf')
 278     uf = u0(ifts:ifte,1,jfts:jfte)
 279     vf = v0(ifts:ifte,1,jfts:jfte)
 280   endif
 281
 282   !  write_fire_wind(filename,frame0_fmw,uf,vf,u_fmw,v_fmw,w_fmw,frame)
 283   if(write_uvw_fmw.gt.0)then
 284     call write_fire_wind(filename,frame0_fmw,uf,vf, &
 285       u(ifts:ifte,kfts:kfte,jfts:jfte),             &
 286       v(ifts:ifte,kfts:kfte,jfts:jfte),             &
 287       w(ifts:ifte,kfts:kfte,jfts:jfte), frame=frame)
 288   else
 289     call write_fire_wind(filename,frame0_fmw,uf,vf,frame=frame)
 290   endif
 291
 292   if(write_debug.gt.0)then
 293     ! write output as is in 3D but with tight dimensions
 294     call write_array(u(ifts:ifte,kfts:kfte,jfts:jfte),'u')
 295     call write_array(v(ifts:ifte,kfts:kfte,jfts:jfte),'v')
 296     call write_array(w(ifts:ifte,kfts:kfte,jfts:jfte),'w')
 297     call write_scalar(rate,'rate')
 298   endif
 299
 300 enddo
 301
 302 call write_fire_wind(filename,frame_terminate,uf,vf)
 303
 304 print *,'femwind end'
 305
 306 contains
 307
 308 subroutine write_average_to_center(F,name)
 309 ! average and return with tight bounds
 310 implicit none
 311 !*** arguments
 312 real, intent(in)::F(ifms:ifme, kfms:kfme, jfms:jfme)
 313 character(len=*)::name
 314 !*** local
 315 real, allocatable::Fm(:,:,:)
 316 integer::ii,kk,jj
 317 real::s
 318 !*** executable
 319 allocate(Fm(ifts:ifte,kfts:kfte,jfts:jfte))
 320 do j=jfts,jfte
 321   do k=kfts,kfte
 322     do i=ifts,ifte
 323       s = 0.
 324       do ii=0,1
 325         do kk=0,1
 326           do jj=0,1
 327             s = s + F(i+ii,k+kk,j+jj)
 328           enddo
 329         enddo
 330       enddo
 331       Fm(i,k,j) = s/8
 332     enddo
 333   enddo
 334 enddo
 335 call write_array(Fm,name)
 336 end subroutine write_average_to_center
 337
 338 end program femwind_wrfout