Merge pull request #21 from pedro-jm/vinterp

More sharable vertical interpolation code

Author: pedro-jm

driver/initialize_mod.F90

@@ -91,10 +91,6 @@ subroutine Init_fire_state (grid, config_flags, wrf)
         if (DEBUG_LOCAL) call Print_message ('    Initializing fire state')
         call grid%Initialization (config_flags, geogrid)
 
-        if (present (wrf)) then
-          if (DEBUG_LOCAL) call Print_message ('    Initializing atmospheric state')
-          call grid%Handle_wrfdata_update (wrf, config_flags)
-        end if
       else
           ! Ideal
         if (DEBUG_LOCAL) call Print_message ('    Initializing fire state')
@@ -103,6 +99,11 @@ subroutine Init_fire_state (grid, config_flags, wrf)
 
       call Init_fire_components (grid, config_flags)
 
+      if (present (wrf)) then
+        if (DEBUG_LOCAL) call Print_message ('    Initializing atmospheric state')
+        call grid%Handle_wrfdata_update (wrf, config_flags)
+      end if
+
       if (DEBUG_LOCAL) then
           ! print lat/lons
         open (newunit = unit_out, file = 'latlons_c.dat')

io/coupling_mod.F90

@@ -1,15 +1,14 @@
   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 interp_mod, only : Interp_horizontal_nearest, Interp_horizontal_bilinear, HINTERP_NEAREST, HINTERP_BILINEAR, Interp_profile
     use stderrout_mod, only : Stop_simulation
 
-
     implicit none
 
     private
 
-    public :: Interp_horizontal
+    public :: Interp_horizontal, Calc_fire_wind
 
   contains
 
@@ -62,4 +61,58 @@ subroutine Interp_horizontal (data_in, proj_data_in, ims, ime, jms, jme, ifms, i
 
     end subroutine Interp_horizontal
 
+    subroutine Calc_fire_wind (u3d, v3d, z_at_w, z0, iims, iime, jims, jime, kims, kime, fire_lsm_zcoupling, &
+        fire_lsm_zcoupling_ref, fire_wind_height, ioms, iome, joms, jome, iops, iope, jops, jope, &
+        u_out, v_out, cap_winds)
+
+      implicit none
+
+      integer, intent (in) :: iims, iime, jims, jime, kims, kime, ioms, iome, joms, jome, iops, iope, jops, jope
+      real, intent (in) :: fire_wind_height, fire_lsm_zcoupling_ref
+      logical, intent (in) :: fire_lsm_zcoupling
+      real, dimension(iims:iime, jims:jime, kims:kime), intent (in) :: u3d, v3d, z_at_w
+      real, dimension(iims:iime, jims:jime), intent (in) :: z0
+      real, dimension(ioms:iome, joms:joms), intent (out) :: u_out, v_out
+      logical, intent (in), optional :: cap_winds
+
+      real :: wspd
+      integer :: i, j
+      logical :: cap_winds_flag
+
+
+!      print *, 'shape u3d = ', shape (u3d)
+!      print *, 'shape v3d = ', shape (v3d)
+!      print *, 'shape z_at_w = ', shape (z_at_w)
+!      print *, 'shape u_out = ', shape (u_out)
+!      print *, 'shape v_out = ', shape (v_out)
+!      print *, 'shape z0 = ', shape (z0)
+
+      if (present (cap_winds)) then
+        cap_winds_flag = cap_winds
+      else
+        cap_winds_flag = .false.
+      end if
+
+      do j = jops, jope
+        do i = iops, iope
+          call Interp_profile (fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, kims, kime, &
+              u3d(i, j, :), v3d(i, j, :), z_at_w(i, j, :), z0(i, j), u_out(i, j), v_out(i, j))
+        end do
+      end do
+
+        ! To avoid arithmatic error
+      if (cap_winds_flag) then
+        do j = jops, jope
+          do i = iops, iope
+            wspd = sqrt (u_out(i, j) ** 2 + v_out(i, j) ** 2)
+            if (wspd < 0.001) then
+              u_out(i, j) = sign (0.001, u_out(i, j))
+              v_out(i, j) = sign (0.001, v_out(i, j))
+            end if
+          end do
+        end do
+      end if
+
+    end subroutine Calc_fire_wind
+
   end module coupling_mod

io/namelist_mod.F90

@@ -49,7 +49,7 @@ module namelist_mod
       integer :: fast_dist_reinit_freq = 600  ! Number of time steps to perform a reinit with fast distance reinit method
 
       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 :: wind_vinterp_opt = 0         ! "mid-flame height wind interpolation option: 0) Interp to specified height, 1) Use WAFs"
       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"
@@ -139,7 +139,7 @@ module namelist_mod
       real :: true_lat_2 = LAT_DEFAULT
 
         ! Atmosphere
-      integer :: kds = 1, kde = 1
+      integer :: kde = 1
     contains
       procedure, public :: Broadcast_nml => Broadcast_nml
       procedure, public :: Check_nml => Check_nml
@@ -411,7 +411,7 @@ subroutine Init_fire_block (this, file_name)
       integer :: fmoist_freq = 0              ! "frequency to run moisture model 0: use fmoist_dt, k>0: every k timesteps" "1"
       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 :: wind_vinterp_opt = 0         ! "mid-flame height wind interpolation option: 0) Interp to specified height, 1) Use WAFs"
       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"