Skip to content

Add Scaling Parameters to Data-Driven Storms #671

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Sep 7, 2025
Merged

Add Scaling Parameters to Data-Driven Storms #671

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7b0cb6e
Add scaling parameters to data based storms
mandli Sep 4, 2025
50215f8
Add read function for surge.data
mandli Sep 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 44 additions & 37 deletions src/2d/shallow/surge/data_storm_module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,12 @@
!
! ==============================================================================
! Copyright (C) Clawpack Developers 2017
! Distributed under the terms of the Berkeley Software Distribution (BSD)
! Distributed under the terms of the Berkeley Software Distribution (BSD)
! license
! http://www.opensource.org/licenses/
! ==============================================================================
module data_storm_module

implicit none
save

Expand Down Expand Up @@ -45,6 +45,9 @@ module data_storm_module

integer, private :: last_storm_index

! Run-time storm modification settings
real(kind=8) :: wind_scaling, pressure_scaling

! Ramping settings
integer :: window_type
real(kind=8) :: window(4), ramp_width
Expand All @@ -63,7 +66,7 @@ module data_storm_module
! Initializes the storm type for an HWRF type data derived storm that is
! saved as a netcdf format
subroutine set_storm(storm_data_path, storm, storm_spec_type, log_unit)

#ifdef NETCDF
use netcdf
#endif
Expand All @@ -86,7 +89,7 @@ subroutine set_storm(storm_data_path, storm, storm_spec_type, log_unit)

! ASCII
real(kind=8) :: sw_lon, sw_lat, dx, dy

! NetCDF
integer :: nc_fid, num_dims, num_vars, var_id
character(len=16) :: name, dim_names(3), var_names(3)
Expand All @@ -105,12 +108,14 @@ subroutine set_storm(storm_data_path, storm, storm_spec_type, log_unit)
read(data_unit, "(a)") storm%time_offset
read(data_unit, "(i2)") file_format
read(data_unit, "(i2)") num_data_files
read(data_unit, *) wind_scaling, pressure_scaling
read(data_unit, "(i2)") window_type
read(data_unit, *) ramp_width

write(log_unit, "('time_offset = ',a)") storm%time_offset
write(log_unit, "('Format = ',i2)") file_format
write(log_unit, "('Num files = ',i2)") num_data_files
write(log_unit, "('format = ',i2)") file_format
write(log_unit, "('num files = ',i2)") num_data_files
write(log_unit, "('scaling = ',2d16.8)") wind_scaling, pressure_scaling
write(log_unit, "('window_type = ',i2)") window_type
write(log_unit, "('ramp_width = ',d16.8)") ramp_width

Expand All @@ -121,7 +126,7 @@ subroutine set_storm(storm_data_path, storm, storm_spec_type, log_unit)
read(data_unit, *)
end if
read(data_unit, *)

if (DEBUG) then
print "('time_offset = ',a)", storm%time_offset
print "('Format = ',i2)", file_format
Expand Down Expand Up @@ -201,13 +206,13 @@ subroutine set_storm(storm_data_path, storm, storm_spec_type, log_unit)
storm, seconds_from_epoch(time(:, 1)))

case(2)
#ifdef NETCDF
#ifdef NETCDF
! Open file and get file ID
! :TODO: Only read in times that are between t0 and tfinal
print *, "Reading storm NetCDF file ", storm%paths(1)
call check_netcdf_error(nf90_open(storm%paths(1), nf90_nowrite, nc_fid))
! Check dim/var number
call check_netcdf_error(nf90_inquire(nc_fid, num_dims, num_vars))
call check_netcdf_error(nf90_inquire(nc_fid, num_dims, num_vars))
if (num_dims /= 3 .and. num_vars /= 3) then
print *, "Invalid number of dimensions/variables."
print *, " num_dims = ", num_dims
Expand All @@ -230,7 +235,7 @@ subroutine set_storm(storm_data_path, storm, storm_spec_type, log_unit)
allocate(storm%longitude(mx))
allocate(storm%latitude(my))
allocate(storm%time(mt))

write(log_unit, *) "Storm header info:"
write(log_unit, "(' mx,mx,mt = ', 3i4)") mx, my, mt

Expand Down Expand Up @@ -300,7 +305,7 @@ end subroutine set_storm
! data structure
subroutine read_OWI_ASCII_header(pressure_file, mx, my, mt, swlon, swlat, &
dy, dx)

use utility_module, only: seconds_from_epoch

implicit none
Expand All @@ -321,20 +326,20 @@ subroutine read_OWI_ASCII_header(pressure_file, mx, my, mt, swlon, swlat, &
"i4,i2,i2,i2,i2)"
character(len=*), parameter :: time_info_format = &
"(t69,i4,i2,i2,i2,i2)"

! Local
integer, parameter :: OWI_unit = 156
integer :: i, time(5, 2), total_time, dt


! Read in start and end dates of file from first header line
open(OWI_unit, file=pressure_file, status='old', action='read')

! Total number of time in seconds of the dataset
read(OWI_unit, header_format) time(1:4, 1), time(1:4, 2)
read(OWI_unit, header_format) time(1:4, 1), time(1:4, 2)
total_time = seconds_from_epoch(time(1:4, 2)) &
- seconds_from_epoch(time(1:4, 1))

! Read second header from file to obtain array sizes and first timestep
read(OWI_unit, full_info_format) my, mx, dx, dy, swlat, swlon, time(:, 1)

Expand All @@ -343,32 +348,32 @@ subroutine read_OWI_ASCII_header(pressure_file, mx, my, mt, swlon, swlat, &
read(OWI_unit, *)
end do
read(OWI_unit, time_info_format) time(:, 2)

close(OWI_unit)

dt = seconds_from_epoch(time(:, 2)) - seconds_from_epoch(time(:, 1))
mt = (total_time / dt) + 1

end subroutine read_OWI_ASCII_header
end subroutine read_OWI_ASCII_header


! === read_OWI_ASCII =======================================================
! reads the data files and fills out the storm object and it's dataarrays
subroutine read_OWI_ASCII(wind_file, pressure_file, mx, my, mt, &
swlat, swlon, &
dx, dy, storm, &
seconds_from_offset)

use utility_module, only: seconds_from_epoch

implicit none

! Input arguments
type(data_storm_type) :: storm
character(len=*), intent(in) :: wind_file, pressure_file
integer, intent(in) :: my, mx, mt, seconds_from_offset
integer, intent(in) :: my, mx, mt, seconds_from_offset
real(kind=8), intent(in) :: swlat, swlon, dx, dy

! Local storage
integer, parameter :: wind_unit = 700, pressure_unit = 800
integer :: i, j, n, time(5)
Expand All @@ -380,17 +385,17 @@ subroutine read_OWI_ASCII(wind_file, pressure_file, mx, my, mt, &
! Skip file headers
read(wind_unit, *)
read(pressure_unit, *)
! Loop over all timesteps
! Loop over all timesteps
do n = 1, mt
! Read each time from the next array
read(wind_unit, '(t69, i4,i2,i2,i2, i2)') time
storm%time(n) = seconds_from_epoch(time) - seconds_from_offset

read(wind_unit, '(8f10.0)') ((storm%wind_u(i,j, n),i=1,mx),j=1,my)
read(wind_unit, '(8f10.0)') ((storm%wind_v(i,j, n),i=1,mx),j=1,my)

read(pressure_unit, *) ! Skip header line since we have it from above
read(pressure_unit, '(8f10.0)') ((storm%pressure(i,j, n),i=1,mx),j=1,my)
read(pressure_unit, '(8f10.0)') ((storm%pressure(i,j, n),i=1,mx),j=1,my)
! Convert from Pa to hPa
storm%pressure(:,:,n) = storm%pressure(:,:,n) * 100.0
end do
Expand All @@ -405,7 +410,7 @@ subroutine read_OWI_ASCII(wind_file, pressure_file, mx, my, mt, &
do j = 1, mx
storm%longitude(j) = swlon + j * dx
end do

end subroutine read_OWI_ASCII


Expand Down Expand Up @@ -575,9 +580,11 @@ subroutine set_data_fields(maux, mbc, mx, my, xlower, ylower, &
ramp = 0.d0
end if
aux(pressure_index,i,j) = &
Pa + (aux(pressure_index,i,j) - Pa) * ramp
Pa + (aux(pressure_index,i,j) - Pa) &
* ramp * pressure_scaling
aux(wind_index:wind_index+1,i,j) = &
aux(wind_index:wind_index+1,i,j) * ramp
aux(wind_index:wind_index+1,i,j) &
* ramp * wind_scaling

else
aux(wind_index:wind_index + 1, i, j) = 0.d0
Expand Down Expand Up @@ -630,7 +637,7 @@ end subroutine set_netcdf_fields
! === interp_time ==========================================================
! Interpolate storm arrays to get current time
subroutine interp_time(storm, t, wind_u, wind_v, pressure)

implicit none
! Subroutine IO
type(data_storm_type), intent(in) :: storm
Expand Down Expand Up @@ -688,9 +695,9 @@ end subroutine interp_time
! === spatial_intrp ========================================================
! Spatially interpolate interp_array onto (x,y) point
pure real(kind=8) function spatial_interp(storm, x, y, interp_array) result(value)

implicit none

! Subroutine I/O
type(data_storm_type), intent(in) :: storm
real(kind=8), intent(in) :: x, y, interp_array(size(storm%longitude), &
Expand Down Expand Up @@ -731,7 +738,7 @@ pure real(kind=8) function spatial_interp(storm, x, y, interp_array) result(valu
endif

value = interp_array(xidx_low, yidx_low)

else
call find_nearest(storm, x - storm_dx, y - storm_dy, llon, llat, xidx_low, yidx_low)
call find_nearest(storm, x + storm_dx, y + storm_dy, ulon, ulat, xidx_high, yidx_high)
Expand All @@ -740,15 +747,15 @@ pure real(kind=8) function spatial_interp(storm, x, y, interp_array) result(valu
q12 = interp_array(xidx_low, yidx_high)
q21 = interp_array(xidx_high, yidx_low)
q22 = interp_array(xidx_high, yidx_high)

! Calculate the value at the center of the box using bilinear interpolation
value = (q11 * (ulon - x) * (ulat - y) + &
q21 * (x - llon) * (ulat - y) + &
q12 * (ulon - x) * (y - llat) + &
q22 * (x - llon) * (y - llat)) / ((ulon - llon) * (ulat - llat) + 0.00)

end if

end function spatial_interp

! === find_nearest =========================================================
Expand All @@ -768,4 +775,4 @@ pure subroutine find_nearest(storm, x, y, lon, lat, xidx, yidx)
lat = storm%latitude(yidx)
end subroutine find_nearest

end module data_storm_module
end module data_storm_module
54 changes: 42 additions & 12 deletions src/python/geoclaw/data.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
"""

import os
from pathlib import Path
import numpy as np
import warnings

Expand Down Expand Up @@ -575,6 +576,43 @@ def __init__(self):
self.add_attribute('storm_specification_type', 0) # Type of parameterized storm
self.add_attribute("storm_file", None) # File containing data

def read(self, path: Path=Path("surge.data"), force: bool=False):
"""Read surge data file"""

with Path(path).open() as data_file:
# Header
data_file.readline()
data_file.readline()
data_file.readline()
data_file.readline()
data_file.readline()
data_file.readline()

self.wind_forcing = bool(data_file.readline())
self.drag_law = int(data_file.readline().split("=:")[0])
self.pressure_forcing = bool(data_file.readline().split("=:")[0])
self.rotation_override = data_file.readline().split("=:")[0]
data_file.readline()

self.wind_index = int(data_file.readline().split("=:")[0]) - 1
self.pressure_index = int(data_file.readline().split("=:")[0]) - 1
self.display_landfall_time = bool(data_file.readline().split("=:")[0])
data_file.readline()

# AMR parameters
self.wind_refine = self._parse_value(data_file.readline())
self.R_refine = self._parse_value(data_file.readline())
data_file.readline()

# Storm specification
self.storm_specification_type = int(data_file.readline().split("=:")[0])
line = data_file.readline().split("=:")[0]
if line[0] == "'":
self.storm_file = line.strip()[1:-1]
else:
raise IOError("Error reading storm file name.")


def write(self, out_file='surge.data', data_source="setrun.py"):
"""Write out the data file to the path given"""

Expand Down Expand Up @@ -697,21 +735,14 @@ def read(self, path="friction.data", force=False):
# Regions
self.friction_regions = []
for n in range(num_regions):
lower = self._convert_line(data_file.readline())
upper = self._convert_line(data_file.readline())
depths = self._convert_line(data_file.readline())
coeff = self._convert_line(data_file.readline())
lower = self._parse_value(data_file.readline())
upper = self._parse_value(data_file.readline())
depths = self._parse_value(data_file.readline())
coeff = self._parse_value(data_file.readline())
self.friction_regions.append([lower, upper, depths, coeff])
data_file.readline()
self.friction_files = [] # Is not supported

def _convert_line(self, line):
values = []
for value in line.split("=:")[0].split(" "):
if len(value) > 1:
values.append(float(value))
return values


def write(self, out_file='friction.data', data_source='setrun.py'):

Expand Down Expand Up @@ -911,4 +942,3 @@ def write(self,out_file='bouss.data',data_source='setrun.py'):
self.data_write('bouss_min_depth')

self.close_data_file()

Loading
Loading