import numpy as np from scipy import constants #''' #* 1) Use only for SW components, better code available soon from Joonas Merikanto or Kalle Nordling! #SW term is futer decomposed using method by taylor et.al (2006) #code used CMIP standard for field names #Original python code Kalle Nordling, 2019 #''' def subs(x1,x2): res = np.zeros(x1.shape) for i in range(0,x1.shape[0]): for j in range(0,x1.shape[1]): for k in range(0,x1.shape[2]): res[i,j,k] = x1[i,j,k]-x2[i,j,k] class RY15: def __init__(self,dataset1,dataset2): self.dataset1 = {} self.dataset2 = {} #self.dataset1 = dataset1 #self.dataset2 = dataset2 for key,val in dataset2.items(): print(key) self.dataset2[key] = val.to_masked_array() for key,val in dataset1.items(): self.dataset1[key] = val.to_masked_array() #input('WAIT') self.results = {} #datastructure to hold results self.proseduce() ''' Main function where calculations are done ''' def proseduce(self): self.addOverCastFields() self.addEpsilonFields() self.results['DeltaT'] = self.calcDelta('tas') self.calcFactorD() self.calcLWterm() self.calcSWterms() self.calcSURFterm() self.calcCONVterms() self.doTaylor2() self.calcERF() #self.results['cloud'] = self.results['Delta_SW_cld']+self.results['LW_cre'] self.results['Dt_tot'] = self.results['LW'] + self.results['SW'] + self.results['SURF'] + self.results['CONV'] self.results['SW_check'] = self.results['dtswin']+self.results['dtAclr']+self.results['dtAcld']+self.results['dtAalfa']+self.results['dtAnl'] #self.calcERF() #self.calcSDCs() # def calcERF(self): # self.dataset1['NET_SW_TOA'] = self.dataset1_fsst['rsdt']-self.dataset1_fsst['rsut'] # self.dataset2['NET_SW_TOA'] = self.dataset2_fsst['rsdt']-self.dataset2_fsst['rsut'] # self.dataset1['rlut_fsst'] = self.dataset1_fsst['rlut'] # self.dataset2['rlut_fsst'] = self.dataset2_fsst['rlut'] # print(self.dataset1['NET_SW_TOA'].shape,self.dataset2['NET_SW_TOA'].shape) # print(self.dataset1_fsst['rsdt'].shape,self.dataset1_fsst['rsut'].shape) # # self.dataset1['NET_SW_TOA_CS'] = self.dataset1_fsst['rsdt']-self.dataset1_fsst['rsutcs'] # self.dataset2['NET_SW_TOA_CS'] = self.dataset2_fsst['rsdt']-self.dataset2_fsst['rsutcs'] # self.dataset1['rlutcs_fsst'] = self.dataset1_fsst['rlutcs'] # self.dataset2['rlutcs_fsst'] = self.dataset2_fsst['rlutcs'] # # self.results['ERF'] = self.calcDelta('NET_SW_TOA')+ self.calcDelta('rlut_fsst') # self.results['ERFCS'] = self.calcDelta('NET_SW_TOA_CS')+ self.calcDelta('rlutcs_fsst') # self.results['ERF_cloud'] = self.results['ERF']-self.results['ERFCS'] def calcERF(self): self.dataset1['NET_SW_TOA'] = self.dataset1['rsdt']-self.dataset1['rsut'] self.dataset2['NET_SW_TOA'] = self.dataset2['rsdt']-self.dataset2['rsut'] print(self.dataset1['NET_SW_TOA'].shape,self.dataset2['NET_SW_TOA'].shape) print(self.dataset1['rsdt'].shape,self.dataset1['rsut'].shape) self.dataset1['NET_SW_TOA_CS'] = self.dataset1['rsdt']-self.dataset1['rsutcs'] self.dataset2['NET_SW_TOA_CS'] = self.dataset2['rsdt']-self.dataset2['rsutcs'] self.results['ERF'] = self.calcDelta('NET_SW_TOA')-self.calcDelta('rlut') self.results['ERFCS'] = self.calcDelta('NET_SW_TOA_CS')- self.calcDelta('rlutcs') self.results['ERF_cloud'] = self.results['ERF']-self.results['ERFCS'] def calcSDC(self,field,deltaT): field_prime = field-np.nanmean(field,axis=0) deltaT_prime =deltaT- np.nanmean(deltaT,axis=0) tmp = field_prime*deltaT_prime tmp2=np.sqrt(deltaT_prime**2) return np.nanmean(tmp,axis=0)/np.nanmean(tmp2,axis=0) def calcSDCs(self): tmp = self.results.copy() for key,value in tmp.items(): self.results[key+"_sdc"] = self.calcSDC(value,tmp['DeltaT']) ''' Function used to calc delta of some fields. If order of calculation is needed to chage it can be done easyly ''' def calcDelta(self,field): try: #res=self.dataset1[field].to_masked_array()-self.dataset2[field].to_masked_array() #print(type(self.dataset1[field])) #print(self.dataset1[field]) #print('calcldelta '+field,self.dataset1[field].shape,self.dataset2[field].shape,res.shape) return self.dataset1[field]-self.dataset2[field] except: print('DELTA calculation faild, cannot find field '+field+' from datasets') ''' function to calculated mean of two fields. Seperate function if method needs to be changes ''' def calcMeanField(self,field): try: return (self.dataset1[field]+self.dataset2[field])*0.5 except: print('MEAN calculation faild, cannot find field '+field+' from datasets') ''' extra helper function to calculate surf term, calculate net surf flux ''' def surfHelp(self,dataset): tmp1 = dataset['hfls']+dataset['hfss'] tmp2 = dataset['rsus']-dataset['rsds'] tmp3 = dataset['rlus']-dataset['rlds'] dataset['NETSURF'] = dataset['rsds']-dataset['rsus']+dataset['rlds']-dataset['rlus']-dataset['hfls']-dataset['hfss'] #dataset['NETSURF'] = tmp1+tmp2+tmp3 ''' extra helpper function for taylor decompostion,epsilon is added to avoid dividing by 0 cases ''' def taylorHelp(self,dataset): epsilon = 1.2345e-7 dataset['Q_hat'] = dataset['rsds']/(dataset['rsdt']+epsilon) dataset['Q_hatcs'] = dataset['rsdscs']/(dataset['rsdt']+epsilon) dataset['Q_hatoc'] = dataset['rsdsoc']/(dataset['rsdt']+epsilon) ''' calculate planetary albedo ''' def calculatePlanetaryAlbedo(self,dataset): epsilon = 1.2345e-7 dataset['A'] = 1-((dataset['rsdt']-dataset['rsut'])/(dataset['rsdt']+epsilon)) dataset['A_cs'] = 1-((dataset['rsdt']-dataset['rsutcs'])/(dataset['rsdt']+epsilon)) dataset['A_oc'] = 1-((dataset['rsdt']-dataset['rsutoc'])/(dataset['rsdt']+epsilon)) ''' calculate mu and gamma parametrs, from taylor 2007 mu = A+Q(1+a), eq. 9, gamma = (mu-Q)/(mu-aQ), eq. 10 ''' def calcMuGamma(self,dataset): # epsilon = 1.2345e-7 epsilon = 1.2345e-10 dataset['mu'] = (1-dataset['alpha'])*dataset['Q_hat']+dataset['A'] dataset['mu_cs'] = (1-dataset['alphacs'])*dataset['Q_hatcs']+dataset['A_cs'] dataset['mu_oc'] = (1-dataset['alphaoc'])*dataset['Q_hatoc']+dataset['A_oc'] dataset['mu_cl'] = dataset['mu_oc']/(dataset['mu_cs']++epsilon) dataset['gamma'] = (dataset['mu']-dataset['Q_hat'])/((dataset['mu']-dataset['alpha']*dataset['Q_hat'])+epsilon) dataset['gamma_cs'] = (dataset['mu_cs']-dataset['Q_hatcs'])/((dataset['mu_cs']-dataset['alphacs']*dataset['Q_hatcs'])+epsilon) dataset['gamma_oc'] = (dataset['mu_oc']-dataset['Q_hatoc'])/((dataset['mu_oc']-dataset['alphaoc']*dataset['Q_hatoc'])+epsilon) dataset['gamma_cl'] = 1-((1-dataset['gamma_oc'])/(1-dataset['gamma_cs']+epsilon)) dataset['A_ocC'] = (dataset['A']-((1-(dataset['clt']))*dataset['A_cs']))/((dataset['clt'])+1.2345e-7) ''' This function treats partly cloudy regions according to Taylor 2007 R = (1-c)*R_clr+c*R_oc eq.3 in Taylors papers, solve R_oc R_oc = (R-(1-c)*R_clr)/(c+epsilon) small epsilon is added to avoid divison by 0 ''' def taylor_overcast(self,R,R_clr,cloud): epsilon = 1.2345e-7 return (R-(cloud-1)*R_clr)/(cloud+epsilon) ''' This function adds calculated overcast fields to dataset structures ''' def calc_overcasted_fileds(self,dataset): dataset['rsutoc']=self.taylor_overcast(dataset['rsut'],dataset['rsutcs'],dataset['clt']) dataset['rsdsoc']=self.taylor_overcast(dataset['rsds'],dataset['rsdscs'],dataset['clt']) dataset['rsusoc']=self.taylor_overcast(dataset['rsus'],dataset['rsuscs'],dataset['clt']) ''' Add overcasted fields to datasets ''' def addOverCastFields(self): self.calc_overcasted_fileds(self.dataset1) self.calc_overcasted_fileds(self.dataset2) ''' calc planetary emissivity LW=outgoing LW radiation (rlut),sigma is stefan-bolzman constant LW = epsilon*sigma*T^4, eq.2 in raisainen paper 2017 , from here solve epsilon epsilon = LW/(sigma*T^4) add result to dataset calculation is done for both all-sky and clear-sky terms ''' def calcEpsilon(self,dataset): dataset['eps'] = dataset['rlut']/(dataset['tas']**4*constants.sigma) dataset['epscs'] = dataset['rlutcs']/(dataset['tas']**4*constants.sigma) #print(dataset['eps'].shape) #input('wait eps shape') ''' function adds epsilon fields to datasets ''' def addEpsilonFields(self): self.calcEpsilon(self.dataset1) self.calcEpsilon(self.dataset2) ''' Function calculates factor D in raisanen paper D=4*sigma[epsilon]*[T]^3 epsilon and T are means of two datasets ''' def calcFactorD(self): mean_T = self.calcMeanField('tas') mean_epsilon = self.calcMeanField('eps') self.results['D'] = 4.0*constants.sigma*mean_epsilon*mean_T**3 ''' calculate surface albedo alpha = Q_s_up/Q_s_down ''' def calcSurfaceAlbedo(self,dataset): epsilon = 1.2345e-7 dataset['alpha'] = dataset['rsus']/(dataset['rsds']+epsilon) dataset['alphacs'] = dataset['rsuscs']/(dataset['rsdscs']+epsilon) dataset['alphaoc'] = dataset['rsusoc']/(dataset['rsdsoc']+epsilon) ''' calc A's eq. 7 in taylor paper ''' def calcA(self,mu,gamma,alpha): epsilon = 1.2345e-7 return mu*gamma+((mu*alpha*(1-gamma)**2)/(1-alpha*gamma+epsilon)) def taylorHepl(self,dataset): eps = 1.2345e-12 dataset['rsutoc'] = (dataset['rsut']+(dataset['c']-1.0)*dataset['rsutcs'])/(dataset['c']+eps) dataset['rsdsoc'] = (dataset['rsds']+(dataset['c']-1.0)*dataset['rsdscs'])/(dataset['c']+eps) dataset['rsusoc'] = (dataset['rsus']+(dataset['c']-1.0)*dataset['rsuscs'])/(dataset['c']+eps) dataset['A'] = dataset['rsut']/(dataset['s']+eps) dataset['Aclr'] = dataset['rsutcs']/(dataset['s']+eps) dataset['Aoc'] = dataset['rsutoc']/(dataset['s']+eps) dataset['alfa']=dataset['rsus']/(dataset['rsds']+eps) dataset['alfaclr']=dataset['rsuscs']/(dataset['rsdscs']+eps) dataset['alfaoc']=dataset['rsusoc']/(dataset['rsdsoc']+eps) ''' for t in range(0,dataset['alfa'].shape[0]): for x in range(0,dataset['alfa'].shape[1]): for y in range(0,dataset['alfa'].shape[2]): for i in ['alfa','alfaclr','alfaoc']: if (dataset[i][t,x,y] > 1): dataset[i][t,x,y] = 1 if (dataset[i][t,x,y] < 0): dataset[i][t,x,y] = 0 ''' dataset['Qrat'] = dataset['rsds']/(dataset['s']+eps) dataset['Qratclr'] = dataset['rsdscs']/(dataset['s']+eps) dataset['Qratoc'] = dataset['rsdsoc']/(dataset['s']+eps) # # dataset['my_ori']=dataset['A']+dataset['Qrat']*(1-dataset['alfa']) # dataset['myclr_ori']=dataset['Aclr'] +dataset['Qratclr']*(1-dataset['alfaclr']) # dataset['myoc_ori']=dataset['Aoc']+dataset['Qratoc']*(1-dataset['alfaoc']) # dataset['my']=np.maximum(dataset['A']+dataset['Qrat']*(1-dataset['alfa']),dataset['Qrat']) dataset['myclr']=np.maximum(dataset['Aclr']+dataset['Qratclr']*(1-dataset['alfaclr']),dataset['Qratclr']) # dataset['myclr']max(=dataset['Aclr'] +dataset['Qratclr']*(1-dataset['alfaclr']),dataset['Qratclr']) dataset['myoc']=np.maximum(dataset['Aoc']+dataset['Qratoc']*(1-dataset['alfaoc']),dataset['Qratoc']) # dataset['gamma']=(dataset['my']-dataset['Qrat'])/(dataset['my']-dataset['alfa']*dataset['Qrat']+eps) dataset['gammaclr']=(dataset['myclr']-dataset['Qratclr'])/(dataset['myclr']-dataset['alfaclr']*dataset['Qratclr']+eps) dataset['gammaoc']=(dataset['myoc']-dataset['Qratoc'])/(dataset['myoc']-dataset['alfaoc']*dataset['Qratoc']+eps) dataset['gammacld']=1-(1-dataset['gammaoc'])/(1-dataset['gammaclr']) dataset['mycld']=(dataset['myoc']+eps)/(dataset['myclr']+eps) def Afrompar(self,xc,xalfaclr,xalfaoc,xmyclr,xmycld,xgammaclr,xgammacld): xgammaoc=1-(1-xgammaclr)*(1-xgammacld) xmyoc=xmyclr*xmycld xAclr=xmyclr*xgammaclr+xmyclr*xalfaclr*(1-xgammaclr)**2/(1-xalfaclr*xgammaclr) xAoc=xmyoc*xgammaoc+xmyoc*xalfaoc*(1-xgammaoc)**2/(1-xalfaoc*xgammaoc) return (1-xc)*xAclr+xc*xAoc def doTaylor2(self): self.dataset1['s'] = self.dataset1['rsdt'] self.dataset2['s'] = self.dataset2['rsdt'] self.dataset1['c'] = self.dataset1['clt'] self.dataset2['c'] = self.dataset2['clt'] self.taylorHepl(self.dataset1) self.taylorHepl(self.dataset2) a=self.calcMeanField('A') s=self.calcMeanField('s') c=self.calcMeanField('c') alfaclr=self.calcMeanField('alfaclr') alfaoc=self.calcMeanField('alfaoc') myclr=self.calcMeanField('myclr') mycld=self.calcMeanField('mycld') gammaclr=self.calcMeanField('gammaclr') gammacld=self.calcMeanField('gammacld') c1=self.dataset1['c'] c2=self.dataset2['c'] myclr1=self.dataset1['myclr'] myclr2=self.dataset2['myclr'] # myclr1=max(self.dataset1['myclr'],self.dataset1['Qrad']) # myclr2=max(self.dataset2['myclr'],self.dataset2['Qrad']) alfaclr1=self.dataset1['alfaclr'] alfaoc1=self.dataset1['alfaoc'] mycld1=self.dataset1['mycld'] gammaclr1=self.dataset1['gammaclr'] gammacld1=self.dataset1['gammacld'] alfaclr2=self.dataset2['alfaclr'] alfaoc2=self.dataset2['alfaoc'] mycld2=self.dataset2['mycld'] gammaclr2=self.dataset2['gammaclr'] gammacld2=self.dataset2['gammacld'] self.dataset1['apar'] = self.Afrompar(c1,alfaclr1,alfaoc1,myclr1,mycld1,gammaclr1,gammacld1) self.dataset2['apar'] = self.Afrompar(c2,alfaclr2,alfaoc2,myclr2,mycld2,gammaclr2,gammacld2) self.dataset1['aalfaclr'] = self.Afrompar(c,alfaclr1,alfaoc,myclr,mycld,gammaclr,gammacld) self.dataset2['aalfaclr'] = self.Afrompar(c,alfaclr2,alfaoc,myclr,mycld,gammaclr,gammacld) self.dataset1['aalfaoc'] = self.Afrompar(c,alfaclr,alfaoc1,myclr,mycld,gammaclr,gammacld) self.dataset2['aalfaoc'] = self.Afrompar(c,alfaclr,alfaoc2,myclr,mycld,gammaclr,gammacld) self.dataset1['amycld'] = self.Afrompar(c,alfaclr,alfaoc,myclr,mycld1,gammaclr,gammacld) self.dataset2['amycld'] = self.Afrompar(c,alfaclr,alfaoc,myclr,mycld2,gammaclr,gammacld) self.dataset1['agammacld'] = self.Afrompar(c,alfaclr,alfaoc,myclr,mycld,gammaclr,gammacld1) self.dataset2['agammacld'] = self.Afrompar(c,alfaclr,alfaoc,myclr,mycld,gammaclr,gammacld2) self.dataset1['ac'] = self.Afrompar(c1,alfaclr,alfaoc,myclr,mycld,gammaclr,gammacld) self.dataset2['ac'] = self.Afrompar(c2,alfaclr,alfaoc,myclr,mycld,gammaclr,gammacld) self.dataset1['amyclr'] = self.Afrompar(c,alfaclr,alfaoc,myclr1,mycld,gammaclr,gammacld) self.dataset2['amyclr'] = self.Afrompar(c,alfaclr,alfaoc,myclr2,mycld,gammaclr,gammacld) self.dataset1['agammaclr'] = self.Afrompar(c,alfaclr,alfaoc,myclr,mycld,gammaclr1,gammacld) self.dataset2['agammaclr'] = self.Afrompar(c,alfaclr,alfaoc,myclr,mycld,gammaclr2,gammacld) self.dataset1['acld'] = self.Afrompar(c1,alfaclr,alfaoc,myclr,mycld1,gammaclr,gammacld1) self.dataset2['acld'] = self.Afrompar(c2,alfaclr,alfaoc,myclr,mycld2,gammaclr,gammacld2) self.dataset1['aalfa'] = self.Afrompar(c,alfaclr1,alfaoc1,myclr,mycld,gammaclr,gammacld) self.dataset2['aalfa'] = self.Afrompar(c,alfaclr2,alfaoc2,myclr,mycld,gammaclr,gammacld) self.dataset1['aclr'] = self.Afrompar(c,alfaclr,alfaoc,myclr1,mycld,gammaclr1,gammacld) self.dataset2['aclr'] = self.Afrompar(c,alfaclr,alfaoc,myclr2,mycld,gammaclr2,gammacld) ddaalfa= self.calcDelta('aalfa') #(aalfa2-aalfa1)S ddaclr= self.calcDelta('aclr') #(aclr2-aclr1) ddacld= self.calcDelta('acld')#(acld2-acld1) daalfa= self.calcDelta('aalfaclr')+self.calcDelta('aalfaoc') #(aalfaclr2-aalfaclr1)+(aalfaoc2-aalfaoc1) daalfaclr= self.calcDelta('aalfaclr')#(aalfaclr2-aalfaclr1) daalfaoc= self.calcDelta('aalfaoc')#(aalfaoc2-aalfaoc1) dacld= self.calcDelta('amycld')+self.calcDelta('agammacld')+self.calcDelta('ac')#(amycld2-amycld1)+(agammacld2-agammacld1)+(ac2-ac1) dacldmy= self.calcDelta('amycld')#(amycld2-amycld1) dacldgamma=self.calcDelta('agammacld') #(agammacld2-agammacld1) dacldc= self.calcDelta('ac')#(ac2-ac1) daclr= self.calcDelta('amyclr')+self.calcDelta('agammaclr')#(amyclr2-amyclr1)+(agammaclr2-agammaclr1) daclrmy= self.calcDelta('amyclr')#(amyclr2-amyclr1) daclrgamma= self.calcDelta('agammaclr')#(agammaclr2-agammaclr1) da=self.calcDelta('A') danl=da-(daalfa+dacld+daclr) ddanl=da-(ddaalfa+ddacld+ddaclr) fs=self.calcDelta('s')*(1-a) fA=-s*da fAclr=-s*daclr fAcld=-s*dacld fAalfa=-s*daalfa fAnl=-s*danl fAcldc=-s*dacldc sensitivity=1.0/self.results['D'] self.results['dtswin']=sensitivity*fs self.results['dtAclr']=sensitivity*fAclr self.results['dtAcld']=sensitivity*fAcld self.results['dtAalfa']=sensitivity*fAalfa self.results['dtAnl']=sensitivity*fAnl self.results['dtAcldc']=sensitivity*fAcldc self.results['daclrmy']=daclrmy self.results['daclrgamma']=daclrgamma self.results['dacldmy']=dacldmy self.results['myclr1']=myclr1 self.results['gammaclr1']=gammaclr1 self.results['amyclr1']=self.dataset1['amyclr'] self.results['myclr2']=myclr2 self.results['gammaclr2']=gammaclr2 self.results['amyclr2']=self.dataset2['amyclr'] self.results['agammaclr1']=self.dataset1['agammaclr'] self.results['agammaclr2']=self.dataset2['agammaclr'] self.results['dacldgamma']=dacldgamma self.results['dacldc']=self.calcDelta('ac') ''' do taylor decomposition ''' def doTaylor(self): self.calcSurfaceAlbedo(self.dataset1) self.calcSurfaceAlbedo(self.dataset2) #calc mean alphas brac_alpha = self.calcMeanField('alpha') brac_alphacs = self.calcMeanField('alphacs') brac_alphaoc = self.calcMeanField('alphaoc') #calc Q_hat terms self.taylorHelp(self.dataset1) self.taylorHelp(self.dataset2) #calc planetary albedo (A term, in taylor paper) self.calculatePlanetaryAlbedo(self.dataset1) self.calculatePlanetaryAlbedo(self.dataset2) brac_a = self.calcMeanField('A') brac_acs = self.calcMeanField('A_cs') brac_aoc = self.calcMeanField('A_oc') #calc mu and gamma self.calcMuGamma(self.dataset1) self.calcMuGamma(self.dataset2) brac_mu = self.calcMeanField('mu') brac_mucs = self.calcMeanField('mu_cs') brac_muoc = self.calcMeanField('mu_oc') brac_mucl = self.calcMeanField('mu_cl') brac_gamma = self.calcMeanField('gamma') brac_gammacs = self.calcMeanField('gamma_cs') brac_gammaoc = self.calcMeanField('gamma_oc') brac_gammacl = self.calcMeanField('gamma_cl') brac_clt = self.calcMeanField('clt') #calc deltaA's #delta alpha gamma_oc = 1.0-(1.0-brac_gammacl)*(1.0-brac_gammacs) mu_oc = brac_mucs*brac_mucl self.dataset1['A_clr1'] = self.calcA(brac_mucs,brac_gammacs,self.dataset1['alphacs']) self.dataset1['A_oc1'] = self.calcA(brac_muoc,brac_gammaoc,self.dataset1['alphaoc']) self.dataset2['A_clr1'] = self.calcA(brac_mucs,brac_gammacs,self.dataset2['alphacs']) self.dataset2['A_oc1'] = self.calcA(brac_muoc,brac_gammaoc,self.dataset2['alphaoc']) delta_A_alpha_clr = (1-brac_clt)*self.calcDelta('A_clr1') delta_A_alpha_oc = brac_clt*self.calcDelta('A_oc1') #delta mu self.dataset1['A_clr1'] = self.calcA(self.dataset1['mu_cs'],brac_gammacs,brac_alpha) self.dataset1['A_oc1'] = self.calcA(brac_mucs*self.dataset1['mu_oc'],brac_gammaoc,brac_alphaoc) self.dataset2['A_clr1'] = self.calcA(self.dataset2['mu_cs'],brac_gammacs,brac_alpha) self.dataset2['A_oc1'] = self.calcA(brac_mucs*self.dataset2['mu_oc'],brac_gammaoc,brac_alphaoc) delta_A_mu_clr = (1-brac_clt)*self.calcDelta('A_clr1') delta_A_mu_oc = brac_clt*self.calcDelta('A_oc1') #delta gamma gamma_oc = 1-(1-self.dataset1['gamma_cl'])*(1-brac_gammacs) self.dataset1['A_clr1'] = self.calcA(brac_mucs,self.dataset1['gamma_cs'],brac_alphacs) self.dataset1['A_oc1'] = self.calcA(brac_mucs,gamma_oc,brac_alphaoc) gamma_oc = 1-(1-self.dataset2['gamma_cl'])*(1-brac_gammacs) self.dataset2['A_clr1'] = self.calcA(brac_mucs,self.dataset2['gamma_cs'],brac_alphacs) self.dataset2['A_oc1'] = self.calcA(brac_mucs,gamma_oc,brac_alphaoc) delta_A_gamma_clr = (1-brac_clt)*self.calcDelta('A_clr1') delta_A_gamma_oc = brac_clt*self.calcDelta('A_oc1') #delta C self.dataset1['A_clr1'] = self.calcA(brac_mucs,brac_gammacs,brac_alphacs) self.dataset1['A_oc1'] = self.calcA(brac_muoc,brac_gammaoc,brac_alphaoc) A_clr = self.calcDelta('A_clr1') A_oc = brac_clt*self.calcDelta('A_oc1') self.dataset1['A_1'] = (1-(self.dataset1['clt']))*A_clr+(self.dataset1['clt'])*A_oc self.dataset2['A_1'] = (1-(self.dataset2['clt']))*A_clr+(self.dataset2['clt'])*A_oc delta_A_c = self.calcDelta('A_1') #add to results self.results['Delta_SW_alpha'] = ((delta_A_alpha_clr+delta_A_alpha_oc)*-1.0*self.calcMeanField('rsdt'))/self.results['D'] self.results['Delta_SW_clr'] = ((delta_A_mu_clr+delta_A_gamma_clr)*-1.0*self.calcMeanField('rsdt'))/self.results['D'] self.results['Delta_SW_cld'] = ((delta_A_mu_oc+delta_A_gamma_oc+delta_A_c)*-1.0*self.calcMeanField('rsdt'))/self.results['D'] self.results['SW_check'] = self.results['Delta_SW_alpha']+self.results['Delta_SW_clr']+self.results['Delta_SW_cld'] self.results['Delta_SW_c'] = ((delta_A_c)*-1.0*self.calcMeanField('rsdt'))/self.results['D'] ''' function calculate LW change factor all sky and cloudy components calculation is done according to raisanen paper eq.4 factor 1 LW = - sigma*Delta_Epsilon*T^4 Temperature here is mean temperature LW = LW_clear+LW_cloudy part -> LW_cre = LW-LW_Clearsky ''' def calcLWterm(self): deltaEPS = self.calcDelta('eps') deltaEPSCS = self.calcDelta('epscs') mean_T_power4 = (self.dataset1['tas']**4+self.dataset2['tas']**4)*0.5 self.results['LW'] = (-1.0*constants.sigma*deltaEPS*mean_T_power4)/self.results['D'] print('LW shape',self.results['LW'].shape) print('D shape',self.results['D'].shape) print('temp 4 shape',mean_T_power4.shape) print('EPS shape',deltaEPS.shape) self.results['LW_clearsky'] = (-1.0*constants.sigma*deltaEPSCS*mean_T_power4)/self.results['D'] self.results['LW_cre'] = self.results['LW']-self.results['LW_clearsky'] ''' Calculate total SW term without decomposition ''' def calcSWterms(self): self.dataset1['SW'] = self.dataset1['rsdt']-self.dataset1['rsut'] self.dataset2['SW'] = self.dataset2['rsdt']-self.dataset2['rsut'] self.dataset1['SW_cs'] = self.dataset1['rsdt']-self.dataset1['rsutcs'] self.dataset2['SW_cs'] = self.dataset2['rsdt']-self.dataset2['rsutcs'] self.results['SW'] = self.calcDelta('SW')/self.results['D'] self.results['SW_cs'] = self.calcDelta('SW_cs')/self.results['D'] self.results['SW_cloud'] = self.results['SW']-self.results['SW_cs'] #self.results['SW_clearsky'] = self.calcDelta('rsutcs')/self.results['D'] #self.results['SW_cre'] = self.results['SW_total']-self.results['SW_clearsky'] ''' calculate change in surface fluxes ''' def calcSURFterm(self): self.surfHelp(self.dataset1) self.surfHelp(self.dataset2) self.results['SURF'] = -self.calcDelta('NETSURF')/self.results['D'] self.results['NETSURF1'] = self.dataset1['NETSURF'] self.results['NETSURF2'] = self.dataset2['NETSURF'] ''' calculate convergence term (horizontal energy transport) Calculated as residual from other terms ''' def calcCONVterm(self,dataset): dataset['CONV'] = dataset['rlut']+dataset['NETSURF']-dataset['SW'] def calcCONVterms(self): self.calcCONVterm(self.dataset1) self.calcCONVterm(self.dataset2) self.results['CONV'] = self.calcDelta('CONV')/self.results['D']