Mercury

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:

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)

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.")

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)

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)