Merge pull request #671 from mandli/add-scaling-storm

Add Scaling Parameters to Data-Driven Storms

Author: rjleveque

src/2d/shallow/surge/data_storm_module.f90

@@ -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
 
@@ -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
@@ -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
@@ -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)
@@ -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
 
@@ -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
@@ -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
@@ -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
 
@@ -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
@@ -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)
 
@@ -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)
@@ -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
@@ -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
 
 
@@ -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
@@ -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
@@ -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),  &
@@ -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)
@@ -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 =========================================================
@@ -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

src/python/geoclaw/data.py

@@ -29,6 +29,7 @@
 """
 
 import os
+from pathlib import Path
 import numpy as np
 import warnings
 
@@ -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"""
 
@@ -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'):
 
@@ -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()
-