diff --git a/CMakeLists.txt b/CMakeLists.txt index 8b17730..818b257 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/io/coupling_mod.F90 b/io/coupling_mod.F90 new file mode 100644 index 0000000..bd912f4 --- /dev/null +++ b/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 diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index b91d2b1..cfd628b 100644 --- a/io/namelist_mod.F90 +++ b/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 diff --git a/io/wrf_mod.F90 b/io/wrf_mod.F90 index 757db4e..05217d5 100644 --- a/io/wrf_mod.F90 +++ b/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) diff --git a/share/interp_mod.F90 b/share/interp_mod.F90 index 57932b3..f653d50 100644 --- a/share/interp_mod.F90 +++ b/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) diff --git a/state/state_mod.F90 b/state/state_mod.F90 index b928fc0..d73323e 100644 --- a/state/state_mod.F90 +++ b/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