def roms_to_swan_bathy_curv(hisfile,outfld):
  ''' 
  Generate a SWAN bathymetry file from either a ROMS history or bathymetry input
  file. 
  
  roms_to_swan_bathy_curv(hisfile,outfld)
  
  Parameters
  ----------
  hisfile  : ROMS history or bathymetry input netCDF file
  outfld   : Folder to save the output files
  
  Returns
  -------
  Two text files (swan_bathy.bot, swan_coord.grd) that contain the bathymetry 
  and coordinates of the grid for SWAN input. 
  
  Notes
  -----
    
    
  '''  
   
  # Load variables of interest from the ocean_his.nc file
  ncfile = netCDF4.Dataset(hisfile,'r')  
  h = ncfile.variables['h'][:]
  x_rho = ncfile.variables['x_rho'][:]
  y_rho = ncfile.variables['y_rho'][:]  
  ncfile.close()
  
   
  # Print text file with extended and interpolated bathymetry
  fid = open(outfld+'/swan_bathy.bot', 'w')  
  for aa in range(h.shape[0]):
      for bb in range(h.shape[1]):
          fid.write('%12.4f' % h[aa,bb])
      fid.write('\n')
  fid.close()
  
  # Print text file with extended and interpolated bathymetry
  fid = open(outfld+'/swan_coord.bot', 'w')  
  for aa in range(x_rho.shape[0]):
      for bb in range(x_rho.shape[1]):
          fid.write('%12.6f' % x_rho[aa,bb])
          fid.write('\n')
  for aa in range(y_rho.shape[0]):
      for bb in range(y_rho.shape[1]):
          fid.write('%12.6f' % y_rho[aa,bb])
          fid.write('\n')          
  fid.close()  
  
  #---------------------------------------------------------- Output for swan.in
  print ' '
  print "========================================================"
  print "Created swan_coord.grd and swan_bathy.bot"
  print ('CGRID CURVILINEAR ' + np.str(h.shape[1]-1) + ' ' + 
         np.str(h.shape[0]-1) + ' CIRCLE ...')
  print ('INPGRID BOTTTOM CURVILINEAR 0 0 ' + np.str(h.shape[1]-1) + ' ' + 
         np.str(h.shape[0]-1) + ' EXC ...')
  print "========================================================"  

 def update_task(self, filename):
     #set scheduled
     with open(filename) as f:
         d = dict(filter(None, csv.reader(f,  delimiter=' ', skipinitialspace=True))) 
         taskid = d['taskID']
         AST= d['AST']                 
         f.close()
 
     #update task status in emoncms
     h = httplib2.Http("/tmp/emoncms/.cache")
     minutes=np.int(AST)%3600
     minutes=minutes/60
     hours=np.int(AST)/3600 
          
     request = "{'status':1,'AST':'"+np.str(hours)+":"+np.str(minutes)+"'}"
     h.request("http://localhost/emoncms/mas/update.json?id="+taskid+"&json="+request+"&apikey="+self.apikey, "GET")
     sys.stderr.write("http://localhost/emoncms/mas/update.json?id="+taskid+"&json="+request+"&apikey="+self.apikey)
     #delay should be not 20 but AST-Time.time
     now = datetime.datetime.now()
     midnight = now.replace(hour=0, minute=0, second=0, microsecond=0)
     seconds = (now - midnight).seconds
     countdown = np.int(AST)-seconds
     #countdown =20
     if(countdown > 0):
         t=Timer(countdown, self.taskexec, [taskid])
         t.start()
         self.scheduled[taskid]=0

def read_input_spectra(fname, numbvals):
    # Read grid of model spectra from the <fname> file
    #
    # :rtype: list of tuples containing model spectra.
    #
    # :param fname: name of the file containing the grid of model spectra.
    # :param numbvals: list with total grid number of model parameters.

    d = np.loadtxt(fname)
    parameter_combinations = np.prod(numbvals)
    numspectra = len(d.T) - 1
    # spectra in columns, 1st column is energy grid, skip it
    if numspectra != parameter_combinations:
        raise NameError(
            fname
            + ": No. of spectra "
            + np.str(numspectra)
            + ", different from declared param combinations "
            + np.str(parameter_combinations)
            + "\n"
        )
    input_spectra = []
    for i in range(len(d.T) - 1):
        # d[:-1] - match no. of energy bins;
        # T[1:] - skip the 1st column with energy vector
        input_spectra.append(tuple(d[:-1].T[1:][i]))

    return input_spectra

def load_fortranfile(T):

	if T < 100:
		file1 = '../outputs/psi1/psi1_0'+np.str(T)+'.bin'
		file2 = '../outputs/psi2/psi2_0'+np.str(T)+'.bin'

	else:
		file1 = '../outputs/psi1/psi1_'+np.str(T)+'.bin'
		file2 = '../outputs/psi2/psi2_'+np.str(T)+'.bin'		


	PSI1 = fortranfiles.FortranFile(file1) 
	PSI1 = PSI1.readReals()


	PSI2 = fortranfiles.FortranFile(file2) 
	PSI2 = PSI2.readReals()



	# Redimensionalizing files
	PSI1 = PSI1.reshape((257,512))
	PSI2 = PSI2.reshape((257,512))



	return PSI1,PSI2

def main():
    if len(sys.argv) < 4:
        print_usage("Not enough arguments!")
    alpha = numpy.float64(sys.argv[1])
    filename= numpy.str(sys.argv[2])
    k_limit = 12
    species = []
    for i, s in enumerate(sys.argv[3:]):
        species.append(numpy.str(s))
    species.sort()
    print("species",str(species))
    print("Reading %s" % filename)
    data = atpy_csv(filename)

    if species == 'Overall' or species == 'overall':
        data = filter_by_keyvalue(data, 'alpha', alpha)

    metrs = []
    for s in species:
        print("Calculating %s at %s" % ('mean', s))
        metrs.append(calc_metrics(data, k_limit, alpha, s))
        print(metrs[-1])

    all_metrics = dict(zip(species, metrs))

    metr_fname = "./figures/%s-ALL-metrics-%d-%s.pdf" % (str(species).replace(']','').replace('[','').replace("'",""), k_limit, str(numpy.around(alpha, decimals=2)))
    plot_metrics_per_k(species, alpha, all_metrics, 'Algorithm: ' + filename.split('.')[0], metr_fname)

    pvr_fname = "./figures/%s-ALL-pvr-%d-%s.pdf" % (str(species).replace(']','').replace('[','').replace("'",""), k_limit, str(numpy.around(alpha, decimals=2)))
    print(pvr_fname)
    plot_precision_v_recall(species, alpha, all_metrics, 'Algorithm: ' + filename.split('.')[0], pvr_fname)

    return 0

    def calc_test_press(self, path_const='T'):

        TOL = 1e-3

        Nsamp = 10001
        eos_mod = self.load_eos(path_const=path_const)

        V0, = eos_mod.get_param_values(param_names='V0')
        V0 += -.137
        eos_mod.set_param_values(V0,param_names='V0')

        V0get, = eos_mod.get_param_values(param_names='V0')

        assert V0 == V0get, 'Must be able to store and retrieve non-integer values'

        assert np.abs(eos_mod.press(V0))<TOL/100,(
            'pressure at V0 must be zero by definition'
        )

        Vmod_a = np.linspace(.7,1.2,Nsamp)*V0
        dV = Vmod_a[1] - Vmod_a[0]

        press_a = eos_mod.press(Vmod_a)
        energy_a = eos_mod.energy(Vmod_a)

        abs_err, rel_err, range_err = self.numerical_deriv(
            Vmod_a, energy_a, press_a, scale=-core.CONSTS['PV_ratio'])

        assert range_err < TOL, 'range error in Press, ' + np.str(range_err) + \
            ', must be less than TOL, ' + np.str(TOL)

    def test_press_simple(self, kind_compress='Vinet',
                          compress_path_const='T',
                          kind_gamma='GammaFiniteStrain',
                          kind_RTpoly='V', RTpoly_order=5, natom=1,
                          kind_electronic='CvPowLaw', apply_electronic=True):

        TOL = 1e-3
        Nsamp = 10001
        eos_mod = self.load_eos(kind_compress=kind_compress,
                                compress_path_const=compress_path_const,
                                kind_gamma=kind_gamma, kind_RTpoly=kind_RTpoly,
                                RTpoly_order=RTpoly_order, natom=natom,
                                kind_electronic=kind_electronic,
                                apply_electronic=apply_electronic)

        refstate_calc = eos_mod.calculators['refstate']
        T0 = refstate_calc.ref_temp()
        V0 = refstate_calc.ref_volume()
        S0 = refstate_calc.ref_entropy()
        # V0, T0, S0 = eos_mod.get_param_values(param_names=['V0','T0','S0'])

        Vmod_a = np.linspace(.7,1.2,Nsamp)*V0
        T = 8000
        dV = Vmod_a[1] - Vmod_a[0]

        P_a = eos_mod.press(Vmod_a, T)
        F_a = eos_mod.helmholtz_energy(Vmod_a, T)
        abs_err, rel_err, range_err = self.numerical_deriv(
              Vmod_a, F_a, P_a, scale=-core.CONSTS['PV_ratio'])

        S_a = eos_mod.entropy(Vmod_a, T)


        assert abs_err < TOL, ('abs error in Press, ' + np.str(abs_err) +
                                 ', must be less than TOL, ' + np.str(TOL))

def freq_spec_1d(eta,dt=1,verbose=True):
    """
    Computes the frequency spectrum from a given time series.
    
    freq,spec = freq_spec_1d(eta,dt,verbose)
    
    PARAMETERS:
    -----------
    eta      : Time series of water surface elevation [m]
    dt       : Time step [s]
    verbose  : Display computed bulk parameters to the screen
    
    RETURNS:
    --------
    freq     : Frequency vector
    spec     : Variance spectrum (Power spectrum)
    
    NOTES:
    ------
    This is really a copy of gsignal.psdraw. If results differ, trust gsignal
      this code will not be updated.
    """
    
    # Remove mean
    eta -= eta.mean()

    # Compute record length
    N = eta.shape[0]
    
    # Compute fourier frequencies
    fj = np.fft.fftfreq(N,dt)

    # Compute power spectral density (Cooley-Tukey Method)
    yf = np.fft.fft(eta)/N
    psd = N*dt*yf*np.conjugate(yf)

    # One sided psd from dft
    if np.mod(N,2) == 0:
        sf = np.concatenate((np.array([psd[0]]),2.0*psd[1:N/2],
                             np.array([psd[N/2]])))
        freq_amp = np.abs(np.concatenate((np.array([fj[0]]),fj[1:N/2],
                                          np.array([fj[N/2]]))))
    else:
        sf = np.concatenate((np.array([psd[0]]),2.0*psd[1:(N+1)/2]))
        freq_amp = np.abs(np.concatenate((np.array([fj[0]]),fj[1:(N+1)/2])))

    sf = sf.real
    
    if verbose:
        print("===============================================")
        print("Bulk Wave Parameters:")
        print("Hs = " + np.str(4.004*np.sqrt(np.trapz(sf,freq_amp))) + "m")
        print("H1 = "+np.str(4.004*np.sqrt(np.trapz(sf,freq_amp))*2.0/3.0)+"m")
        print("Spectral Parameters:")
        print("Nyquist Frequency = " + np.str(1.0/(2.0*dt)) + "Hz")
        print("Frequency interval = " + np.str(1.0/(N*dt)) + "Hz")
        print("===============================================")

    # End of function
    return freq_amp,sf

    def _calc_test_heat_capacity(self, kind_compress='Vinet',
                                 compress_path_const='T',
                                 kind_gamma='GammaFiniteStrain',
                                 kind_RTpoly='V', RTpoly_order=5, natom=1,
                                 kind_electronic='None',
                                 apply_electronic=False):

        TOL = 1e-3
        Nsamp = 10001

        eos_mod = self.load_eos(kind_compress=kind_compress,
                                compress_path_const=compress_path_const,
                                kind_gamma=kind_gamma, kind_RTpoly=kind_RTpoly,
                                RTpoly_order=RTpoly_order, natom=natom,
                                kind_electronic=kind_electronic,
                                apply_electronic=apply_electronic)

        Tmod_a = np.linspace(3000.0, 8000.0, Nsamp)

        V0, = eos_mod.get_param_values(param_names=['V0'])
        # Vmod = V0*(0.6+.5*np.random.rand(Nsamp))
        Vmod = V0*0.7

        thermal_energy_a = eos_mod.thermal_energy(Vmod, Tmod_a)
        heat_capacity_a = eos_mod.heat_capacity(Vmod, Tmod_a)

        abs_err, rel_err, range_err = self.numerical_deriv(
            Tmod_a, thermal_energy_a, heat_capacity_a, scale=1)

        Cvlimfac = eos_mod.calculators['thermal']._get_Cv_limit()
        assert rel_err < TOL, 'rel-error in Cv, ' + np.str(rel_err) + \
            ', must be less than TOL, ' + np.str(TOL)

    def calc_test_RTcoefs(self, kind_compress='Vinet',
                           compress_path_const='T',
                           kind_gamma='GammaFiniteStrain', kind_RTpoly='V',
                           RTpoly_order=5, natom=1, kind_electronic='None',
                           apply_electronic=False):

        TOL = 1e-3

        Nsamp = 10001
        eos_mod = self.load_eos(kind_compress=kind_compress,
                                compress_path_const=compress_path_const,
                                kind_gamma=kind_gamma, kind_RTpoly=kind_RTpoly,
                                RTpoly_order=RTpoly_order, natom=natom,
                                kind_electronic=kind_electronic,
                                apply_electronic=apply_electronic)

        V0, = eos_mod.get_param_values(param_names='V0')
        Vmod_a = np.linspace(.5,1.2,Nsamp)*V0
        dV = Vmod_a[1] - Vmod_a[0]

        bcoef_a = eos_mod.calc_RTcoefs(Vmod_a)
        bcoef_deriv_a = eos_mod.calc_RTcoefs_deriv(Vmod_a)

        b_abs_err, b_rel_err, b_range_err = self.numerical_deriv(
            Vmod_a, bcoef_a, bcoef_deriv_a, scale=1)

        assert b_range_err < TOL, 'range error in bcoef, ' + \
            np.str(b_range_err) + ', must be less than TOL, ' + np.str(TOL)

    def test_RTcoefs(self, kind_compress='Vinet', compress_order=3,
                     compress_path_const='T', kind_RTpoly='V',
                     RTpoly_order=5, natom=1):

        TOL = 1e-3

        Nsamp = 10001
        eos_mod = self.load_eos(kind_compress=kind_compress,
                                compress_order=compress_order,
                                compress_path_const=compress_path_const,
                                kind_RTpoly=kind_RTpoly,
                                RTpoly_order=RTpoly_order,
                                natom=natom)

        V0, = eos_mod.get_param_values(param_names='V0')
        Vmod_a = np.linspace(.5,1.2,Nsamp)*V0
        dV = Vmod_a[1] - Vmod_a[0]

        acoef_a, bcoef_a = eos_mod.calc_RTcoefs(Vmod_a)
        acoef_deriv_a, bcoef_deriv_a = eos_mod.calc_RTcoefs_deriv(Vmod_a)

        a_abs_err, a_rel_err, a_range_err = self.numerical_deriv(
            Vmod_a, acoef_a, acoef_deriv_a, scale=1)

        b_abs_err, b_rel_err, b_range_err = self.numerical_deriv(
            Vmod_a, bcoef_a, bcoef_deriv_a, scale=1)


        assert a_range_err < TOL, 'range error in acoef, ' + \
            np.str(a_range_err) + ', must be less than TOL, ' + np.str(TOL)

        assert b_range_err < TOL, 'range error in bcoef, ' + \
            np.str(b_range_err) + ', must be less than TOL, ' + np.str(TOL)

def laserScan2D(width, height, delta, peakArea, fileName):
    # Scans an area with given width and height at a
    # step rate of delta. peakArea is int value for
    # the location of peak with a scale from 0 to 10,000

    startTime = time.clock()
    h = 0
    w = 0
    n = 0
    m = 0

    x = np.arange(0, width + delta, delta)
    y = np.arange(0, height + delta, delta)
    Y, X = np.meshgrid(y, x)
    tValues = np.zeros((np.size(x), np.size(y)))
    vValues = np.zeros((np.size(x), np.size(y)))

    # set up scope
    scanRange = 1000
    scope = vi.instrument("TCPIP::138.67.12.235::INSTR")
    sRead.setParam(scope, peakArea - scanRange, peakArea + scanRange)

    # get motor and zero location
    motor = mC.setupMotor()

    while w <= width:
        h = 0
        m = 0
        while h <= height:
            mC.moveTo(motor, w, h)
            time.sleep(0.5)
            x, y = sRead.getData(scope, peakArea - scanRange, peakArea + scanRange)
            t, v = findPeak(x, y)
            tValues[n, m] = t
            vValues[n, m] = v
            h = h + delta
            m = m + 1
        w = w + delta
        n = n + 1

        # Estimates Time Left
        timeLeft = (width - w) / w * (time.clock() - startTime) / 60
        print "Est. Time Left " + np.str(timeLeft) + "min"
    mC.moveTo(motor, 0, 0)

    # Contour Plot of Time
    makePlot2D(X, Y, tValues, fileName + " Time")

    # Contour Plot of Voltage
    makePlot2D(X, Y, vValues, fileName + " Voltage")

    # File Output
    np.savez(fileName + ".npz", X=X, Y=Y, tValues=tValues, vValues=vValues)

    # Time Taken Calc
    timeTaken = (time.clock() - startTime) / 60  # in min
    print "Time Taken " + np.str(timeTaken)
    motor.close()
    scope.close()
    return timeTaken, tValues

def wave_length(period, h, verbose=True):
    """
    Compute wave length using linear wave theory

    Parameters
    ----------
    period   : wave period [s]
    h        : water depth [m]

    Results
    -------
    wl_int   : real wave length [m]

    Screen output
    -------------
    wl_deep  : deep water wave length [m]
    wl_sha   : shallow water wave length [m]

    """

    wl_deep = 9.81 * period ** 2 / 2.0 / np.pi
    wl_sha = period * np.sqrt(9.81 * h)
    k = dispersion(period, h)
    wl_int = 9.81 / 2.0 / np.pi * period ** 2 * np.tanh(k * h)

    if verbose:
        print(" ")
        print("---------------------------------------------------------")
        print("Wave Length deep water approx      = " + np.str(wl_deep) + " m")
        print("Wave Length shallow water approx   = " + np.str(wl_sha) + " m")
        print("Wave Length linear wave theory     = " + np.str(wl_int) + " m")
        print("---------------------------------------------------------")
        print(" ")

    return wl_int

def wave_length(period,h):
    '''
    Compute wave length using linear wave theory
    
    Parameters
    ----------
    period   : wave period [s]
    h        : water depth [m]
    
    Results
    -------
    wl_int   : real wave length [m]
    
    Screen output
    -------------
    wl_deep  : deep water wave length [m]
    wl_sha   : shallow water wave length [m]
    
    '''
    
    wl_deep = 9.81 * period**2 / 2.0 / np.pi
    wl_sha = period * np.sqrt(9.81 * h)
    k = dispersion(period,h)
    wl_int = 9.81 / 2.0 / np.pi * period**2 * np.tanh(k*h)
    
    print(' ')
    print('---------------------------------------------------------')
    print('Wave Length deep water approx      = ' + np.str(wl_deep) + ' m')
    print('Wave Length shallow water approx   = ' + np.str(wl_sha) + ' m')
    print('Wave Length linear wave theory     = ' + np.str(wl_int) + ' m')
    print('---------------------------------------------------------')
    print(' ')
    
    return wl_int

def save_orbit( x, y, z, filename ):
	ff = open( filename + '.3d', 'w' )
	for i in range(len(x)):
		ff.write( np.str(x[i]) + "," +
				  np.str(y[i]) + "," +
				  np.str(z[i]) + "\n" )
	ff.close()
	return

def iter_to_str(iteration, maximum):
    """ Converts an iteration number to string.

    Uses the maximum as second input to guarantee equal length for all.

    """
    cur_trial_len = len(np.str(iteration))
    return ((len(np.str(np.int(maximum)+1))-cur_trial_len) * '0') + np.str(iteration)

def time_lag(eta,ot,lags=None):
    """
    Function to compute average time lag between the wave staffs
    
    USAGE:
    ------
    ot_lag = time_lag(eta,ot,lags)
    
    PARAMETERS:
    -----------
    eta    : Numpy array of water surface elevation and time
             eta.shape = (time,points)
    ot     : Time vector (numpy array)
    lags   : Number of lags to compute
    
    RETURNS:
    --------
    ot_lag : Numpy array of the same dimensions as eta with the time lagged
             arrays.
    
    DEPENDENCIES:
    -------------
    gsignal.cross_corr
    
    """
    
    # Verify the requested lags
    if not lags:
        lags = np.floor(ot.shape[0]/2)
    
    # Cumulative lag time 
    cum_lag_time = np.zeros((eta.shape[1],))
    
    # Time interval
    dt = ot[2] - ot[1]
    
    # Loop over points
    for aa in range(1,cum_lag_time.shape[0]):
        
        # Find the time lagged cross-correlation to adjust the time series
        rho,stats = gsignal.cross_corr(eta[:,aa-1],eta[:,aa],lags)
        
        # Identify the maximum auto correlation
        if np.max(rho) < 0.8:
            print('Warning: Correlation is less than 0.8')
            print('  aa = ' + np.str(aa))
            print('  r = ' + np.str(np.max(rho)))
            
        # Compute cumulative lag time
        cum_lag_time[aa] = cum_lag_time[aa-1] + stats[np.argmax(rho),0] * dt
    
    # Create output array based on lag time
    ot_lag = np.zeros_like(eta)
    for aa in range(cum_lag_time.shape[0]):
        ot_lag[:,aa] = ot - cum_lag_time[aa]
        
    # Exit function
    return ot_lag

def makerunfake(rundir, base, param_file, nstart, nruns):
	for i in range(nstart, nstart+nruns):
		fakeparam = param_file+".fake_"+np.str(i)
		outfile = "runfake"+np.str(i)
		f = open(outfile, 'w')
		f.write("cd " + rundir+"\n")
		f.write("dolphot " + base+"_"+np.str(i)+ " -p" + fakeparam + " >> fake.log_"+np.str(i))
		f.close()
		subprocess.call("chmod +x " + outfile, shell=True)

    def computeDFF(self,secsWindow=5,quantilMin=8,method='only_baseline',order='C'):
        """
        compute the DFF of the movie or remove baseline
        In order to compute the baseline frames are binned according to the window length parameter
        and then the intermediate values are interpolated.
        Parameters
        ----------
        secsWindow: length of the windows used to compute the quantile
        quantilMin : value of the quantile
        method='only_baseline','delta_f_over_f','delta_f_over_sqrt_f'

        Returns
        -----------
        self: DF or DF/F or DF/sqrt(F) movies
        movBL=baseline movie
        """

        print("computing minimum ..."); sys.stdout.flush()
        minmov=np.min(self)

        if np.min(self)<=0 and method != 'only_baseline':
            raise ValueError("All pixels must be positive")

        numFrames,linePerFrame,pixPerLine=np.shape(self)
        downsampfact=int(secsWindow*self.fr);
        elm_missing=int(np.ceil(numFrames*1.0/downsampfact)*downsampfact-numFrames)
        padbefore=int(np.floor(old_div(elm_missing,2.0)))
        padafter=int(np.ceil(old_div(elm_missing,2.0)))

        print(('Inizial Size Image:' + np.str(np.shape(self)))); sys.stdout.flush()
        self=movie(np.pad(self,((padbefore,padafter),(0,0),(0,0)),mode='reflect'),**self.__dict__)
        numFramesNew,linePerFrame,pixPerLine=np.shape(self)

        #% compute baseline quickly
        print("binning data ..."); sys.stdout.flush()
#        import pdb
#        pdb.set_trace()
        movBL=np.reshape(self,(downsampfact,int(old_div(numFramesNew,downsampfact)),linePerFrame,pixPerLine),order=order);
        movBL=np.percentile(movBL,quantilMin,axis=0);
        print("interpolating data ..."); sys.stdout.flush()
        print((movBL.shape))

        movBL=scipy.ndimage.zoom(np.array(movBL,dtype=np.float32),[downsampfact ,1, 1],order=0, mode='constant', cval=0.0, prefilter=False)

        #% compute DF/F
        if method == 'delta_f_over_sqrt_f':
            self=old_div((self-movBL),np.sqrt(movBL))
        elif method == 'delta_f_over_f':
            self=old_div((self-movBL),movBL)
        elif method  =='only_baseline':
            self=(self-movBL)
        else:
            raise Exception('Unknown method')

        self=self[padbefore:len(movBL)-padafter,:,:];
        print(('Final Size Movie:' +  np.str(self.shape)))
        return self,movie(movBL,fr=self.fr,start_time=self.start_time,meta_data=self.meta_data,file_name=self.file_name)

def roms_to_swan_bathy_curv(hisfile, outfld):
    """ 
  Generate a SWAN bathymetry file from either a ROMS history or bathymetry input
  file. 
  
  roms_to_swan_bathy_curv(hisfile,outfld)
  
  Parameters
  ----------
  hisfile  : ROMS history or bathymetry input netCDF file
  outfld   : Folder to save the output files
  
  Returns
  -------
  Two text files (swan_bathy.bot, swan_coord.grd) that contain the bathymetry 
  and coordinates of the grid for SWAN input. 
  
  Notes
  -----
    
    
  """

    # Load variables of interest from the ocean_his.nc file
    ncfile = netCDF4.Dataset(hisfile, "r")
    h = ncfile.variables["h"][:]
    x_rho = ncfile.variables["x_rho"][:]
    y_rho = ncfile.variables["y_rho"][:]
    ncfile.close()

    # Print text file with extended and interpolated bathymetry
    fid = open(outfld + "/swan_bathy.bot", "w")
    for aa in range(h.shape[0]):
        for bb in range(h.shape[1]):
            fid.write("%12.4f" % h[aa, bb])
        fid.write("\n")
    fid.close()

    # Print text file with extended and interpolated bathymetry
    fid = open(outfld + "/swan_coord.bot", "w")
    for aa in range(x_rho.shape[0]):
        for bb in range(x_rho.shape[1]):
            fid.write("%12.6f" % x_rho[aa, bb])
            fid.write("\n")
    for aa in range(y_rho.shape[0]):
        for bb in range(y_rho.shape[1]):
            fid.write("%12.6f" % y_rho[aa, bb])
            fid.write("\n")
    fid.close()

    # ---------------------------------------------------------- Output for swan.in
    print(" ")
    print("=====================================================================")
    print("Created swan_coord.grd and swan_bathy.bot")
    print("CGRID CURVILINEAR " + np.str(h.shape[1] - 1) + " " + np.str(h.shape[0] - 1) + " CIRCLE ...")
    print("INPGRID BOTTTOM CURVILINEAR 0 0 " + np.str(h.shape[1] - 1) + " " + np.str(h.shape[0] - 1) + " EXC ...")
    print("=====================================================================")