diff --git a/BC/VPnutrients_CO2Hg.nc b/BC/VPnutrients_CO2Hg.nc new file mode 100644 index 0000000..5d0d392 Binary files /dev/null and b/BC/VPnutrients_CO2Hg.nc differ diff --git a/BC/atlantic/atlantic_generator.py b/BC/atlantic/atlantic_generator.py index b29332f..42f6b8a 100644 --- a/BC/atlantic/atlantic_generator.py +++ b/BC/atlantic/atlantic_generator.py @@ -45,7 +45,7 @@ def argument(): os.system("mkdir -p " + OUTPUTDIR) -OpenMask=Mask(MASKFILE) +OpenMask=Mask.from_file(MASKFILE) nav_lon=OpenMask.xlevels B=OpenMask.bathymetry_in_cells() nav_lev=OpenMask.zlevels @@ -69,7 +69,7 @@ def argument(): mydtype=[(var,np.float32) for var in ALLVARS] BOUNDARY_CONCENTRATION= np.zeros((jpk,),dtype=mydtype) -VARLIST=['N1p', 'N3n', 'N5s','O2o','O3c','O3h'] +VARLIST=['N1p', 'N3n', 'N5s','O2o','O3c','O3h','Hg2','Hg0','MHg','DHg','P1h','P2h','P3h','P4h','Z6h','Z5h','Z4h','Z3h'] for var in ALLVARS: if var in VARLIST: diff --git a/BC/atm_gen.py b/BC/atm_gen.py index 2fe721d..870aca7 100644 --- a/BC/atm_gen.py +++ b/BC/atm_gen.py @@ -6,7 +6,7 @@ filename="atmDep_16sub_openAndCoast_winAndSum.txt" maskfile="/Users/gbolzon/eclipse-workspace/preproc/IC/meshmask.nc" -TheMask=Mask(maskfile) +TheMask=Mask.from_file(maskfile) jpk, jpj, jpi = TheMask.shape Mask0 = TheMask.cut_at_level(0) mask200= TheMask.mask_at_level(200) @@ -21,7 +21,7 @@ dtype = [(sub.name, bool) for sub in OGS.Pred] SUB = np.zeros((jpj,jpi),dtype=dtype) for sub in OGS.Pred: - sbmask = SubMask(sub,maskobject=Mask0).mask + sbmask = SubMask(sub, Mask0).mask SUB[sub.name] = sbmask[0,:,:] N3n_win = np.zeros((jpj,jpi),np.float32) diff --git a/BC/bclib/atmosphere.py b/BC/bclib/atmosphere.py index 0d151dd..42094c3 100644 --- a/BC/bclib/atmosphere.py +++ b/BC/bclib/atmosphere.py @@ -19,32 +19,42 @@ def __init__(self,mask,conf): logging.info("Atmosphere start calculation") _,jpj,jpi = mask.shape - mask0 = mask.cut_at_level(0) EAS=ComposedBasin('eas',[V2.eas,V2.adr1,V2.adr2,V2.aeg],'East with marginal seas') # Here we use basin objet features to run faster EAS.region.border_longitudes = [ 9,36] EAS.region.border_latitudes = [30,46] V2.wes.region.border_longitudes = [-6,18] V2.wes.region.border_latitudes = [32,45] + V2.med.region.border_longitudes = [-6,36] + V2.med.region.border_latitudes = [30,46] #--------------------------------------------- print("wes") - wes = SubMask(V2.wes,maskobject=mask0).mask_at_level(0) + wes = SubMask(V2.wes, mask) print("eas") - eas = SubMask(EAS, maskobject=mask0).mask_at_level(0) + eas = SubMask(EAS, mask) + # print("med") + # med = SubMask(V2.med,mask0).mask_at_level(0) print("done") - e3t = mask0.dz - Volume = mask0.area*e3t[0] - Nwes = Volume[wes].sum() - Neas = Volume[eas].sum() + e3t = mask.dz + Volume = mask.area*e3t[0] + Nwes = Volume[wes[0]].sum() + Neas = Volume[eas[0]].sum() + # Nmed = Volume[med].sum() self.nitrate = np.zeros((jpj,jpi),np.float32); self.phosphate = np.zeros((jpj,jpi),np.float32); + self.HgII = np.zeros((jpj,jpi),np.float32); + self.MMHg = np.zeros((jpj,jpi),np.float32); + self.Hg0atm = np.zeros((jpj,jpi),np.float32); - self.nitrate[wes] = conf.n3n_wes/Nwes - self.phosphate[wes]= conf.po4_wes/Nwes - self.nitrate[eas] = conf.n3n_eas/Neas - self.phosphate[eas]= conf.po4_eas/Neas + self.nitrate[wes[0]] = conf.n3n_wes/Nwes + self.phosphate[wes[0]]= conf.po4_wes/Nwes + self.nitrate[eas[0]] = conf.n3n_eas/Neas + self.phosphate[eas[0]]= conf.po4_eas/Neas + self.HgII[eas[0]] = conf.HgII_eas/Neas + self.MMHg[eas[0]] = conf.MMHg_eas/Neas + self.Hg0atm[eas[0]] = conf.Hg0atm_eas/Neas logging.info("Atmosphere finish calculation") @@ -54,19 +64,17 @@ def check(self,mask): _,jpj,jpi = mask.shape n = 1./14 p = 1./31 + h=1./200.59 + vol_layer0 = mask.area * mask.dz[0] totP = .0 totN = .0 + totHgII = .0 + totMMHg = .0 totN_KTy = .0 totP_KTy = .0 - for jj in range(jpj): - for ji in range(jpi): - VolCell1=mask.area[jj,ji]*mask.dz[0] - totN = totN + self.nitrate[jj,ji] *VolCell1 - totP = totP + self.phosphate[jj,ji]*VolCell1 - totN_KTy = totN_KTy + self.nitrate[jj,ji] *(1.e-3/n)*VolCell1 - totP_KTy = totP_KTy + self.phosphate[jj,ji]*(1.e-3/p)*VolCell1 - return totN,totP, totN_KTy, totP_KTy - + totHgII_KTy=.0 + totMMHg_KTy=.0 + return totN,totP,totHgII,totMMHg, totN_KTy, totP_KTy,totHgII_KTy, totMMHg_KTy def write_netcdf(self,mask,outdir,concentration_dimension="Volume"): @@ -79,27 +87,34 @@ def write_netcdf(self,mask,outdir,concentration_dimension="Volume"): assert concentration_dimension=="Volume" or concentration_dimension =="Area" # Check is done using volume concentrations - totN,totP, totN_KTy, totP_KTy = self.check(mask) + totN,totP,totHgII, totMMHg, totN_KTy, totP_KTy,totHgII_KTy, totMMHg_KTy = self.check(mask) #conversion from Mmol/y to mmol/s + # Hg conversion from Mmol/y to nmol/s w = 1.0e+09 + w_hg = 1.0e+15 t = 1./(365 * 86400) cn = w*t + cn_hg=w_hg*t self.nitrate = self.nitrate*cn self.phosphate = self.phosphate*cn + self.HgII = self.HgII*cn_hg + self.MMHg = self.MMHg*cn_hg if concentration_dimension == "Volume": units ="mmol/m3" + hunits ="nmol/m3" comment="transport model can use these atmospherical contributions as they are" if concentration_dimension =="Area": units = "mmol/m2" + hunits = "nmol/m2" comment = "transport model can use thes atmospherical contributions by dividing them by the depth of the first cell" self.nitrate = self.nitrate * mask.dz[0] self.phosphate = self.phosphate * mask.dz[0] - - + self.HgII = self.HgII * mask.dz[0] + self.MMHg = self.MMHg * mask.dz[0] # now is ready to dump - fileOUT = outdir + "ATM_yyyy0630-00:00:00.nc" + fileOUT = outdir + "ATM_yyyy0630-00:00:00_LOWHg.nc" _,jpj,jpi = mask.shape ncfile = nc.Dataset(fileOUT, "w", format="NETCDF4") @@ -111,11 +126,34 @@ def write_netcdf(self,mask,outdir,concentration_dimension="Volume"): p = ncfile.createVariable('atm_N1p','f', ('lat','lon')) p[:] = self.phosphate[:] setattr(p,"units", units) + h1 = ncfile.createVariable('atm_Hg2','f', ('lat','lon')) + h1[:] = self.HgII[:] + setattr(h1,"units", hunits) + h2 = ncfile.createVariable('atm_MMHg','f', ('lat','lon')) + h2[:] = self.MMHg[:] + setattr(h2,"units", hunits) setattr(ncfile, 'ATM_P_MassBalance_kTON_y', totP_KTy) setattr(ncfile, 'ATM_N_MassBalance_kTON_y', totN_KTy) + setattr(ncfile, 'ATM_HgII_MassBalance_kTON_y', totHgII_KTy) + setattr(ncfile, 'ATM_MMHg_MassBalance_kTON_y', totMMHg_KTy) setattr(ncfile, 'ATM_P_MassBalance_Mmol_y', totP) setattr(ncfile, 'ATM_N_MassBalance_Mmol_y', totN) + setattr(ncfile, 'ATM_HgII_MassBalance_Mmol_y', totHgII) + setattr(ncfile, 'ATM_MMHg_MassBalance_Mmol_y', totMMHg) setattr(ncfile, 'comments', comment) ncfile.close() logging.info("Atmosphere Netcdf written.") + + fileOUT = outdir + "atm_Hg0_yyyy0107-00:00:00.nc" + _,jpj,jpi = mask.shape + + ncfile = nc.Dataset(fileOUT, "w", format="NETCDF4") + ncfile.createDimension('lon', jpi) + ncfile.createDimension('lat', jpj) + h = ncfile.createVariable('atm_Hg0','f', ('lat','lon')) + h[:] = self.Hg0atm[:] + setattr(h,"units", "ug m-3") + ncfile.close() + + logging.info("Hg0 at mNetcdf written.") diff --git a/BC/bclib/bounmask.py b/BC/bclib/bounmask.py index 23b4021..3527c0f 100644 --- a/BC/bclib/bounmask.py +++ b/BC/bclib/bounmask.py @@ -121,7 +121,7 @@ def write_netcdf(self,mask, all_variables=True): navlev_wnc = ncfile.createVariable('nav_lev', 'f', 'z') navlev_wnc[:] = mask.zlevels for jn in range(nudg): - corrf =[1.,1.,1.,1.,1.01,1.01,1.]; + corrf =[1.,1.,1.,1.,1.01,1.01,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.]; aux= self.resto[jn,:,:,:]*corrf[jn]; aux[~mask.mask] = 1.e+20 resto_wnc = ncfile.createVariable("re" + vnudg[jn][0], 'f4', ('time','z','y','x')) @@ -145,9 +145,9 @@ def write_netcdf(self,mask, all_variables=True): if __name__ == '__main__': from bitsea.commons.mask import Mask import config as conf - conf.file_mask="../meshmask_872.nc" +# conf.file_mask="../meshmask_872.nc" conf.file_bmask="bounmask.nc" - TheMask = Mask(conf.file_mask) + TheMask = Mask.from_file(conf.file_mask) B=bounmask(conf) B.generate(TheMask) B.write_netcdf(TheMask) diff --git a/BC/bclib/gib.py b/BC/bclib/gib.py index 96ae8d5..89eafea 100644 --- a/BC/bclib/gib.py +++ b/BC/bclib/gib.py @@ -27,6 +27,18 @@ def nutrient_dataset_by_index(self,jn): if jn==6 : size_nutrients = self.gibilterra.phos.shape return np.ones(size_nutrients,np.float64)*0.0025 # G. Cossarini estimate + if jn==7 : return self.gibilterra.me0 + if jn==8 : return self.gibilterra.me2 + if jn==9 : return self.gibilterra.mmh + if jn==10: return self.gibilterra.dmh + if jn==11: return self.gibilterra.pl1 + if jn==12: return self.gibilterra.pl2 + if jn==13: return self.gibilterra.pl3 + if jn==14: return self.gibilterra.pl4 + if jn==15: return self.gibilterra.zo6 + if jn==16: return self.gibilterra.zo5 + if jn==17: return self.gibilterra.zo4 + if jn==18: return self.gibilterra.zo3 def generate(self, mask, bounmask_obj, all_variables=True): @@ -63,6 +75,7 @@ def generate(self, mask, bounmask_obj, all_variables=True): GIB_matrix = self.nutrient_dataset_by_index(jn) aux = bounmask_obj.resto[jn, ...] isNudg = (aux != 0) & (mask.mask) + print (jn, nudg, vnudg[jn][0]) GIB_matrix[time, ~isNudg] = -1. GIB_matrix[time, ~mask.mask] = 1.e+20 @@ -94,9 +107,10 @@ def generate(self, mask, bounmask_obj, all_variables=True): from bitsea.commons.mask import Mask import config as conf conf.file_nutrients = "../"+ conf.file_nutrients - conf.file_mask="../masks/meshmask_872.nc" +# conf.file_mask="../masks/meshmask_872.nc" + conf.file_mask="../meshmask.nc" conf.dir_out = "../out" - TheMask = Mask(conf.file_mask) + TheMask = Mask.from_file(conf.file_mask) GIB = gib(conf,TheMask) from bounmask import bounmask BOUN=bounmask(conf) diff --git a/BC/bclib/lateral.py b/BC/bclib/lateral.py index 1c4fa59..53aa2d7 100644 --- a/BC/bclib/lateral.py +++ b/BC/bclib/lateral.py @@ -25,6 +25,8 @@ def _extract_information(self): def _convert_information(self,mask): jpt = len(self.season) + print('self.season',self.season) + print('--------------------') jpk, jpj, jpi = mask.shape size_nutrients = (jpt,jpk,jpj, jpi) self.phos = np.zeros(size_nutrients, np.float64) @@ -33,15 +35,42 @@ def _convert_information(self,mask): self.sica = np.zeros(size_nutrients, np.float64) self.dic = np.zeros(size_nutrients, np.float64) self.alk = np.zeros(size_nutrients, np.float64) - tmask4d = np.zeros(size_nutrients, np.bool) + self.me0 = np.zeros(size_nutrients, np.float64) + self.me2 = np.zeros(size_nutrients, np.float64) + self.mmh = np.zeros(size_nutrients, np.float64) + self.dmh = np.zeros(size_nutrients, np.float64) + self.pl1 = np.zeros(size_nutrients, np.float64) + self.pl2 = np.zeros(size_nutrients, np.float64) + self.pl3 = np.zeros(size_nutrients, np.float64) + self.pl4 = np.zeros(size_nutrients, np.float64) + self.zo6 = np.zeros(size_nutrients, np.float64) + self.zo5 = np.zeros(size_nutrients, np.float64) + self.zo4 = np.zeros(size_nutrients, np.float64) + self.zo3 = np.zeros(size_nutrients, np.float64) + tmask4d = np.zeros(size_nutrients, np.bool_) + # tmask4d = np.zeros(size_nutrients, np.bool) nav_lev = mask.zlevels n_lev = jpk + print('n_lev',n_lev) + #n_lev2 = jpk vp_phos = np.zeros((jpt,n_lev)); vp_ntra = np.zeros((jpt,n_lev)); vp_sica = np.zeros((jpt,n_lev)); vp_dox = np.zeros((jpt,n_lev)); vp_dic = np.zeros((jpt,n_lev)); vp_alk = np.zeros((jpt,n_lev)); + vp_Hg0 = np.zeros((jpt,n_lev)); + vp_Hg2 = np.zeros((jpt,n_lev)); + vp_MHg = np.zeros((jpt,n_lev)); + vp_DHg = np.zeros((jpt,n_lev)); + vp_P1h = np.zeros((jpt,n_lev)); + vp_P2h = np.zeros((jpt,n_lev)); + vp_P3h = np.zeros((jpt,n_lev)); + vp_P4h = np.zeros((jpt,n_lev)); + vp_Z6h = np.zeros((jpt,n_lev)); + vp_Z5h = np.zeros((jpt,n_lev)); + vp_Z4h = np.zeros((jpt,n_lev)); + vp_Z3h = np.zeros((jpt,n_lev)); nav_lev_in=np.concatenate(([0],self.lev1,[5500])) @@ -49,17 +78,38 @@ def _convert_information(self,mask): vp_ntra_in=np.column_stack((self.N3n[:,0], self.N3n, self.N3n[:,-2])) vp_sica_in=np.column_stack((self.N5s[:,0], self.N5s, self.N5s[:,-2])) vp_dox_in =np.column_stack((self.O2o[:,0], self.O2o, self.O2o[:,-2])) - end = len(self.ALK) +####### vp_Hg0_in =np.column_stack((self.Hg0[:,0], self.Hg0, self.Hg0[:,-2])) +#### end of two dimensional vars - 1D vars + end = len(self.ALK) + print('END',end) nav_lev_in2 = np.append([0],self.lev2.T) nav_lev_in2 = np.append(nav_lev_in2,[5500]) - - vp_dic_in = np.append(self.DIC[1],self.DIC[:].T) + print('NAVin2',nav_lev_in2 ) + vp_dic_in = np.append(self.DIC[0],self.DIC[:].T) vp_dic_in = np.append(vp_dic_in,self.DIC[end-1]) - vp_alk_in = np.append(self.ALK[1],self.ALK[:].T) + vp_alk_in = np.append(self.ALK[0],self.ALK[:].T) vp_alk_in = np.append(vp_alk_in,self.ALK[end-1]); + print('vp_alk_in0',vp_alk_in) + + end2=len(self.Hg0) + print('end for Hg0',end2) + nav_lev_in3 = self.lev3.T + print('NAVin3',nav_lev_in3) + + vp_Hg0_in = self.Hg0[:].T + # vp_Hg0 = np.zeros((jpt,n_lev)) + print('self.Hg0',self.Hg0) + # print('vp_Hg0_in1',vp_Hg0_in) + # print('self.Hg0[0]',self.Hg0[0]) + # print('self.Hg0[0].T',self.Hg0[0].T) + + vp_Hg2_in = self.Hg2[:].T + vp_MHg_in = self.MMHg[:].T + vp_DHg_in = self.DMHg[:].T + for i in range(jpt): jj= (~ np.isnan(vp_phos_in[i,:])) | ( vp_phos_in[i,:] < 1e19 ) vp_phos[i,:] = np.interp(nav_lev, nav_lev_in[jj], vp_phos_in[i,jj]) @@ -74,7 +124,26 @@ def _convert_information(self,mask): jj=~ np.isnan(vp_alk_in[:]) | ( vp_alk_in[:] < 1e19) vp_alk[i,:] = np.interp(nav_lev,nav_lev_in2[jj],vp_alk_in[jj]) -# + jj=~ np.isnan(vp_Hg0_in[:]) | ( vp_Hg0_in[:] < 1e19) + vp_Hg0[i,:] = np.interp(nav_lev,nav_lev_in3[jj],vp_Hg0_in[jj]) + jj=~ np.isnan(vp_Hg2_in[:]) | ( vp_Hg2_in[:] < 1e19) + vp_Hg2[i,:] = np.interp(nav_lev,nav_lev_in3[jj],vp_Hg2_in[jj]) + jj=~ np.isnan(vp_MHg_in[:]) | ( vp_MHg_in[:] < 1e19) + vp_MHg[i,:] = np.interp(nav_lev,nav_lev_in3[jj],vp_MHg_in[jj]) + jj=~ np.isnan(vp_DHg_in[:]) | ( vp_DHg_in[:] < 1e19) + vp_DHg[i,:] = np.interp(nav_lev,nav_lev_in3[jj],vp_DHg_in[jj]) + + print('vp_Hg0_in[:]',vp_Hg0_in[:].shape) + print('vp_alk_in[:]',vp_alk_in[:].shape) + + print('self.phos',self.phos.shape) + print('self.phos',self.phos.shape) + print('jpts',jpt) + print('n_levs',n_lev) + + print('vp.phos',vp_phos.shape) + print('vp_Hg0',vp_Hg0.shape) + print('tmask4d',tmask4d.shape) # %Loop on time seasonal for jt in range(jpt): for jk in range(n_lev): @@ -84,9 +153,21 @@ def _convert_information(self,mask): self.sica[jt,jk,:,:] = vp_sica[jt,jk] self.dic[ jt,jk,:,:] = vp_dic[ jt,jk] self.alk[ jt,jk,:,:] = vp_alk[ jt,jk] - - - for jt in range(jpt): tmask4d[jt,:,:,:]=mask.mask + self.me0[ jt,jk,:,:] = vp_Hg0[ jt,jk] + self.me2[ jt,jk,:,:] = vp_Hg2[ jt,jk] + self.mmh[ jt,jk,:,:] = vp_MHg[ jt,jk] + self.dmh[ jt,jk,:,:] = vp_DHg[ jt,jk] + self.pl1[ jt,jk,:,:] = vp_P1h[ jt,jk] + self.pl2[ jt,jk,:,:] = vp_P2h[ jt,jk] + self.pl3[ jt,jk,:,:] = vp_P3h[ jt,jk] + self.pl4[ jt,jk,:,:] = vp_P4h[ jt,jk] + self.zo6[ jt,jk,:,:] = vp_Z6h[ jt,jk] + self.zo5[ jt,jk,:,:] = vp_Z5h[ jt,jk] + self.zo4[ jt,jk,:,:] = vp_Z4h[ jt,jk] + self.zo3[ jt,jk,:,:] = vp_Z3h[ jt,jk] + + + for jt in range(jpt): tmask4d[jt,:,:,:]=mask[:] self.phos[~tmask4d]=1.e+20 self.ntra[~tmask4d]=1.e+20 @@ -94,3 +175,7 @@ def _convert_information(self,mask): self.sica[~tmask4d]=1.e+20 self.dic[ ~tmask4d]=1.e+20 self.alk[ ~tmask4d]=1.e+20 + self.me0[ ~tmask4d]=1.e+20 + self.me2[ ~tmask4d]=1.e+20 + self.mmh[ ~tmask4d]=1.e+20 + self.dmh[ ~tmask4d]=1.e+20 diff --git a/BC/bclib/river.py b/BC/bclib/river.py index 0eb2062..dfc7c54 100644 --- a/BC/bclib/river.py +++ b/BC/bclib/river.py @@ -10,18 +10,22 @@ def conversion(var): ''' from mmol or mg to KTONS(N1p,N3n,N5s,O3c,O3h) or Gmol (O2o) + from nmol to KTONS (Hg, MMHg) ''' Conversion={} w= 1.0e-12 n = 14 p = 31 s = 28 + h = 200590000 #/200.59*1000000. Conversion['N1p'] =w*p Conversion['N3n'] =w*n Conversion['N5s'] =w*s Conversion['O3h'] =w Conversion['O3c'] =w Conversion['O2o'] =w + Conversion['Hg2'] =w*h + Conversion['MMHg'] =w*h return Conversion[var] class river(): @@ -254,7 +258,9 @@ def total(self,year): totS = self.river_data["DIS_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; totA = self.river_data["ALK_GmolperYR_NOBLS"][str(year)].sum(axis=1)/12; totD = self.river_data["DIC_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; - return totN,totP,totS,totA,totD + totH = self.river_data["HgII_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; + totM = self.river_data["MMHg_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; + return totN,totP,totS,totA,totD,totH,totM def spread_array(self,array): @@ -283,11 +289,13 @@ def get_monthly_data(self,yearstr,month): A = self.river_data["ALK_GmolperYR_NOBLS"][yearstr][:,month-1] D = self.river_data["DIC_KTperYR_NOBLS" ][yearstr][:,month-1] O = self.river_data["O2o_GmolperYR_NOBLS"][yearstr][:,month-1] + H = self.river_data["HgII_KTperYR_NOBLS"][yearstr][:,month-1] + M = self.river_data["MMHg_KTperYR_NOBLS"][yearstr][:,month-1] DOC = self.river_data["DOC_KTperYR_NOBLS"][yearstr][:,month-1] CDOM = self.river_data["CDOM_KTperYR_NOBLS"][yearstr][:,month-1]*4.0 if (self.nspread==1): - return N,P,S,A,D,O,DOC,CDOM + return N,P,S,A,D,O,H,M,DOC,CDOM else: N_long = self.spread_array(N) P_long = self.spread_array(P) @@ -295,17 +303,19 @@ def get_monthly_data(self,yearstr,month): A_long = self.spread_array(A) D_long = self.spread_array(D) O_long = self.spread_array(O) + H_long = self.spread_array(H) + M_long = self.spread_array(M) DOC_long = self.spread_array(DOC) CDOM_long = self.spread_array(CDOM) - return N_long, P_long, S_long, A_long, D_long, O_long, DOC_long, CDOM_long + return N_long, P_long, S_long, A_long, D_long, O_long,H_long,M_long,DOC_long, CDOM_long - def conversion(self,N,P,S,A,D,O,DOC,CDOM): + def conversion(self,N,P,S,A,D,O,H,M,DOC,CDOM): ''' Performs conversion of all variables - from KT/y to mmol/s or mg/s + from KT/y to mmol/s or mg/s or nmol/s Arguments: Nitrate, Phosphate, Silicate, Alcalinity, Dissoved Inorganic Carbon, Dissolved_Oxygen @@ -321,21 +331,26 @@ def conversion(self,N,P,S,A,D,O,DOC,CDOM): * O * in mmol/s * DOC * mg/s * CDOM * mg/s + * Hg * nmol/s + * MeHg * nmol/s as (nRivers,) numpy array Data are not ready for model, they have to be divided by cell area or cell volume. ''' w= 1.0e+12; + w2= 1.0e+6; t = 1./(365 * 86400) n = 1./14; p = 1./31; s = 1./28; + h = 1./200.59; cn = w*t*n cp = w*t*p cs = w*t*s ca = w*t cc = w*t - return N*cn, P*cp, S*cs,A*ca, D*cc, O*cc, DOC*cc, CDOM*cc + ch = w*w2*h*t + return N*cn, P*cp, S*cs,A*ca, D*cc, O*cc,H*ch,M*ch, DOC*cc, CDOM*cc def generate_monthly_files(self,conf,mask):#, idxt_riv, positions): ''' @@ -363,9 +378,9 @@ def generate_monthly_files(self,conf,mask):#, idxt_riv, positions): for year in range(start_year,end___year): for month in range(1,13): filename = conf.dir_out + "TIN_%d%02d15-00:00:00.nc" %(year, month) - N,P,S,A,D,O,DOC,CDOM = self.get_monthly_data(str(year), month) - N,P,S,A,D,O, DOC, CDOM = self.conversion(N, P, S, A, D, O, DOC, CDOM) - self.dump_file(filename, N/Area, P/Area, S/Area, A/Area, D/Area, O/Area, DOC/Area, CDOM/Area, mask) + N,P,S,A,D,O,H,M,DOC,CDOM = self.get_monthly_data(str(year), month) + N,P,S,A,D,O,H,M,DOC, CDOM = self.conversion(N, P, S, A, D, O,H,M,DOC, CDOM) + self.dump_file(filename, N/Area, P/Area, S/Area, A/Area, D/Area, O/Area,H/Area, M/Area, DOC/Area, CDOM/Area, mask) logging.info("Non climatological TIN file generation : done") def generate_climatological_monthly_files(self,conf,mask): #idxt_riv, positions): @@ -391,13 +406,13 @@ def generate_climatological_monthly_files(self,conf,mask): #idxt_riv, positions) year="yyyy" for month in range(1,13): filename = conf.dir_out+"/TIN_yyyy%02d15-00:00:00.nc" %(month) - N,P,S,A,D,O,DOC,CDOM = self.get_monthly_data(str(year), month) - N,P,S,A,D,O,DOC,CDOM = self.conversion(N, P, S, A, D, O, DOC,CDOM) - self.dump_file(filename, N/Area, P/Area, S/Area, A/Area, D/Area, O/Area, DOC/Area, CDOM/Area, mask) + N,P,S,A,D,O,H,M,DOC,CDOM = self.get_monthly_data(str(year), month) + N,P,S,A,D,O,H,M,DOC,CDOM = self.conversion(N, P, S, A, D, O,H,M, DOC,CDOM) + self.dump_file(filename, N/Area, P/Area, S/Area, A/Area, D/Area, O/Area, H/Area,M/Area,DOC/Area, CDOM/Area, mask) logging.info("Climatological TIN file generation : done") - def dump_file_old(self,filename,N,P,S,A,D,O, idxt_riv,positions): + def dump_file_old(self,filename,N,P,S,A,D,O,H,M, idxt_riv,positions): ''' Writes the single TIN file Variables are dumped as they are, all but positions (incremented by one) @@ -415,6 +430,8 @@ def dump_file_old(self,filename,N,P,S,A,D,O, idxt_riv,positions): riv_a_o3c = ncfile.createVariable('riv_O3c', 'f4', ('riv_idxt',)) riv_a_o3h = ncfile.createVariable('riv_O3h', 'f4', ('riv_idxt',)) riv_a_O2o = ncfile.createVariable('riv_O2o', 'f4', ('riv_idxt',)) + riv_a_Hg2 = ncfile.createVariable('riv_Hg2', 'f4', ('riv_idxt',)) + riv_a_MMHg = ncfile.createVariable('riv_MMHg', 'f4', ('riv_idxt',)) riv_idxt_riv[:] = idxt_riv[:] riv_pos[:,:] = positions+1 @@ -424,9 +441,11 @@ def dump_file_old(self,filename,N,P,S,A,D,O, idxt_riv,positions): riv_a_o3c[:] = D riv_a_o3h[:] = A riv_a_O2o[:] = O + riv_a_Hg2[:] = H + riv_a_MMHg[:] = M ncfile.close() - def dump_file(self,filename,N,P,S,A,D,O,DOC, CDOM, mask): + def dump_file(self,filename,N,P,S,A,D,O,H,M,DOC, CDOM, mask): ''' Writes the single TIN file Variables are dumped as they are, all but positions (incremented by one) @@ -441,6 +460,8 @@ def dump_file(self,filename,N,P,S,A,D,O,DOC, CDOM, mask): riv_a_o3c = ncfile.createVariable('riv_O3c', 'f4', ('lat','lon')) riv_a_o3h = ncfile.createVariable('riv_O3h', 'f4', ('lat','lon')) riv_a_O2o = ncfile.createVariable('riv_O2o', 'f4', ('lat','lon')) + riv_a_Hg2 = ncfile.createVariable('riv_Hg2', 'f4', ('lat','lon')) + riv_a_MMHg = ncfile.createVariable('riv_MMHg', 'f4', ('lat','lon')) riv_a_R3c = ncfile.createVariable('riv_R3c', 'f4', ('lat','lon')) riv_a_R3l = ncfile.createVariable('riv_R3l', 'f4', ('lat','lon')) @@ -450,6 +471,8 @@ def dump_file(self,filename,N,P,S,A,D,O,DOC, CDOM, mask): riv_a_o3c[:] = self.get_map_from_1d_array(D, mask) riv_a_o3h[:] = self.get_map_from_1d_array(A, mask) riv_a_O2o[:] = self.get_map_from_1d_array(O, mask) + riv_a_Hg2[:] = self.get_map_from_1d_array(H, mask) + riv_a_MMHg[:] = self.get_map_from_1d_array(M, mask) riv_a_R3c[:] = self.get_map_from_1d_array(DOC, mask) riv_a_R3l[:] = self.get_map_from_1d_array(CDOM,mask) @@ -459,6 +482,8 @@ def dump_file(self,filename,N,P,S,A,D,O,DOC, CDOM, mask): setattr(riv_a_o3c,'missing_value',np.float32(1.e+20)) setattr(riv_a_o3h,'missing_value',np.float32(1.e+20)) setattr(riv_a_O2o,'missing_value',np.float32(1.e+20)) + setattr(riv_a_Hg2,'missing_value',np.float32(1.e+20)) + setattr(riv_a_MMHg,'missing_value',np.float32(1.e+20)) setattr(riv_a_R3c,'missing_value',np.float32(1.e+20)) setattr(riv_a_R3l,'missing_value',np.float32(1.e+20)) ncfile.close() @@ -507,7 +532,7 @@ def load_from_file(self,filename,var): if __name__=="__main__": from bitsea.commons.mask import Mask import config as conf - TheMask = Mask(conf.file_mask) + TheMask = Mask.from_file(conf.file_mask) R = river(conf) R.modularize(conf) R.gen_map_indexes(TheMask) diff --git a/BC/bclib/river_Hg.py b/BC/bclib/river_Hg.py new file mode 100644 index 0000000..aafec2f --- /dev/null +++ b/BC/bclib/river_Hg.py @@ -0,0 +1,549 @@ +# -*- coding: utf-8 -*- + +# See /gss/gss_work/DRES_OGS_BiGe/gcossarini/River_Generation to know about excel generation + +import numpy as np +from bclib import excel_reader +import logging +import netCDF4 + +def conversion(var): + ''' + from mmol or mg to KTONS(N1p,N3n,N5s,O3c,O3h) or Gmol (O2o) + ''' + Conversion={} + w= 1.0e-12 + n = 14 + p = 31 + s = 28 + h = 200590000 #/200.59*1000000. + Conversion['N1p'] =w*p + Conversion['N3n'] =w*n + Conversion['N5s'] =w*s + Conversion['O3h'] =w + Conversion['O3c'] =w + Conversion['O2o'] =w + Conversion['Hg2'] =h + Conversion['MMHg'] =h + return Conversion[var] + +class river(): + + def __init__(self,conf): + ''' + + Reads data from excel in + xls_data, a dict of dicts + + Example + import config as conf + R=river(conf) + a = R.xls_data['DIP_KTperYR_NOBLS']['2004'] + ''' + self.georef = None + logging.info("Reading from xls file...") + sheet_list=conf.river_data_sheet + river_excel_file = excel_reader.xlsx(conf.file_river) + river_spreadsheet = {} + #read xlsx into a dict + river_spreadsheet["monthly"] = river_excel_file.read_spreadsheet_all("monthly") + for sheet in sheet_list: + river_spreadsheet[sheet] = river_excel_file.read_spreadsheet_all(sheet) + + #extract data + self.lon = river_spreadsheet["monthly"][1:,1].astype(np.float32) + self.lat = river_spreadsheet["monthly"][1:,2].astype(np.float32) + self.forced_coord = river_spreadsheet["monthly"][1:,3:5].astype(np.int32) + self.monthly_mod = river_spreadsheet["monthly"][1:,9:21].astype(np.float32) + + self.nrivers = len(self.lon) + self.xls_data = {} + for sheet in sheet_list: + river_sheet_collected_data = {} + self.available_years = river_spreadsheet[sheet][0][9:] + for iyear, y in enumerate(self.available_years): + river_sheet_collected_data[str(int(y))] = river_spreadsheet[sheet][1:,iyear+9].astype(np.float64) + self.xls_data[sheet] = river_sheet_collected_data.copy() + self.xls_data['monthly'] = river_spreadsheet["monthly"] + logging.info("Done") + + + + def _coast_line_mask(self, mask): + ''' + mask is a 2d array of booleans + Returns: + coast: a 2d array of booleans + ''' + + x,y = mask.shape + coast=np.zeros((x,y),dtype=bool) + for r in range(1,x-1): + for c in range(1,y-1): + expr = (mask[r,c-1] == 0 or mask[r,c+1] == 0 or mask[r-1,c] == 0 or mask[r+1,c] == 0) + expr = expr or (mask[r-1,c-1] == 0 or mask[r+1,c+1] == 0 or mask[r-1,c+1] == 0 or mask[r+1,c-1] == 0) + if mask[r,c] == 1 and expr : + coast[r,c] = True + + return coast + + def gen_map_indexes(self,mask, nspread=1, delta=5): + """ + Arguments: + * mask * a mask object + * nspread * integer, the number of horizontal cells in which a river will be placed + * the nearest "nspread" cells to the corresponding river cell will be choosen + * delta * integer, the dimension in cells (the half side) + of a square box to search for nspread cells. + + Generates the "georef" field + a (nRivers,) numpy array with indLon, indLat, lonmesh, latmesh + or the "georef_spread" field + with the same meaning + + Warning: + * for forced coordinates indLon and indLat are starting from one + * the others are not tested + """ + logging.info("River position calculation: start ") + self.nspread=nspread + there_are_free_points=np.any(self.forced_coord==-1) + if there_are_free_points: + mask1 = mask.mask_at_level(0) + coast = self._coast_line_mask(mask1) + + loncm = mask.xlevels[coast] + latcm = mask.ylevels[coast] + coastline_row_index,coastline_col_index = np.nonzero(coast) + + + ### river contributes + georef = np.zeros((self.nrivers,), dtype=[('indLon',int),('indLat',int),('lonmesh',np.float32),('latmesh',np.float32)]) + georef_spread=np.zeros((self.nrivers*self.nspread),dtype=[('indLon',int),('indLat',int),('lonmesh',np.float32),('latmesh',np.float32)]) + for jr in range (self.nrivers): + if self.forced_coord[jr,0] != -1 and self.forced_coord[jr,1] != -1: + georef['indLon'][jr]=self.forced_coord[jr,0] + georef['indLat'][jr]=self.forced_coord[jr,1] + #if(mask1[georef[jr,1]-1,georef[jr,2]-1] == 0): + # print("RIVER ON THE LAND") + # print(georef[jr,:]) + else: + dist = (loncm-self.lon[jr])**2 + (latcm-self.lat[jr])**2 + ind = np.argmin(dist) + indlon = coastline_col_index[ind]+1 + indlat = coastline_row_index[ind]+1 + if (self.nspread==1): + georef['lonmesh'][jr]=loncm[ind] + georef['latmesh'][jr]=latcm[ind] + + georef['indLon'][jr] = coastline_col_index[ind]+1 + georef['indLat'][jr] = coastline_row_index[ind]+1 + else: + lonmesh=loncm[ind] + latmesh=latcm[ind] + _, jpj, jpi=mask.shape + I,J=np.meshgrid(np.arange(jpi), np.arange(jpj)) + localmask =mask1[indlat-delta:indlat+delta,indlon-delta:indlon+delta] + local_X=mask.xlevels[indlat-delta:indlat+delta,indlon-delta:indlon+delta] + local_Y=mask.ylevels[indlat-delta:indlat+delta,indlon-delta:indlon+delta] + local_I= I[indlat-delta:indlat+delta,indlon-delta:indlon+delta] + local_J= J[indlat-delta:indlat+delta,indlon-delta:indlon+delta] + nPoints=localmask.sum() + NEAREST=np.zeros((nPoints,),dtype=[('I',np.int32), ('J',np.int32),('lon',np.float32), ('lat',np.float32), ('dist',np.float32)]) + NEAREST['I' ]=local_I[localmask] + NEAREST['J' ]=local_J[localmask] + NEAREST['lon' ]=local_X[localmask] + NEAREST['lat' ]=local_Y[localmask] + NEAREST['dist'] = (NEAREST['lon']-lonmesh)**2 + (NEAREST['lat']-latmesh)**2 + NEAREST_ORDERED=np.sort(NEAREST,order='dist') + for k_spread in range(self.nspread): + ind_spread= jr*self.nspread + k_spread + georef_spread['indLon' ][ind_spread]=NEAREST_ORDERED['I' ][k_spread]+1 + georef_spread['indLat' ][ind_spread]=NEAREST_ORDERED['J' ][k_spread]+1 + georef_spread['lonmesh'][ind_spread]=NEAREST_ORDERED['lon'][k_spread] + georef_spread['latmesh'][ind_spread]=NEAREST_ORDERED['lat'][k_spread] + if (self.nspread==1): + self.georef = georef + else: + self.georef_spread = georef_spread + logging.info("River position calculation: done ") + + + def modularize(self,conf): + + ''' + Applyies monthly modulation to yearly data + Generates a new field river_data, + it is a dict of dicts + river_data={'sheet_name':{'2009': (nRivers,12) np.array } ; ... } + Modularize factors are read in "monthly" sheet, + they are pure factors, not percent. + Modularization takes in account weight of eanch month. + ''' + + logging.info("River Modularization") + days_in_month=np.array([31,28,31,30,31,30,31,31,30,31,30,31]) + m=np.zeros((self.nrivers,12),np.float32) + river_data={} + sheet_list = conf.river_data_sheet + for sheet in sheet_list : + years_data={} + for year in self.available_years: + yearstr=str(int(year)) + for r in range(self.nrivers): + ry = self.xls_data[sheet][yearstr][r] + m[r,:] = (self.monthly_mod[r,:]*365/days_in_month)*ry + years_data[yearstr]=m.copy() + river_data[sheet]=years_data.copy() + + self.river_data = river_data + + def get_coordinates(self,mask): + ''' + Arguments: + * mask * a mask object + Returns: + * LAT * numpy ndarray + * LON * numpy ndarray + ''' + LAT = np.zeros((self.nrivers), np.float32) + LON = np.zeros((self.nrivers), np.float32) + for jr in range(self.nrivers): + jj = self.georef[jr]['indLat']-1 + ji = self.georef[jr]['indLon']-1 + LAT[jr] = mask.ylevels[jj,ji] + LON[jr] = mask.xlevels[jj,ji] + + return LAT, LON + + + + def gen_boun_indexes(self,boun_indexes): + ''' + For each river, generates its bounmask index + Arguments + * boun_indexes * is a 3d array of integers, the "index" array of bounmask.nc + Returns + * idxt * integer 1d array (nrivers) with rivers index in bounmask + indexes start from one, is ok for fortran + * positions * integer 2d array (nrivers,3) with river position jk,jj,ji in meshmask + indexes start from zero, is ok for python + + External conditions: + * bounmask.nc['index'] values are supposed start from one + * excel data i,j are indexes starting from one + ''' + + position = np.zeros((self.nrivers,3), dtype = int); + index_riv_a = np.zeros((self.nrivers,) , dtype = int); + for jr in range(self.nrivers): + jj = self.georef[jr]['indLat']-1; + ji = self.georef[jr]['indLon']-1; + index_riv_a[jr] = boun_indexes[0,jj,ji]; + assert index_riv_a[jr] > 0 + position[jr] = [0,jj,ji] + idxt_riv = index_riv_a; #no sort at the moment + return idxt_riv, position + + def total(self,year): + ''' + Returns: + * nitrate * kT/yr + TOTALMENTE INUTILE ADESSO + ''' + + totN = self.river_data["DIN_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; + totP = self.river_data["DIP_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; + totS = self.river_data["DIS_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; + totA = self.river_data["ALK_GmolperYR_NOBLS"][str(year)].sum(axis=1)/12; + totD = self.river_data["DIC_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; + totH = self.river_data["HgII_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; + totM = self.river_data["MMHg_KTperYR_NOBLS" ][str(year)].sum(axis=1)/12; + return totN,totP,totS,totA,totD,totH,totM + + + def spread_array(self,array): + assert len(array)==self.nrivers + long_array=np.zeros((self.nrivers*self.nspread)) + for jr in range(self.nrivers): + for k in range(self.nspread): + ind = jr*self.nspread + k + long_array[ind] = array[jr]/self.nspread + return long_array + + def get_monthly_data(self,yearstr,month): + ''' + Returns monthly data, read from excel and modularized + Arguments: + *yearstr* year as strig + *month* integer from 1 to 12 + + Returns: + N,P,S,A,D : (nRivers, ) numpy arrays + + ''' + # N = self.river_data["DIN_KTperYR_NOBLS" ][yearstr][:,month-1] + # P = self.river_data["DIP_KTperYR_NOBLS" ][yearstr][:,month-1] + # S = self.river_data["DIS_KTperYR_NOBLS" ][yearstr][:,month-1] + # A = self.river_data["ALK_GmolperYR_NOBLS"][yearstr][:,month-1] + # D = self.river_data["DIC_KTperYR_NOBLS" ][yearstr][:,month-1] + # O = self.river_data["O2o_GmolperYR_NOBLS"][yearstr][:,month-1] + H = self.river_data["HgII_KTperYR_NOBLS"][yearstr][:,month-1] + M = self.river_data["MMHg_KTperYR_NOBLS"][yearstr][:,month-1] + # DOC = self.river_data["DOC_KTperYR_NOBLS"][yearstr][:,month-1] + # CDOM = self.river_data["CDOM_KTperYR_NOBLS"][yearstr][:,month-1]*4.0 + + if (self.nspread==1): + # return N,P,S,A,D,O,H,M,DOC,CDOM + return H,M + else: + # N_long = self.spread_array(N) + # P_long = self.spread_array(P) + # S_long = self.spread_array(S) + # A_long = self.spread_array(A) + # D_long = self.spread_array(D) + # O_long = self.spread_array(O) + H_long = self.spread_array(H) + M_long = self.spread_array(M) + # DOC_long = self.spread_array(DOC) + # CDOM_long = self.spread_array(CDOM) + # return N_long, P_long, S_long, A_long, D_long, O_long,H_long,M_long,DOC_long, CDOM_long + return H_long,M_long + + + + + #def conversion(self,N,P,S,A,D,O,H,M,DOC,CDOM): + def conversion(self,H,M): + ''' + Performs conversion of all variables + from KT/y to mmol/s or mg/s + + Arguments: + Nitrate, Phosphate, Silicate, Alcalinity, Dissoved Inorganic Carbon, Dissolved_Oxygen + Dissolved Organic Carbon, CDOM + arrays + + Returns + * N * in mmol/s + * P * in mmol/s + * S * in mmol/s + * A * in mg/s + * D * in mg/s + * O * in mmol/s + * DOC * mg/s + * CDOM * mg/s + * Hg * nmol/s + * MeHg * nmol/s + as (nRivers,) numpy array + + Data are not ready for model, they have to be divided by cell area or cell volume. + ''' + w= 1.0e+12; + w2= 1.0e+6; + t = 1./(365 * 86400) + n = 1./14; + p = 1./31; + s = 1./28; + h = 1./200.59; + # cn = w*t*n + # cp = w*t*p + # cs = w*t*s + # ca = w*t + # cc = w*t + ch = w*w2*t*h #check conversion for mercury + # return N*cn, P*cp, S*cs,A*ca, D*cc, O*cc,H*ch,M*ch, DOC*cc, CDOM*cc + return H*ch,M*ch + + def generate_monthly_files(self,conf,mask):#, idxt_riv, positions): + ''' + Generates TIN files for every year and every month. + + positions is the same of georef, then it starts from zero + ''' + logging.info("Non climatological TIN file generation : start ") + + if (self.nspread==1): + Area=np.zeros((self.nrivers,),float) + for jr in range(self.nrivers): + ji = self.georef['indLon'][jr]-1 + jj = self.georef['indLat'][jr]-1 + Area[jr] = mask.area[jj,ji] + else: + Area = np.zeros((self.nrivers*self.nspread,),float) + for jr in range(self.nrivers*self.nspread): + ji = self.georef_spread['indLon'][jr]-1 + jj = self.georef_spread['indLat'][jr]-1 + Area[jr] = mask.area[jj,ji] + + start_year=conf.simulation_start_time + end___year=conf.simulation_end_time + for year in range(start_year,end___year): + for month in range(1,13): + filename = conf.dir_out + "TIN_%d%02d15-00:00:00.nc" %(year, month) + # N,P,S,A,D,O,H,M,DOC,CDOM = self.get_monthly_data(str(year), month) + # N,P,S,A,D,O,H,M,DOC, CDOM = self.conversion(N, P, S, A, D, O,H,M,DOC, CDOM) + # self.dump_file(filename, N/Area, P/Area, S/Area, A/Area, D/Area, O/Area,H/Area, M/Area, DOC/Area, CDOM/Area, mask) + H, M = self.get_monthly_data(str(year), month) + H,M = self.conversion(H,M) + self.dump_file(filename, H/Area, M/Area) + logging.info("Non climatological TIN file generation : done") + + def generate_climatological_monthly_files(self,conf,mask): #idxt_riv, positions): + ''' + Generates 12 TIN_yyyy*nc files + ''' + logging.info("Climatological TIN file generation : start") + Area=np.zeros((self.nrivers,),np.float32) + + if (self.nspread==1): + Area=np.zeros((self.nrivers,),np.float32) + for jr in range(self.nrivers): + ji = self.georef['indLon'][jr]-1 + jj = self.georef['indLat'][jr]-1 + Area[jr] = mask.area[jj,ji] + else: + Area = np.zeros((self.nrivers*self.nspread,),np.float32) + for jr in range(self.nrivers*self.nspread): + ji = self.georef_spread['indLon'][jr]-1 + jj = self.georef_spread['indLat'][jr]-1 + Area[jr] = mask.area[jj,ji] + + year="yyyy" + for month in range(1,13): + filename = conf.dir_out+"/TIN_yyyy%02d15-00:00:00.nc" %(month) + #N,P,S,A,D,O,H,M,DOC,CDOM = self.get_monthly_data(str(year), month) + # N,P,S,A,D,O,H,M,DOC,CDOM = self.conversion(N, P, S, A, D, O,H,M, DOC,CDOM) + # self.dump_file(filename, N/Area, P/Area, S/Area, A/Area, D/Area, O/Area, H/Area,M/Area,DOC/Area, CDOM/Area, mask) + H, M = self.get_monthly_data(str(year), month) + H,M = self.conversion(H,M) + self.dump_file(filename, H/Area, M/Area,mask) + logging.info("Climatological TIN file generation : done") + + + def dump_file_old(self,filename,N,P,S,A,D,O,H,M, idxt_riv,positions): + ''' + Writes the single TIN file + Variables are dumped as they are, all but positions (incremented by one) + ''' + ncfile = netCDF4.Dataset(filename, 'w') + ncfile.createDimension("riv_idxt",self.nrivers) + ncfile.createDimension("coords",3) + riv_idxt_riv = ncfile.createVariable('riv_idxt', 'i4', ('riv_idxt',)) + riv_pos = ncfile.createVariable('position', 'i4', ('riv_idxt','coords')) + setattr(riv_pos,'order',"k,j,i") + setattr(ncfile,'Units','mmol or mg /(s*m2) ') + # riv_a_n3n = ncfile.createVariable('riv_N3n', 'f4', ('riv_idxt',)) + # riv_a_n1p = ncfile.createVariable('riv_N1p', 'f4', ('riv_idxt',)) + # riv_a_n5s = ncfile.createVariable('riv_N5s', 'f4', ('riv_idxt',)) + # riv_a_o3c = ncfile.createVariable('riv_O3c', 'f4', ('riv_idxt',)) + # riv_a_o3h = ncfile.createVariable('riv_O3h', 'f4', ('riv_idxt',)) + # riv_a_O2o = ncfile.createVariable('riv_O2o', 'f4', ('riv_idxt',)) + riv_a_Hg2 = ncfile.createVariable('riv_Hg2', 'f4', ('riv_idxt',)) + riv_a_MMHg = ncfile.createVariable('riv_MMHg', 'f4', ('riv_idxt',)) + + riv_idxt_riv[:] = idxt_riv[:] + riv_pos[:,:] = positions+1 + # riv_a_n3n[:] = N + # riv_a_n1p[:] = P + # riv_a_n5s[:] = S + # riv_a_o3c[:] = D + # riv_a_o3h[:] = A + # riv_a_O2o[:] = O + riv_a_Hg2[:] = H + riv_a_MMHg[:] = M + ncfile.close() + + #def dump_file(self,filename,N,P,S,A,D,O,H,M,DOC, CDOM, mask): + def dump_file(self,filename,H,M, mask): + ''' + Writes the single TIN file + Variables are dumped as they are, all but positions (incremented by one) + ''' + _,jpj, jpi = mask.shape + ncfile = netCDF4.Dataset(filename, 'w') + ncfile.createDimension('lon',jpi) + ncfile.createDimension('lat',jpj) + # riv_a_n3n = ncfile.createVariable('riv_N3n', 'f4', ('lat','lon')) + # riv_a_n1p = ncfile.createVariable('riv_N1p', 'f4', ('lat','lon')) + # riv_a_n5s = ncfile.createVariable('riv_N5s', 'f4', ('lat','lon')) + # riv_a_o3c = ncfile.createVariable('riv_O3c', 'f4', ('lat','lon')) + # riv_a_o3h = ncfile.createVariable('riv_O3h', 'f4', ('lat','lon')) + # riv_a_O2o = ncfile.createVariable('riv_O2o', 'f4', ('lat','lon')) + riv_a_Hg2 = ncfile.createVariable('riv_Hg2', 'f4', ('lat','lon')) + riv_a_MMHg = ncfile.createVariable('riv_MMHg', 'f4', ('lat','lon')) + # riv_a_R3c = ncfile.createVariable('riv_R3c', 'f4', ('lat','lon')) + # riv_a_R3l = ncfile.createVariable('riv_R3l', 'f4', ('lat','lon')) + + # riv_a_n3n[:] = self.get_map_from_1d_array(N, mask) + # riv_a_n1p[:] = self.get_map_from_1d_array(P, mask) + # riv_a_n5s[:] = self.get_map_from_1d_array(S, mask) + # riv_a_o3c[:] = self.get_map_from_1d_array(D, mask) + # riv_a_o3h[:] = self.get_map_from_1d_array(A, mask) + # riv_a_O2o[:] = self.get_map_from_1d_array(O, mask) + riv_a_Hg2[:] = self.get_map_from_1d_array(H, mask) + riv_a_MMHg[:] = self.get_map_from_1d_array(M, mask) + # riv_a_R3c[:] = self.get_map_from_1d_array(DOC, mask) + # riv_a_R3l[:] = self.get_map_from_1d_array(CDOM,mask) + + # setattr(riv_a_n3n,'missing_value',np.float32(1.e+20)) + # setattr(riv_a_n1p,'missing_value',np.float32(1.e+20)) + # setattr(riv_a_n5s,'missing_value',np.float32(1.e+20)) + # setattr(riv_a_o3c,'missing_value',np.float32(1.e+20)) + # setattr(riv_a_o3h,'missing_value',np.float32(1.e+20)) + # setattr(riv_a_O2o,'missing_value',np.float32(1.e+20)) + setattr(riv_a_Hg2,'missing_value',np.float32(1.e+20)) + setattr(riv_a_MMHg,'missing_value',np.float32(1.e+20)) + # setattr(riv_a_R3c,'missing_value',np.float32(1.e+20)) + # setattr(riv_a_R3l,'missing_value',np.float32(1.e+20)) + ncfile.close() + return + + + def get_map_from_1d_array(self,array,mask): + ''' + Arguments: + * array * an 1D numpuy array + * mask * Mask object + We sum because there are more rivers that can be assigned to the same cell + + Returns: + * OUT * a 2D map + ''' + _, jpj, jpi = mask.shape + OUT = np.ones((jpj, jpi), np.float32) * 1.e+20 + OUT[mask.mask_at_level(0)] = -1. + if (self.nspread==1): + for jr in range(self.nrivers): + ji = self.georef['indLon'][jr] - 1 + jj = self.georef['indLat'][jr] - 1 + if OUT[jj,ji] < 0.: + OUT[jj,ji] = array[jr] # First time + else: + OUT[jj,ji] += array[jr] # To handle cells with multiple river points + return OUT + else: + for jr in range(self.nrivers*self.nspread): + ji = self.georef_spread['indLon'][jr] - 1 + jj = self.georef_spread['indLat'][jr] - 1 + if OUT[jj,ji] < 0.: + OUT[jj,ji] = array[jr] # First time + else: + OUT[jj,ji] += array[jr] # To handle cells with multiple river points + return OUT + + + def load_from_file(self,filename,var): + ''' Useful for check/debug ''' + dset = netCDF4.Dataset(filename, 'r') + M = np.array(dset[var]) + dset.close() + return M +if __name__=="__main__": + from bitsea.commons.mask import Mask + import config as conf + TheMask = Mask(conf.file_mask) + R = river(conf) + R.modularize(conf) + R.gen_map_indexes(TheMask) + diff --git a/BC/check.py b/BC/check.py index af860ad..0c3fa21 100644 --- a/BC/check.py +++ b/BC/check.py @@ -6,7 +6,7 @@ BOUN = bounmask(conf) index_inv = BOUN.load('index_inv') -TheMask = Mask(conf.file_mask) +TheMask = Mask.from_file(conf.file_mask) filename=conf.dir_out + "TIN_20040515-00:00:00.nc" diff --git a/BC/checkRIV.py b/BC/checkRIV.py index ec165b0..9c0a82b 100644 --- a/BC/checkRIV.py +++ b/BC/checkRIV.py @@ -7,13 +7,13 @@ from bitsea.commons import timerequestors from bclib.river import conversion -TheMask=Mask("/gpfs/scratch/userexternal/gbolzon0/OPEN_BOUNDARY/TEST_02/wrkdir/MODEL/meshmask.nc") +TheMask=Mask("/g100_work/OGS23_PRACE_IT/grosati/NECCTON/wrkdir/MODEL/meshmask.nc") Mask_0 = TheMask.cut_at_level(0) area = TheMask.area jpk, jpj, jpi= TheMask.shape -INPUTDIR="/gpfs/scratch/userexternal/gbolzon0/OPEN_BOUNDARY/TEST_02/wrkdir/MODEL/BC" -TI = TimeInterval("2017","2018","%Y") +INPUTDIR="/g100_work/OGS23_PRACE_IT/grosati/NECCTON/wrkdir/MODEL/BC" +TI = TimeInterval("2015","%Y") TL =TimeList.fromfilenames(TI, INPUTDIR, "TIN_2*", prefix="TIN_") nFrames = TL.nTimes @@ -22,7 +22,7 @@ SUBM = np.zeros((jpj, jpi), dtype=bool_dtype) for isub, sub in enumerate(OGS.Pred): - S=SubMask(sub, maskobject=Mask_0) + S=SubMask(sub, Mask_0) SUBM[sub.name]=S.mask SUBM['med']=(SUBM['med'] | SUBM[sub.name]) @@ -30,14 +30,14 @@ -VARLIST=["N3n","N1p","O2o","N5s","O3h","O3c"] +VARLIST=["N3n","N1p","O2o","N5s","O3h","O3c","Hg2","MMHg"] BALANCE_KT_MONTH={} for var in VARLIST: - print var + print (var) MMOL_MONTH = np.zeros((nFrames,),dtype=sub__dtype) KTON_MONTH = np.zeros((nFrames,),dtype=sub__dtype) for iFrame, filename in enumerate(TL.filelist): - req = timerequestors.Monthly_req(2017,iFrame+1) + req = timerequestors.Monthly_req(2015,iFrame+1) B=netcdf4.readfile(filename, "riv_" + var) good = B > -1 for isub, sub in enumerate(OGS.P): @@ -49,8 +49,8 @@ KTON_MONTH[sub.name][iFrame] = river_on_sub *time * conversion(var) BALANCE_KT_MONTH[var]=KTON_MONTH -print BALANCE_KT_MONTH['N3n']['adr1'] -print BALANCE_KT_MONTH['N3n']['adr1'].sum() +print (BALANCE_KT_MONTH['Hg2']['adr1']) +print (BALANCE_KT_MONTH['Hg2']['adr1'].sum()) nvars = len(VARLIST) TABLE= np.zeros((nSub,nvars),np.float32) diff --git a/BC/config.py b/BC/config.py index 7af1bf4..914d436 100644 --- a/BC/config.py +++ b/BC/config.py @@ -1,20 +1,23 @@ # paths -dir_out = "/g100_work/OGS_prod100/OPA/V9C/RUNS_SETUP/PREPROC/BC/out/" -file_mask = "/g100_work/OGS_prod100/OPA/V9C/RUNS_SETUP/PREPROC/MASK/meshmask.nc" +dir_out = "/g100_work/OGS23_PRACE_IT/grosati/NECCTON/PREPROC_ogstm/ogstm_preproc/BC/out/" +file_mask = "/g100_work/OGS23_PRACE_IT/grosati/NECCTON/PREPROC_ogstm/ogstm_preproc/BC/meshmask.nc" file_co2 = "CMIP5_scenarios_RCP_CO2_mixing_ratio.nc" -file_nutrients = "VPnutrients_CO2.nc" -file_bmask = dir_out + "bounmask.nc" -file_river = "Perseus-4.6_38rivers_mesh24.xlsx" +file_nutrients = "VPnutrients_CO2Hg.nc" +file_bmask = "/g100_scratch/userexternal/camadio0/Neccton_hindcast1999_2022/wrkdir/MODEL/bounmask.nc" +file_river = "Perseus-4.6_40rivers_genericmesh_withDissHg.xlsx" -simulation_start_time = 2018 +simulation_start_time = 2015 simulation_end_time = 2020 # Atmosphere settings in Mmol/y # origin of data should be added here -po4_wes = ( 357. + 697)/2 -po4_eas = ( 379. + 957)/2 -n3n_wes = (10042. + 72825)/2 -n3n_eas = ( 6064. + 73621)/2 +po4_wes = ( 357. + 697)*4./2. +po4_eas = ( 379. + 957)*4./2. +n3n_wes = (10042. + 72825)/2. +n3n_eas = ( 6064. + 73621)/2. +HgII_med = 0.19 #HgII deposition Mmol/y +MMHg_med = HgII_med*0.02 #MMHg deposition Mmol/y +Hg0atm_med = 1.6 # ng/g 64. + 73621)/2. # bounmask resto settings @@ -27,11 +30,16 @@ [u'N5s', -8.0,-6.5], [u'O3c', -8.0,-6.5], [u'O3h', -8.0,-6.5], - [u'N6r', -8.0,-6.5]] + [u'N6r', -8.0,-6.5], + [u'Hg0', -8.0,-6.5], + [u'Hg2', -8.0,-6.5], + [u'MHg',-8.0,-6.5], + [u'DHg',-8.0,-6.5]] gib_season = (["0630-00:00:00"]) -RST_FILES ="/g100_work/OGS_prod100/OPA/V9C/RUNS_SETUP/PREPROC/IC/RST_2018/RST*nc" +#RST_FILES ="../PREPROC/IC/RST_2018/RST*nc" +RST_FILES ="/g100_work/OGS23_PRACE_IT/grosati/NECCTON/PREPROC_ogstm/ogstm_preproc/IC/RST_const/RST*nc" river_data_sheet = ['KM3perYR_NOBLS', @@ -41,6 +49,8 @@ 'DIC_KTperYR_NOBLS', 'ALK_GmolperYR_NOBLS', 'O2o_GmolperYR_NOBLS', + 'HgII_KTperYR_NOBLS', + 'MMHg_KTperYR_NOBLS', 'DOC_KTperYR_NOBLS', 'CDOM_KTperYR_NOBLS'] diff --git a/BC/config_rivOnlyHg.py b/BC/config_rivOnlyHg.py new file mode 100644 index 0000000..ea4c845 --- /dev/null +++ b/BC/config_rivOnlyHg.py @@ -0,0 +1,61 @@ +# paths +dir_out = "/g100_work/OGS23_PRACE_IT/grosati/NECCTON/PREPROC_ogstm/ogstm_preproc/BC/out/" +file_mask = "/g100_work/OGS23_PRACE_IT/grosati/NECCTON/PREPROC_ogstm/ogstm_preproc/BC/meshmask.nc" +file_co2 = "CMIP5_scenarios_RCP_CO2_mixing_ratio.nc" +file_nutrients = "VPnutrients_CO2Hg.nc" +file_bmask = "bounmask.nc" +file_river = "Perseus-4.6_40rivers_genericmesh_withDissHg.xlsx" + +simulation_start_time = 2015 +simulation_end_time = 2020 + +# Atmosphere settings in Mmol/y +# origin of data should be added here +po4_wes = ( 357. + 697)*4./2. +po4_eas = ( 379. + 957)*4./2. +n3n_wes = (10042. + 72825)/2. +n3n_eas = ( 6064. + 73621)/2. +HgII_med = 0.19 #HgII deposition Mmol/y +MMHg_med = HgII_med*0.02 #MMHg deposition Mmol/y +Hg0atm_med = 1.6 # ng/g 64. + 73621)/2. + + +# bounmask resto settings +rdpmin = 1./24 +rdpmax = 90. + +variables=[[u'N1p', -8.0, -6.5], + [u'N3n', -8.0,-6.5], + [u'O2o', -8.0,-7.5], + [u'N5s', -8.0,-6.5], + [u'O3c', -8.0,-6.5], + [u'O3h', -8.0,-6.5], + [u'N6r', -8.0,-6.5], + [u'Hg0', -8.0,-6.5], + [u'Hg2', -8.0,-6.5], + [u'MHg',-8.0,-6.5], + [u'DHg',-8.0,-6.5]] + +gib_season = (["0630-00:00:00"]) + +#RST_FILES ="../PREPROC/IC/RST_2018/RST*nc" +RST_FILES ="/g100_work/OGS23_PRACE_IT/grosati/NECCTON/PREPROC_ogstm/ogstm_preproc/IC/RST_const/RST*nc" + + +river_data_sheet = ['KM3perYR_NOBLS', + 'HgII_KTperYR_NOBLS', + 'MMHg_KTperYR_NOBLS'] + +''' +river_data_sheet = ['KM3perYR_NOBLS', + 'DIP_KTperYR_NOBLS', + 'DIN_KTperYR_NOBLS', + 'DIS_KTperYR_NOBLS', + 'DIC_KTperYR_NOBLS', + 'ALK_GmolperYR_NOBLS', + 'O2o_GmolperYR_NOBLS', + 'HgII_KTperYR_NOBLS', + 'MMHg_KTperYR_NOBLS', + 'DOC_KTperYR_NOBLS', + 'CDOM_KTperYR_NOBLS'] + ''' diff --git a/BC/dardanelles/analysis_config.py b/BC/dardanelles/analysis_config.py index 3cdf039..67e388f 100644 --- a/BC/dardanelles/analysis_config.py +++ b/BC/dardanelles/analysis_config.py @@ -14,5 +14,5 @@ INPUTDIR="/gpfs/scratch/userexternal/gbolzon0/TRANSITION_24/AVE/" TI = TimeInterval("2016","2018","%Y") -VARLIST=['N1p', 'N3n', 'N5s','O3c','O3h'] +VARLIST=['N1p', 'N3n', 'N5s','O3c','O3h','HgII','MMHg'] mydtype=[(var,np.float32) for var in VARLIST] diff --git a/BC/dardanelles/dardanelles_generator.py b/BC/dardanelles/dardanelles_generator.py index 6704956..da78549 100644 --- a/BC/dardanelles/dardanelles_generator.py +++ b/BC/dardanelles/dardanelles_generator.py @@ -7,9 +7,11 @@ BOUNDARY_CONCENTRATION= np.zeros((1,),dtype=mydtype) BOUNDARY_CONCENTRATION['N5s'] = 2.0 BOUNDARY_CONCENTRATION['N1p'] = 0.065 # mmol/m3 -BOUNDARY_CONCENTRATION['N3n']= 1.3 # mol/m3 -BOUNDARY_CONCENTRATION['O3c']= 28700 # mg/m3 -BOUNDARY_CONCENTRATION['O3h']= 2800 # mmol/m3 +BOUNDARY_CONCENTRATION['N3n']= 1.3 # mol/m3 +BOUNDARY_CONCENTRATION['O3c']= 28700 # mg/m3 +BOUNDARY_CONCENTRATION['O3h']= 2800 # mmol/m3 +BOUNDARY_CONCENTRATION['HgII']= 2. # nmol/m3 +BOUNDARY_CONCENTRATION['MMHg']= 0.07 # nmol/m3 def dump_files(OpenMask, OUTPUTDIR): @@ -28,4 +30,4 @@ def dump_files(OpenMask, OUTPUTDIR): - \ No newline at end of file + diff --git a/BC/main.py b/BC/main.py index 36e1f1b..63c6280 100755 --- a/BC/main.py +++ b/BC/main.py @@ -9,24 +9,22 @@ import numpy as np from dardanelles import dardanelles_generator -TheMask = Mask(conf.file_mask) +TheMask = Mask.from_file(conf.file_mask) +#CO2 =co2atm(conf) +#CO2.generate(TheMask, experiment="RCP85") + #ATM=atmosphere.atmosphere(TheMask,conf) +#ATM.write_netcdf(TheMask, conf.dir_out, "Area") -CO2 =co2atm(conf) -CO2.generate(TheMask, experiment="RCP85") -ATM=atmosphere.atmosphere(TheMask,conf) -ATM.write_netcdf(TheMask, conf.dir_out, "Area") - - -BOUN = bounmask(conf) -BOUN.generate(TheMask) # can be commented in case of test -BOUN.write_netcdf(TheMask, all_variables=False) +#BOUN = bounmask(conf) +#BOUN.generate(TheMask) # can be commented in case of test +#BOUN.write_netcdf(TheMask, all_variables=False) G = gib(conf,TheMask) G.generate(TheMask, BOUN, all_variables=False) -R = river(conf) # here excel is read -R.modularize(conf) -R.gen_map_indexes(TheMask) +# R = river(conf) # here excel is read +#R.modularize(conf) +#R.gen_map_indexes(TheMask) climatological=True diff --git a/BC/write_VPnutrients.py b/BC/write_VPnutrients.py index 0d986cb..bf2ad12 100644 --- a/BC/write_VPnutrients.py +++ b/BC/write_VPnutrients.py @@ -20,11 +20,11 @@ def argument(): args = argument() -from bitsea.commons.mask import Mask +from bitsea.commons.mesh import Mesh from bitsea.commons import netcdf4 import netCDF4 from bitsea.commons.utils import addsep -TheMask=Mask(args.maskfile, loadtmask=False) +TheMask=Mesh.from_file(args.maskfile) INPUTDIR=addsep(args.inputdir) jpk,_,_ = TheMask.shape