From fc720f2f8db6c4ff1ae654212173b99d6259a75a Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Wed, 28 Jan 2026 11:45:50 -0700 Subject: [PATCH 1/9] Removing atm wind interpolation from fire state --- io/wrf_mod.F90 | 10 +++- nuopc/fire_behavior_nuopc.F90 | 15 +++--- state/state_mod.F90 | 98 +---------------------------------- 3 files changed, 19 insertions(+), 104 deletions(-) diff --git a/io/wrf_mod.F90 b/io/wrf_mod.F90 index 05217d5..7af7b3f 100644 --- a/io/wrf_mod.F90 +++ b/io/wrf_mod.F90 @@ -1067,12 +1067,13 @@ subroutine Fire_tendency( & end subroutine Fire_tendency - subroutine Update_atm_state (this, datetime_now) + subroutine Update_atm_state (this, datetime_now, config_flags) implicit none class (wrf_t), intent(in out) :: this type (datetime_t), intent (in) :: datetime_now + type (namelist_t), intent (in) :: config_flags integer :: i, j, k @@ -1135,6 +1136,13 @@ subroutine Update_atm_state (this, datetime_now) end do call this%Destroy_phl () + do j = 1, this%jde - 1 + do i = 1, this%ide - 1 + call Interp_profile (config_flags%fire_lsm_zcoupling, config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, this%kds, this%kde, & + this%u3d_stag(i, :, j), this%v3d_stag(i, :, j), this%phl_stag(i, :, j) / G, this%z0_stag(i, j), this%ua(i, j), this%va(i, j)) + end do + end do + end subroutine Update_atm_state end module wrf_mod diff --git a/nuopc/fire_behavior_nuopc.F90 b/nuopc/fire_behavior_nuopc.F90 index 1b3cfd9..86b3313 100644 --- a/nuopc/fire_behavior_nuopc.F90 +++ b/nuopc/fire_behavior_nuopc.F90 @@ -18,6 +18,7 @@ module fire_behavior_nuopc use advance_mod, only : Advance_state use constants_mod, only : G, XLV, CP, FVIRT, R_D use stderrout_mod, only : Stop_simulation + use interp_mod, only : Interp_profile implicit none @@ -844,12 +845,14 @@ subroutine Advance(model, rc) case (0) do j = grid%jfps, grid%jfpe do i = grid%ifps, grid%ifpe - call grid%Interpolate_profile (config_flags, & ! for debug output, <= 0 no output - config_flags%fire_wind_height, & ! interpolation height - grid%kfds, grid%kfde, & ! fire grid dimensions - atm_u3d(i,j,:),atm_v3d(i,j,:), & ! atm grid arrays in - atm_ph(i,j,:), & - grid%uf(i,j),grid%vf(i,j),grid%fz0(i,j)) + call Interp_profile (config_flags%fire_lsm_zcoupling, & + config_flags%fire_lsm_zcoupling_ref, & + config_flags%fire_wind_height, & + grid%kfds, grid%kfde, & + atm_u3d(i, j, :), atm_v3d(i, j, :), & + atm_ph(i, j, :) / 9.81, & + grid%fz0(i, j), & + grid%uf(i, j), grid%vf(i, j)) ! avoid arithmatic error wspd = (grid%uf(i,j) ** 2. + grid%vf(i,j) ** 2.) ** .5 diff --git a/state/state_mod.F90 b/state/state_mod.F90 index d73323e..ff3f040 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -114,7 +114,6 @@ module state_mod procedure :: Init_tiles => Init_tiles procedure :: Init_tiles_in_wrf => Init_tiles_in_wrf procedure :: Interpolate_vars_atm_to_fire => Interpolate_vars_atm_to_fire - procedure, public :: Interpolate_profile => Interpolate_profile procedure, public :: Print => Print_domain ! private procedure, public :: Print_tiles => Print_tiles procedure, public :: Save_state => Save_state @@ -238,7 +237,7 @@ subroutine Handle_wrfdata_update (this, wrf, config_flags) if (DEBUG_LOCAL) call Print_message (' Updating WRF atm state...') if (DEBUG_LOCAL) call this%datetime_now%Print_datetime () - call wrf%Update_atm_state (this%datetime_now) + call wrf%Update_atm_state (this%datetime_now, config_flags) if (DEBUG_LOCAL) call Print_message (' Interpolating WRF vars...') call this%Interpolate_vars_atm_to_fire(wrf, config_flags) @@ -754,13 +753,6 @@ subroutine Interpolate_vars_atm_to_fire (this, wrf, config_flags) 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 - end do - 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) @@ -787,94 +779,6 @@ subroutine Interpolate_vars_atm_to_fire (this, wrf, config_flags) end subroutine Interpolate_vars_atm_to_fire - subroutine Interpolate_profile (this, config_flags, fire_wind_height, kfds, kfde, & - uin, vin, phl, uout, vout,z0f) - - implicit none - - class (state_fire_t), intent (in) :: this - type (namelist_t), intent (in) :: config_flags - real, intent (in) :: fire_wind_height - integer, intent (in) :: kfds, kfde - real, intent (in) :: uin(:), vin(:), phl(:) - real, intent (out) :: uout, vout - real, intent (in) :: z0f - - - real, parameter :: VK_KAPPA = 0.4 - real, dimension (kfds:kfde - 1) :: altw, hgt - integer :: k, kdmax - real :: loght, loglast, logz0, logfwh, ht, r_nan, fire_wind_height_local, z0fc, & - ust_d, wsf, wsf1, uf_temp, vf_temp - - - ! max layer to interpolate from, can be less - kdmax = kfde - 2 - do k = kfds, kdmax + 1 - ! altitude of the bottom w-point - altw(k) = phl(k) / G - end do - - do k = kfds, kdmax - ! height of the mass point above the ground - hgt(k) = 0.5 * (altw(k) + altw(k + 1)) - altw(kfds) - end do - - ! extrapolate mid-flame height from fire_lsm_zcoupling_ref? - if (config_flags%fire_lsm_zcoupling) then - logfwh = log (config_flags%fire_lsm_zcoupling_ref) - fire_wind_height_local = config_flags%fire_lsm_zcoupling_ref - else - logfwh = log (fire_wind_height) - fire_wind_height_local = fire_wind_height - end if - - ! interpolate u - if (fire_wind_height_local > z0f)then - do k = kfds, kdmax - ht = hgt(k) - if (ht >= fire_wind_height_local) then - ! found layer k this point is in - loght = log(ht) - if (k == kfds) then - ! first layer, log linear interpolation from 0 at zr - logz0 = log(z0f) - uout = uin(k) * (logfwh - logz0) / (loght - logz0) - vout = vin(k) * (logfwh - logz0) / (loght - logz0) - else - ! log linear interpolation - loglast = log (hgt(k - 1)) - uout = uin(k - 1) + (uin(k) - uin(k - 1)) * (logfwh - loglast) / (loght - loglast) - vout = vin(k - 1) + (vin(k) - vin(k - 1)) * (logfwh - loglast) / (loght - loglast) - end if - exit - end if - if (k == kdmax) then - ! last layer, still not high enough - uout = uin(k) - vout = vin(k) - end if - end do - else - ! roughness higher than the fire wind height - uout = 0.0 - vout = 0.0 - end if - - ! Extrapol wind to target height - if (config_flags%fire_lsm_zcoupling) then - uf_temp = uout - vf_temp = vout - wsf = max (sqrt (uf_temp ** 2.0 + vf_temp ** 2.0), 0.1) - z0fc = z0f - ust_d = wsf * VK_KAPPA / log(config_flags%fire_lsm_zcoupling_ref / z0fc) - wsf1 = (ust_d / VK_KAPPA) * log((fire_wind_height + z0fc) / z0fc) - uout = wsf1 * uf_temp / wsf - vout = wsf1 * vf_temp / wsf - end if - - end subroutine Interpolate_profile - subroutine Print_domain (this) use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT From e6b5525df16fc6ff69a3af18ea14e60c5aba2d9d Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Wed, 28 Jan 2026 17:14:45 -0700 Subject: [PATCH 2/9] Adding a generic sub to calculate the fire winds --- io/coupling_mod.F90 | 36 +++++++++++++++++++++++++++++++++--- io/wrf_mod.F90 | 41 +++++++---------------------------------- 2 files changed, 40 insertions(+), 37 deletions(-) diff --git a/io/coupling_mod.F90 b/io/coupling_mod.F90 index bd912f4..f685470 100644 --- a/io/coupling_mod.F90 +++ b/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,35 @@ 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, fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, ua, va) + + implicit none + + real, intent (in) :: fire_wind_height, fire_lsm_zcoupling_ref + logical, intent (in) :: fire_lsm_zcoupling + real, dimension(:, :, :), intent (in) :: u3d, v3d, z_at_w + real, dimension(:, :), intent (in) :: z0 + real, dimension(:, :), intent (out) :: ua, va + + integer :: i, j, kds, kde + + +! print *, 'shape u3d = ', shape (u3d) +! print *, 'shape v3d = ', shape (v3d) +! print *, 'shape z_at_w = ', shape (z_at_w) +! print *, 'shape ua = ', shape (ua) +! print *, 'shape va = ', shape (va) +! print *, 'shape z0 = ', shape (z0) + + kds = 1 + kde = size (z_at_w, dim = 3) + do j = 1, size (ua, dim = 2) + do i = 1, size (ua, dim = 1) + call Interp_profile (fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, kds, kde, & + u3d(i, j, :), v3d(i, j, :), z_at_w(i, j, :), z0(i, j), ua(i, j), va(i, j)) + end do + end do + + end subroutine Calc_fire_wind + end module coupling_mod diff --git a/io/wrf_mod.F90 b/io/wrf_mod.F90 index 7af7b3f..f1becbe 100644 --- a/io/wrf_mod.F90 +++ b/io/wrf_mod.F90 @@ -7,7 +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 + use coupling_mod, only : Interp_horizontal, Calc_fire_wind implicit none @@ -1075,8 +1075,6 @@ subroutine Update_atm_state (this, datetime_now, config_flags) type (datetime_t), intent (in) :: datetime_now type (namelist_t), intent (in) :: config_flags - integer :: i, j, k - ! Update t2_stag call this%Get_t2 (datetime_now) @@ -1105,45 +1103,20 @@ subroutine Update_atm_state (this, datetime_now, config_flags) ! Update U3D call this%Get_u3d_stag (datetime_now) - do k = this%kds, this%kde - 1 - do j = this%jds, this%jde - 1 - do i = this%ids, this%ide - this%u3d_stag(i, k, j) = this%u3d(i, j, k) - end do - end do - end do - call this%Destroy_u3d () ! Update V3D call this%Get_v3d_stag (datetime_now) - do k = this%kds, this%kde - 1 - do j = this%jds, this%jde - do i = this%ids, this%ide - 1 - this%v3d_stag(i, k, j) = this%v3d(i, j, k) - end do - end do - end do - call this%Destroy_v3d () ! Update geopotential heights call this%Get_phl (datetime_now) - do k = this%kds, this%kde - do j = this%jds, this%jde - 1 - do i = this%ids, this%ide - 1 - this%phl_stag(i, k, j) = this%phl(i, j, k) - end do - end do - end do - call this%Destroy_phl () - do j = 1, this%jde - 1 - do i = 1, this%ide - 1 - call Interp_profile (config_flags%fire_lsm_zcoupling, config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, this%kds, this%kde, & - this%u3d_stag(i, :, j), this%v3d_stag(i, :, j), this%phl_stag(i, :, j) / G, this%z0_stag(i, j), this%ua(i, j), this%va(i, j)) - end do - end do +! call this%Destroy_u3d () +! call this%Destroy_v3d () +! call this%Destroy_phl () + + call Calc_fire_wind (this%u3d(:this%ide-1, :, :), this%v3d(:, :this%jde-1, :), this%phl / G, this%z0_stag(:this%ide-1,:this%jde-1), config_flags%fire_lsm_zcoupling, & + config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, this%ua(:this%ide-1, :this%jde-1), this%va(:this%ide-1, :this%jde-1)) end subroutine Update_atm_state end module wrf_mod - From 7132edc3aee5ce650a6644072bdb296e27e71246 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Wed, 28 Jan 2026 17:46:33 -0700 Subject: [PATCH 3/9] Use winds at mass points in offline runs --- io/wrf_mod.F90 | 6 ++++-- tests/test7.s | 16 ++++++++++++---- tests/test8.s | 10 ++++++---- 3 files changed, 22 insertions(+), 10 deletions(-) diff --git a/io/wrf_mod.F90 b/io/wrf_mod.F90 index f1becbe..c167d85 100644 --- a/io/wrf_mod.F90 +++ b/io/wrf_mod.F90 @@ -1102,10 +1102,10 @@ subroutine Update_atm_state (this, datetime_now, config_flags) call this%Destroy_z0 () ! Update U3D - call this%Get_u3d_stag (datetime_now) + call this%Get_u3d (datetime_now) ! Update V3D - call this%Get_v3d_stag (datetime_now) + call this%Get_v3d (datetime_now) ! Update geopotential heights call this%Get_phl (datetime_now) @@ -1116,6 +1116,8 @@ subroutine Update_atm_state (this, datetime_now, config_flags) call Calc_fire_wind (this%u3d(:this%ide-1, :, :), this%v3d(:, :this%jde-1, :), this%phl / G, this%z0_stag(:this%ide-1,:this%jde-1), config_flags%fire_lsm_zcoupling, & config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, this%ua(:this%ide-1, :this%jde-1), this%va(:this%ide-1, :this%jde-1)) +! call Calc_fire_wind (this%u3d, this%v3d, this%phl / G, this%z0_stag(:this%ide-1,:this%jde-1), config_flags%fire_lsm_zcoupling, & +! config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, this%ua(:this%ide-1, :this%jde-1), this%va(:this%ide-1, :this%jde-1)) end subroutine Update_atm_state diff --git a/tests/test7.s b/tests/test7.s index 91f078b..8c0d41d 100755 --- a/tests/test7.s +++ b/tests/test7.s @@ -8,6 +8,9 @@ # purge_output=1 # 0) No, 1) yes plot=0 # 0) No, 1) yes +update_nml=1 # 0) No, 1) yes +#export OMP_NUM_THREADS=16 +#export OMP_PROC_BIND=true # ################################################# # @@ -25,7 +28,10 @@ file_exe=../install/bin/fire_behavior.exe file_output=test7_output.txt cp ./test7/wrf.nc . -cp ./test7/namelist.fire . +if [ $update_nml -eq 1 ] +then + cp ./test7/namelist.fire . +fi cp ./test7/geo_em.d01.nc . rm -f ./$file_output @@ -80,7 +86,7 @@ then test=$(diff ./file1.dat ./file2.dat | wc -l) echo $test - if [ $test -eq 10 ] + if [ $test -eq 8 ] then echo ' Ignore this difference:' diff ./file1.dat ./file2.dat @@ -169,9 +175,12 @@ fi # # Purge -rm -f ./namelist.fire.output ./file1.dat ./file2.dat ./wrf_input.dat ./geo_em.d01.nc ./namelist.fire +rm -f ./namelist.fire.output ./file1.dat ./file2.dat ./wrf_input.dat if [ $purge_output -eq 1 ] then + rm -f ./wrf.nc + rm -f ./geo_em.d01.nc + rm -f ./namelist.fire rm -rf ./$file_output rm -f ./fire_output_2012-06-25_18:00:??.nc fi @@ -203,7 +212,6 @@ then fi rm -f ./latlons.dat ./latlons_c.dat ./wrf_latlons_atm.dat ./wrf_latlons_fire.dat -rm -f ./wrf.nc # OUTPUT_LATLON_CHECK is true # Check we are doing the projections right diff --git a/tests/test8.s b/tests/test8.s index fd8439c..b9a304e 100755 --- a/tests/test8.s +++ b/tests/test8.s @@ -81,7 +81,7 @@ then test=$(diff ./file1.dat ./file2.dat | wc -l) # Here we allow one difference since we are not expecting bit4bit results echo $test - if [ $test -eq 12 ] + if [ $test -eq 4 ] then echo ' Ignore this difference:' diff ./file1.dat ./file2.dat @@ -105,7 +105,8 @@ then grep "$var" $file_wrf | awk '{print $2, $7}' > ./file2.dat test=$(diff ./file1.dat ./file2.dat | wc -l) - if [ $test -eq 4 ] + echo $test + if [ $test -eq 8 ] then echo ' Ignore this difference:' diff ./file1.dat ./file2.dat @@ -129,8 +130,9 @@ then grep "$var" $file_wrf | awk '{print $2, $7}' > ./file2.dat test=$(diff ./file1.dat ./file2.dat | wc -l) + echo $test # Here we allow one difference since we are not expecting bit4bit results - if [ $test -eq 8 ] + if [ $test -eq 4 ] then echo ' Ignore this difference:' diff ./file1.dat ./file2.dat @@ -156,7 +158,7 @@ then test=$(diff ./file1.dat ./file2.dat | wc -l) # Here we allow one difference since we are not expecting bit4bit results echo $test - if [ $test -eq 12 ] + if [ $test -eq 8 ] then echo ' Ignore this difference:' diff ./file1.dat ./file2.dat From b91acbf00611992ca1ec9b5fea0af53c264964bf Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Thu, 29 Jan 2026 08:08:23 -0700 Subject: [PATCH 4/9] More generic interpolation of wind profiles to calc mid-flame height winds --- io/coupling_mod.F90 | 48 ++++++++--- io/namelist_mod.F90 | 2 +- io/wrf_mod.F90 | 149 +++++++++++++++------------------- nuopc/fire_behavior_nuopc.F90 | 48 +++++------ nuopc/wrf_nuopc.F90 | 15 ++-- state/state_mod.F90 | 8 +- 6 files changed, 142 insertions(+), 128 deletions(-) diff --git a/io/coupling_mod.F90 b/io/coupling_mod.F90 index f685470..d1f641b 100644 --- a/io/coupling_mod.F90 +++ b/io/coupling_mod.F90 @@ -61,35 +61,57 @@ 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, fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, ua, va) + 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(:, :, :), intent (in) :: u3d, v3d, z_at_w - real, dimension(:, :), intent (in) :: z0 - real, dimension(:, :), intent (out) :: ua, va + 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 - integer :: i, j, kds, kde + 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 ua = ', shape (ua) -! print *, 'shape va = ', shape (va) +! print *, 'shape u_out = ', shape (u_out) +! print *, 'shape v_out = ', shape (v_out) ! print *, 'shape z0 = ', shape (z0) - kds = 1 - kde = size (z_at_w, dim = 3) - do j = 1, size (ua, dim = 2) - do i = 1, size (ua, dim = 1) - call Interp_profile (fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, kds, kde, & - u3d(i, j, :), v3d(i, j, :), z_at_w(i, j, :), z0(i, j), ua(i, j), va(i, j)) + 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 diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index cfd628b..a158a69 100644 --- a/io/namelist_mod.F90 +++ b/io/namelist_mod.F90 @@ -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 diff --git a/io/wrf_mod.F90 b/io/wrf_mod.F90 index c167d85..9622a4e 100644 --- a/io/wrf_mod.F90 +++ b/io/wrf_mod.F90 @@ -20,9 +20,8 @@ module wrf_mod type :: wrf_t character (len = 300) :: file_name - real, dimension(:, :, :), allocatable :: u3d, v3d, phl, u3d_stag, v3d_stag, phl_stag - real, dimension(:, :), allocatable :: lats, lons, lats_c, lons_c, t2, q2, z0, psfc, rain, ua, va, & - t2_stag, q2_stag, z0_stag, psfc_stag, rain_stag + real, dimension(:, :, :), allocatable :: u3d, v3d, phl + real, dimension(:, :), allocatable :: lats, lons, lats_c, lons_c, t2, q2, z0, psfc, rain, ua, va integer :: ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte real :: cen_lat, cen_lon, dx, dy, truelat1, truelat2, stand_lon contains @@ -526,16 +525,18 @@ function Wrf_t_const (file_name, config_flags) result (return_value) ! latlon at corners call return_value%Get_latcloncs () + ! Init domain dimensions return_value%ids = 1 - return_value%jds = 1 call Get_netcdf_dim (trim (file_name), 'west_east_stag', att_int32) return_value%ide = att_int32 + return_value%jds = 1 call Get_netcdf_dim (trim (file_name), 'south_north_stag', att_int32) return_value%jde = att_int32 + return_value%kds = 1 + call Get_netcdf_dim (trim (file_name), 'bottom_top_stag', att_int32) + return_value%kde = att_int32 - return_value%kds = config_flags%kds - return_value%kde = config_flags%kde - + ! Init rest of dimensions return_value%ims = return_value%ids return_value%ime = return_value%ide return_value%kms = return_value%kds @@ -552,37 +553,26 @@ function Wrf_t_const (file_name, config_flags) result (return_value) if (DEBUG_LOCAL) call return_value%Print_domain() - allocate (return_value%phl_stag(return_value%ims:return_value%ime, & - return_value%kms:return_value%kme, return_value%jms:return_value%jme)) - return_value%phl_stag = 0.0 - - allocate (return_value%u3d_stag(return_value%ims:return_value%ime, & - return_value%kms:return_value%kme, return_value%jms:return_value%jme)) - return_value%u3d_stag = 0.0 - - allocate (return_value%v3d_stag(return_value%ims:return_value%ime, & - return_value%kms:return_value%kme, return_value%jms:return_value%jme)) - return_value%v3d_stag = 0.0 + ! Init some vars to default values + allocate (return_value%z0(return_value%ids:return_value%ide - 1, return_value%jds:return_value%jde - 1)) + return_value%z0 = DEFAULT_Z0 - allocate (return_value%z0_stag(return_value%ims:return_value%ime, return_value%jms:return_value%jme)) - return_value%z0_stag = DEFAULT_Z0 + allocate (return_value%rain(return_value%ids:return_value%ide - 1, return_value%jds:return_value%jde - 1)) + return_value%rain = DEFAULT_RAIN - allocate (return_value%rain_stag(return_value%ims:return_value%ime, return_value%jms:return_value%jme)) - return_value%rain_stag = DEFAULT_RAIN + allocate (return_value%t2(return_value%ids:return_value%ide - 1, return_value%jds:return_value%jde - 1)) + return_value%t2 = DEFAULT_T2 - allocate (return_value%t2_stag(return_value%ims:return_value%ime, return_value%jms:return_value%jme)) - return_value%t2_stag = DEFAULT_T2 + allocate (return_value%q2(return_value%ids:return_value%ide - 1, return_value%jds:return_value%jde - 1)) + return_value%q2 = DEFAULT_Q2 - allocate (return_value%q2_stag(return_value%ims:return_value%ime, return_value%jms:return_value%jme)) - return_value%q2_stag = DEFAULT_Q2 + allocate (return_value%psfc(return_value%ids:return_value%ide - 1, return_value%jds:return_value%jde - 1)) + return_value%psfc = DEFAULT_PSFC - allocate (return_value%psfc_stag(return_value%ims:return_value%ime, return_value%jms:return_value%jme)) - return_value%psfc_stag = DEFAULT_PSFC - - allocate (return_value%ua(return_value%ims:return_value%ime, return_value%jms:return_value%jme)) + allocate (return_value%ua(return_value%ids:return_value%ide - 1, return_value%jds:return_value%jde - 1)) return_value%ua = 0.0 - allocate (return_value%va(return_value%ims:return_value%ime, return_value%jms:return_value%jme)) + allocate (return_value%va(return_value%ids:return_value%ide - 1, return_value%jds:return_value%jde - 1)) return_value%va = 0.0 if (DEBUG_LOCAL) Call Print_message ('Leaving wrf_t constructor') @@ -669,25 +659,25 @@ subroutine Interp_var2grid (this, lats_out, lons_out, ifms, ifme, jfms, jfme, & ! Get WRF data select case (var_name) case ('t2') - var_wrf = this%t2_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) + var_wrf = this%t2 case ('q2') - var_wrf = this%q2_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) + var_wrf = this%q2 case ('psfc') - var_wrf = this%psfc_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) + var_wrf = this%psfc case ('rain') - var_wrf = this%rain_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) + var_wrf = this%rain case ('fz0') - var_wrf = this%z0_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) + var_wrf = this%z0 case ('uf') - var_wrf = this%ua(this%ids:this%ide - 1, this%jds:this%jde - 1) + var_wrf = this%ua case ('vf') - var_wrf = this%va(this%ids:this%ide - 1, this%jds:this%jde - 1) + var_wrf = this%va case default call Stop_simulation ('Unknown variable name to interpolate') @@ -699,7 +689,6 @@ subroutine Interp_var2grid (this, lats_out, lons_out, ifms, ifme, jfms, jfme, & 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) @@ -1075,49 +1064,45 @@ subroutine Update_atm_state (this, datetime_now, config_flags) type (datetime_t), intent (in) :: datetime_now type (namelist_t), intent (in) :: config_flags - - ! Update t2_stag - call this%Get_t2 (datetime_now) - this%t2_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) = this%t2(:, :) - call this%Destroy_t2 () - - ! Update q2 - call this%Get_q2 (datetime_now) - this%q2_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) = this%q2(:, :) - call this%Destroy_q2 () - - ! Update psfc - call this%Get_psfc (datetime_now) - this%psfc_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) = this%psfc(:, :) - call this%Destroy_psfc () - - ! Update rain - call this%Get_rain (datetime_now) - this%rain_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) = this%rain(:, :) - call this%Destroy_rain () - - ! Update z0 - call this%Get_z0 (datetime_now) - this%z0_stag(this%ids:this%ide - 1, this%jds:this%jde - 1) = this%z0(:, :) - call this%Destroy_z0 () - - ! Update U3D - call this%Get_u3d (datetime_now) - - ! Update V3D - call this%Get_v3d (datetime_now) - - ! Update geopotential heights - call this%Get_phl (datetime_now) - -! call this%Destroy_u3d () -! call this%Destroy_v3d () -! call this%Destroy_phl () - - call Calc_fire_wind (this%u3d(:this%ide-1, :, :), this%v3d(:, :this%jde-1, :), this%phl / G, this%z0_stag(:this%ide-1,:this%jde-1), config_flags%fire_lsm_zcoupling, & - config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, this%ua(:this%ide-1, :this%jde-1), this%va(:this%ide-1, :this%jde-1)) -! call Calc_fire_wind (this%u3d, this%v3d, this%phl / G, this%z0_stag(:this%ide-1,:this%jde-1), config_flags%fire_lsm_zcoupling, & -! config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, this%ua(:this%ide-1, :this%jde-1), this%va(:this%ide-1, :this%jde-1)) + integer :: iims, iime, jims, jime, kims, kime, ioms, iome, joms, jome, iops, iope, jops, jope + + + call this%Get_t2 (datetime_now) + call this%Get_q2 (datetime_now) + call this%Get_psfc (datetime_now) + call this%Get_rain (datetime_now) + call this%Get_z0 (datetime_now) + call this%Get_u3d (datetime_now) + call this%Get_v3d (datetime_now) + call this%Get_phl (datetime_now) + + ! Set input (i) and output (o) indices + iims = this%ids + iime = this%ide - 1 + jims = this%ids + jime = this%ide - 1 + kims = this%kds + kime = this%kde - 1 + + ioms = this%ids + iome = this%ide - 1 + joms = this%ids + jome = this%ide - 1 + + iops = this%ids + iope = this%ide - 1 + jops = this%ids + jope = this%ide - 1 + ! For compatibility with nuopc couplings + ! pass z_at_w with vertical dim kde - 1 instead of kde + call Calc_fire_wind (this%u3d, this%v3d, this%phl(iims:iime, jims:jime, kims:kime) / G, this%z0, & + iims, iime, jims, jime, kims, kime, config_flags%fire_lsm_zcoupling, config_flags%fire_lsm_zcoupling_ref, & + config_flags%fire_wind_height, ioms, iome, joms, jome, iops, iope, jops, jope, this%ua, this%va) + +! call this%Destroy_z0 () + call this%Destroy_u3d () + call this%Destroy_v3d () + call this%Destroy_phl () end subroutine Update_atm_state diff --git a/nuopc/fire_behavior_nuopc.F90 b/nuopc/fire_behavior_nuopc.F90 index 86b3313..eb5d1d5 100644 --- a/nuopc/fire_behavior_nuopc.F90 +++ b/nuopc/fire_behavior_nuopc.F90 @@ -18,7 +18,7 @@ module fire_behavior_nuopc use advance_mod, only : Advance_state use constants_mod, only : G, XLV, CP, FVIRT, R_D use stderrout_mod, only : Stop_simulation - use interp_mod, only : Interp_profile + use coupling_mod, only : Calc_fire_wind implicit none @@ -771,13 +771,13 @@ subroutine Advance(model, rc) real(ESMF_KIND_R8) :: ts type(ESMF_State) :: importState, exportState integer :: i, j - real :: wspd, q0, rho + real :: q0, rho character(len=160) :: msgString real, dimension(:, :, :), allocatable :: atm_u3d, atm_v3d, atm_ph real, dimension(:, :), allocatable :: atm_lowest_t, atm_lowest_q, atm_lowest_pres real, dimension(:, :), allocatable :: grnhfx_kinematic, grnqfx_kinematic, smoke real :: dtratio - + integer :: iims, iime, jims, jime, kims, kime, ioms, iome, joms, jome, iops, iope, jops, jope rc = ESMF_SUCCESS @@ -843,26 +843,28 @@ subroutine Advance(model, rc) select case (config_flags%wind_vinterp_opt) case (0) - do j = grid%jfps, grid%jfpe - do i = grid%ifps, grid%ifpe - call Interp_profile (config_flags%fire_lsm_zcoupling, & - config_flags%fire_lsm_zcoupling_ref, & - config_flags%fire_wind_height, & - grid%kfds, grid%kfde, & - atm_u3d(i, j, :), atm_v3d(i, j, :), & - atm_ph(i, j, :) / 9.81, & - grid%fz0(i, j), & - grid%uf(i, j), grid%vf(i, j)) - - ! avoid arithmatic error - wspd = (grid%uf(i,j) ** 2. + grid%vf(i,j) ** 2.) ** .5 - if (wspd < 0.001) then - grid%uf(i,j) = sign(0.001, grid%uf(i,j)) - grid%vf(i,j) = sign(0.001, grid%vf(i,j)) - endif - - enddo - enddo + + iims = grid%ifps + iime = grid%ifpe + jims = grid%jfps + jime = grid%jfpe + kims = grid%kfds + kime = grid%kfde - 1 + + ioms = grid%ifms + iome = grid%ifme + joms = grid%jfms + jome = grid%jfme + + iops = grid%ifps + iope = grid%ifpe + jops = grid%jfps + jope = grid%jfpe + + call Calc_fire_wind (atm_u3d, atm_v3d, atm_ph / 9.81, grid%fz0, iims, iime, jims, jime, kims, kime, config_flags%fire_lsm_zcoupling, & + config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, ioms, iome, joms, jome, iops, iope, jops, jope, & + grid%uf, grid%vf, cap_winds = .true.) + case (1) do j = grid%jfps, grid%jfpe do i = grid%ifps, grid%ifpe diff --git a/nuopc/wrf_nuopc.F90 b/nuopc/wrf_nuopc.F90 index 651315d..afb6d89 100644 --- a/nuopc/wrf_nuopc.F90 +++ b/nuopc/wrf_nuopc.F90 @@ -109,11 +109,16 @@ subroutine Advertise(model, rc) call Init_atm_state(state, config_flags) - allocate (state%q2(size(state%lats, dim=1), size(state%lats, dim=2))) - allocate (state%t2(size(state%lats, dim=1), size(state%lats, dim=2))) - allocate (state%z0(size(state%lats, dim=1), size(state%lats, dim=2))) - allocate (state%psfc(size(state%lats, dim=1), size(state%lats, dim=2))) - allocate (state%rain(size(state%lats, dim=1), size(state%lats, dim=2))) + if (.not. allocated (state%q2)) & + allocate (state%q2(size(state%lats, dim=1), size(state%lats, dim=2))) + if (.not. allocated (state%t2)) & + allocate (state%t2(size(state%lats, dim=1), size(state%lats, dim=2))) + if (.not. allocated (state%z0)) & + allocate (state%z0(size(state%lats, dim=1), size(state%lats, dim=2))) + if (.not. allocated (state%psfc)) & + allocate (state%psfc(size(state%lats, dim=1), size(state%lats, dim=2))) + if (.not. allocated (state%rain)) & + allocate (state%rain(size(state%lats, dim=1), size(state%lats, dim=2))) allocate (state%u3d(size(state%lats, dim=1), size(state%lats, dim=2), state%kde - 1)) allocate (state%v3d(size(state%lats, dim=1), size(state%lats, dim=2), state%kde - 1)) allocate (state%phl(size(state%lats, dim=1), size(state%lats, dim=2), state%kde - 1)) diff --git a/state/state_mod.F90 b/state/state_mod.F90 index ff3f040..0d4f238 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -408,13 +408,13 @@ subroutine Init_domain (this, config_flags, geogrid, & this%jfps = jps this%jfpe = jpe - this%kfds = config_flags%kds + this%kfds = 1 this%kfde = config_flags%kde - this%kfms = config_flags%kds + this%kfms = 1 this%kfme = config_flags%kde - this%kfps = config_flags%kds + this%kfps = 1 this%kfpe = config_flags%kde - this%kfts = config_flags%kds + this%kfts = 1 this%kfte = config_flags%kde call this%Init_tiles (config_flags) From ea13aef02662ddf5f3b37033b29b07b2be8bf7d8 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Thu, 29 Jan 2026 15:31:47 -0700 Subject: [PATCH 5/9] Generic wind interpolation + WAFS --- driver/initialize_mod.F90 | 9 ++- io/coupling_mod.F90 | 5 +- io/namelist_mod.F90 | 4 +- io/wrf_mod.F90 | 138 ++++++++++++++++++++++++++-------- nuopc/fire_behavior_nuopc.F90 | 7 +- state/state_mod.F90 | 35 ++++++++- 6 files changed, 154 insertions(+), 44 deletions(-) diff --git a/driver/initialize_mod.F90 b/driver/initialize_mod.F90 index c47ecce..1855b0e 100644 --- a/driver/initialize_mod.F90 +++ b/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') diff --git a/io/coupling_mod.F90 b/io/coupling_mod.F90 index d1f641b..e396fd9 100644 --- a/io/coupling_mod.F90 +++ b/io/coupling_mod.F90 @@ -61,8 +61,9 @@ 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) + 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 diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index a158a69..5ff35b4 100644 --- a/io/namelist_mod.F90 +++ b/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" @@ -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" diff --git a/io/wrf_mod.F90 b/io/wrf_mod.F90 index 9622a4e..222118d 100644 --- a/io/wrf_mod.F90 +++ b/io/wrf_mod.F90 @@ -21,7 +21,7 @@ module wrf_mod type :: wrf_t character (len = 300) :: file_name real, dimension(:, :, :), allocatable :: u3d, v3d, phl - real, dimension(:, :), allocatable :: lats, lons, lats_c, lons_c, t2, q2, z0, psfc, rain, ua, va + real, dimension(:, :), allocatable :: lats, lons, lats_c, lons_c, t2, q2, z0, psfc, rain, ua, va, u10, v10 integer :: ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte real :: cen_lat, cen_lon, dx, dy, truelat1, truelat2, stand_lon contains @@ -30,6 +30,8 @@ module wrf_mod procedure, public :: Destroy_rain => Destroy_rain procedure, public :: Destroy_q2 => Destroy_specific_humidity_2m procedure, public :: Destroy_t2 => Destroy_temperature_2m + procedure, public :: Destroy_u10 => Destroy_zonal_wind10 + procedure, public :: Destroy_v10 => Destroy_meridional_wind10 procedure, public :: Destroy_u3d => Destroy_zonal_wind procedure, public :: Destroy_v3d => Destroy_meridional_wind procedure, public :: Destroy_z0 => Destroy_z0 @@ -42,8 +44,10 @@ module wrf_mod procedure, public :: Get_psfc => Get_surface_pressure procedure, public :: Get_q2 => Get_specific_humidity_2m procedure, public :: Get_t2 => Get_temperature_2m + procedure, public :: Get_u10 => Get_zonal_wind10 procedure, public :: Get_u3d => Get_zonal_wind_3d procedure, public :: Get_u3d_stag => Get_zonal_wind_stag_3d + procedure, public :: Get_v10 => Get_meridional_wind10 procedure, public :: Get_v3d => Get_meridional_wind_3d procedure, public :: Get_v3d_stag => Get_meridional_wind_stag_3d procedure, public :: Get_z0 => Get_z0 @@ -154,6 +158,16 @@ subroutine Destroy_meridional_wind (this) end subroutine Destroy_meridional_wind + subroutine Destroy_meridional_wind10 (this) + + implicit none + + class (wrf_t), intent (in out) :: this + + if (allocated(this%v10)) deallocate (this%v10) + + end subroutine Destroy_meridional_wind10 + subroutine Destroy_rain (this) implicit none @@ -214,6 +228,16 @@ subroutine Destroy_zonal_wind (this) end subroutine Destroy_zonal_wind + subroutine Destroy_zonal_wind10 (this) + + implicit none + + class (wrf_t), intent (in out) :: this + + if (allocated(this%u10)) deallocate (this%u10) + + end subroutine Destroy_zonal_wind10 + function Get_datetime_index (this, datetime) result (return_value) use, intrinsic :: iso_fortran_env, only : ERROR_UNIT @@ -282,6 +306,23 @@ subroutine Get_geopotential_levels (this, datetime) end subroutine Get_geopotential_levels + subroutine Get_meridional_wind10 (this, datetime) + + implicit none + + class (wrf_t), intent (in out) :: this + type (datetime_t), intent (in) :: datetime + + real, dimension(:, :, :), allocatable :: var3d + integer :: nt + + + nt = this%Get_datetime_index (datetime) + call Get_netcdf_var (trim (this%file_name), 'V10', var3d) + this%v10 = var3d(:, :, nt) + + end subroutine Get_meridional_wind10 + subroutine Get_meridional_wind_3d (this, datetime) implicit none @@ -437,6 +478,23 @@ subroutine Get_z0 (this, datetime) end subroutine Get_z0 + subroutine Get_zonal_wind10 (this, datetime) + + implicit none + + class (wrf_t), intent (in out) :: this + type (datetime_t), intent (in) :: datetime + + real, dimension(:, :, :), allocatable :: var3d + integer :: nt + + + nt = this%Get_datetime_index (datetime) + call Get_netcdf_var (trim (this%file_name), 'U10', var3d) + this%u10 = var3d(:, :, nt) + + end subroutine Get_zonal_wind10 + subroutine Get_zonal_wind_3d (this, datetime) implicit none @@ -1072,37 +1130,53 @@ subroutine Update_atm_state (this, datetime_now, config_flags) call this%Get_psfc (datetime_now) call this%Get_rain (datetime_now) call this%Get_z0 (datetime_now) - call this%Get_u3d (datetime_now) - call this%Get_v3d (datetime_now) - call this%Get_phl (datetime_now) - - ! Set input (i) and output (o) indices - iims = this%ids - iime = this%ide - 1 - jims = this%ids - jime = this%ide - 1 - kims = this%kds - kime = this%kde - 1 - - ioms = this%ids - iome = this%ide - 1 - joms = this%ids - jome = this%ide - 1 - - iops = this%ids - iope = this%ide - 1 - jops = this%ids - jope = this%ide - 1 - ! For compatibility with nuopc couplings - ! pass z_at_w with vertical dim kde - 1 instead of kde - call Calc_fire_wind (this%u3d, this%v3d, this%phl(iims:iime, jims:jime, kims:kime) / G, this%z0, & - iims, iime, jims, jime, kims, kime, config_flags%fire_lsm_zcoupling, config_flags%fire_lsm_zcoupling_ref, & - config_flags%fire_wind_height, ioms, iome, joms, jome, iops, iope, jops, jope, this%ua, this%va) - -! call this%Destroy_z0 () - call this%Destroy_u3d () - call this%Destroy_v3d () - call this%Destroy_phl () + + select case (config_flags%wind_vinterp_opt) + case (0) + call this%Get_u3d (datetime_now) + call this%Get_v3d (datetime_now) + call this%Get_phl (datetime_now) + + ! Set input (i) and output (o) indices + iims = this%ids + iime = this%ide - 1 + jims = this%ids + jime = this%ide - 1 + kims = this%kds + kime = this%kde - 1 + + ioms = this%ids + iome = this%ide - 1 + joms = this%ids + jome = this%ide - 1 + + iops = this%ids + iope = this%ide - 1 + jops = this%ids + jope = this%ide - 1 + ! For compatibility with nuopc couplings + ! pass z_at_w with vertical dim kde - 1 instead of kde + call Calc_fire_wind (this%u3d, this%v3d, this%phl(iims:iime, jims:jime, kims:kime) / G, this%z0, & + iims, iime, jims, jime, kims, kime, config_flags%fire_lsm_zcoupling, config_flags%fire_lsm_zcoupling_ref, & + config_flags%fire_wind_height, ioms, iome, joms, jome, iops, & + iope, jops, jope, this%ua, this%va) + + call this%Destroy_u3d () + call this%Destroy_v3d () + call this%Destroy_phl () +! call this%Destroy_z0 () + + case (1) + call this%Get_u10 (datetime_now) + call this%Get_v10 (datetime_now) + + this%ua = this%u10 + this%va = this%v10 + + case default + call Stop_simulation ('Error: wrong wind_vinterp_opt') + + end select end subroutine Update_atm_state diff --git a/nuopc/fire_behavior_nuopc.F90 b/nuopc/fire_behavior_nuopc.F90 index eb5d1d5..bd8674d 100644 --- a/nuopc/fire_behavior_nuopc.F90 +++ b/nuopc/fire_behavior_nuopc.F90 @@ -861,9 +861,9 @@ subroutine Advance(model, rc) jops = grid%jfps jope = grid%jfpe - call Calc_fire_wind (atm_u3d, atm_v3d, atm_ph / 9.81, grid%fz0, iims, iime, jims, jime, kims, kime, config_flags%fire_lsm_zcoupling, & - config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, ioms, iome, joms, jome, iops, iope, jops, jope, & - grid%uf, grid%vf, cap_winds = .true.) + call Calc_fire_wind (atm_u3d, atm_v3d, atm_ph / 9.81, grid%fz0, iims, iime, jims, jime, kims, kime, & + config_flags%fire_lsm_zcoupling, config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, & + ioms, iome, joms, jome, iops, iope, jops, jope, grid%uf, grid%vf, cap_winds = .true.) case (1) do j = grid%jfps, grid%jfpe @@ -874,6 +874,7 @@ subroutine Advance(model, rc) end do case default call Stop_simulation ('Error: wrong wind_vinterp_opt') + end select if (grid%datetime_now == grid%datetime_start) call grid%Save_state () diff --git a/state/state_mod.F90 b/state/state_mod.F90 index 0d4f238..f825ba7 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -104,6 +104,7 @@ module state_mod integer :: px, py ! Number of MPI tasks in X and Y, respectively contains procedure, public :: Allocate_vars => Allocate_vars + procedure, public :: Apply_wafs => Apply_wafs procedure, public :: Convert_sb_to_ander => Convert_scottburgan_to_anderson procedure, public :: Handle_output => Handle_output procedure, public :: Handle_wrfdata_update => Handle_wrfdata_update @@ -175,6 +176,36 @@ subroutine Allocate_vars (this, ifms, ifme, jfms, jfme) end subroutine Allocate_vars + subroutine Apply_wafs (this) + + implicit none + + class (state_fire_t), intent(in out) :: this + + integer :: i, j, ij, ifts, ifte, jfts, jfte + real :: waf + + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte, waf) + do ij = 1, this%num_tiles + ifts = this%i_start(ij) + ifte = this%i_end(ij) + jfts = this%j_start(ij) + jfte = this%j_end(ij) + + do j = jfts, jfte + do i = ifts, ifte + waf = this%fuels%waf(int (this%nfuel_cat(i, j))) + this%uf(i, j) = waf * this%uf(i, j) + this%vf(i, j) = waf * this%vf(i, j) + end do + end do + end do + !$OMP END PARALLEL DO + + end subroutine Apply_wafs + subroutine Convert_scottburgan_to_anderson (this) implicit none @@ -740,7 +771,7 @@ subroutine Interpolate_vars_atm_to_fire (this, wrf, config_flags) implicit none class (state_fire_t), intent(in out) :: this - type (wrf_t), intent(inout) :: wrf + type (wrf_t), intent(in out) :: wrf type (namelist_t), intent (in) :: config_flags integer :: i, j @@ -761,6 +792,8 @@ subroutine Interpolate_vars_atm_to_fire (this, wrf, config_flags) config_flags%num_tiles, this%i_start, this%i_end, this%j_start, this%j_end, & 'vf', config_flags%hinterp_opt, this%vf) + if (config_flags%wind_vinterp_opt == 1) call this%Apply_wafs () + 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) From 89e4c89b85370c0b2d8b55e4c45339cd915e7009 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Thu, 29 Jan 2026 19:05:26 -0700 Subject: [PATCH 6/9] Fixes interpolation of winds from nuopc --- nuopc/fire_behavior_nuopc.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/nuopc/fire_behavior_nuopc.F90 b/nuopc/fire_behavior_nuopc.F90 index bd8674d..3d7633f 100644 --- a/nuopc/fire_behavior_nuopc.F90 +++ b/nuopc/fire_behavior_nuopc.F90 @@ -860,8 +860,9 @@ subroutine Advance(model, rc) iope = grid%ifpe jops = grid%jfps jope = grid%jfpe - - call Calc_fire_wind (atm_u3d, atm_v3d, atm_ph / 9.81, grid%fz0, iims, iime, jims, jime, kims, kime, & + ! pass the z0 array without halos + ! for compatibility with offline sims + call Calc_fire_wind (atm_u3d, atm_v3d, atm_ph / 9.81, grid%fz0(iims:iime, jims:jime), iims, iime, jims, jime, kims, kime, & config_flags%fire_lsm_zcoupling, config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, & ioms, iome, joms, jome, iops, iope, jops, jope, grid%uf, grid%vf, cap_winds = .true.) From 335d8548f21e1df801b731b255a02a44072e0d8f Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 30 Jan 2026 13:30:34 -0700 Subject: [PATCH 7/9] wafs work in the wrfdata atm nuopc --- nuopc/wrf_nuopc.F90 | 49 +++++++++++++++++++++++++++++++++++---------- 1 file changed, 38 insertions(+), 11 deletions(-) diff --git a/nuopc/wrf_nuopc.F90 b/nuopc/wrf_nuopc.F90 index afb6d89..ce9bad6 100644 --- a/nuopc/wrf_nuopc.F90 +++ b/nuopc/wrf_nuopc.F90 @@ -28,6 +28,8 @@ module wrf_nuopc real(ESMF_KIND_R8), pointer :: ptr_psfc(:,:) real(ESMF_KIND_R8), pointer :: ptr_rain(:,:) real(ESMF_KIND_R8), pointer :: ptr_t2(:,:) + real(ESMF_KIND_R8), pointer :: ptr_u10(:,:) + real(ESMF_KIND_R8), pointer :: ptr_v10(:,:) real(ESMF_KIND_R8), pointer :: ptr_u3d(:,:,:) real(ESMF_KIND_R8), pointer :: ptr_v3d(:,:,:) real(ESMF_KIND_R8), pointer :: ptr_phl(:,:,:) @@ -113,6 +115,13 @@ subroutine Advertise(model, rc) allocate (state%q2(size(state%lats, dim=1), size(state%lats, dim=2))) if (.not. allocated (state%t2)) & allocate (state%t2(size(state%lats, dim=1), size(state%lats, dim=2))) + + if (.not. allocated (state%u10)) & + allocate (state%u10(size(state%lats, dim=1), size(state%lats, dim=2))) + + if (.not. allocated (state%v10)) & + allocate (state%v10(size(state%lats, dim=1), size(state%lats, dim=2))) + if (.not. allocated (state%z0)) & allocate (state%z0(size(state%lats, dim=1), size(state%lats, dim=2))) if (.not. allocated (state%psfc)) & @@ -481,7 +490,7 @@ subroutine Realize(model, rc) file=__FILE__)) & return ! bail out -! ! exportable field on Grid: inst_pres_levels + ! exportable field on Grid: inst_pres_levels ! field = ESMF_FieldCreate(name="inst_pres_levels", grid=grid, & ! gridToFieldMap=(/1,2/), ungriddedLBound=(/1/), & ! ungriddedUBound=(/state%kde - 1/), & @@ -599,34 +608,46 @@ subroutine Realize(model, rc) file=__FILE__)) & return ! bail out - ! exportable field on Grid: inst_pres_height_lowest_from_phys - field = ESMF_FieldCreate(name="inst_pres_height_lowest_from_phys", grid=grid, & + ! exportable field on Grid: inst_zonal_wind_height10m + field = ESMF_FieldCreate(name="inst_zonal_wind_height10m", grid=grid, & typekind=ESMF_TYPEKIND_R8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out call NUOPC_Realize(exportState, field=field, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + ! Get Field memory + call ESMF_FieldGet(field, localDe=0, farrayPtr=ptr_u10, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out - ! exportable field on Grid: inst_spec_humid_height_lowest_from_phys - field = ESMF_FieldCreate(name="inst_spec_humid_height_lowest_from_phys", grid=grid, & + ! exportable field on Grid: inst_merid_wind_height10m + field = ESMF_FieldCreate(name="inst_merid_wind_height10m", grid=grid, & typekind=ESMF_TYPEKIND_R8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out call NUOPC_Realize(exportState, field=field, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + ! Get Field memory + call ESMF_FieldGet(field, localDe=0, farrayPtr=ptr_v10, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out - ! exportable field on Grid: inst_temp_height_lowest_from_phys - field = ESMF_FieldCreate(name="inst_temp_height_lowest_from_phys", grid=grid, & + ! exportable field on Grid: inst_pres_height_lowest_from_phys + field = ESMF_FieldCreate(name="inst_pres_height_lowest_from_phys", grid=grid, & typekind=ESMF_TYPEKIND_R8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & @@ -638,8 +659,8 @@ subroutine Realize(model, rc) file=__FILE__)) & return ! bail out - ! exportable field on Grid: inst_zonal_wind_height10m - field = ESMF_FieldCreate(name="inst_zonal_wind_height10m", grid=grid, & + ! exportable field on Grid: inst_spec_humid_height_lowest_from_phys + field = ESMF_FieldCreate(name="inst_spec_humid_height_lowest_from_phys", grid=grid, & typekind=ESMF_TYPEKIND_R8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & @@ -651,8 +672,8 @@ subroutine Realize(model, rc) file=__FILE__)) & return ! bail out - ! exportable field on Grid: inst_merid_wind_height10m - field = ESMF_FieldCreate(name="inst_merid_wind_height10m", grid=grid, & + ! exportable field on Grid: inst_temp_height_lowest_from_phys + field = ESMF_FieldCreate(name="inst_temp_height_lowest_from_phys", grid=grid, & typekind=ESMF_TYPEKIND_R8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & @@ -814,6 +835,8 @@ subroutine Update_atm_state (datetime) call state%Get_q2(datetime) call state%Get_psfc(datetime) call state%Get_rain(datetime) + call state%Get_u10(datetime) + call state%Get_v10(datetime) call state%Get_u3d(datetime) call state%Get_v3d(datetime) call state%Get_phl(datetime) @@ -830,6 +853,10 @@ subroutine Update_atm_state (datetime) state%rain(1:size(state%lats, dim=1),1:size(state%lats, dim=2)) ptr_t2(clb(1):cub(1),clb(2):cub(2))= & state%t2(1:size(state%lats, dim=1),1:size(state%lats, dim=2)) + ptr_u10(clb(1):cub(1),clb(2):cub(2))= & + state%u10(1:size(state%lats, dim=1),1:size(state%lats, dim=2)) + ptr_v10(clb(1):cub(1),clb(2):cub(2))= & + state%v10(1:size(state%lats, dim=1),1:size(state%lats, dim=2)) ptr_u3d(clb3(1):cub3(1),clb3(2):cub3(2),clb3(3):cub3(3))= & state%u3d(1:size(state%lats, dim=1),1:size(state%lats, dim=2), 1:state%kde - 1) From 5542fef9d66fd8a63090ef6a76a439c497a43071 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 30 Jan 2026 13:45:08 -0700 Subject: [PATCH 8/9] free up 10m winds from memory after calc mid-flame winds in offline sims --- state/state_mod.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/state/state_mod.F90 b/state/state_mod.F90 index f825ba7..9e3af27 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -792,7 +792,11 @@ subroutine Interpolate_vars_atm_to_fire (this, wrf, config_flags) config_flags%num_tiles, this%i_start, this%i_end, this%j_start, this%j_end, & 'vf', config_flags%hinterp_opt, this%vf) - if (config_flags%wind_vinterp_opt == 1) call this%Apply_wafs () + if (config_flags%wind_vinterp_opt == 1) then + call this%Apply_wafs () + call wrf%Destroy_u10 () + call wrf%Destroy_v10 () + 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, & From 38d635a555eef3ab1f2b607ecf29439ebccd3ce9 Mon Sep 17 00:00:00 2001 From: "Pedro A. Jimenez" Date: Fri, 30 Jan 2026 13:55:00 -0700 Subject: [PATCH 9/9] Introduce parameters to define the wind interpolation options --- io/wrf_mod.F90 | 6 +++--- nuopc/fire_behavior_nuopc.F90 | 5 +++-- share/interp_mod.F90 | 5 +++-- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/io/wrf_mod.F90 b/io/wrf_mod.F90 index 222118d..9d92d23 100644 --- a/io/wrf_mod.F90 +++ b/io/wrf_mod.F90 @@ -6,7 +6,7 @@ module wrf_mod use netcdf_mod, only : Get_netcdf_var, Get_netcdf_att, Get_netcdf_dim, Is_netcdf_file_present use proj_lc_mod, only : proj_lc_t use stderrout_mod, only : Print_message, Stop_simulation - use interp_mod, only : Interp_profile + use interp_mod, only : Interp_profile, VINTERP_WINDS_FROM_3D_WINDS, VINTERP_WINDS_FROM_10M_WINDS use coupling_mod, only : Interp_horizontal, Calc_fire_wind implicit none @@ -1132,7 +1132,7 @@ subroutine Update_atm_state (this, datetime_now, config_flags) call this%Get_z0 (datetime_now) select case (config_flags%wind_vinterp_opt) - case (0) + case (VINTERP_WINDS_FROM_3D_WINDS) call this%Get_u3d (datetime_now) call this%Get_v3d (datetime_now) call this%Get_phl (datetime_now) @@ -1166,7 +1166,7 @@ subroutine Update_atm_state (this, datetime_now, config_flags) call this%Destroy_phl () ! call this%Destroy_z0 () - case (1) + case (VINTERP_WINDS_FROM_10M_WINDS) call this%Get_u10 (datetime_now) call this%Get_v10 (datetime_now) diff --git a/nuopc/fire_behavior_nuopc.F90 b/nuopc/fire_behavior_nuopc.F90 index 3d7633f..87fda8c 100644 --- a/nuopc/fire_behavior_nuopc.F90 +++ b/nuopc/fire_behavior_nuopc.F90 @@ -19,6 +19,7 @@ module fire_behavior_nuopc use constants_mod, only : G, XLV, CP, FVIRT, R_D use stderrout_mod, only : Stop_simulation use coupling_mod, only : Calc_fire_wind + use interp_mod, only: VINTERP_WINDS_FROM_3D_WINDS, VINTERP_WINDS_FROM_10M_WINDS implicit none @@ -842,7 +843,7 @@ subroutine Advance(model, rc) #endif select case (config_flags%wind_vinterp_opt) - case (0) + case (VINTERP_WINDS_FROM_3D_WINDS) iims = grid%ifps iime = grid%ifpe @@ -866,7 +867,7 @@ subroutine Advance(model, rc) config_flags%fire_lsm_zcoupling, config_flags%fire_lsm_zcoupling_ref, config_flags%fire_wind_height, & ioms, iome, joms, jome, iops, iope, jops, jope, grid%uf, grid%vf, cap_winds = .true.) - case (1) + case (VINTERP_WINDS_FROM_10M_WINDS) do j = grid%jfps, grid%jfpe do i = grid%ifps, grid%ifpe grid%uf(i,j) = grid%fuels%waf(int(grid%nfuel_cat(i,j))) * ptr_u10(i,j) diff --git a/share/interp_mod.F90 b/share/interp_mod.F90 index f653d50..1ab1e0c 100644 --- a/share/interp_mod.F90 +++ b/share/interp_mod.F90 @@ -6,9 +6,10 @@ module interp_mod private - integer, parameter :: HINTERP_NEAREST = 1, HINTERP_BILINEAR = 2 + integer, parameter :: HINTERP_NEAREST = 1, HINTERP_BILINEAR = 2, VINTERP_WINDS_FROM_3D_WINDS = 0, VINTERP_WINDS_FROM_10M_WINDS = 1 - public :: Interp_profile, Interp_horizontal_nearest, Interp_horizontal_bilinear, HINTERP_NEAREST, HINTERP_BILINEAR + public :: Interp_profile, Interp_horizontal_nearest, Interp_horizontal_bilinear, HINTERP_NEAREST, HINTERP_BILINEAR, & + VINTERP_WINDS_FROM_3D_WINDS, VINTERP_WINDS_FROM_10M_WINDS contains