""" Tool box for processing NorCPM1 CMIP6 output. """ # import libraries import numpy as np from netCDF4 import Dataset import os import warnings # surpress warnings warnings.filterwarnings("ignore") def setGlobals(baseVal=os.path.realpath('../..'),modelVal='noresm1-me'): global base, dataRaw, dataTmp, dataObs, model, gridFileOcn, gridFileAtm\ ,testExperiment, testMember, testLonRange, testLatRange, testLevRange\ ,testYearRange, testMemRange, testInputFile\ ,daysInMonth, weightsAnnMean # model selection model = modelVal # default # path settings etc base = baseVal dataRaw = base + '/' + 'data_raw' dataTmp = base + '/' + 'data_tmp' dataObs = base + '/' + 'data_obs' # test settings testExperiment = 'historical' testMember = range(1,2) testLonRange = [-65,10] testLatRange = [48.,65.] testLevRange = [0.,2000.] testYearRange = [1950,1950] testMemRange = [1,3] testInputFile = dataTmp + '/' + 'noresm1-cmip6_historical_19500101_mem01.micom.hm.1950-01.nc' if model == 'noresm1-me': gridFileOcn = dataTmp + '/' + 'grid_ocn_gx1v6.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_1.9x2.5_gx1v6.nc' elif model == 'cesm1' or model == 'cesm2': gridFileOcn = dataTmp + '/' + 'grid_ocn_gx1v6.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_0.9x1.25_tnx1v4.nc' elif model == 'noresm2-lm' or model == 'noresm2-ml': gridFileOcn = dataTmp + '/' + 'grid_ocn_tnx1v4.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_1.9x2.5_tnx1v4.nc' elif model == 'noresm2-mm': gridFileOcn = dataTmp + '/' + 'grid_ocn_tnx1v4.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_0.9x1.25_tnx1v4.nc' testInputFile = dataTmp + '/' + 'noresm2-mm_oda_20mem_19500101_mem001.blom.hm.1950-01.nc' elif model == 'noresm2-mh': gridFileOcn = dataTmp + '/' + 'grid_ocn_tnx0.25v4.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_0.9x1.25_tnx0.25v4.nc' # other global fields daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) weightsAnnMean = daysInMonth / np.sum(daysInMonth) setGlobals() def globalMsss(fieldMod='sst',exp='HIN1',leadRange=[0,0],memRange=[1,10],res=[5,5],nboot=1000,blocklen=5,force=False): if fieldMod == 'sst': fieldObs = 'sst' obsName = 'ERSSTv5' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 273.15 levRange = [0,0] landFill = True elif fieldMod == 'ppint': fieldObs = 'PP' #obsName = 'GlobColour' obsName = 'GlobColour2' obsCoverage = [1998,2019] obsFactor = 1/(12*1000*24*3600) # mod = mol C m-2 s-1, obs = mg m-2 day-1 obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'ppintPerfect': fieldMod = 'ppint' fieldObs = 'ppint' obsName = 'dcppA-assim-i1' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'pco2': fieldObs = 'spco2' obsName = 'SOCCOM' obsCoverage = [1982,2017] obsFactor = 1. # mod = uatm , muatm obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'fgco2': fieldObs = 'fgco2' obsName = 'SOCCOM' obsCoverage = [1982,2017] obsFactor = -12/(1000*365*24*3600) # mod = kg C m-2 s-1 , obs = mol/m2/yr obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'templvl': fieldObs = 'temperature' obsName = 'EN4' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,300] landFill = True elif fieldMod == 'salnlvl': fieldObs = 'salinity' obsName = 'EN4' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,300] landFill = True elif fieldMod == 'sealv': fieldObs = 'zo' obsName = 'ARMOR3D' obsCoverage = [1993,2018] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'TREFHT': fieldObs = 'TREFHT' obsName = 'HadCRUT' obsCoverage = [1950,2019] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'PRECT': fieldObs = 'PRECT' obsName = 'CRUPRE' obsCoverage = [1950,2018] obsFactor = 1/(365/12*24*3600*1000) # mod=m/s obs=mm/month obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'Z500': fieldObs = 'Z500' obsName = 'ERA5' obsCoverage = [1979,2019] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'PSL': fieldObs = 'PSL' obsName = 'NCEP' obsCoverage = [1950,2019] obsFactor = 100. obsOffset = 0. levRange = [0,0] landFill = False # lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) tagRes = '_{:d}x{:d}'.format(res[0],res[1]) tagField = '_' + fieldMod if obsName[0:11] != 'dcppA-assim' else '_' + fieldMod + 'Perfect' tagLead = '_LY{:d}'.format(leadRange[0]+1) if leadRange[0] == leadRange[1] else '_LY{:d}-{:d}'.format(leadRange[0]+1,leadRange[1]+1) tagExp = '_' + exp # extract data if exp[0:3] == 'ANA': modCoverage = [1950,2018] syear1 = np.max((modCoverage[0],obsCoverage[0])) else: modCoverage = [1950,2029] syear1 = np.max((1960,obsCoverage[0])) if not exp == 'PERS' else np.max((1960,obsCoverage[0]+1+leadRange[1]-leadRange[0])) syearn = np.min((modCoverage[1],obsCoverage[1]))-leadRange[1]-1 syears = range(syear1,syearn+1) tagYears = '_s{:d}-{:d}'.format(syears[0],syears[-1]) # filePath = 'data/gmsss10mem' + tagField + tagExp + tagYears + tagLead + tagRes + '.nc' if os.path.isfile(filePath) and not force: print(filePath + ' exists. Read result from file.') nc = Dataset(filePath) msss = nc.msss bmsss = nc.bmsss nc.close() else: print(filePath + ' does not exist. Compute and write result to file.') # obs if obsName == 'GlobColour': obs = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes),axis=1) elif obsName[0:11] == 'dcppA-assim': obs = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10],levRange=levRange,suffix=tagRes) else: obs = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes) obs = (obs + obsOffset)*obsFactor # mod if exp == 'HIST': mod = readHindcastLY(fieldMod,'historical',syears,leadRange,yearRange=[1950,2029],memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'PERS': if obsName[0:11] == 'dcppA-assim': mod = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10], levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) elif obsName == 'GlobColour': mod = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False),axis=1) else: mod = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) mod = (mod + obsOffset)*obsFactor elif exp == 'HIN1': mod = readHindcastLY(fieldMod,'dcppA-hindcast-i1',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'HIN2': mod = readHindcastLY(fieldMod,'dcppA-hindcast-i2',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'ANA1': mod = readHindcastLY(fieldMod,'dcppA-assim-i1',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'ANA2': mod = readHindcastLY(fieldMod,'dcppA-assim-i2',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) # if fieldMod == 'pco2': obs = detrendMulti(obs) mod = detrendMulti(mod) # print('mod ave =',np.nanmean(mod)) print('obs ave =',np.nanmean(obs)) # y1, y2 = mod, obs # add dummy member dimension to x1,x2 if necessary if len(y1.shape) < 4: y1 = y1.reshape((y1.shape[0],1) + y1.shape[1:]) if len(y2.shape) < 4: y2 = y2.reshape((y2.shape[0],1) + y2.shape[1:]) # compute anomalies y1 = y1 - np.nanmean(y1,axis=(0,1),keepdims=True) y2 = y2 - np.nanmean(y2,axis=(0,1),keepdims=True) # compute ensemble average y1a = np.nanmean(y1,axis=1) y2a = np.nanmean(y2,axis=1) # compute msss msss = 1 - np.nansum((y1a-y2a)**2 * np.cos(lat2*np.pi/180)) / np.nansum((y1a*0-y2a)**2 * np.cos(lat2*np.pi/180)) # compute local significance maxlen, nmem = y1.shape[0], y1.shape[1] nblock = np.int(np.ceil(maxlen / blocklen)) np.random.seed(0) indBlockStart = np.trunc(np.random.rand(nboot,nblock)*(maxlen-(blocklen-1))) indBlockMem = np.trunc(np.random.rand(nboot,nmem)*nmem).astype(int) ind = np.zeros(nblock*blocklen,dtype=int) bmsss = np.zeros(nboot) for i in range(nboot): for j in range(nblock): ind[j*blocklen:(j+1)*blocklen] = indBlockStart[i,j] + np.arange(blocklen) bm = y1.take(indBlockMem[i,:],axis=1).mean(axis=1).take(ind[0:maxlen],axis=0) bm = bm - bm.mean(axis=0,keepdims=True) bo = y2a.take(ind[0:maxlen],axis=0) bo = bo - bo.mean(axis=0,keepdims=True) bmsss[i] = 1 - np.nansum((bm - bo)**2 * np.cos(lat2*np.pi/180)) / np.nansum((bm*0 - bo)**2 * np.cos(lat2*np.pi/180)) # write result nc = Dataset(filePath,'w') nc.msss = msss nc.bmsss = bmsss nc.close() print('msss =', msss) return msss, bmsss def globalCorr(fieldMod='sst',exp='HIN1',leadRange=[0,0],memRange=[1,10],res=[5,5],nboot=1000,blocklen=5,force=False,normalize=False): if fieldMod == 'sst': fieldObs = 'sst' obsName = 'ERSSTv5' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 273.15 levRange = [0,0] landFill = True elif fieldMod == 'ppint': fieldObs = 'PP' #obsName = 'GlobColour' obsName = 'GlobColour2' #obsCoverage = [1998,2019] obsCoverage = [1998,2021] obsFactor = 1/(12*1000*24*3600) # mod = mol C m-2 s-1, obs = mg m-2 day-1 obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'ppintPerfect': fieldMod = 'ppint' fieldObs = 'ppint' obsName = 'dcppA-assim-i1' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'pco2': fieldObs = 'spco2' obsName = 'SOCCOM' obsCoverage = [1982,2017] obsFactor = 1. # mod = uatm , muatm obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'fgco2': fieldObs = 'fgco2' obsName = 'SOCCOM' obsCoverage = [1982,2017] obsFactor = -12/(1000*365*24*3600) # mod = kg C m-2 s-1 , obs = mol/m2/yr obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'templvl': fieldObs = 'temperature' obsName = 'EN4' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,300] landFill = True elif fieldMod == 'salnlvl': fieldObs = 'salinity' obsName = 'EN4' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,300] landFill = True elif fieldMod == 'sealv': fieldObs = 'zo' obsName = 'ARMOR3D' obsCoverage = [1993,2018] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'TREFHT': fieldObs = 'TREFHT' obsName = 'HadCRUT' obsCoverage = [1950,2019] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'PRECT': fieldObs = 'PRECT' obsName = 'CRUPRE' obsCoverage = [1950,2018] obsFactor = 1/(365/12*24*3600*1000) # mod=m/s obs=mm/month obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'PSL': fieldObs = 'PSL' obsName = 'NCEP' obsCoverage = [1950,2019] obsFactor = 100. obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'Z500': fieldObs = 'Z500' obsName = 'ERA5' obsCoverage = [1979,2019] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = False # lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) tagRes = '_{:d}x{:d}'.format(res[0],res[1]) tagField = '_' + fieldMod if obsName[0:11] != 'dcppA-assim' else '_' + fieldMod + 'Perfect' tagLead = '_LY{:d}'.format(leadRange[0]+1) if leadRange[0] == leadRange[1] else '_LY{:d}-{:d}'.format(leadRange[0]+1,leadRange[1]+1) tagExp = '_' + exp # extract data if exp[0:3] == 'ANA' and not exp == 'ANAL-MM': modCoverage = [1950,2018] syear1 = np.max((modCoverage[0],obsCoverage[0])) elif exp == 'ANAL-MM': modCoverage = [1998,2020] syear1 = np.max((modCoverage[0],obsCoverage[0])) elif exp == 'NUDGE-MM': modCoverage = [1998,2020] syear1 = np.max((modCoverage[0],obsCoverage[0])) elif exp == 'FOSI-MM': modCoverage = [1990,2018] syear1 = np.max((modCoverage[0],obsCoverage[0])) elif exp == 'HIST-MM': modCoverage = [1998,2022] syear1 = np.max((modCoverage[0],obsCoverage[0])) else: if obsName == 'GlobColour2': modCoverage = [1998,2021] syear1 = np.max((1998,obsCoverage[0])) if not exp == 'PERS' else np.max((1998,obsCoverage[0]+1+leadRange[1]-leadRange[0])) else: modCoverage = [1950,2029] syear1 = np.max((1960,obsCoverage[0])) if not exp == 'PERS' else np.max((1960,obsCoverage[0]+1+leadRange[1]-leadRange[0])) syearn = np.min((modCoverage[1],obsCoverage[1]))-leadRange[1]-1 syears = range(syear1,syearn+1) tagYears = '_s{:d}-{:d}'.format(syears[0],syears[-1]) # if normalize: filePath = 'data/gcorrnorm' + tagField + tagExp + tagYears + tagLead + tagRes + '.nc' else: filePath = 'data/gcorr' + tagField + tagExp + tagYears + tagLead + tagRes + '.nc' if os.path.isfile(filePath) and not force: print(filePath + ' exists. Read result from file.') nc = Dataset(filePath) r = nc.corr br = nc.bcorr nc.close() else: print(filePath + ' does not exist. Compute and write result to file.') # obs if obsName == 'GlobColour': obs = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes),axis=1) elif obsName[0:11] == 'dcppA-assim': obs = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10],levRange=levRange,suffix=tagRes) else: obs = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes) obs = (obs + obsOffset)*obsFactor # mod if exp == 'HIST': mod = readHindcastLY(fieldMod,'historical',syears,leadRange,yearRange=[1950,2029],memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'PERS': if obsName[0:11] == 'dcppA-assim': mod = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10], levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) elif obsName == 'GlobColour': mod = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False),axis=1) else: mod = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) mod = (mod + obsOffset)*obsFactor elif exp == 'HIN1': mod = readHindcastLY(fieldMod,'dcppA-hindcast-i1',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'HIN2': mod = readHindcastLY(fieldMod,'dcppA-hindcast-i2',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'ANA1': mod = readHindcastLY(fieldMod,'dcppA-assim-i1',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'ANA2': mod = readHindcastLY(fieldMod,'dcppA-assim-i2',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'HIST-MM': mod = readHindcastLY(fieldMod,'noresm_ctl_f09_tn14_19700101',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'ANAL-MM': mod = readHindcastLY(fieldMod,'norcpm_ana_f09_tn14',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'NUDGE-MM': mod = readHindcastLY(fieldMod,'NHISTfrc2_f09_tn14_nudgeEra5AnomUVPS6h_s19800101_1',syears,leadRange,yearRange=modCoverage,memRange=[0,0],levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'FOSI-MM': mod = readHindcastLY(fieldMod,'NOIIAJRAOC20TR_TL319_tn14_20190710',syears,leadRange,yearRange=modCoverage,memRange=[0,0],levRange=levRange,suffix=tagRes,ensave=False) # if fieldMod == 'pco2': obs = detrendMulti(obs) mod = detrendMulti(mod) # print('mod ave =',np.nanmean(mod)) print('obs ave =',np.nanmean(obs)) # y1, y2 = mod, obs # add dummy member dimension to x1,x2 if necessary if len(y1.shape) < 4: y1 = y1.reshape((y1.shape[0],1) + y1.shape[1:]) if len(y2.shape) < 4: y2 = y2.reshape((y2.shape[0],1) + y2.shape[1:]) # compute anomalies y1 = y1 - np.nanmean(y1,axis=(0,1),keepdims=True) y2 = y2 - np.nanmean(y2,axis=(0,1),keepdims=True) # normalize if normalize: y1 = y1 / np.nanstd(y1,axis=(0,1),keepdims=True) y2 = y2 / np.nanstd(y2,axis=(0,1),keepdims=True) # compute ensemble average y1a = np.nanmean(y1,axis=1) y2a = np.nanmean(y2,axis=1) # compute correlation r = np.nanmean(y1a*y2a * np.cos(lat2*np.pi/180)) / np.sqrt(np.nanmean((0*y1a**2+y2a**2) * np.cos(lat2*np.pi/180))) / np.sqrt(np.nanmean((y1a**2+0*y2a**2) * np.cos(lat2*np.pi/180))) # compute local significance maxlen, nmem = y1.shape[0], y1.shape[1] nblock = np.int(np.ceil(maxlen / blocklen)) np.random.seed(0) indBlockStart = np.trunc(np.random.rand(nboot,nblock)*(maxlen-(blocklen-1))) indBlockMem = np.trunc(np.random.rand(nboot,nmem)*nmem).astype(int) ind = np.zeros(nblock*blocklen,dtype=int) #print(y1.shape) #print(y2.shape) br = np.zeros(nboot) for i in range(nboot): for j in range(nblock): ind[j*blocklen:(j+1)*blocklen] = indBlockStart[i,j] + np.arange(blocklen) bm = y1.take(indBlockMem[i,:],axis=1).mean(axis=1).take(ind[0:maxlen],axis=0) bm = bm - bm.mean(axis=0,keepdims=True) bo = y2a.take(ind[0:maxlen],axis=0) bo = bo - bo.mean(axis=0,keepdims=True) br[i] = np.nanmean(bm*bo * np.cos(lat2*np.pi/180)) / np.sqrt(np.nanmean((0*bm**2 + bo**2) * np.cos(lat2*np.pi/180))) / np.sqrt(np.nanmean((bm**2 + 0*bo**2) * np.cos(lat2*np.pi/180))) # write result nc = Dataset(filePath,'w') nc.corr = r nc.bcorr = br nc.close() print('r =', r) return r, br def globalMeanCorr(fieldMod='sst',exp='HIN1',leadRange=[0,0],memRange=[1,10],res=[5,5],nboot=1000,blocklen=5,force=False,normalize=False): if fieldMod == 'sst': fieldObs = 'sst' obsName = 'ERSSTv5' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 273.15 levRange = [0,0] landFill = True elif fieldMod == 'ppint': fieldObs = 'PP' #obsName = 'GlobColour' obsName = 'GlobColour2' #obsCoverage = [1998,2019] obsCoverage = [1998,2021] obsFactor = 1/(12*1000*24*3600) # mod = mol C m-2 s-1, obs = mg m-2 day-1 obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'ppintPerfect': fieldMod = 'ppint' fieldObs = 'ppint' obsName = 'dcppA-assim-i1' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'pco2': fieldObs = 'spco2' obsName = 'SOCCOM' obsCoverage = [1982,2017] obsFactor = 1. # mod = uatm , muatm obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'fgco2': fieldObs = 'fgco2' obsName = 'SOCCOM' obsCoverage = [1982,2017] obsFactor = -12/(1000*365*24*3600) # mod = kg C m-2 s-1 , obs = mol/m2/yr obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'templvl': fieldObs = 'temperature' obsName = 'EN4' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,300] landFill = True elif fieldMod == 'salnlvl': fieldObs = 'salinity' obsName = 'EN4' obsCoverage = [1950,2018] obsFactor = 1. obsOffset = 0. levRange = [0,300] landFill = True elif fieldMod == 'sealv': fieldObs = 'zo' obsName = 'ARMOR3D' obsCoverage = [1993,2018] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = True elif fieldMod == 'TREFHT': fieldObs = 'TREFHT' obsName = 'HadCRUT' obsCoverage = [1950,2019] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'PRECT': fieldObs = 'PRECT' obsName = 'CRUPRE' obsCoverage = [1950,2018] obsFactor = 1/(365/12*24*3600*1000) # mod=m/s obs=mm/month obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'PSL': fieldObs = 'PSL' obsName = 'NCEP' obsCoverage = [1950,2019] obsFactor = 100. obsOffset = 0. levRange = [0,0] landFill = False elif fieldMod == 'Z500': fieldObs = 'Z500' obsName = 'ERA5' obsCoverage = [1979,2019] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = False # lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) tagRes = '_{:d}x{:d}'.format(res[0],res[1]) tagField = '_' + fieldMod if obsName[0:11] != 'dcppA-assim' else '_' + fieldMod + 'Perfect' tagLead = '_LY{:d}'.format(leadRange[0]+1) if leadRange[0] == leadRange[1] else '_LY{:d}-{:d}'.format(leadRange[0]+1,leadRange[1]+1) tagExp = '_' + exp # extract data if exp[0:3] == 'ANA' and not exp == 'ANAL-MM': modCoverage = [1950,2018] syear1 = np.max((modCoverage[0],obsCoverage[0])) elif exp == 'ANAL-MM': modCoverage = [1998,2020] syear1 = np.max((modCoverage[0],obsCoverage[0])) elif exp == 'NUDGE-MM': modCoverage = [1998,2020] syear1 = np.max((modCoverage[0],obsCoverage[0])) elif exp == 'FOSI-MM': modCoverage = [1958,2018] syear1 = np.max((modCoverage[0],obsCoverage[0])) elif exp == 'HIST-MM': modCoverage = [1998,2022] syear1 = np.max((modCoverage[0],obsCoverage[0])) else: modCoverage = [1950,2029] syear1 = np.max((1960,obsCoverage[0])) if not exp == 'PERS' else np.max((1960,obsCoverage[0]+1+leadRange[1]-leadRange[0])) syearn = np.min((modCoverage[1],obsCoverage[1]))-leadRange[1]-1 syears = range(syear1,syearn+1) tagYears = '_s{:d}-{:d}'.format(syears[0],syears[-1]) # if normalize: filePath = 'data/gmcorrnorm10mem' + tagField + tagExp + tagYears + tagLead + tagRes + '.nc' else: filePath = 'data/gmcorr10mem' + tagField + tagExp + tagYears + tagLead + tagRes + '.nc' if os.path.isfile(filePath) and not force: print(filePath + ' exists. Read result from file.') nc = Dataset(filePath) r = nc.corr br = nc.bcorr nc.close() else: print(filePath + ' does not exist. Compute and write result to file.') # obs if obsName == 'GlobColour': obs = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes),axis=1) elif obsName[0:11] == 'dcppA-assim': obs = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10],levRange=levRange,suffix=tagRes) else: obs = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes) obs = (obs + obsOffset)*obsFactor # mod if exp == 'HIST': mod = readHindcastLY(fieldMod,'historical',syears,leadRange,yearRange=[1950,2029],memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'PERS': if obsName[0:11] == 'dcppA-assim': mod = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10], levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) elif obsName == 'GlobColour': mod = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False),axis=1) else: mod = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) mod = (mod + obsOffset)*obsFactor elif exp == 'HIN1': mod = readHindcastLY(fieldMod,'dcppA-hindcast-i1',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'HIN2': mod = readHindcastLY(fieldMod,'dcppA-hindcast-i2',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'ANA1': mod = readHindcastLY(fieldMod,'dcppA-assim-i1',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'ANA2': mod = readHindcastLY(fieldMod,'dcppA-assim-i2',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'HIST-MM': mod = readHindcastLY(fieldMod,'noresm_ctl_f09_tn14_19700101',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'ANAL-MM': mod = readHindcastLY(fieldMod,'norcpm_ana_f09_tn14',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'NUDGE-MM': mod = readHindcastLY(fieldMod,'NHISTfrc2_f09_tn14_nudgeEra5AnomUVPS6h_s19800101_1',syears,leadRange,yearRange=modCoverage,memRange=[0,0],levRange=levRange,suffix=tagRes,ensave=False) elif exp == 'FOSI-MM': mod = readHindcastLY(fieldMod,'NOIIAOC20TR_T62_tn14_20190628',syears,leadRange,yearRange=modCoverage,memRange=[0,0],levRange=levRange,suffix=tagRes,ensave=False) # if fieldMod == 'pco2': obs = detrendMulti(obs) mod = detrendMulti(mod) # print('mod ave =',np.nanmean(mod)) print('obs ave =',np.nanmean(obs)) # y1, y2 = mod, obs # add dummy member dimension to x1,x2 if necessary if len(y1.shape) < 4: y1 = y1.reshape((y1.shape[0],1) + y1.shape[1:]) if len(y2.shape) < 4: y2 = y2.reshape((y2.shape[0],1) + y2.shape[1:]) # compute anomalies y1 = y1 - np.nanmean(y1,axis=(0,1),keepdims=True) y2 = y2 - np.nanmean(y2,axis=(0,1),keepdims=True) # normalize if normalize: y1 = y1 / np.nanstd(y1,axis=(0,1),keepdims=True) y2 = y2 / np.nanstd(y2,axis=(0,1),keepdims=True) # compute global average and ensemble average print(((y1+0*y2)*np.cos(lat2*np.pi/180)).shape) print(((1+0*y1+0*y2)*np.cos(lat2*np.pi/180)).shape) y1gm = np.nanmean((y1+0*y2)*np.cos(lat2*np.pi/180),axis=(2,3))/np.nanmean((1+0*y1+0*y2)*np.cos(lat2*np.pi/180),axis=(2,3)) y2gm = np.nanmean((0*y1+y2)*np.cos(lat2*np.pi/180),axis=(2,3))/np.nanmean((1+0*y1+0*y2)*np.cos(lat2*np.pi/180),axis=(2,3)) y1=y1gm y2=y2gm y1a = np.nanmean(y1,axis=1) y2a = np.nanmean(y2,axis=1) # compute correlation r = np.nanmean(y1gm*y2gm) / np.sqrt(np.nanmean(y1gm**2)*np.nanmean(y2gm**2)) # compute significance maxlen, nmem = y1.shape[0], y1.shape[1] nblock = np.int(np.ceil(maxlen / blocklen)) np.random.seed(0) indBlockStart = np.trunc(np.random.rand(nboot,nblock)*(maxlen-(blocklen-1))) indBlockMem = np.trunc(np.random.rand(nboot,nmem)*nmem).astype(int) ind = np.zeros(nblock*blocklen,dtype=int) br = np.zeros(nboot) for i in range(nboot): for j in range(nblock): ind[j*blocklen:(j+1)*blocklen] = indBlockStart[i,j] + np.arange(blocklen) bm = y1.take(indBlockMem[i,:],axis=1).mean(axis=1).take(ind[0:maxlen],axis=0) bm = bm - bm.mean(axis=0,keepdims=True) bo = y2a.take(ind[0:maxlen],axis=0) bo = bo - bo.mean(axis=0,keepdims=True) br[i] = np.nanmean(bm*bo) / np.sqrt(np.nanmean(bo**2)) / np.sqrt(np.nanmean(bm**2)) # write result nc = Dataset(filePath,'w') nc.corr = r nc.bcorr = br nc.close() print('r =', r) return r, br def plotYeagerMaps(inData=['sst',['analysis','hindcast'],[[5,5]],True,True,True]): print(inData) fieldMod = inData[0] products = inData[1] resOptions = inData[2] doACC = inData[3] doMSSS = inData[4] doLabels = inData[5] polar = '' if fieldMod == 'sst': fieldObs = 'sst' obsName = 'ERSSTv5' obsCoverage = [1950,2018] obsFactor = 1. levRange = [0,0] landFill = True if fieldMod in ['aiceNH','aiceSH']: polar = 'NH' if fieldMod == 'aiceNH' else 'SH' fieldMod = 'aice' fieldObs = 'aice' obsName = 'HadISST' obsCoverage = [1950,2019] obsFactor = 1. levRange = [0,0] landFill = False elif fieldMod == 'ppint': fieldObs = 'PP' #obsName = 'GlobColour' obsName = 'GlobColour2' #obsCoverage = [1998,2019] obsCoverage = [1998,2021] obsFactor = 1/(12*1000*24*3600) # mod = mol C m-2 s-1, obs = mg m-2 day-1 levRange = [0,0] landFill = True elif fieldMod == 'ppintPerfect': fieldMod = 'ppint' fieldObs = 'ppint' obsName = 'dcppA-assim-i1' obsCoverage = [1950,2018] obsFactor = 1. levRange = [0,0] landFill = True elif fieldMod == 'pco2': fieldObs = 'spco2' obsName = 'SOCCOM' obsCoverage = [1982,2017] obsFactor = 1. # mod = uatm , muatm levRange = [0,0] landFill = True elif fieldMod == 'fgco2': fieldObs = 'fgco2' obsName = 'SOCCOM' obsCoverage = [1982,2017] obsFactor = 12/(1000*365*24*3600) # mod = kg C m-2 s-1 , obs = mol/m2/yr levRange = [0,0] landFill = True elif fieldMod == 'templvl': fieldObs = 'temperature' obsName = 'EN4' obsCoverage = [1950,2018] obsFactor = 1. levRange = [0,300] landFill = True elif fieldMod == 'salnlvl': fieldObs = 'salinity' obsName = 'EN4' obsCoverage = [1950,2018] obsFactor = 1. levRange = [0,300] landFill = True elif fieldMod == 'sealv': fieldObs = 'zo' obsName = 'ARMOR3D' obsCoverage = [1993,2018] obsFactor = 1. levRange = [0,0] landFill = True elif fieldMod == 'TREFHT': fieldObs = 'TREFHT' obsName = 'HadCRUT' obsCoverage = [1950,2019] obsFactor = 1. levRange = [0,0] landFill = False elif fieldMod == 'PRECT': fieldObs = 'PRECT' obsName = 'CRUPRE' obsCoverage = [1950,2018] obsFactor = 1/(365/12*24*3600*1000) # mod=m/s obs=mm/month levRange = [0,0] landFill = False elif fieldMod == 'PSL': fieldObs = 'PSL' obsName = 'NCEP' obsCoverage = [1950,2019] obsFactor = 100. levRange = [0,0] landFill = False elif fieldMod == 'Z500': fieldObs = 'Z500' obsName = 'ERA5' obsCoverage = [1950,2019] obsFactor = 1. obsOffset = 0. levRange = [0,0] landFill = False for product in products: if product == 'analysis': leadRanges = [[-1,-1]] expOptions = [['ANA1',''],['ANA1','HIST'],['ANA2','ANA1']] memRange = [1,30] else: leadRanges = [[0,0],[1,4],[5,8]] expOptions = [['HIN1','PERS'],['HIN1',''],['HIN1','ANA1'],['HIN2','HIN1'], ['HIN1','HIST']] memRange = [1,10] for leadRange in leadRanges: for res in resOptions: lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) mskfdr = np.where(lat2 < 80, 1, 0) tagRes = '_{:d}x{:d}'.format(res[0],res[1]) tagField = '_' + fieldMod if obsName[0:11] != 'dcppA-assim' else '_' + fieldMod + 'Perfect' tagLead = '_LY{:d}'.format(leadRange[0]+1) if leadRange[0] == leadRange[1] else '_LY{:d}-{:d}'.format(leadRange[0]+1,leadRange[1]+1) # extract data for expOption in expOptions: if expOption[0][0:3] == 'ANA' or expOption[1][0:3] == 'ANA': modCoverage = [1950,2018] else: modCoverage = [1950,2029] if product == 'analysis': syear1 = np.max((modCoverage[0],obsCoverage[0])) else: syear1 = np.max((1960,obsCoverage[0])) if not expOption[1] == 'PERS' else np.max((1960,obsCoverage[0]+1+leadRange[1]-leadRange[0])) syearn = np.min((modCoverage[1],obsCoverage[1]))-leadRange[1]-1 syears = range(syear1,syearn+1) tagYears = '_s{:d}-{:d}'.format(syears[0],syears[-1]) if obsName == 'GlobColour': obs = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes),axis=1) elif obsName[0:11] == 'dcppA-assim': obs = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10],levRange=levRange,suffix=tagRes) else: obs = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes) if fieldObs == 'fgco2': obs = -obs if expOption[0] == 'HIST': fld1 = readHindcastLY(fieldMod,'historical',syears,leadRange,yearRange=[1950,2029],memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif expOption[0] == 'PERS': if obsName[0:11] == 'dcppA-assim': fld1 = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10], levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) elif obsName == 'GlobColour': fld1 = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False),axis=1) else: fld1 = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) elif expOption[0] == 'HIN1': fld1 = readHindcastLY(fieldMod,'dcppA-hindcast-i1',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif expOption[0] == 'HIN2': fld1 = readHindcastLY(fieldMod,'dcppA-hindcast-i2',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif expOption[0] == 'ANA1': fld1 = readHindcastLY(fieldMod,'dcppA-assim-i1',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif expOption[0] == 'ANA2': fld1 = readHindcastLY(fieldMod,'dcppA-assim-i2',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) # if expOption[1] == 'HIST': fld2 = readHindcastLY(fieldMod,'historical',syears,leadRange,yearRange=[1950,2029],memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif expOption[1] == 'PERS': if obsName[0:11] == 'dcppA-assim': fld2 = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10], levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) elif obsName == 'GlobColour': fld2 = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False),axis=1) else: fld2 = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage, levRange=levRange,suffix=tagRes,persistence='mean',ensave=False) elif expOption[1] == 'HIN1': fld2 = readHindcastLY(fieldMod,'dcppA-hindcast-i1',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif expOption[1] == 'HIN2': fld2 = readHindcastLY(fieldMod,'dcppA-hindcast-i2',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif expOption[1] == 'ANA1': fld2 = readHindcastLY(fieldMod,'dcppA-assim-i1',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) elif expOption[1] == 'ANA2': fld2 = readHindcastLY(fieldMod,'dcppA-assim-i2',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False) # if fieldMod == 'pco2': obs = detrendMulti(obs) fld1 = detrendMulti(fld1) if expOption[1] != '': fld2 = detrendMulti(fld2) # print('fld1 ave =',np.nanmean(fld1)) print('obs ave =',np.nanmean(obs)) tagExp = '_' + expOption[0] if expOption[1] == '' else '_{:s}-{:s}'.format(expOption[0],expOption[1]) if doLabels: titleString = expOption[0] if expOption[1] == '' else expOption[0] + ' - ' + expOption[1] title2 = tagLead[1:] if product == 'hindcast' else '' else: titleString = ' ' title2 = ' ' if doACC: filePrefix = 'ACC10mem' + tagField + tagExp + tagYears + tagLead + tagRes print(filePrefix) fld = corrMultiArrayYeager(fld1,obs*obsFactor) if expOption[1] == '' else corrMultiArrayDiffYeager(fld1,obs,fld2) plotACC(lon=lon,lat=lat,fld=fld,filePrefix=filePrefix,lbLabelBarOn=False, title=titleString,title2=title2,landFill=landFill,plottype='ACC',polar=polar) if doMSSS: filePrefix = 'MSSS10mem' + tagField + tagExp + tagYears + tagLead + tagRes print(filePrefix) fld = MSSSyeager(fld1,obs*obsFactor) if expOption[1] == '' else MSSSyeager(fld1,obs,fld2) plotACC(lon=lon,lat=lat,fld=fld,filePrefix=filePrefix,lbLabelBarOn=False, title=titleString,title2=title2,landFill=landFill,plottype='MSSS',polar=polar) def zostoga(temp,saln): return def boxMean3d(field,experiment,yearRange,lonRange=[0,360],latRange=[-90,90],levRange=[0,8000],memRange=[0,0],syear=0): """ Compute global or box mean of 3d field. """ if lonRange == [0,360] and latRange == [-90,90]: suffix = '_gm' globalDomain = True else: suffix = '_' + boxTag(lonRange,latRange) globalDomain = False filePathGlob = outPath(field,experiment,yearRange,memRange=memRange,levRange=levRange,syear=syear,suffix=suffix) if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read global mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate global mean.') compWeightsOcn3d(lonRange,latRange,levRange) if memRange[1] == 0: gm = np.zeros([yearRange[1]-yearRange[0]+1]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): ncin = Dataset(filePathMon(field,experiment,year,month+1,mem,syear=syear)) glob = np.nansum(ncin.variables[field][0,:]*weightsOcn3d) ncin = ncin.close() if memRange[1] == 0: gm[year-yearRange[0]] = gm[year-yearRange[0]] + glob*daysInMonth[month]/365. else: gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + glob*daysInMonth[month]/365. nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('year',yearRange[1]-yearRange[0]+1) if memRange[1] == 0: ncgm = nc.createVariable(field,'f4',['year']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) ncgm = nc.createVariable(field,'f4',['year','member']) print(gm.shape) ncgm[:] = gm nc = nc.close() return gm def boxMean2d(field,experiment,yearRange,memRange=[0,0],month1=1,syear=0,lonRange=[0,360],latRange=[-90,90],levRange=[0,0],force=False,annual=False,ensave=False,singleInputFile=False,debug=False): """ Compute global mean of 2d field. """ if debug: print("boxMean2d field: ", field) print("boxMean2d experiment: ", experiment) print("boxMean2d yearRange: ", yearRange) print("boxMean2d memRange: ", memRange) print("boxMean2d lonRange: ", lonRange) print("boxMean2d latRange: ", latRange) print("boxMean2d levRange: ", levRange) print("boxMean2d force: ", force) print("boxMean2d annual: ", annual) print("boxMean2d singleInputFile: ", singleInputFile) filePathGlob = outPath(field,experiment,yearRange,memRange=memRange,levRange=levRange,month1=month1,syear=syear,annual=annual,ensave=ensave,suffix='_' + boxTag(lonRange,latRange)) print(filePathGlob) if os.path.isfile(filePathGlob) and not(force): print(filePathGlob + ' exists. Read box mean from file.') nc = Dataset(filePathGlob) bm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate box mean.') filePath = outPath(field,experiment,yearRange,memRange=memRange,levRange=levRange,month1=month1,syear=syear,annual=False,ensave=ensave) print(filePath) compWeightsOcn2d(lonRange,latRange) if experiment == 'EN4': nc = Dataset('data_tmp/coords_EN4.nc') lon_EN4 = nc.variables['lon'][:] lat_EN4 = nc.variables['lat'][:] lon2_EN4, lat2_EN4 = np.meshgrid(lon_EN4,lat_EN4) area_EN4 = np.cos(lat2_EN4*np.pi/180) print(latRange,lonRange) for j in range(area_EN4.shape[0]): for i in range(area_EN4.shape[1]): if not (lon2_EN4[j,i] > lonRange[0]+360 and lon2_EN4[j,i] < lonRange[1]+360 and lat2_EN4[j,i] > latRange[0] and lat2_EN4[j,i] < latRange[1]): area_EN4[j,i] = 0. nc.close() nmonth = 1 if annual else 12 if os.path.isfile(filePath) and False: ncin = Dataset(filePath) bm = np.sum(ncin.variables[field][:,:,0,:,:]*areaNorm,axis=(-1,-2)) ncin = ncin.close() else: if memRange[1] == 0 or ensave: bm = np.zeros((yearRange[1]-yearRange[0]+1)*nmonth-month1+1) else: bm = np.zeros([(yearRange[1]-yearRange[0]+1)*nmonth-month1+1,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): monthLast = 1 if annual else 12 monthFirst = month1 if year == yearRange[0] else 1 for month in range(monthFirst,monthLast+1): #print(filePath) if not(os.path.isfile(filePath)): if singleInputFile: fpath = outPath(field,experiment,yearRange,levRange=levRange,memRange=memRange,annual=False,ensave=ensave) else: print(month) fpath = filePathMon(field,experiment,year,month,mem,syear=syear) print(fpath) ncin = Dataset(fpath) if singleInputFile: rec = (year - yearRange[0])*12 + month - 1 fld2d = ncin.variables[field][rec,:] if mem == 0: if experiment == 'EN4': boxAve = np.nansum(fld2d*area_EN4,axis=(0,1))/np.nansum(fld2d*0+area_EN4,axis=(0,1)) else: boxAve = np.nansum(fld2d*weightsOcn2d,axis=(0,1)) else: if experiment == 'EN4': boxAve = np.nansum(fld2d*area_EN4,axis=(0,1))/np.nansum(fld2d*0+area_EN4,axis=(0,1)) else: boxAve = np.nansum(fld2d*weightsOcn2d,axis=(0,1)) #boxAve = np.nansum(ncin.variables[field][mem-1,rec,:]*weightsOcn2d,axis=(0,1)) else: fld2d = ncin.variables[field][0,:] if experiment == 'EN4': boxAve = np.nansum(fld2d*area_EN4,axis=(0,1))/np.nansum(fld2d*0+area_EN4,axis=(0,1)) else: boxAve = np.nansum(fld2d*weightsOcn2d,axis=(0,1)) #boxAve = np.nansum(ncin.variables[field][0,:]*weightsOcn2d,axis=(0,1)) ncin = ncin.close() else: ncin = Dataset(filePath) record = (year-yearRange[0])*12+month-1-month1+1 #print('rec=',record,', mem=',mem) if mem == 0: fld2d = ncin.variables[field][record,:] #boxAve = np.nansum(ncin.variables[field][record,:]*weightsOcn2d,axis=(-1,-2)) else: fld2d = ncin.variables[field][record,mem-1,:] if experiment == 'EN4': boxAve = np.nansum(fld2d*area_EN4,axis=(-1,-2))/np.nansum(fld2d*0+area_EN4,axis=(-1,-2)) else: boxAve = np.nansum(fld2d*weightsOcn2d,axis=(-1,-2)) #boxAve = np.nansum(ncin.variables[field][record,mem-1,:]*weightsOcn2d,axis=(-1,-2)) ncin = ncin.close() if annual: if memRange[1] == 0: bm[year-yearRange[0]] = bm[year-yearRange[0]] +boxAve*daysInMonth[month]/365. else: bm[year-yearRange[0],mem-memRange[0]] = bm[year-yearRange[0],mem-memRange[0]] + boxAve*daysInMonth[month]/365. else: if memRange[1] == 0: bm[(year-yearRange[0])*12+month-month1] = boxAve else: bm[(year-yearRange[0])*12+month-month1,mem-memRange[0]] = boxAve nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('record',(yearRange[1]-yearRange[0]+1)*nmonth-month1+1) if memRange[1] == 0: ncbm = nc.createVariable(field,'f4',['record']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) ncbm = nc.createVariable(field,'f4',['record','member']) ncbm[:] = bm nc = nc.close() return bm def globalMean2d(field,experiment,yearRange,memRange=[0,0],syear=0,lonRange=[0,360],latRange=[-90,90],annual=True): """ Compute global mean of 2d field. """ if lonRange == [0,360] and latRange == [-90,90]: suffix = '_gm' globalDomain = True else: suffix = '_' + boxTag(lonRange,latRange) globalDomain = False filePathGlob = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear,suffix=suffix,annual=annual) if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read global mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate global mean.') if getComponentFromField(field) == 'atm': if not 'areaAtmNorm' in globals(): readGridAtm() areaNorm = areaAtmNorm lon, lat = lonAtm, latAtm lon, lat = np.meshgrid(lon,lat) else: if not 'areaOcnNorm' in globals(): readGridOcn() areaNorm = areaOcnNorm lon, lat = lonOcn, latOcn if experiment in ['ERA5']: ncin = Dataset(outPath(field,experiment,yearRange)) lon, lat = ncin.variables['lon'][:], ncin.variables['lat'][:] lon, lat = np.meshgrid(lon,lat) area = np.cos(lat*np.pi/180.) areaNorm = area/np.sum(area) ncin.close() if not(globalDomain): for j in range(areaNorm.shape[0]): for i in range(areaNorm.shape[1]): if lonRange[0] < lonRange[1]: # box lies within [-180 180] if not (lon[j,i] > lonRange[0] and lon[j,i] < lonRange[1] and lat[j,i] > latRange[0] and lat[j,i] < latRange[1]): areaNorm[j,i] = 0. else: # box crosses 180 if not ((lon[j,i] > lonRange[0] or lon[j,i] < lonRange[1]) and lat[j,i] > latRange[0] and lat[j,i] < latRange[1]): areaNorm[j,i] = 0. areaNorm = areaNorm / np.nansum(areaNorm) filePath = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear,ensave=False,annual=annual) print(filePath) if os.path.isfile(filePath): ncin = Dataset(filePath) dims = ncin.variables[field].dimensions if experiment in ['ERA5']: dims = ncin.variables[field].dimensions nrec = len(ncin.dimensions[dims[0]]) print(nrec) gm = np.zeros((nrec)) for r in range(nrec): gm[r] = np.nansum(ncin.variables[field][r,:]*areaNorm,axis=(-1,-2)) else: dims = ncin.variables[field].dimensions nmem = len(ncin.dimensions[dims[0]]) nrec = len(ncin.dimensions[dims[1]]) print(nmem,nrec) gm = np.zeros((nrec,nmem)) for m in range(nmem): for r in range(nrec): gm[r,m] = np.nansum(ncin.variables[field][m,r,:]*areaNorm,axis=(-1,-2)) ncin = ncin.close() else: if memRange[1] == 0: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1]) else: gm = np.zeros((yearRange[1]-yearRange[0]+1)*12) else: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) else: gm = np.zeros([(yearRange[1]-yearRange[0]+1)*12,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): print(filePathMon(field,experiment,year,month+1,mem,syear=syear)) ncin = Dataset(filePathMon(field,experiment,year,month+1,mem,syear=syear)) data = ncin.variables[field][:] glob = np.nansum(data*areaNorm,axis=(1,2)) ncin = ncin.close() if memRange[1] == 0: if annual: gm[year-yearRange[0]] = gm[year-yearRange[0]] + glob*daysInMonth[month]/365. else: gm[(year-yearRange[0])*12+month] = glob else: if annual: gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + glob*daysInMonth[month]/365. else: gm[(year-yearRange[0])*12+month,mem-memRange[0]] = glob nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') if annual: nc.createDimension('record',yearRange[1]-yearRange[0]+1) else: nc.createDimension('record',(yearRange[1]-yearRange[0]+1)*12) if memRange[1] == 0: ncgm = nc.createVariable(field,'f4',['record']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) ncgm = nc.createVariable(field,'f4',['record','member']) print(gm.shape) ncgm[:] = gm nc = nc.close() return gm def boxHADISST(yearRange,memRange=[0,0],lonRange=[-180,180],latRange=[-90,90],syear=0,annual=True): field = 'sst' experiment = 'HADISST' res=[1.,1.] lon = np.arange(res[0]/2,360,res[0])-180 lat = -np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) for j in range(len(lat)): area[j,:] = area[j,:] if lat[j] > latRange[0] and lat[j] < latRange[1] else 0. if lonRange[0] < lonRange[1]: # box lies within [-180 180] for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] and lon[i] < lonRange[1] else 0. else: for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] or lon[i] < lonRange[1] else 0. areasum_before=np.nansum(area) filePathGlob = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=annual) print(filePathGlob) if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read global mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate .') filePath = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=annual) if os.path.isfile(filePath): ncin = Dataset(filePath) gm = ncin.variables[field][:] ncin = ncin.close() else: if memRange[1] == 0: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,12]) else: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,12,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): fld2d=readHADISST_sst(year,month+1,readCoords=False) area = area + fld2d*0 print(areasum_before,np.nansum(area)) areaNorm = area / np.nansum(area) print(areaNorm) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): fld2d=readHADISST_sst(year,month+1,readCoords=False) if memRange[1] == 0: if annual: gm[year-yearRange[0]] = gm[year-yearRange[0]] + np.nansum(fld2d*areaNorm,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],month] = gm[year-yearRange[0],month] + np.nansum(fld2d*areaNorm,axis=(0,1)) else: if annual: gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + np.nansum(fld2d*area,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],month,mem-memRange[0]] = gm[year-yearRange[0],month,mem-memRange[0]] + np.nansum(fld2d*area,axis=(0,1)) nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('year',yearRange[1]-yearRange[0]+1) if not(annual): nc.createDimension('month',12) if memRange[1] == 0: if annual: ncgm = nc.createVariable(field,'f4',['year']) else: ncgm = nc.createVariable(field,'f4',['year','month']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) if annual: ncgm = nc.createVariable(field,'f4',['year','member']) else: ncgm = nc.createVariable(field,'f4',['year','month','member']) ncgm[:] = gm nc = nc.close() return gm def boxOISST(yearRange,memRange=[0,0],lonRange=[0,360],latRange=[-90,90],syear=0,annual=True,daily=False,hires=True): daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) field = 'sst' experiment = 'OISST' if hires: res=[1./4.,1./4.] lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) else: res=[1.,1.] lon = np.arange(res[0]/2,360,res[0]) lat = np.flipud(np.arange(-90+res[1]/2,90,res[1])) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) for j in range(len(lat)): area[j,:] = area[j,:] if lat[j] > latRange[0] and lat[j] < latRange[1] else 0. if lonRange[0] < lonRange[1]: # box lies within [-180 180] for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] and lon[i] < lonRange[1] else 0. else: for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] or lon[i] < lonRange[1] else 0. areasum_before=np.nansum(area) if daily: filePathGlob = outPathDaily(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear) else: filePathGlob = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=annual) print(filePathGlob) if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read global mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate .') if daily: filePath = outPathDaily(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear) else: filePath = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=annual) if os.path.isfile(filePath): ncin = Dataset(filePath) gm = ncin.variables[field][:] ncin = ncin.close() else: if memRange[1] == 0: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1]) else: if daily: gm = np.zeros([yearRange[1]-yearRange[0]+1,365]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,12]) else: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) else: if daily: gm = np.zeros([yearRange[1]-yearRange[0]+1,12,memRange[1]-memRange[0]+1]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,365,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): if daily: for day in range(daysInMonth[month]): fld2d,lonOISST,latOISST=readOISST(field,year,month+1,day+1,readCoords=True,hires=hires) area = area + fld2d*0 else: fld2d,lonOISST,latOISST=readMonOISST(field,year,month+1,readCoords=True,hires=hires) area = area + fld2d*0 print(areasum_before,np.nansum(area)) areaNorm = area / np.nansum(area) print(areaNorm) #print(lon-lonOISST) #print(lat-latOISST) #print(lat) #print(latOISST) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): jday = - 1 for month in range(12): if daily: for day in range(daysInMonth[month]): fld2d=readOISST(field,year,month+1,day+1,readCoords=False,hires=hires) jday = jday + 1 if memRange[1] == 0: gm[year-yearRange[0],jday] = gm[year-yearRange[0],jday] + np.nansum(fld2d*areaNorm,axis=(0,1)) else: gm[year-yearRange[0],jday,mem-memRange[0]] = gm[year-yearRange[0],jday,mem-memRange[0]] + np.nansum(fld2d*area,axis=(0,1)) else: fld2d=readMonOISST(field,year,month+1,readCoords=False,hires=hires) if memRange[1] == 0: if annual: gm[year-yearRange[0]] = gm[year-yearRange[0]] + np.nansum(fld2d*areaNorm,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],month] = gm[year-yearRange[0],month] + np.nansum(fld2d*areaNorm,axis=(0,1)) else: if annual: gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + np.nansum(fld2d*area,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],month,mem-memRange[0]] = gm[year-yearRange[0],month,mem-memRange[0]] + np.nansum(fld2d*area,axis=(0,1)) print(gm.shape) nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('year',yearRange[1]-yearRange[0]+1) if not(annual): if daily: nc.createDimension('jday',365) else: nc.createDimension('month',12) if memRange[1] == 0: if annual: ncgm = nc.createVariable(field,'f4',['year']) else: if daily: ncgm = nc.createVariable(field,'f4',['year','jday']) else: ncgm = nc.createVariable(field,'f4',['year','month']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) if annual: ncgm = nc.createVariable(field,'f4',['year','member']) else: if daily: ncgm = nc.createVariable(field,'f4',['year','jday','member']) else: ncgm = nc.createVariable(field,'f4',['year','month','member']) ncgm[:] = gm nc = nc.close() return gm def boxAVISO(yearRange,memRange=[0,0],lonRange=[0,360],latRange=[-90,90],syear=0,annual=True): field = 'sla' experiment = 'AVISO' res=[1./4.,1./4.] lon = np.arange(res[0]/2,360,res[0]) print(-90+res[1]/2) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) for j in range(len(lat)): area[j,:] = area[j,:] if lat[j] > latRange[0] and lat[j] < latRange[1] else 0. if lonRange[0] < lonRange[1]: # box lies within [-180 180] for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] and lon[i] < lonRange[1] else 0. else: for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] or lon[i] < lonRange[1] else 0. areasum_before=np.nansum(area) filePathGlob = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=annual) print(filePathGlob) if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read global mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate .') filePath = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=annual) if os.path.isfile(filePath): ncin = Dataset(filePath) gm = ncin.variables[field][:] ncin = ncin.close() else: if memRange[1] == 0: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,12]) else: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,12,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): print(fileAVISO(year=year,month=month+1)) ncin = Dataset(fileAVISO(year=year,month=month+1)) fld2d = ncin.variables[field][0,:] area = area + fld2d*0 ncin = ncin.close() print(np.nansum(area)) print(areasum_before,np.nansum(area)) areaNorm = area / np.nansum(area) print(areaNorm) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): ncin = Dataset(fileAVISO(year=year,month=month+1)) fld2d = ncin.variables[field][0,:] ncin = ncin.close() if memRange[1] == 0: if annual: gm[year-yearRange[0]] = gm[year-yearRange[0]] + np.nansum(fld2d*areaNorm,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],month] = gm[year-yearRange[0],month] + np.nansum(fld2d*areaNorm,axis=(0,1)) else: if annual: gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + np.nansum(fld2d*area,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],month,mem-memRange[0]] = gm[year-yearRange[0],month,mem-memRange[0]] + np.nansum(fld2d*area,axis=(0,1)) print(gm.shape) nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('year',yearRange[1]-yearRange[0]+1) if not(annual): nc.createDimension('month',12) if memRange[1] == 0: if annual: ncgm = nc.createVariable(field,'f4',['year']) else: ncgm = nc.createVariable(field,'f4',['year','month']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) if annual: ncgm = nc.createVariable(field,'f4',['year','member']) else: ncgm = nc.createVariable(field,'f4',['year','month','member']) ncgm[:] = gm nc = nc.close() return gm def globalMeanAVISO(yearRange,memRange=[0,0],syear=0,annual=True): """ Compute global mean of observed SSH. """ field = 'sla' experiment = 'AVISO' res=[1./4.,1./4.] lon = np.arange(res[0]/2,360,res[0]) print(-90+res[1]/2) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) areasum_before=np.nansum(area) filePathGlob = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear,suffix='_gm',annual=annual) print(filePathGlob) if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read global mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate global mean.') filePath = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear) if os.path.isfile(filePath): ncin = Dataset(filePath) gm = ncin.variables[field][:] ncin = ncin.close() else: if memRange[1] == 0: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,12]) else: if annual: gm = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,12,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): print(fileAVISO(year=year,month=month+1)) ncin = Dataset(fileAVISO(year=year,month=month+1)) fld2d = ncin.variables[field][0,:] area = area + fld2d*0 ncin = ncin.close() print(np.nansum(area)) print(areasum_before,np.nansum(area)) areaNorm = area / np.nansum(area) print(areaNorm) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): ncin = Dataset(fileAVISO(year=year,month=month+1)) fld2d = ncin.variables[field][0,:] ncin = ncin.close() if memRange[1] == 0: if annual: gm[year-yearRange[0]] = gm[year-yearRange[0]] + np.nansum(fld2d*areaNorm,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],month] = gm[year-yearRange[0],month] + np.nansum(fld2d*areaNorm,axis=(0,1)) else: if annual: gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + np.nansum(fld2d*areaNorm,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],month,mem-memRange[0]] = gm[year-yearRange[0],month,mem-memRange[0]] + np.nansum(fld2d*areaNorm,axis=(0,1)) print(gm.shape) nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('year',yearRange[1]-yearRange[0]+1) if not(annual): nc.createDimension('month',12) if memRange[1] == 0: if annual: ncgm = nc.createVariable(field,'f4',['year']) else: ncgm = nc.createVariable(field,'f4',['year','month']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) if annual: ncgm = nc.createVariable(field,'f4',['year','member']) else: ncgm = nc.createVariable(field,'f4',['year','month','member']) ncgm[:] = gm nc = nc.close() return gm def globalMeanARMOR3D(yearRange,memRange=[0,0],syear=0,hiRes=False): """ Compute global mean of observed SSH. """ field = 'zo' if hiRes: experiment = 'ARMOR3D025' res=[1./4.,1./4.] lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-82-res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) else: experiment = 'ARMOR3D' res=[1.,1.] lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) filePathGlob = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear,suffix='_gm') print(filePathGlob) if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read global mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate global mean.') filePath = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear,suffix='_gm') if os.path.isfile(filePath): ncin = Dataset(filePath) full = ncin.variables[field][:] gm = np.sum(full*areaNorm,axis=(1,2)) ncin = ncin.close() else: if memRange[1] == 0: gm = np.zeros([yearRange[1]-yearRange[0]+1]) full = np.zeros([yearRange[1]-yearRange[0]+1,180,360]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) full = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1],180,360) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): ncin = Dataset(fileARMOR3D(year=year,month=month+1,hiRes=hiRes)) fld2d = ncin.variables[field][0,0,:] if memRange[1] == 0: full[year-yearRange[0],:,:] = fld2d else: full[year-yearRange[0],mem-memRange[0],:,:] = fld2d area = area + fld2d*0 ncin = ncin.close() areaNorm = area / np.nansum(area) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): if memRange[1] == 0: gm[year-yearRange[0]] = gm[year-yearRange[0]] + np.nansum(full[year-yearRange[0],:,:]*areaNorm,axis=(0,1))*daysInMonth[month]/365. else: gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + np.nansum(full[year-yearRange[0],mem-memRange[0],:,:]*area,axis=(0,1))*daysInMonth[month]/365. print(gm.shape) nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('year',yearRange[1]-yearRange[0]+1) if memRange[1] == 0: ncgm = nc.createVariable(field,'f4',['year']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) ncgm = nc.createVariable(field,'f4',['year','member']) ncgm[:] = gm nc = nc.close() return gm def spgARMOR3D(yearRange,memRange=[0,0],syear=0,hiRes=False): """ Compute SPG mean of observed SSH. """ field = 'zo' if hiRes: experiment = 'ARMOR3D025' res=[1./4.,1./4.] lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-82-res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) else: experiment = 'ARMOR3D' res=[1.,1.] lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) lonRange = [360-60,360-15] latRange = [48,65] for j in range(len(lat)): area[j,:] = area[j,:] if lat[j] > latRange[0] and lat[j] < latRange[1] else 0. if lonRange[0] < lonRange[1]: # box lies within [-180 180] for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] and lon[i] < lonRange[1] else 0. else: for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] or lon[i] < lonRange[1] else 0. filePathGlob = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear,suffix='_spg') print(filePathGlob) yearRange if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read SPG mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate SPG mean.') filePath = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear) if os.path.isfile(filePath): ncin = Dataset(filePath) full = ncin.variables[field][:] gm = np.sum(full*areaNorm,axis=(1,2)) ncin = ncin.close() else: if memRange[1] == 0: gm = np.zeros([yearRange[1]-yearRange[0]+1]) full = np.zeros([yearRange[1]-yearRange[0]+1,180,360]) else: gm = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) full = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1],180,360) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): ncin = Dataset(fileARMOR3D(year=year,month=month+1,hiRes=hiRes)) fld2d = ncin.variables[field][0,0,:] if memRange[1] == 0: full[year-yearRange[0],:,:] = fld2d else: full[year-yearRange[0],mem-memRange[0],:,:] = fld2d area = area + fld2d*0 ncin = ncin.close() areaNorm = area / np.nansum(area) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): if memRange[1] == 0: #print(gm.shape) # print(area.shape) #print(np.nansum(full[year-yearRange[0],:,:]*area,axis=(0,1)).shape) gm[year-yearRange[0]] = gm[year-yearRange[0]] + np.nansum(full[year-yearRange[0],:,:]*areaNorm)*daysInMonth[month]/365. else: gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + np.nansum(full[year-yearRange[0],mem-memRange[0],:,:]*areaNorm,axis=(0,1))*daysInMonth[month]/365. print(gm.shape) print(area) nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('year',yearRange[1]-yearRange[0]+1) if memRange[1] == 0: ncgm = nc.createVariable(field,'f4',['year']) else: nc.createDimension('member',memRange[1]-memRange[0]+1) ncgm = nc.createVariable(field,'f4',['year','member']) ncgm[:] = gm nc = nc.close() return gm def boxARMOR3D(yearRange,memRange=[0,0],lonRange=[0,360],latRange=[-90,90],syear=0,hiRes=False,staticMask=False): """ Compute box mean of observed SSH. """ field = 'zo' if hiRes: experiment = 'ARMOR3D025' res=[1./4.,1./4.] lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-82-res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) else: experiment = 'ARMOR3D' res=[1.,1.] lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180) for j in range(len(lat)): area[j,:] = area[j,:] if lat[j] > latRange[0] and lat[j] < latRange[1] else 0. if lonRange[0] < lonRange[1]: # box lies within [-180 180] for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] and lon[i] < lonRange[1] else 0. else: for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] or lon[i] < lonRange[1] else 0. if staticMask: filePathGlob = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=False,suffix='_staticMask') else: filePathGlob = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=False) print(filePathGlob) yearRange if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Read SPG mean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate SPG mean.') filePath = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear) if os.path.isfile(filePath) and False: ncin = Dataset(filePath) full = ncin.variables[field][:] gm = np.sum(full*areaNorm,axis=(1,2)) ncin = ncin.close() else: if memRange[1] == 0: gm = np.zeros((yearRange[1]-yearRange[0]+1)*12) full = np.zeros([yearRange[1]-yearRange[0]+1,12,len(lat),len(lon)]) else: gm = np.zeros([(yearRange[1]-yearRange[0]+1)*12,memRange[1]-memRange[0]+1]) full = np.zeros([yearRange[1]-yearRange[0]+1,12,memRange[1]-memRange[0]+1,len(lat),len(lon)]) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): ncin = Dataset(fileARMOR3D(year=year,month=month+1,hiRes=hiRes)) print(fileARMOR3D(year=year,month=month+1,hiRes=hiRes)) fld2d = ncin.variables[field][0,0,:] if memRange[1] == 0: full[year-yearRange[0],:,:] = fld2d else: full[year-yearRange[0],mem-memRange[0],:,:] = fld2d area = area + fld2d*0 ncin = ncin.close() print(np.nansum(area)) areaNorm = area / np.nansum(area) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for mem in range(memRange[0],memRange[1]+1): for month in range(12): print(fileARMOR3D(year=year,month=month+1,hiRes=hiRes)) ncin = Dataset(fileARMOR3D(year=year,month=month+1,hiRes=hiRes)) fld2d = ncin.variables[field][0,0,:] if memRange[1] == 0: full[year-yearRange[0],month,:,:] = fld2d else: full[year-yearRange[0],month,mem-memRange[0],:,:] = fld2d ncin = ncin.close() for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): if memRange[1] == 0: #print(gm.shape) # print(area.shape) #print(np.nansum(full[year-yearRange[0],:,:]*area,axis=(0,1)).shape) #gm[year-yearRange[0]] = gm[year-yearRange[0]] + np.nansum(full[year-yearRange[0],:,:]*areaNorm)*daysInMonth[month]/365. fld2d = full[year-yearRange[0],month,:,:] if not(staticMask): areaNorm = np.where(fld2d>=10.,0,area) areaNorm = areaNorm / np.nansum(areaNorm) gm[(year-yearRange[0])*12+month] = np.nansum(fld2d*areaNorm) #gm[(year-yearRange[0])*12+month] = np.nansum(fld2d*area)/np.nansum(0*fld2d+area) #print(np.nansum(fld2d*areaNorm)) else: #gm[year-yearRange[0],mem-memRange[0]] = gm[year-yearRange[0],mem-memRange[0]] + np.nansum(full[year-yearRange[0],mem-memRange[0],:,:]*areaNorm,axis=(0,1))*daysInMonth[month]/365. gm[(year-yearRange[0])*12+month,mem-memRange[0]] = np.nansum(full[year-yearRange[0],month,mem-memRange[0],:,:]*areaNorm,axis=(0,1)) print(gm.shape) print(area) nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') nc.createDimension('month',(yearRange[1]-yearRange[0]+1)*12) if memRange[1] == 0: ncgm = nc.createVariable(field,'f4',['month']) else: nc.createDimension('member',(memRange[1]-memRange[0]+1)*12) ncgm = nc.createVariable(field,'f4',['month','member']) ncgm[:] = gm nc = nc.close() return gm def boxGlobColour2(field,yearRange,memRange=[0,0],lonRange=[0,360],latRange=[-90,90],syear=0,annual=True,res=[0.5,0.5]): experiment = 'GlobColour2' lon = np.arange(res[0]/2,360,res[0]) lat = np.arange(-90+res[1]/2,90,res[1]) lon2, lat2 = np.meshgrid(lon,lat) if res == [0.5,0.5]: tagRes = '_05x05' else: tagRes = '_{:d}x{:d}'.format(res[0],res[1]) area = np.cos(lat2*np.pi/180) for j in range(len(lat)): area[j,:] = area[j,:] if lat[j] > latRange[0] and lat[j] < latRange[1] else 0. if lonRange[0] < lonRange[1]: # box lies within [-180 180] for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] and lon[i] < lonRange[1] else 0. else: for i in range(len(lon)): area[:,i] = area[:,i] if lon[i] > lonRange[0] or lon[i] < lonRange[1] else 0. filePathGlob = outPath(field,experiment,yearRange,lonRange,latRange,memRange=memRange,levRange=[0,0],syear=syear,annual=annual) yearRange if os.path.isfile(filePathGlob): print(filePathGlob + ' exists. Readmean from file.') nc = Dataset(filePathGlob) gm = nc.variables[field][:] nc.close() else: print(filePathGlob + ' does not exist. Generate mean.') filePath = outPath(field,experiment,yearRange,memRange=memRange,levRange=[0,0],syear=syear,annual=annual,suffix=tagRes) print(filePath) ncin = Dataset(filePath) if annual: gm = np.zeros((yearRange[1]-yearRange[0]+1)) else: gm = np.zeros((yearRange[1]-yearRange[0]+1)*12) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) if annual: rec = year - yearRange[0] fld2d = np.array(ncin.variables[field][rec,:,:]) areaNorm = area + fld2d*0 areaNorm = areaNorm / np.nansum(areaNorm) gm[(year-yearRange[0])] = np.nansum(fld2d*areaNorm,axis=(0,1)) ncin = ncin.close() nc = Dataset(filePathGlob, 'w', format='NETCDF4_CLASSIC') if annual: nc.createDimension('year',(yearRange[1]-yearRange[0]+1)) ncgm = nc.createVariable(field,'f4',['year']) ncgm[:] = gm nc.close() return gm def msftbarot(experiment,yearRange,memRange,singleInputFile=False,annual=True): fieldout = 'msftbarot' filePath = outPath(fieldout,experiment,yearRange,memRange=memRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Read barotropic stream function from file.') nc = Dataset(filePath) strmf = nc.variables[fieldout][:] nc.close() else: nc = Dataset(outPath('uflxlvl',experiment,yearRange,memRange=memRange,annual=annual)) umflx = np.nansum(nc['uflxlvl'][0,:,:,:],axis=0) umflx = np.where(umflx<1e20,umflx,0) nc.close() nc = Dataset(outPath('vflxlvl',experiment,yearRange,memRange=memRange,annual=annual)) vmflx = np.nansum(nc['vflxlvl'][0,:,:,:],axis=0) vmflx = np.where(vmflx<1e20,vmflx,0) nc.close() # [jdm, idm] = umflx.shape strmf = np.zeros([jdm, idm]) for j in range(1,jdm): strmf[j,0] = strmf[j-1,0]-umflx[j-1,0] for j in range(jdm-1,-1,-1): for i in range(1,idm): strmf[j,i] = strmf[j,i-1]+vmflx[j,i-1] # # Move reference to America strmf -= strmf[350,25] # Convert mass units to Sv strmf = strmf/1028*1e-6 # Smooth for j in range(0,jdm-1): jp1 = (j+1)%jdm for i in range(0,idm): ip1 = (i+1)%idm strmf[j,i] = 0.25 * (strmf[j,i]+strmf[j,ip1]+strmf[jp1,i]+strmf[jp1,ip1]) # Mask strmf = np.where(np.abs(umflx)+np.abs(vmflx)>0,strmf,np.nan) # nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('record',1) nc.createDimension('x',idm) nc.createDimension('y',jdm) nc.createVariable(fieldout,'f4',['record','y','x'])[0,:,:] = strmf nc = nc.close() return strmf def AMOC(experiment,yearRange,memRange,lat=26.5,useTestInput=False,singleInputFile=False,month1=1,annual=True): """ Read AMOC strength at specified latitude and compute annual mean. """ global model field = 'MOC' if model == "cesm1" or model == 'cesm2' else 'mmflxd' if lat > np.int32(lat): fieldout = 'AMOC{:d}'.format(np.int32(lat*10)) else: fieldout = 'AMOC{:d}'.format(np.int32(lat)) filePath = outPath(fieldout,experiment,yearRange,memRange=memRange,month1=month1,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Read AMOC from file.') nc = Dataset(filePath) amoc265 = nc.variables[fieldout][:] nc = nc.close() else: print(filePath + ' does not exist. Generate raw AMOC.') if annual: amoc265 = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) else: amoc265 = np.zeros([(yearRange[1]-yearRange[0]+1)*12-month1+1,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): #print('year=',year) for mem in range(memRange[0],memRange[1]+1): #monthLast = 1 if annual else 12 monthLast = 12 monthFirst = month1 if year == yearRange[0] else 1 for month in range(monthFirst,monthLast+1): if singleInputFile: fpath = outPath(field,experiment,yearRange,memRange=memRange,month1=month1,annual=False) else: fpath = filePathMon(field,experiment,year,month,mem,syear=0,useTestInput=useTestInput) nc = Dataset(fpath) #print('read lat from ',fpath) if field == 'MOC': latvec = nc.variables['lat_aux_grid'][:] else: latvec = nc.variables['lat'][:] ind26 = np.where(np.amin((latvec-lat)**2) == (latvec-lat)**2) ind26 = ind26[0] #print(ind26) if field == 'MOC': if singleInputFile: rec = (year - yearRange[0])*12 + month - 1 - month1 + 1 stmf26 = np.sum(nc.variables[field][mem-1,rec,1,:,:,ind26],axis=0) else: stmf26 = np.sum(nc.variables[field][0,1,:,:,ind26],axis=0) else: if singleInputFile: rec = (year - yearRange[0])*12 + month - 1 - month1 + 1 stmf26 = (nc.variables[field][mem-1,rec,0,:,ind26]+nc.variables[field][mem-1,rec,0,:,ind26+1])*0.5 else: stmf26 = (nc.variables[field][0,0,:,ind26]+nc.variables[field][0,0,:,ind26+1])*0.5 if annual: amoc265[year-yearRange[0],mem-memRange[0]] = amoc265[year-yearRange[0],mem-memRange[0]] + np.amax(stmf26)*daysInMonth[month-1]/365. else: amoc265[(year-yearRange[0])*12+month-1-month1+1,mem-memRange[0]] = np.amax(stmf26) nc = nc.close() nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('member',memRange[1]-memRange[0]+1) if annual: nc.createDimension('year',yearRange[1]-yearRange[0]+1) ncamoc265 = nc.createVariable(fieldout,'f4',['year','member']) else: nc.createDimension('record',(yearRange[1]-yearRange[0]+1)*12-month1+1) ncamoc265 = nc.createVariable(fieldout,'f4',['record','member']) ncamoc265[:] = amoc265 nc = nc.close() return amoc265 def AMOC265(experiment,yearRange,memRange,useTestInput=False): """ Read AMOC strength at 26.5N and compute annual. """ global ind26, model field = 'MOC' if model == "cesm1" or model == 'cesm2' else 'mmflxd' filePath = outPath('AMOC265',experiment,yearRange,memRange=memRange) if os.path.isfile(filePath): print(filePath + ' exists. Read AMOC from file.') nc = Dataset(filePath) amoc265 = nc.variables['amoc265'][:] nc = nc.close() else: print(filePath + ' does not exist. Generate raw AMOC.') amoc265 = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): nc = Dataset(filePathMon(field,experiment,year,month+1,mem,syear=0,useTestInput=useTestInput)) if not 'ind26' in globals(): print('read lat from ',filePathMon(field,experiment,year,month+1,mem,syear=0,useTestInput=useTestInput)) if field == 'MOC': lat = nc.variables['lat_aux_grid'][:] ind26 = np.where(np.amin((lat-26)**2) == (lat-26)**2) else: lat = nc.variables['lat'][:] ind26 = np.where(np.amin((lat-26.5)**2) == (lat-26.5)**2) ind26 = ind26[0] print(ind26) if field == 'MOC': stmf26 = np.sum(nc.variables[field][0,1,:,:,ind26],axis=0) else: stmf26 = (nc.variables[field][0,0,:,ind26]+nc.variables[field][0,0,:,ind26+1])*0.5 amoc265[year-yearRange[0],mem-memRange[0]] = amoc265[year-yearRange[0],mem-memRange[0]] + np.amax(stmf26)*daysInMonth[month]/365. nc = nc.close() #print(year,mem) nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('member',memRange[1]-memRange[0]+1) nc.createDimension('year',yearRange[1]-yearRange[0]+1) ncamoc265 = nc.createVariable('amoc265','f4',['year','member']) ncamoc265[:] = amoc265 nc = nc.close() return amoc265 def mhflx(experiment,yearRange,memRange,lat=48,useTestInput=False,singleInputFile=False,annual=True): """ Get annual mean of meridional heat flux at given latitude. """ #global indlat, model global model field = 'N_HEAT' if model == "cesm1" or model == 'cesm2' else 'mhflx' filePath = outPath(f'mhflx{lat}',experiment,yearRange,memRange=memRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Read from file.') nc = Dataset(filePath) data = nc.variables[f'mhflx{lat}'][:] nc = nc.close() else: print(filePath + ' does not exist. Generate raw heat transport.') if annual: data = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) else: data = np.zeros([(yearRange[1]-yearRange[0]+1)*12,memRange[1]-memRange[0]+1]) daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): if singleInputFile: fpath = outPath(field,experiment,yearRange,memRange=memRange,annual=False,ensave=False) else: fpath = filePathMon(field,experiment,year,month+1,mem,syear=0,useTestInput=useTestInput) nc = Dataset(fpath) if not 'indlat' in globals(): print('read lat from ',fpath) if model == 'cesm1' or model == 'cesm2': latvec = nc.variables['lat_aux_grid'][:] indlat = np.where(np.amin((latvec-lat)**2) == (latvec-lat)**2) else: latvec = nc.variables['lat'][:] indlat = np.where(np.amin((latvec-lat)**2) == (latvec-lat)**2) indlat = indlat[0] if model == 'cesm1' or model == 'cesm2': if singleInputFile: rec = (year - yearRange[0])*12 + month -1 fld = nc.variables[field][mem-1,rec,1,1,indlat]+ nc.variables[field][mem-1,rec,1,3,indlat]+nc.variables[field][mem-1,rec,1,4,indlat] else: fld = nc.variables[field][0,1,1,indlat]+nc.variables[field][0,1,3,indlat]+nc.variables[field][0,1,4,indlat] else: fld = (nc.variables[field][0,0,indlat]+nc.variables[field][0,0,indlat+1])*0.5 if annual: data[year-yearRange[0],mem-memRange[0]] = data[year-yearRange[0],mem-memRange[0]] + fld*daysInMonth[month]/365. else: data[(year-yearRange[0])*12+month,mem-memRange[0]] = fld print(year,mem,month,fld) nc = nc.close() print(year,mem) nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('member',memRange[1]-memRange[0]+1) if annual: nc.createDimension('year',yearRange[1]-yearRange[0]+1) ncmhflx = nc.createVariable(f'mhflx{lat}','f4',['year','member']) else: nc.createDimension('record',(yearRange[1]-yearRange[0]+1)*12) ncmhflx = nc.createVariable(f'mhflx{lat}','f4',['record','member']) ncmhflx[:] = data nc = nc.close() return data def mhflxatl(experiment,yearRange,memRange,useTestInput=False,singleInputFile=False,annual=True): """ Get annual mean of meridional heat flux at given latitude. """ #global indlat, model global model field = 'N_HEAT' if model == "cesm1" or model == 'cesm2' else 'mhflx' filePath = outPath(f'mhflx',experiment,yearRange,memRange=memRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Read from file.') nc = Dataset(filePath) data = nc.variables[f'mhflxatl'][:] nc = nc.close() else: print(filePath + ' does not exist. Generate raw heat transport.') if annual: data = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) else: data = np.zeros([(yearRange[1]-yearRange[0]+1)*12,memRange[1]-memRange[0]+1]) daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): if singleInputFile: fpath = outPath(field,experiment,yearRange,memRange=memRange,annual=False,ensave=False) else: fpath = filePathMon(field,experiment,year,month+1,mem,syear=0,useTestInput=useTestInput) nc = Dataset(fpath) if not 'indlat' in globals(): print('read lat from ',fpath) if model == 'cesm1' or model == 'cesm2': latvec = nc.variables['lat_aux_grid'][:] else: latvec = nc.variables['lat'][:] if model == 'cesm1' or model == 'cesm2': if singleInputFile: rec = (year - yearRange[0])*12 + month -1 fld = nc.variables[field][mem-1,rec,1,1,:]+ nc.variables[field][mem-1,rec,1,3,:]+nc.variables[field][mem-1,rec,1,4,indlat] else: fld = nc.variables[field][0,1,1,:]+nc.variables[field][0,1,3,:]+nc.variables[field][0,1,4,:] else: fld = nc.variables[field][0,0,:] if annual: if year == yearRange[0] and month == 0 and mem == memRange[0]: data = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1,len(latvec)]) else: data[year-yearRange[0],mem-memRange[0],:] = data[year-yearRange[0],mem-memRange[0],:] + fld*daysInMonth[month]/365. else: if year == yearRange[0] and month == 0 and mem == memRange[0]: data = np.zeros([(yearRange[1]-yearRange[0]+1)*12,memRange[1]-memRange[0]+1,len(latvec)]) else: data[(year-yearRange[0])*12+month,mem-memRange[0],:] = fld print(experiment,year,mem) nc = nc.close() print(year,mem) nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('member',memRange[1]-memRange[0]+1) nc.createDimension('lat',len(latvec)) if annual: nc.createDimension('year',yearRange[1]-yearRange[0]+1) ncmhflx = nc.createVariable(f'mhflxatl','f4',['year','member','lat']) else: nc.createDimension('record',(yearRange[1]-yearRange[0]+1)*12) ncmhflx = nc.createVariable(f'mhflxatl','f4',['record','member','lat']) ncmhflx[:] = data nc = nc.close() return data def msflx(experiment,yearRange,memRange,lat=48,useTestInput=False,singleInputFile=False,annual=True): """ Get annual mean of meridional heat flux at given latitude. """ #global indlat, model global model field = 'N_SALT' if model == "cesm1" or model == 'cesm2' else 'msflx' filePath = outPath(f'msflx{lat}',experiment,yearRange,memRange=memRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Read from file.') nc = Dataset(filePath) data = nc.variables[f'msflx{lat}'][:] nc = nc.close() else: print(filePath + ' does not exist. Generate raw heat transport.') if annual: data = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) else: data = np.zeros([(yearRange[1]-yearRange[0]+1)*12,memRange[1]-memRange[0]+1]) daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): if singleInputFile: fpath = outPath(field,experiment,yearRange,memRange=memRange,annual=False,ensave=False) else: fpath = filePathMon(field,experiment,year,month+1,mem,syear=0,useTestInput=useTestInput) nc = Dataset(fpath) if not 'indlat' in globals(): print('read lat from ',fpath) if model == 'cesm1' or model == 'cesm2': latvec = nc.variables['lat_aux_grid'][:] indlat = np.where(np.amin((latvec-lat)**2) == (latvec-lat)**2) else: latvec = nc.variables['lat'][:] indlat = np.where(np.amin((latvec-lat)**2) == (latvec-lat)**2) indlat = indlat[0] if model == 'cesm1' or model == 'cesm2': if singleInputFile: rec = (year - yearRange[0])*12 + month -1 fld = nc.variables[field][mem-1,rec,1,1,indlat]+ nc.variables[field][mem-1,rec,1,3,indlat]+nc.variables[field][mem-1,rec,1,4,indlat] else: fld = nc.variables[field][0,1,1,indlat]+nc.variables[field][0,1,3,indlat]+nc.variables[field][0,1,4,indlat] else: fld = (nc.variables[field][0,0,indlat]+nc.variables[field][0,0,indlat+1])*0.5 if annual: data[year-yearRange[0],mem-memRange[0]] = data[year-yearRange[0],mem-memRange[0]] + fld*daysInMonth[month]/365. else: data[(year-yearRange[0])*12+month,mem-memRange[0]] = fld print(year,mem,month,fld) nc = nc.close() print(year,mem) nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('member',memRange[1]-memRange[0]+1) if annual: nc.createDimension('year',yearRange[1]-yearRange[0]+1) ncmhflx = nc.createVariable(f'msflx{lat}','f4',['year','member']) else: nc.createDimension('record',(yearRange[1]-yearRange[0]+1)*12) ncmhflx = nc.createVariable(f'msflx{lat}','f4',['record','member']) ncmhflx[:] = data nc = nc.close() return data def ann2mon(fld,experiment,yearRange,memRange,res): return #SHF_LENS_192001-208012_mem01-40_05x05.nc def climatology(fld,experiment,yearRange,memRange,annual=False): filePath = outPath(fld,experiment,yearRange,memRange=memRange,suffix='_clim') if os.path.isfile(filePath): print(filePath + ' exists. Read climatology from file.') nc = Dataset(filePath,mode='r') clim = nc.variables[fld][:] nc.close() else: print(filePath + ' does not exist. Compute climatology from model output.') first = True nmonth = 1 if annual else 12 for year in range(yearRange[0],yearRange[1]+1): print(year) for mem in range(memRange[0],memRange[1]+1): for month in range(nmonth): if annual: nc = Dataset(filePathMon(fld,experiment,year,0,mem)) w = 1 else: nc = Dataset(filePathMon(fld,experiment,year,month+1,mem)) w = weightsAnnMean[month] data = nc.variables[fld][0,:] nc.close() if first: clim = np.zeros(data.shape) first = False clim = clim + data/(yearRange[1]-yearRange[0]+1)/(memRange[1]-memRange[0]+1)*w nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') dims = [] if len(clim.shape) >= 4: nc.createDimension('l',clim.shape[-4]) dims.append('l') if len(clim.shape) >= 3: nc.createDimension('k',clim.shape[-3]) dims.append('k') if len(clim.shape) >= 2: nc.createDimension('j',clim.shape[-2]) dims.append('j') nc.createDimension('i',clim.shape[-1]) dims.append('i') varid = nc.createVariable(fld,'f4',dims) varid[:] = clim nc = nc.close() return clim def merHeatTrans(experiment,yearRange,memRange,annualOcn=False): filePath = outPath('merHeatTrans',experiment,yearRange,memRange=memRange) atmFlds = ['lat_atm','mhflx_atm','mhflx_atm_n','mhflx_atm_s'] ocnFlds = ['lat_ocn','mhflx_ocn','mhflx_atl','mhflx_pac'] if os.path.isfile(filePath): print(filePath + ' exists. Read meridional heat transport from file.') nc = Dataset(filePath,mode='r') out={} for fld in atmFlds + ocnFlds: out[fld] = nc.variables[fld][:] nc.close() else: print(filePath + ' does not exist. Extract meridional heat transport from model output.') readGridOcn() readGridAtm() out={} # # compute atmospheric transport with residual method res = np.zeros(areaAtm.shape) for year in range(yearRange[0],yearRange[1]+1): print(year) for mem in range(memRange[0],memRange[1]+1): for month in range(12): nc = Dataset(filePathMon('FLNT',experiment,year,month+1,mem)) atmIn = {} for fld in ['FLNT','FSNT','FLNS','FSNS','LHFLX','SHFLX']: atmIn[fld] = nc.variables[fld][0,:] #PREC=(ncread(fname,'PRECSC')+ncread(fname,'PRECSL'))*3.337e8; atmIn['PREC'] = 0 nc.close() res = res - (-atmIn['FSNT']+atmIn['FLNT']+atmIn['FSNS']-atmIn['FLNS']-atmIn['LHFLX']-atmIn['SHFLX']-atmIn['PREC'])*weightsAnnMean[month]/(yearRange[1]-yearRange[0]+1)/(memRange[1]-memRange[0]+1) resIntX = np.sum(res*areaAtm*6.37122e6**2*1e-15,axis=1) nlat = len(resIntX) # out['lat_atm'] = np.zeros(nlat+1) out['lat_atm'][0], out['lat_atm'][-1] = -90, 90 for j in range(nlat-1): out['lat_atm'][j+1] = 0.5*(latAtm[j]+latAtm[j+1]) # out['mhflx_atm_n'] = np.zeros(nlat+1) for j in reversed(range(nlat)): out['mhflx_atm_n'][j] = out['mhflx_atm_n'][j+1] - resIntX[j] out['mhflx_atm_s'] = np.zeros(nlat+1) for j in range(nlat): out['mhflx_atm_s'][j+1] = out['mhflx_atm_s'][j] + resIntX[j] out['mhflx_atm'] = np.zeros(nlat+1) for j in range(nlat+1): out['mhflx_atm'][j] = (1-j/nlat)*out['mhflx_atm_s'][j] + j/nlat*out['mhflx_atm_n'][j] # # extract ocean transports first = True nmonth = 1 if annualOcn else 12 for year in range(yearRange[0],yearRange[1]+1): print(year) for mem in range(memRange[0],memRange[1]+1): for month in range(nmonth): if annualOcn: nc = Dataset(filePathMon('mhflx',experiment,year,0,mem)) else: nc = Dataset(filePathMon('mhflx',experiment,year,month+1,mem)) if first: out['lat_ocn'] = nc.variables['lat'][:] nlat = len(out['lat_ocn']) mhflx = np.zeros(nlat) first = False w = 1 if annualOcn else weightsAnnMean[month] mhflx = mhflx + 1e-15*nc.variables['mhflx'][0,:]/(yearRange[1]-yearRange[0]+1)/(memRange[1]-memRange[0]+1)*w nc.close() out['mhflx_atl'] = mhflx[0,:] out['mhflx_pac'] = mhflx[1,:] out['mhflx_ocn'] = mhflx[2,:] nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('lat_atm',len(out['lat_atm'])) nc.createDimension('lat_ocn',len(out['lat_ocn'])) for fld in atmFlds: varid = nc.createVariable(fld,'f4',['lat_atm']) varid[:] = out[fld] for fld in ocnFlds: varid = nc.createVariable(fld,'f4',['lat_ocn']) varid[:] = out[fld] nc = nc.close() return out def merHeatFlux(experiment,yearRange,memRange,lat1,useTestInput=False): """ Extract monthly Atlantic meridional heat transport at selected latitude and compute annual mean. Units are in TW. """ global latMer,lat1Last,indMer,weightsMer filePath = outPath("merHeatFlux{:.0f}N".format(lat1),experiment,yearRange,memRange=memRange) if os.path.isfile(filePath): print(filePath + ' exists. Read meridional heat transport from file.') nc = Dataset(filePath) merHeatFlux = nc.variables['merHeatFlux'][:] nc = nc.close() else: print(filePath + ' does not exist. Extract meridional heat transport from model output.') merHeatFlux = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in range(12): nc = Dataset(filePathMon('mhflx',experiment,year,month+1,mem,syear=0,useTestInput=useTestInput)) if not 'latMer' in globals(): latMer = nc.variables['lat'][:] if not 'lat1Last' in globals() or lat1Last != lat1: indMer = [0,0] for i in range(latMer.size): if latMer[i] <= lat1: indMer[0] = i indMer[1] = min(indMer[0]+1,latMer.size) weightsMer = np.zeros(2) weightsMer[0] = (latMer[indMer[1]]-lat1)/(latMer[indMer[1]]-latMer[indMer[0]]) weightsMer[1] = 1 - weightsMer[0] merHeatFlux[year-yearRange[0],mem-memRange[0]] = merHeatFlux[year-yearRange[0],mem-memRange[0]] + daysInMonth[month]/365.*1e-12*np.dot(weightsMer,nc.variables['mhflx'][0,0,indMer]) nc.close() print(year,mem,merHeatFlux[year-yearRange[0],mem-memRange[0]]) nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('member',memRange[1]-memRange[0]+1) nc.createDimension('year',yearRange[1]-yearRange[0]+1) ncmerHeatFlux = nc.createVariable('merHeatFlux','f4',['year','member']) ncmerHeatFlux[:] = merHeatFlux nc = nc.close() return merHeatFlux def plotACC(lon,lat,fld,msk=False,imageFormat='png',filePrefix='test',lbLabelBarOn=False,title='',title2='',plottype='',box=False,landFill=True,oceanFill=False,lakeFill=False,polar='',clipCorner='',printFDR=True,printLOC=True,printGM=False): import Ngl # use difference if two fields are provided pfdr = '' if len(fld) >= 4: if fld[3] == 0: pfdr = 'pFDR = NA' pfdr = 'pFDR = NA' else: pfdr = "pFDR = " + "{:.3f}".format(fld[3]) pfdr = "pFDR = " + "{:.3f}".format(fld[3]) cfac = 1 if type(fld) == tuple: msk = fld[1] if len(fld) > 2: mskg = fld[2] else: mskg = np.zeros(msk.shape) fld = fld[0] msk = np.where(np.isnan(fld),0,msk) msk = np.where(np.absolute(fld)<0.1,1,msk) mskg = np.where(np.isnan(fld),0,mskg) mskg = np.where(np.absolute(fld)<0.1,0,mskg) elif type(fld) == list: fld = fld[0] - fld[1] cfac = 1 # 0.5 fld = np.where(fld < -10.,-10.,fld) if printGM: lon2, lat2 = np.meshgrid(lon,lat) area = np.cos(lat2*np.pi/180.) + 0*fld area = area/np.nansum(area) gm = np.nansum(area*fld) # dublicate first row lon = np.ma.concatenate((lon,lon[0:1]+360)) fld = np.ma.concatenate((fld,fld[:,0:1]),axis=1) msk = np.ma.concatenate((msk,msk[:,0:1]),axis=1) mskg = np.ma.concatenate((mskg,mskg[:,0:1]),axis=1) # extend to poles if polar in ['NH','SH']: lat = np.ma.concatenate(([-90],lat,[90])) fld = np.ma.concatenate((fld[0:1,:],fld,fld[-1:,:]),axis=0) msk = np.ma.concatenate((msk[0:1,:],msk,msk[-1:,:]),axis=0) mskg = np.ma.concatenate((mskg[0:1,:],mskg,mskg[-1:,:]),axis=0) # create workstation suffix = '_' + polar if polar in ['SH','NH'] else '' wks = Ngl.open_wks(imageFormat,'Figs/' + filePrefix + suffix) res = res if 'res' in globals() else Ngl.Resources() # Contour resources res.cnFillOn = True res.cnFillMode = "RasterFill" #res.cnFillPalette = '(/(/.8,.8,.8/), (/.047,.392,.498/), (/.074,.525,.635/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.741,.890,.890/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.929,.223/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.898,.152,.176/), (/.580,.066,.105/), (/.8,.8,.8/)/)' #res.cnFillPalette = '(/(/.9,.9,.9/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/), (/.047,.392,.498/), (/.074,.525,.635/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.741,.890,.890/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.929,.223/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.898,.152,.176/), (/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/), (/.9,.9,.9/)/)' res.cnFillPalette = '(/(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/), (/.047,.392,.498/), (/.074,.525,.635/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.741,.890,.890/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.929,.223/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.898,.152,.176/), (/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/), (/0.8,0.8,0.8/)/)' #res.cnFillPalette = '(/(/.047,.392,.498/), (/.074,.525,.635/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.992,.992,.992/), (/.992,.992,.992/),(/.992,.992,.992/), (/.992,.992,.992/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.898,.152,.176/), (/.580,.066,.105/)/)' # cnres.cnFillPalette = '(/(/.074,.525,.635/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.741,.890,.890/), (/.992,.992,.992/),(/.992,.992,.992/), (/.992,.992,.992/), (/.992,.992,.992/),(/.996,.929,.223/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.898,.152,.176/),/)' # cnres.cnFillPalette = '(/ (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.992,.992,.992/), (/.992,.992,.992/),(/.992,.992,.992/), (/.992,.992,.992/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/)/)' # cnres.cnFillPalette = '(/(/.047,.392,.498/), (/.109,.686,.776/), (/.113,.698,.776/), (/.113,.721,.694/), (/.525,.815,.811/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.882,.196/), (/.964,.588,.192/), (/.941,.403,.184/), (/.929,.282,.184/), (/.580,.066,.105/)/)' #cnres.cnFillPalette = '(/(/.047,.392,.498/), (/.109,.686,.776/), (/.113,.698,.776/), (/.113,.721,.694/), (/.992,.992,.992/), (/.992,.992,.992/),(/.992,.992,.992/), (/.992,.992,.992/), (/.964,.588,.192/), (/.941,.403,.184/), (/.929,.282,.184/), (/.580,.066,.105/)/)' #res.cnFillPalette = '(/ (/.8,.8,.8/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.741,.890,.890/), (/.992,.992,.992/),(/.992,.992,.992/), (/.992,.992,.992/),(/.992,.992,.992/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.929,.223/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.8,.8,.8/) /)' #res.cnFillPalette = "BlueWhiteOrangeRed" #"BlueYellowRed" # New in PyNGL 1.5.0 res.cnLinesOn = False res.cnLineLabelsOn = False res.cnLevelSelectionMode = "ExplicitLevels" #-- use explicit levels res.cnLevels = np.arange(-2,2.1,0.1)*cfac #res.cnLevels = np.arange(-1,1.1,0.1)*cfac #res.cnLevels = np.arange(-0.9,1,0.1)*cfac #cnres.cnLevels = np.arange(-0.8,1,0.2)*cfac # Labelbar resource res.lbOrientation = "horizontal" #res.lbOrientation = "vertical" res.lbLabelBarOn = lbLabelBarOn res.lbLabelFontHeightF = 0.01 # Scalar field resources res.sfXArray = lon res.sfYArray = lat #res.sfYCStartV = -30. # Map resources res.mpOutlineOn = (not (landFill or oceanFill)) res.mpGeophysicalLineThicknessF = 2 res.mpGridLonSpacingF = 30. res.mpGridLatSpacingF = 30. #res.mpGridAndLimbDrawOrder = "PreDraw" res.mpCenterLonF = 0# -60. #-90. # -150. res.mpFillOn = (landFill or oceanFill) res.mpFillDrawOrder = "PostDraw" res.mpLandFillColor = "Black" if landFill else "Transparent" res.mpOceanFillColor = "White" if oceanFill else "Transparent" res.mpInlandWaterFillColor = "Black" if lakeFill else "Transparent" res.pmTickMarkDisplayMode = "NoCreate" if polar in ['NH','SH']: res.mpProjection = 'Stereographic' #-- set projection res.mpEllipticalBoundary = True #-- map projection area is limited to an ellipse #-- inscribed within the normal rectangular #-- perimeter of the viewport #res.mpDataSetName = 'Earth..4' #-- change map data set #res.mpDataBaseVersion = 'MediumRes' #-- choose higher map resolution res.mpLimitMode = 'LatLon' if polar == 'NH': res.mpMaxLatF = 90. res.mpMinLatF = 50. #-- minimum latitude res.mpCenterLatF = 90. #-- center latitude else: res.mpMaxLatF = -50. res.mpMinLatF = -90. #-- minimum latitude res.mpCenterLatF = -90. #-- center latitude #res.mpLimitMode = "LatLon" # Limit the map view. #res.mpMinLonF = float(-120) #res.mpMaxLonF = float(0) #res.mpMinLatF = float(-45) #res.mpMaxLatF = float(90) #res.mpPerimOn = True # Turn on map perimeter. # Title #res.tiMainString = title # plot to file #if type(box) != bool: res.nglDraw = False res.nglFrame = False map = Ngl.map(wks,res) plot = Ngl.contour(wks,fld,res) Ngl.overlay(map,plot) # mark areas with / that are not locally significant if printLOC: res.cnMonoFillColor = True res.cnLevels = [0.5] res.cnFillPatterns = [-1,3] res.cnMonoFillScale = True res.cnFillScaleF = 0.5 res.cnMonoFillPattern = False res.cnFillMode = "AreaFill" plotmask = Ngl.contour(wks,msk,res) Ngl.overlay(map,plotmask) # mark areas with . that are field significant if printFDR: npoint = np.sum(mskg).astype(int) lonp = np.zeros(npoint) latp = np.zeros(npoint) count = 0 for j in range(len(lat)): for i in range(len(lon)): if mskg[j,i] > 0.5 and fld[j,i] < 1.99: lonp[count] = lon[i] latp[count] = lat[j] count = count + 1 res.gsMarkerIndex = 1 res.gsMarkerColor = 'black' Ngl.add_polymarker(wks,map,lonp,latp,res) # plot labels #res.txFontHeightF = 0.024 #Ngl.add_text(wks,map,plottype,90,55,res) #Ngl.add_text(wks,map,plottype,-90,-84,res) #Ngl.add_text(wks,map,title,0,-84,res) #Ngl.add_text(wks,map,pfdr,90,-85,res) #res.txFontHeightF = 0.018 ; Ngl.add_text(wks,map,pfdr,80,-80,res) res.txFontColor = 'white' if not landFill: res.txFontColor = 'black' # res.txFont = 4 # res.txFontThicknessF = 8 if polar == 'NH': res.txFontHeightF = 0.028 ; Ngl.add_text(wks,map,title2,100,60,res) res.txFontHeightF = 0.030 ; Ngl.add_text(wks,map,title,90,60,res) if printFDR: res.txFontHeightF = 0.018 ; Ngl.add_text(wks,map,pfdr,80,60,res) elif polar == 'SH': res.txFontHeightF = 0.028 ; Ngl.add_text(wks,map,title2,0,-90,res) res.txFontHeightF = 0.030 ; Ngl.add_text(wks,map,title,0,-85,res) if printFDR: res.txFontHeightF = 0.018 ; Ngl.add_text(wks,map,pfdr,0,-90,res) else: res.txFontHeightF = 0.030 ; Ngl.add_text(wks,map,title2,10,22,res) res.txFontHeightF = 0.035 ; Ngl.add_text(wks,map,title,85,55,res) if printFDR: #res.txFontHeightF = 0.018 ; Ngl.add_text(wks,map,pfdr,85,35,res) res.txFontHeightF = 0.030 ; Ngl.add_text(wks,map,pfdr,60,-82,res) if printGM: res.txFont = 4 res.txFontThicknessF = 6 res.txFontHeightF = 0.025 ; Ngl.add_text(wks,map,'r_global = ' + '{:.2f}'.format(gm),60,-82,res) #if not landFill: # res.txFontColor = 'black' # res.txFont = 4 # res.txFontThicknessF = 13 # res.txFontHeightF = 0.030 ; Ngl.add_text(wks,map,title2,10,22,res) # res.txFontHeightF = 0.035 ; Ngl.add_text(wks,map,title,85,55,res) # res.txFontHeightF = 0.018 ; Ngl.add_text(wks,map,pfdr,85,35,res) if type(box) != bool: lres = Ngl.Resources() lres.gsLineColor = 'purple' lres.gsLineThicknessF = 5. Ngl.add_polyline(wks, map, [box[0][0],box[0][1],box[0][1],box[0][0],box[0][0]], [box[1][0],box[1][0],box[1][1],box[1][1],box[1][0]], lres) #res.gsLineColor = 'Purple' #Ngl.add_polyline(wks,plotId,np.array[-20,-20,10,10,-20],[65,75,75,65,65],cnres) Ngl.draw(map) Ngl.frame(wks) Ngl.delete_wks(wks) def plotBias(lon,lat,fld,msk=False,cfac=3,rasterSmoothing=False,imageFormat='png',filePrefix='test',lbLabelBarOn=False,title='',title2='',plottype='',box=False,landFill=True): # use difference if two fields are provided if len(fld) >= 4: if fld[3] == 0: pfdr = 'pFDR = NA' pfdr = 'pFDR = NA' else: pfdr = "pFDR = " + "{:.3f}".format(fld[3]) pfdr = "pFDR = " + "{:.3f}".format(fld[3]) if type(fld) == tuple: msk = fld[1] if len(fld) > 2: mskg = fld[2] else: mskg = np.zeros(msk.shape) fld = fld[0] msk = np.where(np.isnan(fld),0,msk) #msk = np.where(np.absolute(fld)<0.1,1,msk) mskg = np.where(np.isnan(fld),0,mskg) mskg = np.where(np.absolute(fld)<0.1,0,mskg) elif type(fld) == list: fld = fld[0] - fld[1] fld = np.where(fld < -10.,-10.,fld) # dublicate first row lon = np.ma.concatenate((lon,lon[0:1]+360)) fld = np.ma.concatenate((fld,fld[:,0:1]),axis=1) msk = np.ma.concatenate((msk,msk[:,0:1]),axis=1) mskg = np.ma.concatenate((mskg,mskg[:,0:1]),axis=1) # create workstation wks = Ngl.open_wks(imageFormat,'Figs/' + filePrefix) res = res if 'res' in globals() else Ngl.Resources() # Contour resources res.cnFillOn = True res.cnFillMode = "RasterFill" #res.cnFillMode = "AreaFill" if rasterSmoothing: res.cnRasterSmoothingOn = True res.cnFillPalette = '(/(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/), (/.047,.392,.498/), (/.074,.525,.635/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.741,.890,.890/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.929,.223/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.898,.152,.176/), (/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/), (/0.8,0.8,0.8/)/)' #res.cnFillPalette = '(/(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/),(/.047,.392,.498/), (/.047,.392,.498/), (/.074,.525,.635/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.741,.890,.890/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.929,.223/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.898,.152,.176/),(//) (/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/),(/.580,.066,.105/), (/0.8,0.8,0.8/)/)' #res.cnFillPalette = '(/(/.047,.392,.498/), (/.047,.392,.498/), (/.074,.525,.635/), (/.109,.686,.776/), (/.113,.698,.776/),(/.113,.717,.772/), (/.113,.721,.694/), (/.333,.760,.749/), (/.525,.815,.811/), (/.741,.890,.890/), (/.992,.992,.992/),(/.992,.992,.992/), (/.996,.929,.223/), (/.996,.882,.196/), (/.984,.709,.200/), (/.964,.588,.192/), (/.949,.490,.188/), (/.941,.403,.184/), (/.929,.282,.184/), (/.898,.152,.176/), (/.580,.066,.105/), (/.580,.066,.105/)/)' res.cnLinesOn = False res.cnLineLabelsOn = False res.cnLevelSelectionMode = "ExplicitLevels" #-- use explicit levels res.cnLevels = np.arange(-2,2.1,0.1)*cfac #res.cnLevels = np.arange(-1,1.1,0.1)*cfac # Labelbar resource res.lbOrientation = "horizontal" #res.lbOrientation = "vertical" #res.lbBoxEndCapStyle = "TriangleBothEnds" res.lbLabelBarOn = lbLabelBarOn res.lbLabelFontHeightF = 0.01 # Scalar field resources res.sfXArray = lon res.sfYArray = lat #res.sfYCStartV = -30. # Map resources res.mpOutlineOn = (not landFill) res.mpGeophysicalLineThicknessF = 2 res.mpGridLonSpacingF = 30. res.mpGridLatSpacingF = 30. res.mpCenterLonF = 0# -60. #-90. # -150. res.mpFillOn = landFill res.mpFillDrawOrder = "PostDraw" res.mpLandFillColor = "Black" res.mpOceanFillColor = "Transparent" res.mpInlandWaterFillColor = "Transparent" res.pmTickMarkDisplayMode = "NoCreate" res.nglDraw = False res.nglFrame = False map = Ngl.map(wks,res) plot = Ngl.contour(wks,fld,res) Ngl.overlay(map,plot) # mark areas with / that are not locally significant res.cnMonoFillColor = True res.cnLevels = [0.5] res.cnFillPatterns = [-1,3] res.cnMonoFillScale = True res.cnFillScaleF = 0.5 res.cnMonoFillPattern = False res.cnFillMode = "AreaFill" #plotmask = Ngl.contour(wks,msk,res) #Ngl.overlay(map,plotmask) # mark areas with . that are field significant npoint = np.sum(mskg).astype(int) lonp = np.zeros(npoint) latp = np.zeros(npoint) count = 0 for j in range(len(lat)): for i in range(len(lon)): if mskg[j,i] > 0.5 and fld[j,i] < 1.99: lonp[count] = lon[i] latp[count] = lat[j] count = count + 1 res.gsMarkerIndex = 1 res.gsMarkerColor = 'black' Ngl.add_polymarker(wks,map,lonp,latp,res) res.txFontColor = 'white' if not landFill: res.txFontColor = 'black' res.txFontHeightF = 0.030 ; Ngl.add_text(wks,map,title2,10,22,res) res.txFontHeightF = 0.035 ; Ngl.add_text(wks,map,title,85,55,res) #res.txFontHeightF = 0.018 ; Ngl.add_text(wks,map,pfdr,85,35,res) if type(box) != bool: lres = Ngl.Resources() lres.gsLineColor = 'purple' lres.gsLineThicknessF = 5. Ngl.add_polyline(wks, map, [box[0][0],box[0][1],box[0][1],box[0][0],box[0][0]], [box[1][0],box[1][0],box[1][1],box[1][1],box[1][0]], lres) Ngl.draw(map) Ngl.frame(wks) Ngl.delete_wks(wks) def readIceext(experiment,syearRange,smonth=10,leadRange=[1,24],memRange=[0,0]): global wn, ws, wnobs, wsobs filePath = outPathIceext(experiment,syearRange,smonth=smonth,leadRange=leadRange,memRange=memRange) if os.path.isfile(filePath): print(filePath + ' exists. Read ice extent from file.') nc = Dataset(filePath) iceext = nc.variables['iceext'][:] nc.close() else: nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('hemisphere', 2) nc.createDimension('leadmonth', leadRange[1]-leadRange[0]+1) nc.createDimension('member', memRange[1]-memRange[0]+1) nc.createDimension('startyear', syearRange[1]-syearRange[0]+1) ncvar = nc.createVariable('iceext','f4',['hemisphere','leadmonth','member','startyear']) iceext = np.zeros([2,leadRange[1]-leadRange[0]+1,memRange[1]-memRange[0]+1,syearRange[1]-syearRange[0]+1]) for syear in range(syearRange[0],syearRange[1]+1): print(syear) for leadMonth in range(leadRange[0],leadRange[1]+1): if experiment == 'HadISST-persistence': year = syear month = smonth else: year = syear + int((smonth + leadMonth)/12) month = ((smonth + leadMonth - 1) % 12 ) + 1 for member in range(memRange[0],memRange[1]+1): if experiment[0:7] == 'HadISST': filePathIn = dataObs + '/regridded/HadISST/HadISST_ice.nc' ncIn = Dataset(filePathIn) if not 'wnobs' in globals(): lon = ncIn.variables['longitude'][:] lat = ncIn.variables['latitude'][:] dlat = np.abs(lat[1]-lat[0]) * np.pi / 180 dlon = np.abs(lon[1]-lon[0]) * np.pi / 180 rearth = 6.37122e6 wnobs = np.zeros([lat.size,lon.size]) wsobs = np.zeros([lat.size,lon.size]) for j in range(lat.size): wnobs[j,:] = dlon * dlat * np.cos(lat[j]*np.pi/180) * rearth**2 if lat[j]>0 else 0 wsobs[j,:] = dlon * dlat * np.cos(lat[j]*np.pi/180) * rearth**2 if lat[j]<0 else 0 data = ncIn.variables['sic'][(year-1870)*12+month-1,:,:] ncIn.close() iceext[0,leadMonth-leadRange[0],member-memRange[0],syear-syearRange[0]] = np.nansum(data * wnobs) iceext[1,leadMonth-leadRange[0],member-memRange[0],syear-syearRange[0]] = np.nansum(data * wsobs) else: if experiment[0:7] == 'dcppA-h': filePathIn = filePathMon('aice',experiment,year,month,member=member,syear=syear,useTestInput=False) else: filePathIn = filePathMon('aice',experiment,year,month,member=member,useTestInput=False) ncIn = Dataset(filePathIn) if not 'wn' in globals(): TLAT = ncIn.variables['TLAT'][:] tarea = ncIn.variables['tarea'][:] wn = np.where(TLAT>0,tarea,0) ws = np.where(TLAT<0,tarea,0) data = ncIn.variables['aice'][:] ncIn.close() iceext[0,leadMonth-leadRange[0],member-memRange[0],syear-syearRange[0]] = np.nansum(data * wn) * 0.01 iceext[1,leadMonth-leadRange[0],member-memRange[0],syear-syearRange[0]] = np.nansum(data * ws) * 0.01 ncvar[:,:,:,:]= iceext nc.close() return iceext def corrIceext(y1,y2,detrend=True): # compute anomalies y1 = y1 - np.nanmean(y1,axis=-1,keepdims=True) y2 = y2 - np.nanmean(y2,axis=-1,keepdims=True) # compute ACC for ensemble average y1a = np.nanmean(y1,axis=-2) y2a = np.nanmean(y2,axis=-2) if detrend: x = np.arange(y1a.shape[-1]) ; x = x - np.mean(x) s = np.nanmean(x*y1a,axis=-1,keepdims=True)/np.nanmean(x*x) y1a = y1a - x * s x = np.arange(y2a.shape[-1]) ; x = x - np.mean(x) s = np.nanmean(x*y2a,axis=-1,keepdims=True)/np.nanmean(x*x) y2a = y2a - x * s r = np.nanmean(y1a*y2a,axis=-1)/np.sqrt(np.nanmean(y1a*y1a,axis=-1)*np.nanmean(y2a*y2a,axis=-1)) return r def runmean(x,wlen,axis=0): yshape = list(x.shape) yshape[axis] = yshape[axis] - wlen + 1 y = np.zeros(yshape) for n in range(yshape[axis]): if axis==0: y[n,] = np.mean(x[n:n+wlen,],axis=axis) elif axis==1: y[:,n,] = np.mean(x[:,n:n+wlen,],axis=axis) elif axis==2: y[:,:,n,] = np.mean(x[:,:,n:n+wlen,],axis=axis) elif axis==3: y[:,:,:,n,] = np.mean(x[:,:,:,n:n+wlen,],axis=axis) elif axis==4: y[:,:,:,:,n,] = np.mean(x[:,:,:,:,n:n+wlen,],axis=axis) elif axis==5: y[:,:,:,:,:,n,] = np.mean(x[:,:,:,:,:,n:n+wlen,],axis=axis) return y def readHindcastLY(field,experiment,syears,leadRange,yearRange=[0,0], memRange=[0,0],levRange=[0,0],suffix='',persistence='',ensave=True): for syear in syears: if yearRange[0] == 0: # hindcasts filePath = outPath(field,experiment,[syear+1,syear+10],memRange=memRange, levRange=levRange,syear=syear,suffix=suffix,ensave=ensave) nc = Dataset(filePath) data = nc.variables[field][:] nc.close() if syear == syears[0]: outData = np.zeros((len(syears),) + data.shape[1:]) outData[syear-syears[0],:] = np.mean(data[leadRange[0]:leadRange[1]+1,:],axis=0) else: #print(field) #print(experiment) #print(yearRange) #print(memRange) #print(levRange) #print(suffix) #print(ensave) filePath = outPath(field,experiment,yearRange,memRange=memRange, levRange=levRange,suffix=suffix,ensave=ensave) #print(filePath) if syear == syears[0]: nc = Dataset(filePath) data = nc.variables[field][:] nc.close() outData = np.zeros((len(syears),) + data.shape[1:]) offset = syear - yearRange[0] + 1 # offset is greater equal 1 if persistence == 'latest': outData[syear-syears[0],:] = data[offset-1,:] elif persistence == 'mean': outData[syear-syears[0],:] = np.mean( data[offset-1-leadRange[1]+leadRange[0]:offset,:],axis=0) else: #print(offset+leadRange[0],offset+leadRange[1]+1) outData[syear-syears[0],:] = np.mean( data[offset+leadRange[0]:offset+leadRange[1]+1,:],axis=0) #print(np.nanmax(outData[syear-syears[0],:])) outData = np.where(outData<=-999.,np.nan,outData) outData = np.where(outData==0.0,np.nan,outData) return np.squeeze(outData) def readHindcastLM(field,experiment,syears,leadRange,month1=1,yearRange=[0,0], memRange=[0,0],levRange=[0,0],suffix='',persistence='',ensave=True): for syear in syears: if yearRange[0] == 0: # hindcasts filePath = outPath(field,experiment,[syear,syear+10],memRange=memRange, levRange=levRange,month1=month1,syear=syear,annual=False,suffix=suffix,ensave=ensave) nc = Dataset(filePath) data = nc.variables[field][:] nc.close() if syear == syears[0]: outData = np.zeros((len(syears),) + data.shape[1:]) outData[syear-syears[0],:] = np.mean(data[leadRange[0]:leadRange[1]+1,:],axis=0) else: filePath = outPath(field,experiment,yearRange,memRange=memRange, levRange=levRange,month1=1,annual=False,suffix=suffix,ensave=ensave) if syear == syears[0]: nc = Dataset(filePath) data = nc.variables[field][:] nc.close() outData = np.zeros((len(syears),) + data.shape[1:]) offset = (syear - yearRange[0])*12 + month1 - 1 # offset is greater equal 1 if persistence == 'latest': outData[syear-syears[0],:] = data[offset,:] elif persistence == 'mean': outData[syear-syears[0],:] = np.mean( data[offset-leadRange[1]+leadRange[0]:offset+1,:],axis=0) else: outData[syear-syears[0],:] = np.mean( data[offset+leadRange[0]:offset+leadRange[1]+1,:],axis=0) return np.squeeze(outData) def readHindcast(field,experiment,syears,yearRange=[0,0],memRange=[0,0],levRange=[0,0],suffix=''): for syear in syears: if yearRange[0] == 0: filePath = outPath(field,experiment,[syear+1,syear+10],memRange=memRange, levRange=levRange,syear=syear,suffix=suffix) nc = Dataset(filePath) data = nc.variables[field][:] nc.close() if syear == syears[0]: outData = np.zeros((len(syears),) + data.shape) outData[syear-syears[0],:] = data else: filePath = outPath(field,experiment,yearRange,memRange=memRange, levRange=levRange,suffix=suffix) if syear == syears[0]: nc = Dataset(filePath) data = nc.variables[field][:] nc.close() if persistence: outData = np.zeros((len(syears),1) + data.shape[1:]) else: outData = np.zeros((len(syears),10) + data.shape[1:]) offset = syear - yearRange[0] + 1 if persistence: outData[syear-syears[0],:] = data[offset-1:offset,:] else: outData[syear-syears[0],:] = data[offset:offset+10,:] return outData def readHindcastBox(field,experiment,lonRange,latRange,syears,yearRange=[0,0],memRange=[0,0],levRange=[0,0],suffix='',persistence=False): for syear in syears: if yearRange[0] == 0: filePath = outPath(field,experiment,[syear+1,syear+10],memRange=memRange, levRange=levRange,syear=syear,suffix=suffix) nc = Dataset(filePath) data = nc.variables[field][:] nc.close() if syear == syears[0]: outData = np.zeros((len(syears),) + data.shape[0]) outData[syear-syears[0],:] = compBox2d1x1(lonRange,latRange,data) else: filePath = outPath(field,experiment,yearRange,memRange=memRange, levRange=levRange,suffix=suffix) if syear == syears[0]: nc = Dataset(filePath) data = nc.variables[field][:] nc.close() if persistence: outData = np.zeros((len(syears),1)) else: outData = np.zeros((len(syears),11)) offset = syear - yearRange[0] + 1 if persistence: outData[syear-syears[0],0] = compBox2d1x1(lonRange,latRange,data[offset-1,:]) for i in range(1,10): outData[syear-syears[0],i] = outData[syear-syears[0],0] else: outData[syear-syears[0],:] = compBox2d1x1(lonRange,latRange,data[offset-1:offset+10,:]) return outData def regress(x,y,pctl=(5,95),nboot=1000): xave, yave = np.mean(x), np.mean(y) xano, yano = x - xave, y - yave xvar = np.mean(xano*xano) s = np.mean(xano*yano) / xvar i = yave - xave * s if nboot > 0: np.random.seed(0) nmax = x.shape[0] dof = nmax ind = np.trunc(np.random.rand(nboot,dof)*nmax) bs = np.zeros(nboot) bi = np.zeros(nboot) for n in range(nboot): bx = np.take(x,ind[n,:].astype(int),axis=0) by = np.take(y,ind[n,:].astype(int),axis=0) bxave, byave = np.mean(bx), np.mean(by) bxano, byano = bx - bxave, by - byave bxvar = np.mean(bxano*bxano) bs[n] = np.mean(bxano*byano) / bxvar bi[n] = byave - bxave * bs[n] spctl = np.percentile(bs,pctl) ipctl = np.percentile(bi,pctl) return s, i, spctl, ipctl else: return s, i def ecsGregory(x,y,pctl=(5,95),nboot=1000): xave, yave = np.mean(x), np.mean(y) xano, yano = x - xave, y - yave xvar = np.mean(xano*xano) s = np.mean(xano*yano) / xvar i = yave - xave * s ecs = - i/s /2 if nboot > 0: np.random.seed(0) nmax = x.shape[0] dof = nmax ind = np.trunc(np.random.rand(nboot,dof)*nmax) bs = np.zeros(nboot) bi = np.zeros(nboot) becs = np.zeros(nboot) for n in range(nboot): bx = np.take(x,ind[n,:].astype(int),axis=0) by = np.take(y,ind[n,:].astype(int),axis=0) bxave, byave = np.mean(bx), np.mean(by) bxano, byano = bx - bxave, by - byave bxvar = np.mean(bxano*bxano) bs[n] = np.mean(bxano*byano) / bxvar bi[n] = byave - bxave * bs[n] becs[n] = - bi[n]/bs[n] /2 pctl = np.percentile(becs,pctl) return ecs, pctl else: return ecs def corrSelf(y): nmem=y.shape[1] r = np.zeros([nmem,y.shape[2],y.shape[3]]) for m in range(nmem): y1 = np.take(y,m,axis=1) y2 = np.mean(np.take(y,np.delete(np.arange(nmem),m),axis=1),axis=1) y1 = y1 - np.mean(y1,axis=0,keepdims=True) y2 = y2 - np.mean(y2,axis=0,keepdims=True) r[m,:,:] = np.sum(y1*y2,axis=0)/np.sqrt(np.sum(y1*y1,axis=0)*np.sum(y2*y2,axis=0)) return r def corrMultiArray(y1,y2,axis=0,nboot=1000,dof=10,pval=0.10,ensave=True): # compute correlation y1 = y1 - np.mean(y1,axis=axis,keepdims=True) y2 = y2 - np.mean(y2,axis=axis,keepdims=True) if ensave: r = np.sum(y1*y2,axis=axis)/np.sqrt(np.sum(y1*y1,axis=axis)*np.sum(y2*y2,axis=axis)) msk = np.zeros(r.shape) else: y1ave = np.mean(y1,axis=1) r = np.sum(y1ave*y2,axis=axis)/np.sqrt(np.sum(y1ave*y1ave,axis=axis)*np.sum(y2*y2,axis=axis)) # compute significance msk = np.zeros(r.shape) if ensave and nboot > 0: np.random.seed(0) nmax = y1.shape[axis] ind = np.trunc(np.random.rand(nboot,dof)*nmax) count = np.zeros(r.shape) for i in range(nboot): boot1 = np.take(y1,ind[i,:].astype(int),axis=axis) boot1 = boot1 - np.mean(boot1,axis=axis,keepdims=True) boot2 = np.take(y2,ind[i,:].astype(int),axis=axis) boot2 = boot2 - np.mean(boot2,axis=axis,keepdims=True) rboot = np.sum(boot1*boot2,axis=axis)/np.sqrt(np.sum(boot1*boot1,axis=0)*np.sum(boot2*boot2,axis=axis)) count = np.where(rboot * np.sign(r) < 0, count+1/nboot, count) msk = np.where(count0: return r,msk else: return r def detrendMulti(y): # add dummy member dimension if necessary y = y.reshape((y.shape[0],1) + (y.shape[1:])) if len(y.shape) == 3 else y # compute anomalies y = y - np.nanmean(y,axis=0,keepdims=True) # compute linear trend N = y.shape[0] yEnsAve = np.nanmean(y,axis=1,keepdims=True) x = np.ones(yEnsAve.shape) * (np.arange(0,N)-(N-1)/2).reshape((N,1,1,1)) slope = np.nanmean(x*yEnsAve,axis=0)/np.nanmean(x*x,axis=0) trend = slope * x return np.squeeze(y - trend) def bias(y1,y2): """ Compute bias. """ y1Clim = np.nanmean(y1,axis=(0,1)) y2Clim = np.nanmean(y2,axis=0) b = y1Clim - y2Clim return b def bootstrapMean(y,prcnt,nboot=10000): np.random.seed(0) indMem=np.trunc(np.random.rand(nboot,y.shape[-1])*y.shape[-1]).astype(int) ymean=np.zeros([nboot,y.shape[0]]) for i in range(nboot): ymean[i,:]=np.mean(y.take(indMem[i,:],axis=1),axis=1) return np.percentile(ymean,prcnt,axis=0) def corrMultiArrayYeager(y1,y2,nboot=4000,blocklen=5,pval=0.10): """ Same as corrMultiArray but using block bootstrapping for significance testing. The algorithm estimates ACC 5-95% confidence bounds by resampling mod,obs paired blocks. Each bootstrap uses a different combination of members for estimating the ensemble mean. We regard the correlation estimate significantely different from zero if both confidence bounds have the same sign. p - probability to by chance obtain a correlation value that is as large or larger than r """ # compute anomalies y1 = y1 - np.nanmean(y1,axis=0,keepdims=True) y2 = y2 - np.nanmean(y2,axis=0,keepdims=True) # compute ACC for ensemble average y1a = np.nanmean(y1,axis=1) r = np.nanmean(y1a*y2,axis=0)/np.sqrt(np.nanmean(y1a*y1a,axis=0)*np.nanmean(y2*y2,axis=0)) # compute local significance maxlen, nmem = y1.shape[0], y1.shape[1] nblock = np.int32(np.ceil(maxlen / blocklen)) np.random.seed(0) indBlockStart = np.trunc(np.random.rand(nboot,nblock)*(maxlen-(blocklen-1))) indBlockMem = np.trunc(np.random.rand(nboot,nmem)*nmem).astype(int) ind = np.zeros(nblock*blocklen,dtype=int) p = np.zeros(r.shape) for i in range(nboot): for j in range(nblock): ind[j*blocklen:(j+1)*blocklen] = indBlockStart[i,j] + np.arange(blocklen) y1a = np.nanmean(y1.take(indBlockMem[i,:],axis=1),axis=1) boot1 = y1a.take(ind[0:maxlen],axis=0) boot1 = boot1 - np.nanmean(boot1,axis=0,keepdims=True) boot2 = y2.take(ind[0:maxlen],axis=0) boot2 = boot2 - np.nanmean(boot2,axis=0,keepdims=True) rboot = np.nanmean(boot1*boot2,axis=0)/np.sqrt(np.nanmean(boot1*boot1,axis=0)*np.nanmean(boot2*boot2,axis=0)) p = np.where(rboot * np.sign(r) < 0, p+1/nboot, p) mskl = np.where(p2d}'.format(year)) monthLast = 1 if annual else 12 for month in range(1,monthLast+1): rec = year-yearRange[0] if annual else (year-yearRange[0])*12 + month - 1 if ensave: mem1, memn = 1, 1 else: mem1, memn = memRange[0], memRange[1] for mem in range(mem1,memn+1): # read input data dataIn = ncIn.variables[field][rec,:] if ensave else ncIn.variables[field][rec,mem-1,:] # compute weights if year == yearRange[0] and month == 1 and mem==1: if type(levBounds) == int: if cesm: ncdz = Dataset(dataTmp + '/cesm1_dz.nc') dz = ncdz.variables['dz'][:]*1e-2 ncdz.close() levBounds = np.zeros((dz.shape[0],2)) for k in range(dz.shape[0]): levBounds[k,1] = levBounds[k,0] + dz[k] if k < dz.shape[0] - 1: levBounds[k+1,0] = levBounds[k,1] else: #levBounds = ncIn.variables['depth_bnds'][:] if field == 'templvl' or field == 'salnlvl': if dataIn.shape[0] == 35: ncdz = Dataset(dataTmp + '/coords_noresm1_L35.nc') elif dataIn.shape[0] == 70: ncdz = Dataset(dataTmp + '/coords_noresm1_L70.nc') else: ncdz = Dataset(dataTmp + '/coords_EN4.nc') levBounds = ncdz.variables['depth_bnds'][:] ncdz.close() weights = np.zeros((levBounds.shape[0],dataIn.shape[-2],dataIn.shape[-1])) if type(bath) == int: maxDepth = np.zeros((dataIn.shape[-2],dataIn.shape[-1])) else: maxDepth = bath print(maxDepth.shape) print(bath) for k in range(levBounds.shape[0]): print(dataIn.mask.shape) if type(bath) == int: maxDepth = np.where(dataIn.mask[k,:,:],maxDepth,levBounds[k,1]) zTop = np.minimum(max(levRange[0],levBounds[k,0]),maxDepth) zBot = np.minimum(min(levRange[1],levBounds[k,1]),maxDepth) print(zTop.shape,zBot.shape) weights[k,:,:] = np.maximum(zBot - zTop,1e-10) weights = weights / np.sum(weights,axis=0,keepdims=True) # compute vertical average data = np.nansum(dataIn * weights, axis=-3) # prepare output file if year == yearRange[0] and month == 1 and mem==1: nc.createDimension('time', None) if ensave: dims=['time'] else: nc.createDimension('member',memRange[1]-memRange[0]+1) dims=['time','member'] if len(data.shape) > 3: nc.createDimension('l',data.shape[len(data.shape)-4]) dims.extend('l') if len(data.shape) > 2: nc.createDimension('k',data.shape[len(data.shape)-3]) dims.extend('k') if len(data.shape) > 1: nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') if len(data.shape) > 0: nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') ncvar = nc.createVariable(field,'f4',dims,zlib=True) ncdepth = nc.createVariable('depth','f4',['j','i'],zlib=True) #print(dims) #print(data.shape) #print(ensave) if rec == 0: ncdepth[:,:] = maxDepth if ensave: ncvar[rec,:] = data else: ncvar[rec,mem-1,:] = data ncIn.close() nc.close() if clim: yearRange[1] = yearRange1_old def regrid2reg(field,experiment,yearRange,memRange=[0,0],syear=0,res=[1.,1.],lonBounds=[0.,360.], latBounds=[-90.,90.],levRange=[0,0],month1=1,annual=True,clim=False,ensave=True,isMasked=True,hasDepth=False,suffix=''): if res[0] == 0.25: osuffix = suffix + '_025x025' elif res[0] == 0.5: osuffix = suffix + '_05x05' elif res[0] == 1.25: osuffix = suffix + '_f09' elif res[0] >= 1: osuffix = suffix + '_' + '{:d}'.format(res[0]) + 'x' + '{:d}'.format(res[1]) filePathIn = outPath(field,experiment,yearRange,syear=syear, memRange=memRange,levRange=levRange,month1=month1,annual=annual,suffix=suffix,ensave=ensave) filePath = outPath(field,experiment,yearRange,month1=month1,syear=syear,memRange=memRange,levRange=levRange, annual=annual,suffix=osuffix,ensave=ensave) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') ncIn = Dataset(filePathIn) if hasDepth: print(filePathIn) depth = ncIn.variables['depth'][:,:] nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') if clim: yearRange1_old = yearRange[1] yearRange[1] = yearRange[0] for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) monthLast = 1 if annual else 12 monthFirst = month1 if year == yearRange[0] else 1 for month in range(monthFirst,monthLast+1): rec = year-yearRange[0] if annual else (year-yearRange[0])*12 + month - 1 - (month1 - 1) # read input data #print(rec,ncIn.variables[field][:].shape) dataIn = np.squeeze(ncIn.variables[field][rec,:]) print(dataIn.shape) # regrid data if experiment == 'EN4': data, lon, lat = reg2reg(dataIn,resIn=1.,lonBoundsIn=[0.5,360.5],latBoundsIn=[-83.5,89.5], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized') elif experiment == 'ERSSTv5': data, lon, lat = reg2reg(dataIn,resIn=2.,lonBoundsIn=[-1.0,359.0],latBoundsIn=[-89.0,89.0], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized') elif experiment == 'ERA5': data, lon, lat = reg2reg(dataIn,resIn=1./4.,lonBoundsIn=[-0.25,359.75],latBoundsIn=[-90.125,90.125], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized') elif experiment == 'HadISST': data, lon, lat = reg2reg(dataIn,resIn=1.,lonBoundsIn=[-180.,180.0],latBoundsIn=[-90.0,90.0], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized') elif experiment == 'ARMOR3D': data, lon, lat = reg2reg(dataIn,resIn=1.,lonBoundsIn=[-0.5,359.5],latBoundsIn=[-90.0,90.0], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized',dynamicMask=True) elif experiment == 'OISST': data, lon, lat = reg2reg(dataIn,resIn=1./4.,lonBoundsIn=[0.,360.],latBoundsIn=[-90.0,90.0], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized',dynamicMask=False) elif experiment == 'ARMOR3D025': data, lon, lat = reg2reg(dataIn,resIn=1./4.,lonBoundsIn=[-0.25,359.75],latBoundsIn=[-82.25,90.0], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized',dynamicMask=True) elif experiment == 'AVISO': data, lon, lat = reg2reg(dataIn,resIn=1./4.,lonBoundsIn=[-0.25,359.75],latBoundsIn=[-90.0,90.0], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized',dynamicMask=True) elif experiment == 'ESACCIOC': data, lon, lat = reg2reg(dataIn,resIn=1./24.,lonBoundsIn=[-180.,180.-1e-3],latBoundsIn=[-90.0,90.0-1e-3], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized',dynamicMask=True) elif experiment == 'GlobColour2': data, lon, lat = reg2reg(dataIn,resIn=1./24.,lonBoundsIn=[-180.,180.-1e-3],latBoundsIn=[-90.0,90.0-1e-3], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized',dynamicMask=True) elif experiment == 'SOCCOM': data, lon, lat = reg2reg(dataIn,resIn=1.,lonBoundsIn=[-180.,180.],latBoundsIn=[-90.0,90.0], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized',dynamicMask=True) elif experiment == 'CRUPRE': data, lon, lat = reg2reg(dataIn,resIn=0.5,lonBoundsIn=[-180.,180.],latBoundsIn=[-90.0,90.0], resOut=res,lonBoundsOut=lonBounds,latBoundsOut=latBounds, oceanMaskType='maximized',dynamicMask=True) elif field == 'PRECT' or field == 'TREFHT' or field == 'PSL' or field == 'Z500': if model[:7] == 'noresm1' or model=='noresm2-lm' or model=='noresm2-ml': N = 73 if experiment == 'NCEP' else 96*1-1 dlon = 2.5 else: N = 73 if experiment == 'NCEP' else 96*2-1 dlon = 1.25 dlat = 180/N if len(dataIn.shape) == 3: dataIn = dataIn[:,1:-1,:] else: dataIn = dataIn[1:-1,:] data, lon, lat = reg2reg(dataIn,resIn=[dlon,dlat],lonBoundsIn=[-dlon/2,360-dlon/2], latBoundsIn=[-90+dlat/2,90-dlat/2],resOut=res,lonBoundsOut=lonBounds, latBoundsOut=latBounds,oceanMaskType='maximized') else: data, lon, lat = micom2reg(dataIn,gridFileOcn,res,lonBounds=lonBounds, latBounds=latBounds,oceanMaskType='maximized',isMasked=isMasked) if rec == 0 and hasDepth: depth, lon, lat = micom2reg(depth,gridFileOcn,res,lonBounds=lonBounds, latBounds=latBounds,oceanMaskType='maximized',isMasked=isMasked) # write result if year == yearRange[0] and month == month1: nc.createDimension('time', None) if ensave: dims=['time'] if len(data.shape) > 3: nc.createDimension('l',data.shape[len(data.shape)-4]) dims.extend('l') if len(data.shape) > 2: nc.createDimension('k',data.shape[len(data.shape)-3]) dims.extend('k') if len(data.shape) > 1: nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') if len(data.shape) > 0: nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') else: nc.createDimension('member',memRange[1]-memRange[0]+1) dims=['time','member'] if len(data.shape) > 4: nc.createDimension('l',data.shape[len(data.shape)-4]) dims.extend('l') if len(data.shape) > 3: nc.createDimension('k',data.shape[len(data.shape)-3]) dims.extend('k') if len(data.shape) > 2: nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') if len(data.shape) > 1: nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') ncvar = nc.createVariable(field,'f4',dims,zlib=True) if hasDepth: ncdepth = nc.createVariable('depth','f4',['j','i'],zlib=True) ncdepth[:,:] = depth ncvar[rec,:] = data nc.sync() ncIn.close() nc.close() if clim: yearRange[1] = yearRange1_old def heatContent(field,experiment,yearRange,yearRangeS=[0,0],memRange=[0,0],syear=0,levRange=[0,10000],month1=1,annual=False,clim=False,ensave=True,useBath=False,suffix=''): global bathOcn if field in {'TEMP', 'templvl', 'temperature'}: fieldo='heatContent' elif field in {'SALT', 'salnlvl', 'salinity'}: fieldo='saltContent' else: print('do not recognise field ',field) annualIn = True if clim else False fileIn = outPath(field,experiment,yearRange,syear=syear, memRange=memRange,month1=month1,annual=annualIn,suffix=suffix,ensave=ensave) fileOut = outPath(fieldo,experiment,yearRange,month1=month1,syear=syear,memRange=memRange,levRange=levRange, annual=annual,suffix=suffix,ensave=ensave) if os.path.isfile(fileOut): print(fileOut + ' exists. Do nothing...') return else: print(fileOut + ' does not exist. Creating...') if os.path.isfile(fileIn): ncIn = Dataset(fileIn) print(fileIn) print(fileIn,ncIn.variables[field].shape) else: print(fileIn + ' does not exist. Abort') return # determine input shape shapeIn = ncIn.variables[field].shape ndimsIn = len(shapeIn) isEnsemble = True if ndimsIn >= 5 else False nmem = memRange[1] - memRange[0] + 1 nrec = (yearRange[1] - yearRange[0] + 1)*12 #shapeIn[0] nrecout = int(nrec/12) if annual else nrec if clim: nrec, nrecout = 1, 1 klen = shapeIn[-3] jlen = shapeIn[-2] ilen = shapeIn[-1] print(nrec,nmem,klen,jlen,ilen) # get levels if field in {'templvl', 'temperature', 'salnlvl', 'salinity'}: if field in {'templvl', 'salnlvl'}: fileCoords = dataTmp + '/coords_noresm1_L{:0>2d}.nc'.format(klen) else: fileCoords = dataTmp + '/coords_EN4.nc'.format(klen) fieldCoords = 'depth_bnds' if os.path.isfile(fileCoords): ncCoords = Dataset(fileCoords) depth_bnds = ncCoords.variables[fieldCoords][:,:] elif field in {'TEMP', 'SALT'}: fileCoords = dataTmp + '/coords_cesm1.nc' if os.path.isfile(fileCoords): ncCoords = Dataset(fileCoords) dz = ncCoords.variables['dz'][:]/100. # cm -> m #z_w_top = ncCoords.variables['z_w_top'][:]/100. # cm -> m #z_w_bot = ncCoords.variables['z_w_bot'][:]/100. # cm -> m #depth_bnds = np.concatenate((z_w_top[:,np.newaxis],z_w_bot[:,np.newaxis]),axis=1) depth_bnds = np.zeros((dz.shape[0],2)) for k in range(dz.shape[0]): depth_bnds[k,0] = depth_bnds[max(0,k-1),1] depth_bnds[k,1] = depth_bnds[k,0] + dz[k] ncCoords.close() if useBath: if not 'bathOcn' in globals(): readGridOcn() bath = bathOcn else: bath = 10000. recout_last = -1 for rec in range(nrec): recout = int(rec/12) if annual else rec print(rec,recout) if annual: print(yearRange[0]+recout) else: print(yearRange[0]+int(rec/12),rec%12+1) if rec == 0: nc = Dataset(fileOut, 'w', format='NETCDF4_CLASSIC') nc.createDimension('time', None) if ensave or nmem==1: dims=['time','j','i'] data = np.zeros([nrecout,jlen,ilen]) else: dims=['time','member','j','i'] nc.createDimension('member',nmem) data = np.zeros([nrecout,nmem,jlen,ilen],dtype=np.float32) nc.createDimension('j',jlen) nc.createDimension('i',ilen) ncvar = nc.createVariable(fieldo,'f8',dims,zlib=True) for mem in range(memRange[0],memRange[1]+1): if not ensave or mem == memRange[0]: z = np.zeros([jlen,ilen]) for k in range(klen): bot = np.minimum(min([levRange[1],depth_bnds[k,1]]),bath) top = np.minimum(bot,max([levRange[0],depth_bnds[k,0]])) #print(k,top,bot) if isEnsemble: print(shapeIn,rec,mem-1,k) if ndimsIn == 5: if nmem == shapeIn[0]: fld = ncIn.variables[field][mem-1,rec,k,:,:] else: fld = ncIn.variables[field][rec,mem-1,k,:,:] else: fld = ncIn.variables[field][rec,mem-1,0,k,:,:] else: fld = ncIn.variables[field][rec,k,:,:] #dz = fld * (bot - top) / 100 * 1.030 * 0.41868 # #fldH = fldT * (np.minimum(depth,levRange[1]) - np.minimum(depth,levRange[0]))/100 * 1.030 * 0.41868 # degC -> GJ/m2 dz = fld * (bot - top) * 1026 * 3996 * 1e-9 # #fldH = fldT * (np.minimum(depth,levRange[1]) - np.minimum(depth,levRange[0]))/100 * 1.030 * 0.41868 # degC -> GJ/m2 #dz = eosben07_p_alpha(bot,top,T,S)*1e3 if annual: z = np.where(np.isnan(dz),z,z+dz/12) # correct? else: z = np.where(np.isnan(dz),z,z+dz) z = np.where(z == 0, np.nan, z) if ensave: if annual and recout == recout_last: data[recout,:,:] = data[recout,:,:] + z ncvar[recout,:,:] = data[recout,:,:] + z else: data[recout,:,:] = z ncvar[recout,:,:] = z else: if annual and recout == recout_last: data[recout,mem-1,:,:] = data[recout,mem-1,:,:] + z ncvar[recout,mem-1,:,:] = data[recout,mem-1,:,:] + z else: data[recout,mem-1,:,:] = z ncvar[recout,mem-1,:,:] = z nc.sync() recout_last = recout ncIn.close() nc.close() return data def steric(fieldT,fieldS,experiment,yearRange,yearRangeS=[0,0],memRange=[0,0],syear=0,levRange=[0,10000],month1=1,annual=False,ensave=True,useBath=False,suffix='',sal35psu=False,tem4degC=False): global bathOcn if sal35psu: field='thermosteric' elif tem4degC: field='halosteric' else: field='steric' fileInS = outPath(fieldS,experiment,yearRange,syear=syear, memRange=memRange,month1=month1,annual=False,suffix=suffix,ensave=ensave) fileInT = outPath(fieldT,experiment,yearRange,syear=syear, memRange=memRange,month1=month1,annual=False,suffix=suffix,ensave=ensave) fileOut = outPath(field,experiment,yearRange,month1=month1,syear=syear,memRange=memRange,levRange=levRange, annual=annual,suffix=suffix,ensave=ensave) if os.path.isfile(fileOut): print(fileOut + ' exists. Do nothing...') return else: print(fileOut + ' does not exist. Creating...') if os.path.isfile(fileInS): ncInS = Dataset(fileInS) print(fileInS,ncInS.variables[fieldS].shape) else: print(fileInS + ' does not exist. Abort') return if os.path.isfile(fileInT): ncInT = Dataset(fileInT) print(fileInT) print(fileInT,ncInT.variables[fieldT].shape) else: print(fileInT + ' does not exist. Abort') return # determine input shape shapeIn = ncInT.variables[fieldT].shape ndimsIn = len(shapeIn) isEnsemble = True if ndimsIn >= 5 else False nmem = memRange[1] - memRange[0] + 1 nrec = shapeIn[0] nrecout = int(nrec/12) if annual else nrec klen = shapeIn[-3] jlen = shapeIn[-2] ilen = shapeIn[-1] print(nrec,nmem,klen,jlen,ilen) # get levels if fieldT == 'templvl' or fieldT == 'temperature': if fieldT == 'templvl': fileCoords = dataTmp + '/coords_noresm1_L{:0>2d}.nc'.format(klen) else: fileCoords = dataTmp + '/coords_EN4.nc'.format(klen) fieldCoords = 'depth_bnds' if os.path.isfile(fileCoords): ncCoords = Dataset(fileCoords) depth_bnds = ncCoords.variables[fieldCoords][:,:] elif fieldT == 'TEMP': fileCoords = dataTmp + '/coords_cesm1.nc' if os.path.isfile(fileCoords): ncCoords = Dataset(fileCoords) z_w_top = ncCoords.variables['z_w_top'][:]/100. z_w_bot = ncCoords.variables['z_w_bot'][:]/100. depth_bnds = np.concatenate((z_w_top[:,np.newaxis],z_w_bot[:,np.newaxis]),axis=1) ncCoords.close() if useBath: if not 'bathOcn' in globals(): readGridOcn() bath = bathOcn else: bath = 10000. recout_last = -1 for rec in range(nrec): recout = int(rec/12) if annual else rec print(rec,recout) if annual: print(yearRange[0]+recout) else: print(yearRange[0]+int(rec/12),rec%12+1) if rec == 0: nc = Dataset(fileOut, 'w', format='NETCDF4_CLASSIC') nc.createDimension('time', None) if ensave or nmem==1: dims=['time','j','i'] data = np.zeros([nrecout,jlen,ilen]) else: dims=['time','member','j','i'] nc.createDimension('member',nmem) data = np.zeros([nrecout,nmem,jlen,ilen],dtype=single) nc.createDimension('j',jlen) nc.createDimension('i',ilen) ncvar = nc.createVariable(field,'f8',dims,zlib=True) for mem in range(memRange[0],memRange[1]+1): if not ensave or mem == memRange[0]: z = np.zeros([jlen,ilen]) for k in range(klen): bot = np.minimum(min([levRange[1],depth_bnds[k,1]]),bath) top = np.minimum(bot,max([levRange[0],depth_bnds[k,0]])) #print(k,top,bot) if isEnsemble: if not sal35psu: S = ncInS.variables[fieldS][rec,mem-1,k,:,:] if not tem4degC: T = ncInT.variables[fieldT][rec,mem-1,k,:,:] else: if not sal35psu: S = ncInS.variables[fieldS][rec,k,:,:] if not tem4degC: T = ncInT.variables[fieldT][rec,k,:,:] if sal35psu: S = T * 0 + 35 if tem4degC: T = S * 0 + 4 dz = eosben07_p_alpha(bot,top,T,S)*1e3 if annual: z = np.where(np.isnan(dz),z,z+dz*(1/12)) else: z = np.where(np.isnan(dz),z,z+dz) z = np.where(z == 0, np.nan, z) if ensave: if annual and recout == recout_last: data[recout,:,:] = data[recout,:,:] + z ncvar[recout,:,:] = data[recout,:,:] + z else: data[recout,:,:] = z ncvar[recout,:,:] = z else: if annual and recout == recout_last: data[recout,mem-1,:,:] = data[recout,mem-1,:,:] + z ncvar[recout,mem-1,:,:] = data[recout,mem-1,:,:] + z else: data[recout,mem-1,:,:] = z ncvar[recout,mem-1,:,:] = z nc.sync() recout_last = recout ncInS.close() ncInT.close() nc.close() return data def writeEns(field,experiment,yearRange,memRange,month1=1,syear=0,annual=True,useTestInput=False,zlib=False): filePath = outPath(field,experiment,yearRange,month1=month1,syear=syear,memRange=memRange,annual=annual,ensave=False) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) monthLast = 1 if annual else 12 monthFirst = month1 if year == yearRange[0] else 1 for month in range(monthFirst,monthLast+1): rec = year-yearRange[0] if annual else (year-yearRange[0])*12 + month - 1 - month1 + 1 for member in range(memRange[0],memRange[1]+1): if annual: if field == 'fgco2': data = readAnn('co2fxd',experiment,year,member,syear,useTestInput) - readAnn('co2fxu',experiment,year,member,syear,useTestInput) elif model[:7] == 'noresm2' and field == 'PRECT': data = readAnn('PRECC',experiment,year,member,syear,useTestInput) + readAnn('PRECL',experiment,year,member,syear,useTestInput) else: data = readAnn(field,experiment,year,member,syear,useTestInput) else: if field == 'fgco2': data = readAnn('co2fxd',experiment,year,month,member,syear,useTestInput) - readAnn('co2fxu',experiment,year,month,member,syear,useTestInput) elif model[:7] == 'noresm2' and field == 'PRECT': data = readAnn('PRECT',experiment,year,month,member,syear,useTestInput) + readAnn('PRECL',experiment,year,month,member,syear,useTestInput) else: data = readMon(field,experiment,year,month,member,syear,useTestInput) if member == memRange[0] and year == yearRange[0] and month == month1: nc.createDimension('time', None) nc.createDimension('member', memRange[1]-memRange[0]+1) dims=['time','member'] if len(data.shape) > 3: nc.createDimension('l',data.shape[len(data.shape)-4]) dims.extend('l') if len(data.shape) > 2: nc.createDimension('k',data.shape[len(data.shape)-3]) dims.extend('k') if len(data.shape) > 1: nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') if len(data.shape) > 0: nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') ncvar = nc.createVariable(field,'f4',dims,zlib=zlib) ncvar[rec,member-1,:] = data nc.close() def writeEnsAve(field,experiment,yearRange,memRange,month1=1,syear=0,annual=True,useTestInput=False): filePath = outPath(field,experiment,yearRange,month1=month1,syear=syear,memRange=memRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') w = 1. / (memRange[1] - memRange[0] + 1) nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) monthLast = 1 if annual else 12 monthFirst = month1 if year == yearRange[0] else 1 for month in range(monthFirst,monthLast+1): for member in range(memRange[0],memRange[1]+1): if annual: if field == 'fgco2': data = readAnn('co2fxd',experiment,year,member,syear,useTestInput) - readAnn('co2fxu',experiment,year,member,syear,useTestInput) else: data = readAnn(field,experiment,year,member,syear,useTestInput) else: if field == 'fgco2': data = readAnn('co2fxd',experiment,year,month,member,syear,useTestInput) - readAnn('co2fxu',experiment,year,month,member,syear,useTestInput) else: data = readMon(field,experiment,year,month,member,syear,useTestInput) out = w * data if member == memRange[0] else out + w * data if member == memRange[0] and year == yearRange[0] and month == month1: time = nc.createDimension('time', None) dims=['time'] if len(data.shape) > 4: nc.createDimension('l',data.shape[len(data.shape)-4]) dims.extend('l') if len(data.shape) > 3: nc.createDimension('k',data.shape[len(data.shape)-3]) dims.extend('k') if len(data.shape) > 2: nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') if len(data.shape) > 1: nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') ncvar = nc.createVariable(field,'f4',dims,zlib=True) rec = year-yearRange[0] if annual else (year-yearRange[0])*12 + month - 1 - month1 + 1 ncvar[rec,:] = out nc.close() def concatERSST(yearRange,annual=True): field = 'sst' filePath = outPath(field,'ERSSTv5',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): ncIn = Dataset(fileERSST(year=year,month=month)) data = ncIn.variables[field][0,0,:,:] if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = (data - 273.15) * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + (data - 273.15) * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data - 273.15 ncIn.close() nc.close() def concatHadSLP(yearRange,annual=True): filePathIn = dataObs + '/regridded/PRES/hadslp2r.asc' filePath = outPath('PSL','HadSLP2r',yearRange,annual=annual,suffix='_5x5') if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') # read data file = open(filePathIn,'r') years = [1850,2019] idm = 72 jdm = 37 kdm = (years[1]-years[0]+1)*12 psl = np.zeros((kdm,jdm,idm)) for line in file: if len(line) < 100: year = int(line[0:7]) month = int(line[7:14]) if month == 1: print(year) j = 0 else: k = (year-years[0])*12 + month - 1 for i in range(idm): psl[k,j,i] = float(line[i*8:(i+1)*8]) j = j + 1 file.close() # regrid psl, lon, lat = reg2reg(np.fliplr(psl[:,1:36,:]),resIn=5., lonBoundsIn=[-182.5,177.5],latBoundsIn=[-87.5,87.5], resOut=5.,lonBoundsOut=[0,360],latBoundsOut=[-90,90], oceanMaskType='maximized',dynamicMask=True) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): data = psl[(year-1850)*12+month-1] if year == yearRange[0] and month == 1: lon = lon[0,:] lat = lat[:,0] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable('PSL','f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data nc.close() def concatNCEPSLP(yearRange,annual=True): field = 'slp' filePathIn = dataObs + '/regridded/PRES/ncep.slp.mon.mean.nc' filePath = outPath('PSL','NCEP',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') ncIn = Dataset(filePathIn) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): data = np.flipud(ncIn.variables[field][(year-1948)*12+month-1,:,:]) if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = np.flipud(ncIn.variables['lat'][:]) nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable('PSL','f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def concatHadISSTice(yearRange,annual=True): field = 'aice' filePathIn = dataObs + '/regridded/HadISST/HadISST_ice.nc' filePath = outPath(field,'HadISST',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') filePathIn = dataObs + '/regridded/HadISST/HadISST_ice.nc' ncIn = Dataset(filePathIn) lon = ncIn.variables['longitude'][:] lat = np.flip(ncIn.variables['latitude'][:]) nc.createDimension('time', None) dims=['time'] nc.createDimension('j',len(lat)) dims.extend('j') nc.createDimension('i',len(lon)) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): data = np.flip(ncIn.variables['sic'][(year-1870)*12+month-1,:,:],axis=0) if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def concatHadCRUT(yearRange,annual=True): field = 'temperature_anomaly' filePathIn = dataObs + '/regridded/SAT/HadCRUT.4.6.0.0.median.nc' filePath = outPath('TREFHT','HadCRUT',yearRange,annual=annual,suffix='_5x5') if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') ncIn = Dataset(filePathIn) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): data = ncIn.variables[field][(year-1850)*12+month-1,:,:] data = np.ma.concatenate((data[:,36:],data[:,:36]),axis=1) if year == yearRange[0] and month == 1: lon = ncIn.variables['longitude'][:] lon = np.concatenate((lon[36:],lon[:36]+360)) lat = ncIn.variables['latitude'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable('TREFHT','f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def concatERA5T2M(yearRange,annual=True): field = 't2m' filePathIn = dataObs + '/regridded/SAT/ERA5.t2m.197901-202009.nc' filePath = outPath('TREFHT','ERA5',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') ncIn = Dataset(filePathIn) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): data = np.flip(ncIn.variables[field][(year-1979)*12+month-1,0,:,:],axis=0) if year == yearRange[0] and month == 1: lon = ncIn.variables['longitude'][:] lat = np.flip(ncIn.variables['latitude'][:]) nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable('TREFHT','f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def concatERA5PR(yearRange,annual=True): field = 'tp' filePathIn = dataObs + '/regridded/PRECIP/ERA5/ERA5.tp.1940-2023.nc' filePath = outPath('PRECT','ERA5',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') ncIn = Dataset(filePathIn) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): data = np.flip(ncIn.variables[field][(year-1940)*12+month-1,0,:,:],axis=0) print(data.shape) if year == yearRange[0] and month == 1: lon = ncIn.variables['longitude'][:] lat = np.flip(ncIn.variables['latitude'][:]) nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable('PRECT','f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def concatERA5Z500(yearRange,annual=True): field = 'z' filePath = outPath('Z500','ERA5',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) if year < 1979: filePathIn = dataObs + '/regridded/Z500/ERA5.z.1950-1978.nc' year1 = 1950 else: filePathIn = dataObs + '/regridded/Z500/ERA5.z.1979-2019.nc' year1 = 1979 ncIn = Dataset(filePathIn) for month in range(1,12+1): data = np.flip(ncIn.variables[field][(year-year1)*12+month-1,:,:],axis=0) if year == yearRange[0] and month == 1: lon = ncIn.variables['longitude'][:] lat = np.flip(ncIn.variables['latitude'][:]) nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable('Z500','f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def concatCRUPRE(yearRange,annual=True): field = 'pre' filePathIn = dataObs + '/regridded/PRECIP/cru_ts4.03.1901.2018.pre.dat.nc4' filePath = outPath('PRECT','CRUPRE',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') ncIn = Dataset(filePathIn) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): data = ncIn.variables[field][(year-1901)*12+month-1,:,:] if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable('PRECT','f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def readERA5day(field,year,month,day,readCoords=False): rec = 0 for m in range(month-0): rec = rec + daysInMonth[m] rec = rec + day - 1 print(year,month,day,rec) nc = Dataset(fileERA5day(year=year)) fld = nc.variables[field][rec,0,:,:] if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() if readCoords: return fld, lon, lat else: return fld def readERA5sixhourly(field,year,month,day,record,readCoords=False): rec = (day - 1)*4 + record FPATH=fileERA5sixhourly(field,year,month) if field == 't850': field = 't' print(field,FPATH,year,month,day,rec) nc = Dataset(FPATH) fld = nc.variables[field][rec,:,:] if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() if readCoords: return fld, lon, lat else: return fld def readHADISST_sst(year,month,readCoords=False): fileHADISST = '/nird/projects/NS9039K/shared/obsdata/regridded/HadISST/HadISST_sst.nc' rec = (year - 1870)*12 + month - 1 nc = Dataset(fileHADISST) fld = nc.variables['sst'][rec,:,:] fld = np.where(np.abs(np.array(fld)) < 1001, np.where(np.array(fld)>-999.,np.array(fld),-1.8), np.nan) if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() if readCoords: return fld, lon, lat else: return fld def readOSTIA_sst(year,month,readCoords=False): fileOSTIA = '/nird/projects/NS9039K/shared/obsdata/regridded/OSTIA/OSTIA_198110-202508.nc' rec = (year - 1981)*12 + month - 10 print(year,month,rec) nc = Dataset(fileOSTIA) fld = nc.variables['sst'][rec,:,:] msk = nc.variables['sst'][rec,:,:] fld = np.where(msk==2,np.nan,fld) if readCoords: lon = nc.variables['lon'][:] lat = nc.variables['lat'][:] nc.close() if readCoords: return fld, lon, lat else: return fld def readOISST(field,year,month,day,readCoords=False,hires=True): nc = Dataset(fileOISST(year=year,month=month,day=day,hires=hires)) fld = nc.variables[field][0,0,:,:] if hires else nc.variables[field][0,:,:] fld = np.where(np.abs(np.array(fld)) < 50, np.array(fld), np.nan) if readCoords: lon = nc.variables['lon'][:] lat = nc.variables['lat'][:] nc.close() if readCoords: return fld, lon, lat else: return fld def readMonOISST(field,year,month,readCoords=False,hires=True): daysInMonth = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) for day in range(1,daysInMonth[month-1]+1): w = 1. / daysInMonth[month-1] if readCoords and day == 1: fldday, lon, lat = readOISST(field,year,month,day,readCoords=readCoords,hires=hires) else: fldday = readOISST(field,year,month,day,readCoords=False,hires=hires) if field == 'ice': sst = readOISST('sst',year,month,day,readCoords=False,hires=hires) fldday = np.where(fldday<1e3, fldday, np.where(np.isnan(sst),np.nan,0)) if day == 1: fld = fldday * w else: fld = fld + fldday * w if readCoords: return fld, lon, lat else: return fld def readMonARMOR3D(year,month,hiRes=False,readCoords=False): nc = Dataset(fileARMOR3D(year=year,month=month,hiRes=hiRes)) fld = nc.variables['zo'][0,0,:,:] fld = np.where(np.array(fld) < 10, np.array(fld), np.nan) if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() if readCoords: return fld, lon, lat else: return fld def readMonAVISO(year,month,readCoords=False,anomaly=True): nc = Dataset(fileAVISO(year=year,month=month)) fld = nc.variables['sla'][0,:,:] fld = np.where(np.abs(np.array(fld)) < 10, np.array(fld), np.nan) if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() if anomaly: eraRange = [1993, 2022] gm = globalMeanAVISO(eraRange,annual=False) fld = fld - gm[year-eraRange[0],month-1] if readCoords: return fld, lon, lat else: return fld def climMonARMOR3D(year1,yearn,hiRes=False,readCoords=False): # check if output file already exists field = 'zo' if hiRes: filePath = outPath(field,'ARMOR3D025',[year1,yearn],suffix='_clim',annual=False) else: filePath = outPath(field,'ARMOR3D',[year1,yearn],suffix='_clim',annual=False) if os.path.isfile(filePath): print(filePath + ' exists. Read climatology from file...') nc = Dataset(filePath,'r') clim = nc.variables[field][:,:,:] if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() if readCoords: return clim, lon, lat else: return clim else: # read coordinate variables and determine dimensions nc = Dataset(fileARMOR3D(year=year1,month=1,hiRes=hiRes),'r') if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() idm = len(lon) jdm = len(lat) clim = np.zeros([12,jdm,idm]) mask = np.zeros([12,jdm,idm]) for year in range(year1,yearn+1): print(year) for m in range(12): fld = readMonARMOR3D(year,m+1,hiRes=hiRes,readCoords=False) clim[m,:,:] = np.where(np.isnan(fld), clim[m,:,:], clim[m,:,:] + fld) mask[m,:,:] = np.where(np.isnan(fld), mask[m,:,:], mask[m,:,:] + 1) for m in range(12): clim[m,:,:] = np.where(mask[m,:,:] == 0, np.nan, clim[m,:,:]/mask[m,:,:]) mask = mask / (yearn - year1 + 1) # write result nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('time', 12) nc.createDimension('longitude',idm) nc.createDimension('latitude',jdm) nclon = nc.createVariable('longitude','f8',('longitude')) nclat = nc.createVariable('latitude','f8',('latitude')) ncfld = nc.createVariable(field,'f4',('time','latitude','longitude'),zlib=True) ncmsk = nc.createVariable('mask','f4',('time','latitude','longitude'),zlib=True) nclon[:] = lon nclat[:] = lat ncfld[:,:,:] = clim ncmsk[:,:,:] = mask nc.close() if readCoords: return clim, lon, lat else: return clim def climMonAVISO(year1,yearn,readCoords=False): # check if output file already exists field = 'sla' filePath = outPath(field,'AVISO',[year1,yearn],suffix='_clim',annual=False) if os.path.isfile(filePath): print(filePath + ' exists. Read climatology from file...') nc = Dataset(filePath,'r') clim = nc.variables[field][:,:,:] if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() if readCoords: return clim, lon, lat else: return clim else: # read coordinate variables and determine dimensions nc = Dataset(fileAVISO(year=year1,month=1),'r') if readCoords: lon = nc.variables['longitude'][:] lat = nc.variables['latitude'][:] nc.close() idm = len(lon) jdm = len(lat) clim = np.zeros([12,jdm,idm]) mask = np.zeros([12,jdm,idm]) for year in range(year1,yearn+1): print(year) for m in range(12): fld = readMonAVISO(year,m+1,readCoords=False) clim[m,:,:] = np.where(np.isnan(fld), clim[m,:,:], clim[m,:,:] + fld) mask[m,:,:] = np.where(np.isnan(fld), mask[m,:,:], mask[m,:,:] + 1) for m in range(12): clim[m,:,:] = np.where(mask[m,:,:] == 0, np.nan, clim[m,:,:]/mask[m,:,:]) mask = mask / (yearn - year1 + 1) # write result nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('time', 12) nc.createDimension('longitude',idm) nc.createDimension('latitude',jdm) nclon = nc.createVariable('longitude','f8',('longitude')) nclat = nc.createVariable('latitude','f8',('latitude')) ncfld = nc.createVariable(field,'f4',('time','latitude','longitude'),zlib=True) ncmsk = nc.createVariable('mask','f4',('time','latitude','longitude'),zlib=True) nclon[:] = lon nclat[:] = lat ncfld[:,:,:] = clim ncmsk[:,:,:] = mask nc.close() if readCoords: return clim, lon, lat else: return clim def concatAVISO(yearRange,annual=True): field = 'sla' filePath = outPath(field,'AVISO',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): ncIn = Dataset(fileAVISO(year=year,month=month)) data = ncIn.variables[field][0,:,:] if year == yearRange[0] and month == 1: lon = ncIn.variables['longitude'][:] lat = ncIn.variables['latitude'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def concatARMOR3D(yearRange,annual=True,hiRes=False): field = 'zo' if hiRes: filePath = outPath(field,'ARMOR3D025',yearRange,annual=annual) else: filePath = outPath(field,'ARMOR3D',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): ncIn = Dataset(fileARMOR3D(year=year,month=month,hiRes=hiRes)) data = ncIn.variables[field][0,0,:,:] if year == yearRange[0] and month == 1: if hiRes: lon = ncIn.variables['longitude'][:] lat = ncIn.variables['latitude'][:] else: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = data * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + data * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() nc.close() def concatESACCIOC(yearRange,annual=True): field = 'chlor_a' filePath = outPath(field,'ESACCIOC',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): ncIn = Dataset(fileESACCIOC(year=year,month=month)) data = ncIn.variables[field][0,:,:] if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: msk=np.ma.getmask(data) if month == 1: outmsk=np.where(msk,0,weightsAnnMean[month-1]) outann=np.where(msk,0,data*weightsAnnMean[month-1]) else: outmsk=np.where(msk,outmsk,outmsk+weightsAnnMean[month-1]) outann=np.where(msk,outann,outann+data*weightsAnnMean[month-1]) else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() if annual: rec = year-yearRange[0] ncvar[rec,:] = np.where(outmsk == 0, np.nan, outann/outmsk) nc.close() def concatGlobColour(yearRange,annual=True): field = 'PP' filePath = outPath(field,'GlobColour',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): ncIn = Dataset(fileGlobColour(year=year,month=month)) data = ncIn.variables[field][0,:,:] if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: msk=np.ma.getmask(data) if month == 1: outmsk=np.where(msk,0,weightsAnnMean[month-1]) outann=np.where(msk,0,data*weightsAnnMean[month-1]) else: outmsk=np.where(msk,outmsk,outmsk+weightsAnnMean[month-1]) outann=np.where(msk,outann,outann+data*weightsAnnMean[month-1]) else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() if annual: rec = year-yearRange[0] ncvar[rec,:] = np.where(outmsk == 0, np.nan, outann/outmsk) nc.close() def concatGlobColour2(field,yearRange,annual=True,fill=False): if fill: filePath = outPath(field,'GlobColour2fill',yearRange,annual=annual) else: filePath = outPath(field,'GlobColour2',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') if fill: filePathClim = outPath(field,'GlobColour2clim',yearRange,annual=annual) if not(os.path.isfile(filePathClim)): print(filePathClim + ' does not exist. Creating...') nc = Dataset(filePathClim, 'w', format='NETCDF4_CLASSIC') for month in range(1,12+1): for year in range(yearRange[0],yearRange[1]+1): print('{:0>4d}-{:0>2d}'.format(year,month)) ncIn = Dataset(fileGlobColour2(field,year=year,month=month)) data = np.flip(ncIn.variables[field][0,:,:],axis=0) if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] nc.createDimension('time', 12) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat ncIn.close() mskins=np.ma.getmask(data) if year == yearRange[0]: msk=np.where(mskins,0,1) clim=np.where(mskins,0,data) else: msk=np.where(mskins,msk,msk+1) clim=np.where(mskins,clim,clim+data) ncvar[month-1,:] = np.where(msk == 0, np.nan, clim/msk) nc.close() print('Read monthly climatology from file ' + filePathClim) nc = Dataset(filePathClim,'r') clim = nc.variables[field][:,:,:] clim = np.where(np.isnan(clim),0,clim) nc.close() nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): ncIn = Dataset(fileGlobColour2(field,year=year,month=month)) data = np.flip(ncIn.variables[field][0,:,:],axis=0) if fill: data = np.where(np.ma.getmask(data),clim[month-1,:,:],data) if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: msk=np.ma.getmask(data) if month == 1: outmsk=np.where(msk,0,weightsAnnMean[month-1]) outann=np.where(msk,0,data*weightsAnnMean[month-1]) else: outmsk=np.where(msk,outmsk,outmsk+weightsAnnMean[month-1]) outann=np.where(msk,outann,outann+data*weightsAnnMean[month-1]) else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() if annual: rec = year-yearRange[0] ncvar[rec,:] = np.where(outmsk == 0, np.nan, outann/outmsk) nc.close() def concatMILAGPV(yearRange,annual=True): field = 'MLD' filePath = outPath(field,'MILAGPV',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): ncIn = Dataset(fileMILAGPV(year=year,month=month)) data = ncIn.variables[field][0,:,:] if year == yearRange[0] and month == 1: lon = ncIn.variables['LONGITUDE'][:] lat = ncIn.variables['LATITUDE'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: msk=np.ma.getmask(data) if month == 1: outmsk=np.where(msk,0,weightsAnnMean[month-1]) outann=np.where(msk,0,data*weightsAnnMean[month-1]) else: outmsk=np.where(msk,outmsk,outmsk+weightsAnnMean[month-1]) outann=np.where(msk,outann,outann+data*weightsAnnMean[month-1]) else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data ncIn.close() if annual: rec = year-yearRange[0] ncvar[rec,:] = np.where(outmsk == 0, np.nan, outann/outmsk) nc.close() def concatSOCCOM(yearRange,field = 'spco2',annual=True): # field = 'fgco2' filePathIn = dataObs + '/regridded/fgco2/MPI-SOM_FFN_SOCCOMv2018.nc' filePath = outPath(field,'SOCCOM',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') ncIn = Dataset(filePathIn) for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): data = ncIn.variables[field][(year-1982)*12+month-1,:,:] data = np.where(np.abs(data) > 1e10, np.nan, data) if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] nc.createDimension('time', None) dims=['time'] nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat if annual: msk=np.ma.getmask(data) if month == 1: outmsk=np.where(msk,0,weightsAnnMean[month-1]) outann=np.where(msk,0,data*weightsAnnMean[month-1]) else: outmsk=np.where(msk,outmsk,outmsk+weightsAnnMean[month-1]) outann=np.where(msk,outann,outann+data*weightsAnnMean[month-1]) else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data if annual: rec = year-yearRange[0] ncvar[rec,:] = np.where(outmsk == 0, np.nan, outann/outmsk) ncIn.close() nc.close() def fileERSST(year,month): PREFIX = dataObs + '/regridded/ERSSTv5/ersst.v5.' SUFFIX = '.nc' return PREFIX + '{:0>4d}{:0>2d}'.format(year,month) + SUFFIX def fileERA5day(year): PREFIX = '/scratch/ingo/data/ERA5_daily/era5_daily_' SUFFIX = '.nc' return PREFIX + '{:0>4d}'.format(year) + SUFFIX def fileERA5sixhourly(field,year,month): if field == 't850': PREFIX = '/scratch/ingo/data/ERA5_sixhourly/t850_era5_sixhourly_' else: PREFIX = '/scratch/ingo/data/ERA5_sixhourly/era5_sixhourly_' SUFFIX = '.nc' return PREFIX + '{:0>4d}{:0>2d}'.format(year,month) + SUFFIX def fileOISST(year,month,day,hires=True): #if year == 2023 and month >= 11: # fill missing days # month = 11 # day = 1 if hires: PREFIX = dataObs + '/satellite/oisst_hires/oisst-avhrr-v02r01.' SUFFIX = '.nc' file = PREFIX + '{:0>4d}{:0>2d}{:0>2d}'.format(year,month,day) + SUFFIX else: PREFIX = dataObs + '/satellite/oisst_1x1/' SUFFIX = '.nc' file = PREFIX + '{:0>4d}_{:0>2d}_{:0>2d}'.format(year,month,day) + SUFFIX return file def fileAVISO(year,month): PREFIX = dataObs + '/regridded/AVISO/' SUFFIX = '.nc' return PREFIX + '{:0>4d}_{:0>2d}'.format(year,month) + SUFFIX def fileARMOR3D(year,month,hiRes=False): if hiRes: PREFIX = dataObs + '/regridded/ARMOR3D/025x025/' SUFFIX = '.nc' return PREFIX + 'armor3d_zo_{:0>4d}{:0>2d}'.format(year,month) + SUFFIX else: PREFIX = dataObs + '/regridded/ARMOR3D/' SUFFIX = '.nc' return PREFIX + '{:0>4d}_{:0>2d}'.format(year,month) + SUFFIX def fileESACCIOC(year,month): PREFIX = dataObs + '/satellite/ESACCI-OC/ESACCI-OC-L3S-CHLOR_A-MERGED-1M_MONTHLY_4km_GEO_PML_OCx-' #PREFIX = '/Data/skd/work/ibe062/ESACCI-OC/ESACCI-OC-L3S-CHLOR_A-MERGED-1M_MONTHLY_4km_GEO_PML_OCx-' SUFFIX = '-fv4.2.nc' return PREFIX + '{:0>4d}{:0>2d}'.format(year,month) + SUFFIX def fileGlobColourOld(year,month): PREFIX = dataObs + '/satellite/GlobColour/oc-glo-bio-multi-l4-pp_4km_monthly-rep_' #PREFIX = dataTmp + '/Data/skd/work/ibe062/GlobColour/oc-glo-bio-multi-l4-pp_4km_monthly-rep_' SUFFIX = '.nc' return PREFIX + '{:0>4d}-{:0>2d}'.format(year,month) + SUFFIX def fileGlobColour2(field,year,month): #/projects/NS9039K/data/external/observation/CMEMS/GlobColour/original #20090801-20090831_cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M.nc def isLeapYear(y): return y%4 == 0 and y%100 != 0 or y%400 == 0 def daysInMonth(y,m): daysInMonthNoLeap = [31,28,31,30,31,30,31,31,30,31,30,31] daysInMonthLeap = [31,29,31,30,31,30,31,31,30,31,30,31] if isLeapYear(int(y)): days = daysInMonthLeap[int(m)-1] else: days = daysInMonthNoLeap[int(m)-1] return days PREFIX = '/projects/NS9039K/data/external/observation/CMEMS/GlobColour/original' if field == 'PP': SUFFIX = '_cmems_obs-oc_glo_bgc-pp_my_l4-multi-4km_P1M.nc' elif field == 'CHL': SUFFIX = '_cmems_obs-oc_glo_bgc-plankton_my_l4-multi-4km_P1M.nc' else: print('No ' + field + ' in ' + PREFIX) quit() return PREFIX + '/{:0>4d}{:0>2d}01-{:0>4d}{:0>2d}{:0>2d}'.format(year,month,year,month,daysInMonth(year,month)) + SUFFIX def fileMILAGPV(year,month): PREFIX = dataObs + '/regridded/MILAGPV/mila_mldmld_' SUFFIX = '_010.nc' return PREFIX + '{:0>4d}{:0>2d}00'.format(year,month) + SUFFIX def concatEn4(field,yearRange,annual=True): filePath = outPath(field,'EN4',yearRange,annual=annual) if os.path.isfile(filePath): print(filePath + ' exists. Do nothing...') return else: print(filePath + ' does not exist. Creating...') nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') for year in range(yearRange[0],yearRange[1]+1): print('year = ' + '{:0>2d}'.format(year)) for month in range(1,12+1): ncIn = Dataset(fileEn4(year=year,month=month)) data = ncIn.variables[field][:] if year == yearRange[0] and month == 1: lon = ncIn.variables['lon'][:] lat = ncIn.variables['lat'][:] depth = ncIn.variables['depth'][:] depth_bnds = ncIn.variables['depth_bnds'][:] nc.createDimension('time', None) dims=['time'] if len(data.shape) > 4: nc.createDimension('l',data.shape[len(data.shape)-4]) dims.extend('l') if len(data.shape) > 3: nc.createDimension('k',data.shape[len(data.shape)-3]) dims.extend('k') if len(data.shape) > 2: nc.createDimension('j',data.shape[len(data.shape)-2]) dims.extend('j') if len(data.shape) > 1: nc.createDimension('i',data.shape[len(data.shape)-1]) dims.extend('i') nc.createDimension('bnds', 2) nclon = nc.createVariable('lon','f8',(dims[-1])) nclat = nc.createVariable('lat','f8',(dims[-2])) ncdepth = nc.createVariable('depth','f8',(dims[-3])) ncdepth_bnds = nc.createVariable('depth_bnds','f8',(dims[-3],'bnds')) ncvar = nc.createVariable(field,'f4',dims,zlib=True) nclon[:] = lon nclat[:] = lat ncdepth[:] = depth ncdepth_bnds[:] = depth_bnds if annual: rec = year-yearRange[0] if month == 1: ncvar[rec,:] = (data - 273.15) * weightsAnnMean[0] else: ncvar[rec,:] = ncvar[rec,:] + (data - 273.15) * weightsAnnMean[month-1] else: rec = (year-yearRange[0])*12 + month - 1 ncvar[rec,:] = data - 273.15 ncIn.close() nc.close() def fileEn4(year,month): #PREFIX = dataObs + '/regridded/EN.4.2.1/EN.4.2.1.f.analysis.g10.' PREFIX = dataObs + '/regridded/EN.4.2.2/EN.4.2.2.f.analysis.g10.' SUFFIX = '.nc' return PREFIX + '{:0>4d}'.format(year) + '{:0>2d}'.format(month) + SUFFIX def readAnn(field,experiment,year,member=0,syear=0,useTestInput=False): """ Compute annual avarages from monthly output. """ for month in range(12): data = readMon(field,experiment,year,month+1,member,syear,useTestInput) out = data * weightsAnnMean[0] if month == 0 else out + data * weightsAnnMean[month] return out def readMon(field,experiment,year,month,member=0,syear=0,useTestInput=False): """ Read field from monthly output. """ nc = Dataset(filePathMon(field,experiment,year,month,member,syear,useTestInput)) out = nc.variables[field][:].squeeze() nc = nc.close() if field in ['tempinslvl','salninslvl']: monthLast = ((month - 1) % 12) + 1 yearLast = year -1 if month == 1 else year filePathLast=filePathMon(field,experiment,yearLast,monthLast,member,syear,useTestInput) if os.path.isfile(filePathLast): nc = Dataset(filePathLast) outLast = nc.variables[field][:] nc = nc.close() else: outLast = out out = (out - outLast) / daysInMonth[month-1] return out def NAO(experiment,yearRange,memRange,syear=0,useTestInput=False): """ Read pressure from Reyjavik (64.128288, -21.827774) and Lisbon (38.736946, -9.142685) from monthly output and compute DJFM means. """ global weightsR, weightsL filePath = outPath('NAO',experiment,yearRange,memRange=memRange,syear=syear) if os.path.isfile(filePath): print(filePath + ' exists. Read raw NAO from file.') nc = Dataset(filePath) slpR = nc.variables['slpR'][:] slpL = nc.variables['slpL'][:] nc = nc.close() return slpR, slpL else: print(filePath + ' does not exist. Generate raw NAO.') lonR, latR = 360-21.827774, 64.128288 lonL, latL = 360-9.142685, 38.736946 delta = 2.5 slpR = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) slpL = np.zeros([yearRange[1]-yearRange[0]+1,memRange[1]-memRange[0]+1]) for year in range(yearRange[0],yearRange[1]+1): for mem in range(memRange[0],memRange[1]+1): for month in (np.arange(4)+11) % 12 + 1: if month == 12: nc = Dataset(filePathMon('PSL',experiment,year-1,month,mem,syear,useTestInput)) else: nc = Dataset(filePathMon('PSL',experiment,year,month,mem,syear,useTestInput)) if not 'weightsR' in globals(): lon = nc.variables['lon'][:] lat = nc.variables['lat'][:] weightsR, weightsL = np.zeros([lat.size,lon.size]), np.zeros([lat.size,lon.size]) for j in range(lat.size): for i in range(lon.size): weightsR[j,i] = np.cos(lat[j]*np.pi/180.) if lon[i]>=lonR-delta and lon[i]<=lonR+delta and lat[j]>=latR-delta and lat[j]<=latR+delta else 0 weightsL[j,i] = np.cos(lat[j]*np.pi/180.) if lon[i]>=lonL-delta and lon[i]<=lonL+delta and lat[j]>=latL-delta and lat[j]<=latL+delta else 0 weightsR, weightsL = weightsR/np.sum(weightsR), weightsL/np.sum(weightsL) out = nc.variables['PSL'][:] nc = nc.close() slpR[year-yearRange[0],mem-memRange[0]] = slpR[year-yearRange[0],mem-memRange[0]] + np.sum(out * weightsR) * daysInMonth[month-1]/121. slpL[year-yearRange[0],mem-memRange[0]] = slpL[year-yearRange[0],mem-memRange[0]] + np.sum(out * weightsL) * daysInMonth[month-1]/121. nc = Dataset(filePath, 'w', format='NETCDF4_CLASSIC') nc.createDimension('member',memRange[1]-memRange[0]+1) nc.createDimension('year',yearRange[1]-yearRange[0]+1) ncslpR = nc.createVariable('slpR','f4',['year','member']) ncslpL = nc.createVariable('slpL','f4',['year','member']) ncslpR[:], ncslpL[:] = slpR, slpL nc = nc.close() return slpR, slpL def readGridOcn(): """Read ocean grid information.""" global lonOcn, latOcn, areaOcn, maskOcn, bathOcn, areaOcnNorm, depth_bnds if model == 'noresm1-me': gridFileOcn = dataTmp + '/' + 'grid_ocn_gx1v6.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_1.9x2.5_gx1v6.nc' elif model == 'noresm2-lm': gridFileOcn = dataTmp + '/' + 'grid_ocn_tnx1v4.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_1.9x2.5_tnx1v4.nc' elif model == 'noresm2-mm': gridFileOcn = dataTmp + '/' + 'grid_ocn_tnx1v4.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_0.9x1.25_tnx1v4.nc' elif model == 'cesm1' or model == 'cesm2': gridFileOcn = dataTmp + '/' + 'grid_ocn_gx1v6.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_0.9x1.25_tnx1v4.nc' nc = Dataset(gridFileOcn) lonOcn = nc.variables['plon'][:] latOcn = nc.variables['plat'][:] maskOcn = nc.variables['pmask'][:] if model == 'cesm1' or model == 'cesm2': nc2 = Dataset(dataTmp + '/' + 'coords_cesm1.nc') areaOcn = nc2.variables['TAREA'][:]*1e-4 nc2.close() else: areaOcn = nc.variables['parea'][:] areaOcnNorm = areaOcn / np.nansum(areaOcn*maskOcn) bathOcn = nc.variables['pdepth'][:] nc.close() nc = Dataset(testInputFile) depth_bnds = nc.variables['depth_bnds'][:] nc.close() def readGridAtm(): """Read atmospheric grid information.""" global lonAtm, latAtm, areaAtm, areaAtmNorm, gridFileOcn, gridFileAtm if model == 'noresm1-me': gridFileOcn = dataTmp + '/' + 'grid_ocn_gx1v6.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_1.9x2.5_gx1v6.nc' elif model == 'noresm2-lm': gridFileOcn = dataTmp + '/' + 'grid_ocn_tnx1v4.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_1.9x2.5_tnx1v4.nc' elif model == 'noresm2-mm': gridFileOcn = dataTmp + '/' + 'grid_ocn_tnx1v4.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_0.9x1.25_tnx1v4.nc' elif model == 'cesm1' or model == 'cesm2': gridFileOcn = dataTmp + '/' + 'grid_ocn_gx1v6.nc' gridFileAtm = dataTmp + '/' + 'grid_atm_0.9x1.25_tnx1v4.nc' nc = Dataset(gridFileAtm) lonAtm = nc.variables['lon'][:] latAtm = nc.variables['lat'][:] areaAtm = nc.variables['area'][:] areaAtmNorm = areaAtm / np.sum(areaAtm) nc.close() def filePathMon(field,experiment,year,month,member=0,syear=0,useTestInput=False): """ Create relative path to NorCPM1 input file. Valid components: - atm - lnd - ice - ocn - bgc Valid experiment names: - piControl - historical - 1pctCO2 - abrupt4XCO2 - amip - dcppA-assim-i1 - dcppA-assim-i2 - dcppA-hindcast-i1 - dcppA-hindcast-i2 For month==1 annual ocean output is read. """ if useTestInput: filePath = testInputFile else: component = getComponentFromField(field) if month == 0: (componentTag,relPath) = componentInfo(component,annual=True) else: (componentTag,relPath) = componentInfo(component,annual=False) memTag = '_mem' + '{:0>2d}'.format(member) if experiment == 'piControl': casePrefix = 'noresm1-cmip6_preindustrial_00010101' elif experiment == 'piControlCMIP5': casePrefix = 'N1850AERCNOC_f19_g16_CTRL_02' elif experiment == 'piControl_aerosolbug': casePrefix = 'noresm1-cmip6_preindustrial_areasolbug_00010101' elif experiment == 'piControl_cmip5atm': casePrefix = 'noresm1-cmip6_preindustrial_cmip5atm_00010101' elif experiment == 'piControl_cmip5ice': casePrefix = 'noresm1-cmip6_preindustrial_cmip5ice_00010101' elif experiment == 'piControl_cmip5lnd': casePrefix = 'noresm1-cmip6_preindustrial_cmip5lnd_00010101' elif experiment == 'piControl_cmip5ozone': casePrefix = 'noresm1-cmip6_preindustrial_cmip5ozone_00010101' elif experiment == 'piControl_cmip5volc': casePrefix = 'noresm1-cmip6_preindustrial_cmip5volc_00010101' elif experiment == 'historical_aerosolbug': casePrefix = 'noresm1-cmip6_historical_aerosolbug_18500101' elif experiment == 'historical_cmip5atm': casePrefix = 'noresm1-cmip6_historical_cmip5atm_18500101' elif experiment == 'historical_cmip5ice': casePrefix = 'noresm1-cmip6_historical_cmip5ice_18500101' elif experiment == 'historical_cmip5lnd': casePrefix = 'noresm1-cmip6_historical_cmip5lnd_18500101' elif experiment == 'historical_cmip5ozone': casePrefix = 'noresm1-cmip6_historical_cmip5ozone_18500101' elif experiment == 'historical_cmip5volc': casePrefix = 'noresm1-cmip6_historical_cmip5volc_18500101' elif experiment == 'historical': if year < 1870: startDate = '18500101' elif year < 1950: startDate = '18700101' elif year < 2015: startDate = '19500101' else: startDate = '20150101' casePrefix = 'noresm1-cmip6_historical_' + startDate elif experiment == 'historicalCMIP5': casePrefix = 'NorESM1-ME_historicalExt_noAssim' elif experiment == 'amip': casePrefix = 'noresm1-cmip6_amip_19500101' elif experiment == '1pctCO2': casePrefix = 'noresm1-me_1pctCO2' elif experiment == 'abrupt4XCO2': casePrefix = 'noresm1-me_abrupt4XCO2' elif experiment == 'abrupt4XCO2CMIP5': casePrefix = 'N18504XAERCNOC_f19_g16_02' elif experiment == 'norcpm-cf-system1_assim_19921015': casePrefix = 'norcpm-cf-system1_assim_19921015' elif experiment == 'dcppA-assim-i1': casePrefix = 'norcpm-cmip6_analysis_19500115' elif experiment == 'dcppA-assim-i2': casePrefix = 'noresm1-cmip6_analysis_19500115' elif experiment == 'dcppA-assim-i2b': casePrefix = 'norcpm1_assim-i2_19500115' elif experiment == 'dcppA-assim-i3': casePrefix = 'norcpm1_assim-i3_19500115' elif experiment == 'dcppA-assim-i4': casePrefix = 'norcpm1_assim-i4_19500115' elif experiment == 'dcppA-assim-i5': casePrefix = 'norcpm1_assim-i5_19500115' elif experiment == 'dcppA-hindcast-i1': casePrefix = 'norcpm-cmip6_hindcast' elif experiment == 'dcppA-hindcast-i2': casePrefix = 'noresm1-cmip6_hindcast' elif experiment == 'noresm2-hist1': casePrefix = 'noresm_ctl_f09_tn14_19700101' elif experiment == 'noresm2-assim1': casePrefix = 'norcpm_ana_f09_tn14' elif experiment == 'noresm1-tau': casePrefix = 'N20TREXTAERCN_f19_g16_nudgeEra5AnomUVPS6h_s19500101_4' memTag = '_' + '{:0>1d}'.format(member) elif experiment == 'cesm1-tau': casePrefix = 'B20TRLENS_nudgeEra5AnomUVPS6h_s19500101_05' memTag = '_' + '{:0>2d}'.format(member) elif experiment == 'cesm1-gmw': casePrefix = 'B20TRLENS_nudgeEra5AnomUVPS6hGMW_s19900101_new' memTag = '_' + '{:0>2d}'.format(member) elif experiment == 'cesm1-nogmw': casePrefix = 'B20TRLENS_nudgeEra5AnomUVPS6hNOGMW_s19900101_new' memTag = '_' + '{:0>2d}'.format(member) elif experiment in {'NorESM2-LM_free','noresm2-lm_odamon', 'noresm2-lm_odaday','noresm2-lmesm_assim_10mem','noresm2-lmesm_free_10mem','noresm2-lm_odadaymon_20mem','noresm2-lm_odadaymon_40mem','noresm2-lm_hosing01Sv_19500101','noresm2-lm_hosing-01Sv_19900101','noresm2-lm_hosing-02Sv_19900101','noresm2-lm_histssp370_19900101','noresm2-lm_hosing-01Sv-2_19900101_19900101','noresm2-lm_hosing-01Sv-2_19900101_19900101','noresm2-lm_hosing-02Sv-2_19900101_19900101','noresm2-lm_PIhosing-01Sv_19900101','noresm2-lm_PIhosing-02Sv_19900101','noresm2-lm_dcpp-c-ref','noresm2-lm_dcpp-c-alstephfx','noresm2-lm_dcpp-c-alstephfxtau','noresm2-lm_dcpp-c-naostephfx','noresm2-lm_dcpp-c-naostephfxtau','noresm2-lm_dcpp-c-samsteptau','noresm2-lm_dcpp-c-histnaohfxtau_640pes','noresm2-lm_dcpp-c-histnaohfx','noresm2-lm_dcpp-c-histnaohfx_640pes','noresm2-lm_dcpp-c-histnudge'}: casePrefix = experiment memTag = '_mem' + '{:0>3d}'.format(member) elif experiment in {'NHISTfrc2_f19_tn14_LESFMIPhist-all'}: casePrefix = experiment memTag = '_{:0>3d}'.format(member) elif experiment in {'noresm2-mm_oda_20mem_merged'}: casePrefix = experiment memTag = '_mem{:0>3d}'.format(member) elif experiment in {'noresm2-mm-seaclim_hindcast'}: casePrefix = experiment memTag = '_mem{:0>3d}'.format(member) #print(experiment,componentTag) if componentTag == '.blom.hm.': componentTag = '.blom.hmphyglb.' else: casePrefix = experiment #print(memTag) if month == 0: if member == 0: filePath = dataRaw + '/' + casePrefix + '/' + relPath + '/' + casePrefix + componentTag + '{:0>4d}'.format(year) + '.nc' else: if syear == 0: caseName = casePrefix + memTag filePath = dataRaw + '/' + casePrefix + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '.nc' if experiment == 'noresm1-tau' or experiment == 'cesm1-tau': filePath = dataRaw + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '.nc' else: ensPrefix = casePrefix + '_' + '{:0>4d}'.format(syear) + '1015' caseName = ensPrefix + memTag filePath = dataRaw + '/' + casePrefix + '/' + ensPrefix + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '.nc' else: if member == 0: filePath = dataRaw + '/' + casePrefix + '/' + relPath + '/' + casePrefix + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' else: if syear == 0: caseName = casePrefix + memTag if experiment in {'NorESM2-LM_free','noresm2-lm_odamon', 'noresm2-lm_odaday'}: caseName = casePrefix + '_19960101' + memTag filePath = dataRaw + '/' + casePrefix + '/' + casePrefix + '_19960101' + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' elif experiment in {'noresm2-lm_odadaymon_20mem','noresm2-lm_odadaymon_40mem'}: caseName = casePrefix + '_19820101' + memTag filePath = dataRaw + '/' + casePrefix + '/' + casePrefix + '_19820101' + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' elif experiment in {'noresm2-lm_dcpp-c-ref','noresm2-lm_dcpp-c-alstephfx','noresm2-lm_dcpp-c-alstephfxtau','noresm2-lm_dcpp-c-naostephfx','noresm2-lm_dcpp-c-naostephfxtau','noresm2-lm_dcpp-c-samsteptau','noresm2-lm_dcpp-c-histnaohfxtau_640pes','noresm2-lm_dcpp-c-histnaohfx_640pes','noresm2-lm_dcpp-c-histnaohfx','noresm2-lm_dcpp-c-histnudge'}: caseName = casePrefix + '_19500101' + memTag filePath = dataRaw + '/' + casePrefix + '/' + casePrefix + '_19500101' + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' else: filePath = dataRaw + '/' + casePrefix + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' if experiment == 'noresm1-tau' or experiment == 'cesm1-tau': filePath = dataRaw + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' else: if experiment in {'noresm2-mm-seaclim_hindcast'}: ensPrefix = casePrefix + '_' + '{:0>4d}'.format(syear) + '1101' caseName = ensPrefix + memTag filePath = dataRaw + '/' + casePrefix + '/' + ensPrefix + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' elif casePrefix == 'norcpm-cmip6_hindcast' and syear >= 2022: ensPrefix = 'norcpm1_hindcast-i1' + '_' + '{:0>4d}'.format(syear) + '1015' caseName = ensPrefix + memTag filePath = dataRaw + '/' + casePrefix + '/' + ensPrefix + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' else: ensPrefix = casePrefix + '_' + '{:0>4d}'.format(syear) + '1015' caseName = ensPrefix + memTag filePath = dataRaw + '/' + casePrefix + '/' + ensPrefix + '/' + caseName + '/' + relPath + '/' + caseName + componentTag + '{:0>4d}'.format(year) + '-' + '{:0>2d}'.format(month) + '.nc' #print(filePath,syear,experiment,memTag,member) return filePath def componentInfo(component,annual=False): """Maps components atm, lnd, ice, ocn and bgc to relative paths atm/hist, lnd/hist, ice/hist, ocn/hist and file tags cam2.h0, clm2.h0, cice.h, micom.hm and micom.hbgcm.""" if model[:7] == 'noresm1': ocncomp = 'micom' atmcomp = 'cam2' elif model[:4] == 'cesm': ocncomp = 'pop' atmcomp = 'cam' else: ocncomp = 'blom' #ocncomp = 'micom' atmcomp = 'cam' if component == 'atm': tag = '.' + atmcomp + '.h0.' rpath = 'atm/hist' elif component == 'lnd': tag = '.clm2.h0.' rpath = 'lnd/hist' elif component == 'ice': tag = '.cice.h.' rpath = 'ice/hist' elif component == 'ocn': if annual: tag = '.' + ocncomp + '.hy.' else: if model[:5] == 'cesm1': tag = '.' + ocncomp + '.h.' else: tag = '.' + ocncomp + '.hm.' rpath = 'ocn/hist' elif component == 'bgc': if annual: tag = '.' + ocncomp + '.hbgcy.' else: tag = '.' + ocncomp + '.hbgcm.' rpath = 'ocn/hist' return tag, rpath def getComponentFromField(field): """Returns component of input argument field.""" if field in ['TREFHT','PSL','FLNT','FSNT','PRECT','PRECC','PRECL','Z500']: component = 'atm' elif field in ['templvl','tempinslvl','salnlvl','salninslvl','sst','sss','SST','tos','mmflxd','mmflxl','mhflx','sealv','N_HEAT','N_SALT','MOC','SSH','SALT','TEMP','SHF','SHF_QSW','LWDN_F','LWUP_F','SENH_F','PREC_F','EVAP_F','nsf','swa','hflx','heatContent','saltContent','QFLUX','MELTH_F','FW','MELT_F','IOFF_F','ROFF_F','SNOW_F','SFWF','ztx','mty','taux','tauy','dcpptaux','dcpptauy','dcpphfx']: component = 'ocn' elif field in ['phyc','srfphyc','phyclvl','phycint','srfpp','ppint','pco2','fgco2','co2fxd','co2fxu']: component = 'bgc' elif field in ['aice']: component = 'ice' else: print('Cannot find component for field ' + field) return component def boxTag(lonRange,latRange,niceFormat=False): """Create tag for lon/lat box.""" tag = '' if niceFormat: if latRange[0] >= 0: tag = tag + "{:.0f}".format(latRange[0]) + '\u00B0N\u2013' else: tag = tag + "{:.0f}".format(-latRange[0]) + '\u00B0S\u2013' if latRange[1] >= 0: tag = tag + "{:.0f}".format(latRange[1]) + '\u00B0N' else: tag = tag + "{:.0f}".format(-latRange[1]) + '\u00B0S' tag = tag + ', ' if lonRange[0] >= 0: tag = tag + "{:.0f}".format(lonRange[0]) + '\u00B0E\u2013' else: tag = tag + "{:.0f}".format(-lonRange[0]) + '\u00B0W\u2013' if lonRange[1] >= 0: tag = tag + "{:.0f}".format(lonRange[1]) + '\u00B0E' else: tag = tag + "{:.0f}".format(-lonRange[1]) + '\u00B0W' else: if latRange[0] >= 0: tag = tag + "{:.0f}".format(latRange[0]) + 'N' else: tag = tag + "{:.0f}".format(-latRange[0]) + 'S' if latRange[1] >= 0: tag = tag + "{:.0f}".format(latRange[1]) + 'N' else: tag = tag + "{:.0f}".format(-latRange[1]) + 'S' if lonRange[0] >= 0: tag = tag + "{:.0f}".format(lonRange[0]) + 'E' else: tag = tag + "{:.0f}".format(-lonRange[0]) + 'W' if lonRange[1] >= 0: tag = tag + "{:.0f}".format(lonRange[1]) + 'E' else: tag = tag + "{:.0f}".format(-lonRange[1]) + 'W' return tag def compBox2d1x1(lonRange,latRange,data): """Compute weights for horizontal averaging of ocean output.""" global weights, mask, maskOld if len(data.shape) > 2: mask = np.where(np.isnan(data[0,:,:]),0.,1.) else: mask = np.where(np.isnan(data),0.,1.) if not 'maskOld' in globals() or np.sum(mask-maskOld) != 0 : maskOld = mask res=[1,1] lon=np.arange(res[0]/2,360,res[0]) lat=np.arange(-90+res[1]/2,90,res[1]) weights = mask for j in range(len(lat)): weights[j,:] = weights[j,:] * np.cos(lat[j]*np.pi/180.) if lat[j] > latRange[0] and lat[j] < latRange[1] else 0. if lonRange[0] < lonRange[1]: # box lies within [-180 180] for i in range(len(lon)): weights[:,i] = weights[:,i] if lon[i] > lonRange[0] and lon[i] < lonRange[1] else 0. else: for i in range(len(lon)): weights[:,i] = weights[:,i] if lon[i] > lonRange[0] or lon[i] < lonRange[1] else 0. weights = weights / np.sum(weights) if len(data.shape) > 2: dataBox = np.zeros(data.shape[0]) for i in range(data.shape[0]): dataBox[i] = np.nansum(data[i,:]*weights) else: dataBox = np.nansum(data*weights) return dataBox def compWeightsOcn2d(lonRange,latRange): """Compute weights for horizontal averaging of ocean output.""" global lonOcn, latOcn, areaOcn, maskOcn, areatotOcn, weightsOcn2d if not 'areaOcn' in globals(): readGridOcn() weightsOcn2d = areaOcn * maskOcn for j in range(maskOcn.shape[0]): for i in range(maskOcn.shape[1]): if lonRange[0] < lonRange[1]: # box lies within [-180 180] if not (lonOcn[j,i] > lonRange[0] and lonOcn[j,i] < lonRange[1] and latOcn[j,i] > latRange[0] and latOcn[j,i] < latRange[1]): weightsOcn2d[j,i] = 0. else: # box crosses 180 if not ((lonOcn[j,i] > lonRange[0] or lonOcn[j,i] < lonRange[1]) and latOcn[j,i] > latRange[0] and latOcn[j,i] < latRange[1]): weightsOcn2d[j,i] = 0. areatotOcn = np.nansum(weightsOcn2d) weightsOcn2d = weightsOcn2d / areatotOcn def compWeightsOcn3d(lonRange,latRange,levRange): """Compute weights for 3d averaging of ocean output.""" global bath, depth_bnds, areatotOcn, voltot, weightsOcn2d, weightsOcn3d readGridOcn() bath = bathOcn weightsOcn3d = np.zeros((depth_bnds.shape[0],bath.shape[0],bath.shape[1])) compWeightsOcn2d(lonRange,latRange) for k in range(depth_bnds.shape[0]): bot = np.minimum(min([levRange[1],depth_bnds[k,1]]),bath) top = np.minimum(bot,max([levRange[0],depth_bnds[k,0]])) weightsOcn3d[k,:,:] = areatotOcn * weightsOcn2d * (bot - top) voltot = np.nansum(weightsOcn3d) weightsOcn3d = weightsOcn3d / voltot def outPath(field,experiment,yearRange,lonRange=[0,0],latRange=[0,0],memRange=[0,0],levRange=[0,0],month1=1,syear=0,annual=True,suffix='',ensave=True): prefix = '' if dataTmp == '' else dataTmp + '/' regTag = '' if lonRange[0] == lonRange[1] else '_' + boxTag(lonRange,latRange) if memRange[1] == 0: memTag = '' else: if ensave: memTag = '_mem' + '{:0>2d}'.format(memRange[0]) + '-' + '{:0>2d}'.format(memRange[1]) else: memTag = '_ens' + '{:0>2d}'.format(memRange[0]) + '-' + '{:0>2d}'.format(memRange[1]) levTag = '' if levRange[1] == 0 else '_' + '{:.0f}'.format(levRange[0]) + '-' + '{:.0f}'.format(levRange[1]) + 'm' startTag = '' if syear == 0 else '_s' + '{:0>4d}'.format(syear) yearTag = '_' + str(yearRange[0]) + '-' + str(yearRange[1]) if annual else '_' + str(yearRange[0]) + '{:0>2d}'.format(month1) + '-' + str(yearRange[1]) + '12' return prefix + field + '_' + experiment + startTag + yearTag + regTag + levTag + memTag + suffix + '.nc' def outPathDaily(field,experiment,yearRange,lonRange=[0,0],latRange=[0,0],memRange=[0,0],levRange=[0,0],month1=1,syear=0,suffix='',ensave=True): prefix = '' if dataTmp == '' else dataTmp + '/' regTag = '' if lonRange[0] == lonRange[1] else '_' + boxTag(lonRange,latRange) if memRange[1] == 0: memTag = '' else: if ensave: memTag = '_mem' + '{:0>2d}'.format(memRange[0]) + '-' + '{:0>2d}'.format(memRange[1]) else: memTag = '_ens' + '{:0>2d}'.format(memRange[0]) + '-' + '{:0>2d}'.format(memRange[1]) levTag = '' if levRange[1] == 0 else '_' + '{:.0f}'.format(levRange[0]) + '-' + '{:.0f}'.format(levRange[1]) + 'm' startTag = '' if syear == 0 else '_s' + '{:0>4d}'.format(syear) yearTag = '_' + str(yearRange[0]) + '{:0>2d}01'.format(month1) + '-' + str(yearRange[1]) + '1231' return prefix + field + '_' + experiment + startTag + yearTag + regTag + levTag + memTag + suffix + '.nc' def outPathIceext(experiment,syearRange,smonth=10,leadRange=[0,24],memRange=[0,0],suffix=''): field = 'iceext' prefix = '' if dataTmp == '' else dataTmp + '/' if memRange[1] == 0: memTag = '' else: memTag = '_mem' + '{:0>2d}'.format(memRange[0]) + '-' + '{:0>2d}'.format(memRange[1]) syearTag = '_startyear{:0>4d}'.format(syearRange[0]) if syearRange[0] == syearRange[1] else '_startyear{:0>4d}-{:0>4d}'.format(syearRange[0],syearRange[1]) smonthTag = '_startmonth{:0>2d}'.format(smonth) leadTag = '_leadmonth{:d}'.format(leadRange[0]) if leadRange[0] == leadRange[1] else '_leadmonth{:d}-{:d}'.format(leadRange[0],leadRange[1]) return prefix + field + '_' + experiment + syearTag + smonthTag + leadTag + memTag + suffix + '.nc' def reg2reg(fldIn, resIn, lonBoundsIn, latBoundsIn, resOut, lonBoundsOut=(0,360), latBoundsOut=(-90,90), method='conservative', oceanMaskType='best', dynamicMask=True): """ Regrids from regular lon/lat grid to a regular lon/lat grid. Parameters: ----------- Mandatory fldIn : input field on MICOM grid; the innermost/rightmost two dimensions must be y,x and correspond to the global dimensions in MICOM's grid file; additional dimensions (e.g. vertical and time) are allowed gridPath : full path to MICOM's grid file resIn : input resolution in deg; if skalar is specified the value is used for both longitude and latitude; if tuple is specified the first value is used for longitude and the second for latitude lonBoundsIn : tuple that specifies the western and eastern domain bounds; note that the first longitude value will be shifted by res/2 latBoundsIn : tuple that specifies the southern and northern domain bounds; note that the first longitude value will be shifted by res/2 resOut : output resolution Optional lonBoundsOut : output longitude bounds specification latBoundsOut : output latitude bounds specification method : interpolation method; must be one of 'conservative', 'bilinear', 'patch', 'nearest_s2d' and 'nearest_d2s' oceanMaskType : must be one of 'conservative', 'best' and 'maximized'; 'conservative' sets points to nan if they contain any land, 'best' if more than 50% land and 'maximized' if 100% land Returns: -------- fldInterp, lon2d, lat2d : output field, longitude matrix, latitude matrix History: -------- 2019.10.30 created (Ingo.Bethke@uib.no) """ import xesmf as xe global regridder, gridOut, maskOutInv, shapeInOld, resInOld, resOutOld, \ methodOld, oceanMaskTypeOld # # make sure that input is a masked array; convert missing to mask fldIn = np.ma.masked_invalid(fldIn,copy=False) # # prepare weights if necessary if (not 'regridder' in globals() or fldIn.shape != shapeInOld or resIn != resInOld or resOut != resOutOld or method !=methodOld or oceanMaskType != oceanMaskTypeOld or dynamicMask): shapeInOld, resInOld, resOutOld, methodOld, oceanMaskTypeOld = \ fldIn.shape, resIn, resOut, method, oceanMaskType # # prepare input coordinate information periodic=True if type(resIn) == int or type(resIn) == float: resLon, resLat = resIn, resIn else: resLon, resLat = resIn[0], resIn[1] if lonBoundsIn[0] % 360 != lonBoundsIn[1] % 360: periodic=False gridIn = xe.util.grid_2d(lonBoundsIn[0],lonBoundsIn[1],resLon, latBoundsIn[0],latBoundsIn[1],resLat) # # prepare output coordinate information periodic=True if type(resOut) == int or type(resOut) == float: resLon, resLat = resOut, resOut else: resLon, resLat = resOut[0], resOut[1] if lonBoundsOut[0] % 360 != lonBoundsOut[1] % 360: periodic=False gridOut = xe.util.grid_2d(lonBoundsOut[0],lonBoundsOut[1],resLon, latBoundsOut[0],latBoundsOut[1],resLat) # # compute weights if necessary weightsFile = '{:s}_{:d}x{:d}_{:d}x{:d}.nc'.format(method,gridIn.y.size, \ gridIn.x.size,gridOut.y.size,gridOut.x.size) if os.path.isfile(weightsFile): regridder = xe.Regridder(gridIn, gridOut, method, periodic=periodic, reuse_weights=True,filename=weightsFile) else: regridder = xe.Regridder(gridIn, gridOut, method, periodic=periodic, reuse_weights=False) # # create ocean mask from data and interpolate to output grid maskIn = np.where(fldIn.mask, np.zeros(fldIn.shape), np.ones(fldIn.shape)) maskOut = regridder(maskIn) if oceanMaskType == 'conservative': maskOut = np.where(maskOut>1.-1e-10, maskOut, np.nan) elif oceanMaskType == 'maximized': maskOut = np.where(maskOut>0.+1e-10, maskOut, np.nan) else: if oceanMaskType != 'best': sys.exit('Unknown value for oceanMaskType. Will use "best".') maskOut = np.where(maskOut>0.5, maskOut, np.nan) maskOutInv = 1. / maskOut # # interpolate, normalize and mask # NOTE: To allow interpolation if one or more of contributing input cells # contain missing values, all missing values are set to zero in input. # This has the same effect as setting the grid cell area to zero, so # that cell effectively does not contribute. To ensure conservation, # the same is done in the interpolation of the ocean mask which is # used to normalize the output field. fldOut = regridder(np.where(fldIn.mask,0.,fldIn)) * maskOutInv return fldOut, np.array(gridOut.lon), np.array(gridOut.lat) def micom2reg(fldIn, gridPath, res, lonBounds=(0,360), latBounds=(-90,90), method='conservative', oceanMaskType='best', isMasked=True): """ Regrids MICOM output to a regular lon/lat grid. Parameters: ----------- Mandatory fldIn : input field on MICOM grid; the innermost/rightmost two dimensions must be y,x and correspond to the global dimensions in MICOM's grid file; additional dimensions (e.g. vertical and time) are allowed gridPath : full path to MICOM's grid file res : resolution in deg; if skalar is specified the value is used for both longitude and latitude; if tuple is specified the first value is used for longitude and the second for latitude Optional lonBounds : tuple that specifies the western and eastern domain bounds; note that the first longitude value will be shifted by res/2 latBounds : tuple that specifies the southern and northern domain bounds; note that the first longitude value will be shifted by res/2 method : interpolation method; must be one of 'conservative', 'bilinear', 'patch', 'nearest_s2d' and 'nearest_d2s' oceanMaskType : must be one of 'conservative', 'best' and 'maximized'; 'conservative' sets points to nan if they contain any land, 'best' if more than 50% land and 'maximized' if 100% land Returns: -------- fldInterp, lon2d, lat2d : output field, longitude matrix, latitude matrix Example: -------- # import libraries import matplotlib.pyplot as plt from netCDF4 import Dataset import micom2reg as mr # read gx1 sample data (3d temperature) nc = Dataset('datasample_gx1.nc') T_gx1 = nc.variables['templvl'][:] depth = nc.variables['depth'][:] nc.close() # interpolate data to 1x1 deg lonlat grid T, lon, lat = mr.micom2reg(T_gx1,'gridsample_gx1.nc',1) # plot map of surface temperature plt.pcolormesh(lon,lat,T[0,0,:,:]) # horizontal plot of SST # plot meridional temperature section through Atlantic plt.pcolormesh(lat[:,0],-depth,T[0,:,:,330]) History: -------- 2019.10.30 created (Ingo.Bethke@uib.no) """ import xesmf as xe global regridder, gridOut, maskOutInv, shapeInOld, gridPathOld, resOld, \ methodOld, oceanMaskTypeOld # # make sure that input is a masked array; convert missing to mask fldIn = np.ma.masked_invalid(fldIn,copy=False) # # prepare weights if necessary if (not 'regridder' in globals() or fldIn.shape != shapeInOld or gridPath != gridPathOld or res != resOld or method !=methodOld or oceanMaskType != oceanMaskTypeOld): shapeInOld, gridPathOld, resOld, methodOld, oceanMaskTypeOld = \ fldIn.shape, gridPath, res, method, oceanMaskType # # prepare MICOM coordinate information nc = Dataset(gridPath) lon, lat = nc.variables['plon'][:], nc.variables['plat'][:] clon, clat = nc.variables['pclon'][:], nc.variables['pclat'][:] if not isMasked: pmask = nc.variables['pmask'][:] nc.close() ydm, xdm = lon.shape lonC, latC = np.zeros((ydm+1,xdm+1)), np.zeros((ydm+1,xdm+1)) lonC[0:ydm,0:xdm], latC[0:ydm,0:xdm] = clon[0,:,:], clat[0,:,:] lonC[ydm,0:xdm], latC[ydm,0:xdm] = clon[3,ydm-1,:], clat[3,ydm-1,:] lonC[0:ydm,xdm], latC[0:ydm,xdm] = clon[1,:,xdm-1], clat[1,:,xdm-1] lonC[ydm,xdm], latC[ydm,xdm] = clon[2,ydm-1,xdm-1], clat[2,ydm-1,xdm-1] gridIn = {'lon': lon, 'lon_b': lonC, 'lat': lat, 'lat_b': latC} # # prepare regular grid coordinate information periodic=True if type(res) == int or type(res) == float: resLon, resLat = res, res else: resLon, resLat = res[0], res[1] if lonBounds[0] % 360 != lonBounds[1] % 360: periodic=False gridOut = xe.util.grid_2d(lonBounds[0],lonBounds[1],resLon, latBounds[0],latBounds[1],resLat) # # compute weights weightsFile = '{:s}_{:d}x{:d}_{:d}x{:d}.nc'.format(method,ydm,xdm, \ gridOut.y.size,gridOut.x.size) if os.path.isfile(weightsFile): regridder = xe.Regridder(gridIn, gridOut, method, periodic=periodic, reuse_weights=True,filename=weightsFile) else: regridder = xe.Regridder(gridIn, gridOut, method, periodic=periodic, reuse_weights=False) # # create ocean mask from data and interpolate to output grid if not isMasked: maskIn = np.where(pmask < 0.5, np.zeros(fldIn.shape), np.ones(fldIn.shape)) else: maskIn = np.where(fldIn.mask, np.zeros(fldIn.shape), np.ones(fldIn.shape)) maskOut = regridder(maskIn) if oceanMaskType == 'conservative': maskOut = np.where(maskOut>1.-1e-10, maskOut, np.nan) elif oceanMaskType == 'maximized': maskOut = np.where(maskOut>0.+1e-10, maskOut, np.nan) else: if oceanMaskType != 'best': sys.exit('Unknown value for oceanMaskType. Will use "best".') maskOut = np.where(maskOut>0.5, maskOut, np.nan) maskOutInv = 1. / maskOut # # interpolate, normalize and mask # NOTE: To allow interpolation if one or more of contributing input cells # contain missing values, all missing values are set to zero in input. # This has the same effect as setting the grid cell area to zero, so # that cell effectively does not contribute. To ensure conservation, # the same is done in the interpolation of the ocean mask which is # used to normalize the output field. fldOut = regridder(np.where(fldIn.mask,0.,fldIn)) * maskOutInv return fldOut, np.array(gridOut.lon), np.array(gridOut.lat) def eosben07_p_alpha(p,p0,th,s): """ eosben07_p_alpha Integration of specific volume with respect to pressure. eosben07_p_alpha(P, P0, TH, S) integrates seawater specific volume from pressure P0 to P [dbar] with potential temperature (referenced to 0 dbar), TH [degC], and salinity, S [g kg-1]. Units m3/kg. """ a11= 9.9985372432159340e+02; a12= 1.0380621928183473e+01; a13= 1.7073577195684715e+00; a14=-3.6570490496333680e-02; a15=-7.3677944503527477e-03; a16=-3.5529175999643348e-03; a21= 1.0; a22= 1.0316374535350838e-02; a23= 8.9521792365142522e-04; a24=-2.8438341552142710e-05; a25=-1.1887778959461776e-05; a26=-4.0163964812921489e-06; b11= 1.7083494994335439e-02; b12= 7.1567921402953455e-05; b13= 1.2821026080049485e-05; b21= 1.1995545126831476e-05; b22= 5.5234008384648383e-08; b23= 8.4310335919950873e-09; b1i=1/(b11+b12*th+b13*s); a1=(a11+th*(a12+a14*th+a15*s)+s*(a13+a16*s))*b1i; a2=(a21+th*(a22+a24*th+a25*s)+s*(a23+a26*s))*b1i; b2=(b21+b22*th+b23*s)*b1i; r=b2*(p-p0)+(a2-a1*b2)*np.log((a1+p)/(a1+p0)); return r def prepare_aice_hindcast(syear): field = 'aice' annual = False year1 = syear + 0 yearn = syear + 10 writeEnsAve(field,'dcppA-hindcast-i1',[year1,yearn],[1,10],month1=10,syear=syear,annual=annual) writeEnsAve(field,'dcppA-hindcast-i2',[year1,yearn],[1,10],month1=10,syear=syear,annual=annual) for ensave in [True]: for res in [[1,1],[2,2],[5,5]]: for field in ['aice']: for experiment in ['dcppA-hindcast-i1','dcppA-hindcast-i2']: regrid2reg(field,experiment,[year1,yearn],res=res,memRange=[1,10],syear=syear,month1=10,annual=annual,ensave=ensave) def prepare_aice_historical(yearRange=[1950,2029],memRange=[1,30],annual=False): field = 'aice' writeEnsAve(field,'historical',yearRange,memRange,annual=annual) for res in [[1,1],[2,2],[5,5]]: regrid2reg(field,'historical',yearRange,res=res,memRange=memRange,annual=annual,ensave=True) def prepare_sst_monthly(annual=False): field = 'sst' writeEnsAve(field,'historical',[1950,2029],[1,30],annual=annual) writeEnsAve(field,'dcppA-assim-i1',[1950,2018],[1,30],annual=annual) writeEnsAve(field,'dcppA-assim-i2',[1950,2018],[1,30],annual=annual) for res in [[1,1],[2,2],[5,5]]: regrid2reg(field,'historical',[1950,2029],res=res,memRange=[1,30],annual=annual,ensave=True) regrid2reg(field,'dcppA-assim-i1',[1950,2018],res=res,memRange=[1,30],annual=annual,ensave=True) regrid2reg(field,'dcppA-assim-i2',[1950,2018],res=res,memRange=[1,30],annual=annual,ensave=True) for syear in range(1979,2016+1): year1 = syear + 0 yearn = syear + 10 for experiment in ['dcppA-hindcast-i1','dcppA-hindcast-i2']: writeEnsAve(field,experiment,[year1,yearn],[1,10],month1=10,syear=syear,annual=annual) for res in [[1,1],[2,2],[5,5]]: regrid2reg(field,experiment,[year1,yearn],res=res,memRange=[1,10],syear=syear,month1=10,annual=annual,ensave=True)