- Timestamp:
- Nov 12, 2024, 5:52:49 PM (12 days ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/utilities/photochemistry/photochem_postproc.py
r3431 r3511 15 15 warnings.filterwarnings("ignore", message="The following kwargs were not used by contour: 'lat'") 16 16 17 M = {'co2':44, # Molar masses18 'o':16, # TODO: automatic parser19 'o1d':16,20 'o2':32,21 'o3':48,22 'h':1,23 'h2':2,24 'oh':17,25 'h2o_vap':18,26 'ho2':33,27 'h2o2':34,28 'co':28,29 'cho':29,30 'ch2o':30}31 32 17 background = 'co2' # background gas of the atmosphere 33 18 … … 46 31 Methods 47 32 ------- 48 get_ profile(field,**kw)49 Get p rofile of a field (either local or averaged)33 get_subset(field,**kw) 34 Get pa subset at fixed given coordinate of the data 50 35 plot_meridional_slice(field,logcb,labelcb,**kw) 51 36 Plot a meridional slice of a field 52 37 plot_time_evolution(field,logcb,labelcb,**kw) 53 38 Plot time evolution of a field (as a time vs altitude contour plot) 39 plot_time_series(field,lat,lon,alt,logy) 40 Plot time series of a field (at a specific location (lon, lat, alt), averged, or a combination thereof) 41 plot_atlas(field,t,alt,logcb,labelcb) 42 Plot atlas of a field 54 43 plot_profile(field,**kw) 55 44 Plot a profile of a field … … 80 69 self.data[field] = value 81 70 82 def get_profile(self,field,**kw): 83 """ Get profile of a field (either local or averaged) 84 85 Parameters 86 ---------- 87 field : str 88 Field name 89 t : float (optional) 90 Time at which to select (if nothing specified use time-average) 91 lat : float (optional) 92 Latitude at which to select (if nothing specified use area-weighted meridional average) 93 lon : float (optional) 94 Longitude at which to select (if nothing specified use zonal average) 95 96 """ 97 98 if self['latitude'].size == 1: 99 # 1D 100 return self[field][:,0,0] 101 else: 102 # 3D 103 if 'lat' in kw: 104 if 'lon' in kw: 105 return self[field].sel(latitude=kw['lat'],method='nearest').sel(longitude=kw['lon'],method='nearest') 106 else: 107 return self[field].sel(latitude=kw['lat'],method='nearest').mean(dim='longitude') 108 else: 109 # Latitude-averaged profile: need to weight by the grid cell surface area 110 if 'lon' in kw: 111 return (self['aire']*self[field]/self['aire'].mean(dim='latitude')).mean(dim='latitude').sel(longitude=kw['lon'],method='nearest') 112 else: 113 return (self['aire']*self[field]/self['aire'].mean(dim='latitude').mean(dim='longitude')).mean(dim='latitude').mean(dim='longitude') 71 def __area_weight__(self,data_array): 72 return self['aire']*data_array/self['aire'].mean('latitude') 114 73 115 74 def get_subset(self,field='all',**kw): 116 75 """ Get a subset at fixed given coordinate of the data 76 77 Can also average over a given dimension. In this case, 78 the meridional average is area-weighted. 117 79 118 80 Parameters … … 133 95 Altitude of the slice. If nothing 134 96 specified, use time-average 135 136 Raise 137 ----- 138 Slice direction not provided 139 """ 140 if len(kw) == 0: 141 raise Exception('Slice direction not provided') 142 97 """ 143 98 if field == 'all': 144 99 data_subset = self.data … … 147 102 148 103 if 't' in kw: 149 data_subset = data_subset.sel(Time=kw['t'],method='nearest') 104 if kw['t'] == 'avg': 105 data_subset = data_subset.mean(dim='Time') 106 else: 107 data_subset = data_subset.sel(Time=kw['t'],method='nearest') 108 if 'lat' in kw: 109 if kw['lat'] == 'avg': 110 data_subset = self.__area_weight__(data_subset).mean(dim='latitude') 111 else: 112 data_subset = data_subset.sel(latitude=kw['lat'],method='nearest') 150 113 if 'lon' in kw: 151 data_subset = data_subset.sel(longitude=kw['lon'],method='nearest') 152 if 'lat' in kw: 153 data_subset = data_subset.sel(latitude=kw['lat'],method='nearest') 114 if kw['lon'] == 'avg': 115 data_subset = data_subset.mean(dim='longitude') 116 else: 117 data_subset = data_subset.sel(longitude=kw['lon'],method='nearest') 154 118 if 'alt' in kw: 155 data_subset = data_subset.sel(altitude=kw['alt'],method='nearest') 119 if kw['alt'] == 'avg': 120 data_subset = data_subset.mean(dim='altitude') 121 else: 122 data_subset = data_subset.sel(altitude=kw['alt'],method='nearest') 156 123 157 124 return data_subset … … 178 145 raise Exception('Trying to plot a meridional slice of a 1D simulation') 179 146 180 meridional_slice = self [field]181 182 if t == 'avg':183 meridional_slice = meridional_slice.mean(dim='Time')147 meridional_slice = self.get_subset(field,t=t,lon=lon) 148 149 if logcb: 150 plt.contourf(meridional_slice['latitude'],meridional_slice['altitude'],meridional_slice,locator=tk.LogLocator(),**plt_kw) 184 151 else: 185 meridional_slice = meridional_slice.sel(Time=t,method='nearest') 186 if lon == 'avg': 187 meridional_slice = meridional_slice.mean(dim='longitude') 188 else: 189 meridional_slice = meridional_slice.sel(longitude=lon,method='nearest') 190 191 if logcb: 192 plt.contourf(self['latitude'],self['altitude'],meridional_slice,locator=tk.LogLocator(),**plt_kw) 193 else: 194 plt.contourf(self['latitude'],self['altitude'],meridional_slice,**plt_kw) 152 plt.contourf(meridional_slice['latitude'],meridional_slice['altitude'],meridional_slice,**plt_kw) 195 153 196 154 plt.colorbar(label=field if labelcb==None else labelcb) … … 216 174 """ 217 175 218 time_evolution = self [field]219 220 if l at == 'avg':221 time_evolution = (self['aire']*time_evolution/self['aire'].mean(dim='latitude')).mean(dim='latitude')176 time_evolution = self.get_subset(field,lon=lon,lat=lat) 177 178 if logcb: 179 plt.contourf(time_evolution['Time'],time_evolution['altitude'],time_evolution.T,locator=tk.LogLocator(),**plt_kw) 222 180 else: 223 time_evolution = time_evolution.sel(latitude=lat,method='nearest') 224 if lon == 'avg': 225 time_evolution = time_evolution.mean(dim='longitude') 226 else: 227 time_evolution = time_evolution.sel(longitude=lon,method='nearest') 228 229 if logcb: 230 plt.contourf(self['Time'],self['altitude'],time_evolution.T,locator=tk.LogLocator(),**plt_kw) 231 else: 232 plt.contourf(self['Time'],self['altitude'],time_evolution.T,**plt_kw) 181 plt.contourf(time_evolution['Time'],time_evolution['altitude'],time_evolution.T,**plt_kw) 233 182 234 183 plt.colorbar(label=field if labelcb==None else labelcb) 235 184 plt.xlabel('time [day]') 236 185 plt.ylabel('altitude [km]') 186 187 def plot_time_series(self,field,lat='avg',lon='avg',alt='avg',logy=False,**plt_kw): 188 """ Plot time series of a field (at a specific location (lon, lat, alt), averged, or a combination thereof) 189 190 Parameters 191 ---------- 192 field : str 193 Field name to plot 194 lat : float (optional) 195 Latitude at which to plot (if nothing specified use area-weighted meridional average) 196 lon : float (optional) 197 Longitude at which to plot (if nothing specified use zonal average) 198 logy : bool (optional) 199 Use logarithmic y-axis 200 matplotlib plot keyword arguments 201 """ 202 203 time_series = self.get_subset(field,lon=lon,lat=lat,alt=alt) 204 205 if not 'label' in plt_kw: 206 plt_kw['label'] = self[field].units 207 if logy: 208 plt.semilogy(time_series['Time'],time_series,**plt_kw) 209 else: 210 plt.plot(time_series['Time'],time_series,**plt_kw) 211 212 plt.xlabel('time [day]') 213 plt.ylabel(field+' ['+self[field].units+']') 237 214 238 215 def plot_atlas(self,field,t='avg',alt='avg',logcb=False,labelcb=None,**plt_kw): … … 254 231 """ 255 232 256 atlas = self[field] 257 258 if t == 'avg': 259 atlas = atlas.mean(dim='Time') 233 if 'altitude' in self[field].dims: 234 atlas = self.get_subset(field,t=t,alt=alt) 260 235 else: 261 atlas = atlas.sel(Time=t,method='nearest') 262 263 if 'altitude' in atlas.coords: 264 if alt == 'avg': 265 atlas = atlas.mean(dim='altitude') 266 else: 267 atlas = atlas.sel(altitude=alt,method='nearest') 236 atlas = self.get_subset(field,t=t) 268 237 269 238 if logcb: 270 plt.contourf( self['longitude'],self['latitude'],atlas,locator=tk.LogLocator(),**plt_kw)239 plt.contourf(atlas['longitude'],atlas['latitude'],atlas,locator=tk.LogLocator(),**plt_kw) 271 240 else: 272 plt.contourf( self['longitude'],self['latitude'],atlas,**plt_kw)241 plt.contourf(atlas['longitude'],atlas['latitude'],atlas,**plt_kw) 273 242 274 243 plt.colorbar(label=field if labelcb==None else labelcb) … … 294 263 """ 295 264 296 profile = self [field]265 profile = self.get_subset(field,t=t,lon=lon,lat=lat) 297 266 298 if t == 'avg':299 p rofile = profile.mean(dim='Time')267 if logx: 268 plt.semilogx(profile,profile['altitude'],**plt_kw) 300 269 else: 301 profile = profile.sel(Time=t,method='nearest') 302 303 if self['latitude'].size > 1: 304 305 if lat == 'avg': 306 profile = (self['aire']*profile/self['aire'].mean(dim='latitude')).mean(dim='latitude') 307 else: 308 profile = profile.sel(latitude=lat,method='nearest') 309 if lon == 'avg': 310 profile = profile.mean(dim='longitude') 311 else: 312 profile = profile.sel(longitude=lon,method='nearest') 313 314 if logx: 315 plt.semilogx(profile,self['altitude'],**plt_kw) 316 else: 317 plt.plot(profile,self['altitude'],**plt_kw) 270 plt.plot(profile,profile['altitude'],**plt_kw) 318 271 319 272 plt.xlabel(field+' ['+self[field].units+']') 320 273 plt.ylabel('altitude [km]') 274 275 def to_chempath(self,t,dt,filename_suffix='_chempath',lon='avg',lat='avg',alt='avg'): 276 """ Create files redable by the chemical path analyzer chempath (DOI: 10.5194/gmd-2024-163) 277 278 Parameters 279 ---------- 280 t : float 281 Initial time 282 dt : float 283 Timestep 284 filename_suffix : str (optional) 285 Suffix to append to the files being created (species.txt, 286 reactions.txt, model_time.dat, rates.dat, concentrations.dat) 287 lat : float (optional) 288 Latitude at which to plot (if nothing specified use area-weighted meridional average) 289 lon : float (optional) 290 Longitude at which to plot (if nothing specified use zonal average) 291 alt : float (optional) 292 Altitude of the slice. If nothing 293 matplotlib's plot / semilogx keyword arguments 294 """ 295 296 # Save pecies list in chempath format 297 with open(self.path+'/species'+filename_suffix+'.txt', 'w') as fp: 298 for sp in self.species: 299 fp.write("%s\n" % sp) 300 301 # Save reactions list in chempath format 302 with open(self.path+'/reactions'+filename_suffix+'.txt', 'w') as fp: 303 for r in self.reactions: 304 current_formula = self.reactions[r].formula 305 current_formula = current_formula.replace(' ','') 306 current_formula = current_formula.replace('->','=') 307 fp.write("%s\n" % current_formula) 308 309 # Save rates, concentrations and times in chempath format 310 rates = np.array([],dtype=np.float128) 311 conc = np.array([],dtype=np.float128) 312 conc_dt = np.array([],dtype=np.float128) 313 times = np.array([t,t+dt],dtype=np.float128) 314 315 out = self.get_subset(t=t,lon=lon,lat=lat,alt=alt) 316 for r in self.reactions: 317 rates = np.append(rates,out[f'rate ({r})']) 318 for sp in self.species: 319 conc = np.append(conc,out[f'{sp} vmr']) 320 321 out = self.get_subset(t=t+dt,lon=lon,lat=lat,alt=alt) 322 for sp in self.species: 323 conc_dt = np.append(conc_dt,out[f'{sp} vmr']) 324 conc = np.vstack((conc,conc_dt)) 325 326 times.tofile(self.path+'/model_time'+filename_suffix+'.dat') 327 rates.tofile(self.path+'/rates'+filename_suffix+'.dat') 328 conc.tofile(self.path+'/concentrations'+filename_suffix+'.dat') 321 329 322 330 class reaction: … … 603 611 if reactions == 'read': 604 612 reactions = read_reactfile(s.path) 613 s.tracers = read_traceurs(s.path) 605 614 s.species = [] 606 615 s.reactions = {} # reactions dict will be merged at the end … … 609 618 for r in reactions: 610 619 for sp in reactions[r].reactants: 620 if not sp in s.tracers or not s.tracers[sp].is_chim: 621 raise Warning(sp, 'is not recorded as a chimical tracer') 611 622 if not sp in s.species: 612 623 s.species.append(sp) 613 624 for sp in reactions[r].products: 625 if not sp in s.tracers or not s.tracers[sp].is_chim: 626 raise Warning(sp, 'is not recorded as a chimical tracer') 614 627 if not sp in s.species: 615 628 s.species.append(sp) … … 622 635 for sp in s.species: 623 636 # volume mixing ratios 624 s[sp+' vmr'] = s[sp] * M['co2'] / M[sp]637 s[sp+' vmr'] = s[sp] * s.tracers['co2'].M / s.tracers[sp].M 625 638 s[sp+' vmr'] = s[sp+' vmr'].assign_attrs({'units':'m^3/m^3'}) 626 639 # molecular densities … … 678 691 return s 679 692 693 class trac: 694 """ Tracer class 695 696 A class to store useful informations about simulations tracers 697 698 Attributes 699 ---------- 700 name : string 701 Tracer's name 702 M : float 703 Tracer's molar mass [g/mol] 704 is_chim : bool 705 Is it a photochemical tracer? 706 """ 707 def __init__(self,name,M,is_chim): 708 self.name = name 709 self.M = M 710 self.is_chim = is_chim 711 712 def read_traceurs(path): 713 """ Read the traceurs of a simulation 714 715 Parameters 716 ---------- 717 path : string 718 path to simulation 719 720 Returns 721 ------- 722 dict 723 dictionnaries of all tracers 724 """ 725 tracdict = {} 726 with open(path+'/traceur.def') as tracfile: 727 728 for iline,line in enumerate(tracfile): 729 730 # Empty line 731 if len(line.split()) == 0: 732 continue 733 734 # First line 735 elif iline == 0: 736 if not '#ModernTrac-v1' in line: 737 raise Exception('Can only read modern traceur.def') 738 continue 739 740 # Second line (number of tracers) 741 elif iline == 1: 742 ntrac = int(line) 743 continue 744 745 # Commented line 746 elif line[0] == '!': 747 continue 748 749 # Regular entry 750 else: 751 line = line.split() 752 name = line[0] 753 is_chim = False 754 755 for param in line[1:]: 756 if param[:4] == 'mmol': 757 M = float(param[5:]) 758 elif param[:7] == 'is_chim': 759 is_chim = bool(float(param[8:])) 760 761 tracdict[name] = trac(name,M,is_chim) 762 763 if len(tracdict) != ntrac: 764 raise Exception('Mismatch between announced an read number of tracers') 765 766 return tracdict
Note: See TracChangeset
for help on using the changeset viewer.
