import numpy as np from scipy.io import netcdf # Settings eruption_file = 'nature14565-s6_stratAerosolMass.txt' volc_samples_file = 'volc_samples.nc' ccsm4_forcing_file = 'CCSM4_volcanic_1850-2008_prototype1.nc' year1 = 1 # first year in forcing file (actual forcing data starts one month earlier) yearn = 9999 # last year in forcing file (actual forcing data ends one month later) nmem = 3 # number of synthetic forcing members output_netcdf = True output_text = True # Read peak stratospheric volcanic aerosol mass loads # These are derived from sulphate depositions of Sigl et al. 2015 and # converted to stratospheric aerosol mass using the Gao et al. 2007 scaling volc_data = np.genfromtxt(eruption_file) year = volc_data[:,0] # extract years flag = volc_data[:,1] # extract location flag (1=tropical, 2=NH, 3=SH) mass = volc_data[:,2] # extract stratospheric volcanic aerosol mass (Tg) # Load the CCSM data for horizontal shape ncin = netcdf.netcdf_file(volc_samples_file,'r',mmap=False) month = ncin.variables['month'][:] St = ncin.variables['Tropical_Volcanoe'][:,:] Sh = ncin.variables['Hemisphere_Volcanoe'][:,:] lat = ncin.variables['lat'][:] ncin.close() # Normalise horizontal shapes by time maximum of global integral area_earth = 4*np.pi*6371000**2 weights = np.cos(lat/180.*np.pi) weights = weights/sum(weights)*area_earth St = St/max(sum(np.transpose(St*weights)))*1e9 # convert from Tg to kg Sh = Sh/max(sum(np.transpose(Sh*weights)))*1e9 # convert from Tg to kg # Load the CCSM data for extracting vertical shape and coordinate info ncin = netcdf.netcdf_file(ccsm4_forcing_file,'r',mmap=False) lev = ncin.variables['lev'][:] lat = ncin.variables['lat'][:] ccsm_volc = ncin.variables['MMRVOLC'][:,:,:] ccsm_colmass = ncin.variables['colmass'][:,:] ncin.close() # Create vertical profile profile_lat = np.zeros((64,8)) for i in range(8): zlev = ccsm_volc[:,i,:]/np.maximum(ccsm_colmass,1e-10) profile_lat[:,i] = np.nanmean(zlev, axis=0) # Create time variable time = [year1-1/24*np.sign(year1-1)] # set first record to December 15th of previous year date = [(year1-1)*10000+(12*100+15)*np.sign(year1-1)] for y in range(year1,yearn+1): for m in range(1,13): time.append(y+m/12-1/24) date.append(y*10000+m*100+15) time.append(yearn+1/24*np.sign(yearn)) # set last record to January 15th of next year date.append((yearn+1)*10000+(1*100+15)*np.sign(yearn+1)) time = np.array(time) date = np.array(date) nmonth = time.shape[0] datesec = np.int32(np.zeros(nmonth)) # create CAM forcing by re-sampling eruption events from reconstruction np.random.seed(0) P=1./2500./12. # probability for erupution event to occur in arbitrary month for member in range(nmem): # Create empty col_mass variable new_colmass = np.zeros((nmonth,ccsm_colmass.shape[1])) if output_text: f = open('CCSM4_volcanic_forcing_{:0>4d}-{:0>4d}_mem{:0>2d}.txt'.format(year1,yearn,member+1), 'w') # loop over months for m in range(nmonth): # loop over eruptions in Sigl el al. 2015 for eruption in range(len(mass)): # check if eruption is triggered and add to column mass if yes if ( np.random.rand() < P ): maxlen=min(St.shape[1],nmonth-m) if flag[eruption] == 1: # tropical eruption new_colmass[m:m+maxlen,:] = new_colmass[m:m+maxlen,:] + St[0:maxlen,:] * mass[eruption] elif flag[eruption] == 2: # NH eruption new_colmass[m:m+maxlen,:] = new_colmass[m:m+maxlen,:] + Sh[0:maxlen,:] * mass[eruption] elif flag[eruption] == 3: # SH eruption new_colmass[m:m+maxlen,:] = new_colmass[m:m+maxlen,:] + np.fliplr(Sh[0:maxlen,:]) * mass[eruption] # dump eruption cronology to text file (columns: year, location, mass load, month) if output_text: year_sim = year1 - 1 if m==0 else year1 + int((m-1)/12) # simulation year month_sim = (m-1)%12+1 # simulation month f.write('{:6d} {:1d} {:6.2f} {:2d}\n'.format(year_sim,int(flag[eruption]),mass[eruption],month_sim)) if output_text: f.close() if output_netcdf: # Create new_mmrvolc from the new_colmass as it was created in Explain new_mmrvolc = np.zeros((time.shape[0],lev.shape[0],lat.shape[0])) for i in range(8): prof_temp = profile_lat[:,i] prof_temp = np.tile(prof_temp,(time.shape[0],1)) new_mmrvolc[:,i,:] = new_colmass * prof_temp # Output new NorESM forcing file # Determine size of dimensions ntime = time.shape[0] nlev = lev.shape[0] nlat = lat.shape[0] # open a new netCDF file for writing. fnout = 'CCSM4_volcanic_forcing_{:0>4d}-{:0>4d}_mem{:0>2d}.nc'.format(year1,yearn,member+1) ncfile = netcdf.netcdf_file(fnout,'w') # create the lat and lon dimensions. ncfile.createDimension('time',None) ncfile.createDimension('date',ntime) ncfile.createDimension('lev',nlev) ncfile.createDimension('lat',nlat) # Define the coordinate variables. time_out = ncfile.createVariable('time',np.dtype('float64').char,('time',)) date_out = ncfile.createVariable('date',np.dtype('float64').char,('date',)) lev_out = ncfile.createVariable('lev',np.dtype('float64').char,('lev',)) lat_out = ncfile.createVariable('lat',np.dtype('float64').char,('lat',)) # Assign units attributes to coordinate variable data. date_out.units = 'yyyymmdd' date_out.FillValue = -9999. lev_out.units = 'hPa' lat_out.units = 'degrees_north' # write data to coordinate vars. time_out[:] = time date_out[:] = date lev_out[:] = lev lat_out[:] = lat # create main variables mmrvolc_out = ncfile.createVariable('MMRVOLC',np.dtype('float64').char,('time','lev','lat')) colmass_out = ncfile.createVariable('colmass',np.dtype('float64').char,('time','lat')) datesec_out = ncfile.createVariable('datesec',np.dtype('int32').char,('time',)) # set the units attribute. mmrvolc_out.units = 'kg kg-1' mmrvolc_out.FillValue = -9999. mmrvolc_out.long_name = 'layer volcanic aerosol mass mixing ratio' colmass_out.FillValue = -999. datesec_out.FillValue = -999 # write data to variables along record (unlimited) dimension. mmrvolc_out[:,::] = new_mmrvolc[:,:] colmass_out[:,::] = new_colmass[:,:] datesec_out[:] = datesec[:] # close the file. ncfile.close()