Adding horizontal interpolation option

CMakeLists.txt

@@ -74,6 +74,7 @@ list(APPEND _share_files share/proj_lc_mod.F90
 
 # IO
 list(APPEND _io_files io/namelist_mod.F90
+                      io/coupling_mod.F90
                       io/geogrid_mod.F90
                       io/wrf_mod.F90
 		      io/netcdf_mod.F90

io/coupling_mod.F90

@@ -0,0 +1,65 @@
+  module coupling_mod
+
+    use proj_lc_mod, only : proj_lc_t
+    use interp_mod, only : Interp_horizontal_nearest, Interp_horizontal_bilinear, HINTERP_NEAREST, HINTERP_BILINEAR
+    use stderrout_mod, only : Stop_simulation
+
+
+    implicit none
+
+    private
+
+    public :: Interp_horizontal
+
+  contains
+
+    subroutine Interp_horizontal (data_in, proj_data_in, ims, ime, jms, jme, ifms, ifme, jfms, jfme, &
+        num_tiles, i_start, i_end, j_start, j_end, hinterp_opt, lats_out, lons_out, data_out)
+
+      implicit none
+
+      type (proj_lc_t), intent (in) :: proj_data_in
+      integer, intent (in) :: ifms, ifme, jfms, jfme, ims, ime, jms, jme, num_tiles, hinterp_opt
+      integer, dimension (num_tiles), intent (in) :: i_start, i_end, j_start, j_end
+      real, dimension(ims:ime, jms:jme), intent (in) :: data_in
+      real, dimension (ifms:ifme, jfms:jfme), intent (in) :: lats_out, lons_out
+      real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: data_out
+
+      integer :: ij, ifts, ifte, jfts, jfte
+
+
+      Hinterp: select case (hinterp_opt)
+        case (HINTERP_NEAREST)
+          !$OMP PARALLEL DO   &
+          !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte)
+          do ij = 1, num_tiles
+            ifts = i_start(ij)
+            ifte = i_end(ij)
+            jfts = j_start(ij)
+            jfte = j_end(ij)
+            call Interp_horizontal_nearest (data_in, proj_data_in, ims, ime, jms, jme, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, &
+                lats_out, lons_out, data_out)
+          end do
+          !$OMP END PARALLEL DO
+
+        case (HINTERP_BILINEAR)
+          !$OMP PARALLEL DO   &
+          !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte)
+          do ij = 1, num_tiles
+            ifts = i_start(ij)
+            ifte = i_end(ij)
+            jfts = j_start(ij)
+            jfte = j_end(ij)
+            call Interp_horizontal_bilinear (data_in, proj_data_in, ims, ime, jms, jme, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, &
+                lats_out, lons_out, data_out)
+          end do
+          !$OMP END PARALLEL DO
+
+        case default
+          call Stop_simulation ('The horizontal interpolation option selected does not exist.')
+
+      end select Hinterp
+
+    end subroutine Interp_horizontal
+
+  end module coupling_mod

io/namelist_mod.F90

@@ -50,6 +50,7 @@ module namelist_mod
 
       real :: fire_wind_height = 6.096        ! "height of uah,vah wind in fire spread formula" "m"
       integer :: wind_vinterp_opt = 0         ! "wind (adjustment factor) interpolation option"
+      integer :: hinterp_opt = 1              ! "Horizontal interpolation from atm to fire (offline option): 1) ngp, 2)bi-linear"
       logical :: fire_lsm_zcoupling = .false. ! "flag to activate reference velocity at a different height from fire_wind_height"
       real :: fire_lsm_zcoupling_ref = 50.0   ! "reference height from wich u at fire_wind_hegiht is calculated using a logarithmic profile" "m"
 
@@ -191,6 +192,7 @@ subroutine Broadcast_nml (this)
       call Broadcast_integer (this%ros_opt)
       call Broadcast_integer (this%emis_opt)
       call Broadcast_integer (this%wind_vinterp_opt)
+      call Broadcast_integer (this%hinterp_opt)
 
       call Broadcast_integer (this%fire_num_ignitions)
 
@@ -410,6 +412,7 @@ subroutine Init_fire_block (this, file_name)
       real :: fmoist_dt = 600                 ! "moisture model time step" "s"
       real :: fire_wind_height = 6.096        ! "height of uah,vah wind in fire spread formula" "m"
       integer :: wind_vinterp_opt = 0         ! "wind (adjustment factor) interpolation option"
+      integer :: hinterp_opt = 1              ! "Horizontal interpolation from atm to fire (offline option): 1) ngp, 2) bi-linear"
       logical :: fire_is_real_perim = .false. ! .false. = point/line ignition, .true. = observed perimeter"
       real :: frac_fburnt_to_smoke = 0.02     ! "parts per unit of burned fuel becoming smoke " "g_smoke/kg_air"
       real :: fuelmc_g = 0.08                 ! Fuel moisture content ground (Dead FMC)
@@ -463,6 +466,7 @@ subroutine Init_fire_block (this, file_name)
             ! objects
           fuel_opt, ros_opt, fmc_opt, emis_opt, &
           wind_vinterp_opt, &
+          hinterp_opt, &
             ! Ignitions
           fire_num_ignitions, &
             ! Ignition 1
@@ -527,6 +531,7 @@ subroutine Init_fire_block (this, file_name)
       this%fmc_opt = fmc_opt
       this%emis_opt = emis_opt
       this%wind_vinterp_opt = wind_vinterp_opt
+      this%hinterp_opt = hinterp_opt
 
       this%fire_num_ignitions = fire_num_ignitions
 

io/wrf_mod.F90

@@ -7,6 +7,7 @@ module wrf_mod
     use proj_lc_mod, only : proj_lc_t
     use stderrout_mod, only : Print_message, Stop_simulation
     use interp_mod, only : Interp_profile
+    use coupling_mod, only : Interp_horizontal
 
     implicit none
 
@@ -47,7 +48,7 @@ module wrf_mod
       procedure, public :: Get_v3d => Get_meridional_wind_3d
       procedure, public :: Get_v3d_stag => Get_meridional_wind_stag_3d
       procedure, public :: Get_z0 => Get_z0
-      procedure, public :: Interp_var2grid_nearest => Interp_var2grid_nearest
+      procedure, public :: Interp_var2grid => Interp_var2grid
       procedure, public :: Print_domain => Print_domain
       procedure, public :: Update_atm_state => Update_atm_state
     end type wrf_t
@@ -644,39 +645,28 @@ end subroutine Write_latlon_check
 
     end subroutine Get_latcloncs
 
-    subroutine Interp_var2grid_nearest (this, lats_in, lons_in, var_name, vals_out)
-
-      use, intrinsic :: iso_fortran_env, only : ERROR_UNIT
+    subroutine Interp_var2grid (this, lats_out, lons_out, ifms, ifme, jfms, jfme, &
+        num_tiles, i_start, i_end, j_start, j_end, var_name, hinterp_opt, data_out)
 
       implicit none
 
       class (wrf_t), intent(in) :: this
-      real, dimension(:, :), intent(in) :: lats_in, lons_in
+      integer, intent (in) :: ifms, ifme, jfms, jfme, num_tiles
+      integer, dimension (num_tiles), intent (in) :: i_start, i_end, j_start, j_end
+      real, dimension(ifms:ifme, jfms:jfme), intent(in) :: lats_out, lons_out
       character (len = *), intent (in) :: var_name
-      real, dimension(:, :), allocatable, intent(out) :: vals_out
+      integer, intent (in) :: hinterp_opt
+      real, dimension(ifms:ifme, jfms:jfme), intent(in out) :: data_out
 
-      real, parameter :: DEFAULT_INIT = 0.0
-      real, dimension(:, :), allocatable :: ds, var_wrf
-      real :: d, i_real, j_real
+      real, dimension(:, :), allocatable :: var_wrf
       type (proj_lc_t) :: proj
-      integer :: nx, ny, nx_wrf, ny_wrf, i, j, i_wrf, j_wrf
+      integer :: ims, ime, jms, jme
 
 
         ! Init
-      nx_wrf = this%ide - 1
-      ny_wrf = this%jde - 1
-      nx = size (lats_in, dim = 1)
-      ny = size (lats_in, dim = 2)
-
-      if (allocated (vals_out)) deallocate (vals_out)
-      allocate (vals_out(nx, ny))
-      vals_out = DEFAULT_INIT
-
-      allocate (ds(nx, ny))
-      ds = huge (d)
-
       proj = this%Get_projection ()
 
+        ! Get WRF data
       select case (var_name)
         case ('t2')
           var_wrf = this%t2_stag(this%ids:this%ide - 1, this%jds:this%jde - 1)
@@ -700,25 +690,20 @@ subroutine Interp_var2grid_nearest (this, lats_in, lons_in, var_name, vals_out)
           var_wrf = this%va(this%ids:this%ide - 1, this%jds:this%jde - 1)
 
         case default
-          write (ERROR_UNIT, *) 'Unknown variable name to interpolate'
-          stop
+          call Stop_simulation ('Unknown variable name to interpolate')
+
       end select
 
-        ! Algorithm
-      do j = 1, ny
-        do i = 1, nx
-          call proj%Calc_ij (lats_in(i, j), lons_in(i, j), i_real, j_real)
-          i_wrf = min (max (1, nint (i_real)), nx_wrf)
-          j_wrf = min (max (1, nint (j_real)), ny_wrf)
-          d = (this%lats(i_wrf, j_wrf) - lats_in(i, j)) ** 2 + (this%lons(i_wrf, j_wrf) - lons_in(i, j)) ** 2
-          if (d < ds(i, j)) then
-            vals_out(i, j) = var_wrf(i_wrf, j_wrf)
-            ds(i, j) = d
-          end if
-        end do
-      end do
+        ! Interpolate
+      ims = 1
+      ime = size (var_wrf, dim = 1)
+      jms = 1
+      jme = size (var_wrf, dim = 2)
+
+      call Interp_horizontal (var_wrf, proj, ims, ime, jms, jme, ifms, ifme, jfms, jfme, &
+          num_tiles, i_start, i_end, j_start, j_end, hinterp_opt, lats_out, lons_out, data_out)
 
-    end subroutine Interp_var2grid_nearest
+    end subroutine Interp_var2grid
 
     subroutine Interp_wrf2dvar_to_cfbm (wrfatm2dvar, ims, ime, jms, jme, ifms, ifme, jfms, jfme, ifps, ifpe, jfps, jfpe, &
         lats_in, lons_in, proj, vals_out)

share/interp_mod.F90

@@ -1,13 +1,83 @@
   module interp_mod
 
+    use proj_lc_mod, only : proj_lc_t
+
     implicit none
 
     private
 
-    public :: Interp_profile
+    integer, parameter :: HINTERP_NEAREST = 1, HINTERP_BILINEAR = 2
+
+    public :: Interp_profile, Interp_horizontal_nearest, Interp_horizontal_bilinear, HINTERP_NEAREST, HINTERP_BILINEAR
 
   contains
 
+    subroutine Interp_horizontal_nearest (data_in, proj_data_in, ims, ime, jms, jme, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, &
+        lats_out, lons_out, data_out)
+
+    ! Purpose: Nearest neighbor interpolation/extrapolation
+
+      implicit none
+
+      type (proj_lc_t), intent (in) :: proj_data_in
+      integer, intent (in) :: ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, ims, ime, jms, jme
+      real, dimension(ims:ime, jms:jme), intent (in) :: data_in
+      real, dimension (ifms:ifme, jfms:jfme), intent (in) :: lats_out, lons_out
+      real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: data_out
+
+      integer :: i, j, i_in, j_in
+      real :: i_real, j_real
+
+
+      do j = jfts, jfte
+        do i = ifts, ifte
+          call proj_data_in%Calc_ij (lats_out(i, j), lons_out(i, j), i_real, j_real)
+          i_in = min (max (ims, nint (i_real)), ime)
+          j_in = min (max (jms, nint (j_real)), ime)
+          data_out(i, j) = data_in(i_in, j_in)
+        end do
+      end do
+
+    end subroutine Interp_horizontal_nearest
+
+    subroutine Interp_horizontal_bilinear (data_in, proj_data_in, ims, ime, jms, jme, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, &
+        lats_out, lons_out, data_out)
+
+    ! Purpose: bi-linear interpolation + nearest neighbor extrapolation
+
+      implicit none
+
+      type (proj_lc_t), intent (in) :: proj_data_in
+      integer, intent (in) :: ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, ims, ime, jms, jme
+      real, dimension(ims:ime, jms:jme), intent (in) :: data_in
+      real, dimension (ifms:ifme, jfms:jfme), intent (in) :: lats_out, lons_out
+      real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: data_out
+
+      integer :: i, j, i0, j0, i1, j1
+      real :: i_real, j_real, di, dj
+
+
+      do j = jfts, jfte
+        do i = ifts, ifte
+          call proj_data_in%Calc_ij (lats_out(i, j), lons_out(i, j), i_real, j_real)
+
+          i0 = max (ims, min (ime - 1, int (floor (i_real))))
+          j0 = max (jms, min (jme - 1, int (floor (j_real))))
+          i1 = i0 + 1
+          j1 = j0 + 1
+
+          di = max (0.0, min (1.0, i_real - real (i0)))
+          dj = max (0.0, min (1.0, j_real - real (j0)))
+
+          data_out(i, j) = (1.0 - di) * (1.0 - dj) * data_in(i0, j0) + &
+              di * (1.0 - dj) * data_in(i1, j0) + &
+              (1.0 - di) * dj * data_in(i0, j1) + &
+              di * dj * data_in(i1, j1)
+        end do
+      end do
+
+    end subroutine Interp_horizontal_bilinear
+
     subroutine Interp_profile (fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, kfds, kfde, &
         uin, vin, z_at_w, z0f, uout, vout)
 

state/state_mod.F90

@@ -241,7 +241,7 @@ subroutine Handle_wrfdata_update (this, wrf, config_flags)
         call wrf%Update_atm_state (this%datetime_now)
 
         if (DEBUG_LOCAL) call Print_message ('  Interpolating WRF vars...')
-        call this%interpolate_vars_atm_to_fire(wrf, config_flags)
+        call this%Interpolate_vars_atm_to_fire(wrf, config_flags)
 
         call this%datetime_next_atm_update%Add_seconds (config_flags%interval_atm)
 
@@ -744,52 +744,46 @@ subroutine Interpolate_vars_atm_to_fire (this, wrf, config_flags)
       type (wrf_t), intent(inout) :: wrf
       type (namelist_t), intent (in) :: config_flags
 
-      real, dimension(:, :), allocatable :: var2d
       integer :: i, j
 
 
-        ! We need the fire grid lat/lon
-      if (allocated (this%lats) .and. allocated (this%lons)) then
+      if (.not. allocated (this%lats) .or. .not. allocated (this%lons)) &
+          call Stop_simulation ('Init lats/lons before calling hinterp atm variables')
 
-        If_start: if (this%datetime_now == this%datetime_start) then
-          call wrf%interp_var2grid_nearest (this%lats(this%ifps:this%ifpe, this%jfps:this%jfpe), &
-              this%lons(this%ifps:this%ifpe, this%jfps:this%jfpe), 'fz0', var2d)
-              this%fz0(this%ifps:this%ifpe, this%jfps:this%jfpe) = var2d
-        endif If_start
+      if (this%datetime_now == this%datetime_start) call wrf%Interp_var2grid (this%lats, this%lons, &
+          this%ifms, this%ifme, this%jfms, this%jfme, config_flags%num_tiles, this%i_start, this%i_end, &
+          this%j_start, this%j_end, 'fz0', config_flags%hinterp_opt, this%fz0)
 
-        do j = 1, wrf%jde
-          do i = 1, wrf%ide
-            call this%Interpolate_profile (config_flags, config_flags%fire_wind_height, this%kfds, this%kfde, &
-                wrf%u3d_stag(i,:,j),wrf%v3d_stag(i,:,j), wrf%phl_stag(i,:,j), wrf%ua(i,j),wrf%va(i,j),wrf%z0_stag(i,j))
-          end do
+      do j = 1, wrf%jde
+        do i = 1, wrf%ide
+          call this%Interpolate_profile (config_flags, config_flags%fire_wind_height, this%kfds, this%kfde, &
+              wrf%u3d_stag(i,:,j),wrf%v3d_stag(i,:,j), wrf%phl_stag(i,:,j), wrf%ua(i,j),wrf%va(i,j),wrf%z0_stag(i,j))
         end do
+      end do
 
-        call wrf%interp_var2grid_nearest (this%lats(this%ifps:this%ifpe, this%jfps:this%jfpe), &
-            this%lons(this%ifps:this%ifpe, this%jfps:this%jfpe), 'uf', var2d)
-            this%uf(this%ifps:this%ifpe, this%jfps:this%jfpe) = var2d
-
-        call wrf%interp_var2grid_nearest (this%lats(this%ifps:this%ifpe, this%jfps:this%jfpe), &
-            this%lons(this%ifps:this%ifpe, this%jfps:this%jfpe), 'vf', var2d)
-            this%vf(this%ifps:this%ifpe, this%jfps:this%jfpe) = var2d
+      call wrf%Interp_var2grid (this%lats, this%lons, this%ifms, this%ifme, this%jfms, this%jfme, &
+          config_flags%num_tiles, this%i_start, this%i_end, this%j_start, this%j_end, &
+          'uf', config_flags%hinterp_opt, this%uf)
 
-        call wrf%interp_var2grid_nearest (this%lats(this%ifps:this%ifpe, this%jfps:this%jfpe), &
-            this%lons(this%ifps:this%ifpe, this%jfps:this%jfpe), 't2', var2d)
-            this%fire_t2(this%ifps:this%ifpe, this%jfps:this%jfpe) = var2d
+      call wrf%Interp_var2grid (this%lats, this%lons, this%ifms, this%ifme, this%jfms, this%jfme, &
+          config_flags%num_tiles, this%i_start, this%i_end, this%j_start, this%j_end, &
+          'vf', config_flags%hinterp_opt, this%vf)
 
-        call wrf%interp_var2grid_nearest (this%lats(this%ifps:this%ifpe, this%jfps:this%jfpe), &
-            this%lons(this%ifps:this%ifpe, this%jfps:this%jfpe), 'q2', var2d)
-            this%fire_q2(this%ifps:this%ifpe, this%jfps:this%jfpe) = var2d
+      call wrf%Interp_var2grid (this%lats, this%lons, this%ifms, this%ifme, this%jfms, this%jfme, &
+          config_flags%num_tiles, this%i_start, this%i_end, this%j_start, this%j_end, &
+          't2', config_flags%hinterp_opt, this%fire_t2)
 
-        call wrf%interp_var2grid_nearest (this%lats(this%ifps:this%ifpe, this%jfps:this%jfpe), &
-            this%lons(this%ifps:this%ifpe, this%jfps:this%jfpe), 'psfc', var2d)
-            this%fire_psfc(this%ifps:this%ifpe, this%jfps:this%jfpe) = var2d
+      call wrf%Interp_var2grid (this%lats, this%lons, this%ifms, this%ifme, this%jfms, this%jfme, &
+          config_flags%num_tiles, this%i_start, this%i_end, this%j_start, this%j_end, &
+          'q2', config_flags%hinterp_opt, this%fire_q2)
 
-        call wrf%interp_var2grid_nearest (this%lats(this%ifps:this%ifpe, this%jfps:this%jfpe), &
-            this%lons(this%ifps:this%ifpe, this%jfps:this%jfpe), 'rain', var2d)
-            this%fire_rain(this%ifps:this%ifpe, this%jfps:this%jfpe) = var2d
+      call wrf%Interp_var2grid (this%lats, this%lons, this%ifms, this%ifme, this%jfms, this%jfme, &
+          config_flags%num_tiles, this%i_start, this%i_end, this%j_start, this%j_end, &
+          'psfc', config_flags%hinterp_opt, this%fire_psfc)
 
-        deallocate (var2d)
-      end if
+      call wrf%Interp_var2grid (this%lats, this%lons, this%ifms, this%ifme, this%jfms, this%jfme, &
+          config_flags%num_tiles, this%i_start, this%i_end, this%j_start, this%j_end, &
+          'rain', config_flags%hinterp_opt, this%fire_rain)
 
     end subroutine Interpolate_vars_atm_to_fire