New bc nudging on r3l

BC/bclib/r3l.py

@@ -0,0 +1,32 @@
+import numpy as np
+from bitsea.commons import netcdf4
+from bitsea.commons.mask import Mask
+
+def generate(mask,filename):
+    val = 1./3600
+    jpk,jpj,jpi = mask.shape
+    V = np.ones((jpk,jpj,jpi),np.float32)*val
+    V[~mask.mask] = 1.e+20
+    netcdf4.write_3d_file(V, "reR3l", filename, mask, compression=True)
+
+if __name__=="__main__":
+    import argparse
+    def argument():
+        parser = argparse.ArgumentParser(description = '''
+        Creates R3l restoration file                    
+                                   ''')
+        parser.add_argument(   '--outfile', '-o',
+                                    type = str,
+                                    required = True,
+                                    help="bounmask.nc"
+                                    )
+        parser.add_argument(   '--maskfile', '-m',
+                                    type = str,
+                                    required = True,
+                                    help ='/g100_work/OGS_devC/Benchmark/SETUP/PREPROC/MASK/meshmask.nc'
+                                    )
+        return parser.parse_args()
+    args = argument()
+
+    TheMask = Mask.from_file(args.maskfile)
+    generate(TheMask,args.outfile)

IC/create_2basins_restarts.py

@@ -30,6 +30,7 @@ def argument():
 from bitsea.basins.basin import ComposedBasin
 from IC import RSTwriter
 import numpy as np
+from bitsea.commons import netcdf4
 
 
 TheMask=Mask.from_file(args.maskfile)
@@ -70,4 +71,5 @@ def CDOM_profile(surface_value,bottom_value,depth,depth0):
 #RSTwriter(args.outfile_prefix + ".R2l.nc", "R2l", R2l, TheMask)
 RSTwriter(args.outfile_prefix + ".R3l.nc", "R3l", R3l, TheMask)
 print("done")
+netcdf4.write_3d_file(R3l, "R3l", "R3l_yyyy0630-00:00:00.nc", TheMask, compression=True)
 

IC/create_N-basins_restarts.py

@@ -0,0 +1,107 @@
+from bitsea.commons.mask import Mask
+from bitsea.commons.submask import SubMask
+from bitsea.basins import V2
+from bitsea.basins.basin import ComposedBasin
+from bitsea.commons.interpolators import shift
+from IC import RSTwriter
+import numpy as np
+import pandas as pd
+from bitsea.commons import netcdf4
+import argparse
+
+def argument():
+    parser = argparse.ArgumentParser(description = '''
+    Generates restart and nudging file for R3l
+    ''')
+    parser.add_argument(   '--rst_file', '-r',
+                                type = str,
+                                required = True,
+                                default = "RST.19950101-00:00:00.R3l.nc",
+                                help = "Path of restart file"
+                                )
+    parser.add_argument(   '--ndg_file', '-n',
+                                type = str,
+                                required = True,
+                                default = "R3l.yyyy0630-00:00:00.R3l.nc",
+                                help = 'Path of nudging file'
+                                )
+    parser.add_argument(   '--maskfile','-m',
+                                type = str,
+                                required = True,
+                                help = 'Path of the mask file'
+                                )   
+    parser.add_argument(   '--csvfile','-c',
+                                type = str,
+                                required = True,
+                                help = 'Path of the csv file with surface R3l values'
+                                )
+    return parser.parse_args()
+
+def CDOM_profile(surface_value, bottom_value, depth, depth0):
+    k = 0.05
+    res = surface_value + (bottom_value - surface_value) / (1.0 + np.exp(-k * (depth - depth0)))
+    return res
+
+def smoother3D(mask, RST):
+    #NsmoothingCycles = 20
+    NsmoothingCycles = 2000
+    jpk,jpj,jpi = mask.shape
+    RST[mask ==0 ] = np.nan
+    auxs = np.zeros((5,jpk,jpj,jpi),np.float32)
+    for ss in range(NsmoothingCycles):
+        print(f'{ss}/{NsmoothingCycles}')
+        for j in range(jpk):
+           auxs[0,j,:,:] = RST[j,:,:]
+           auxs[1,j,:,:] = shift(RST[j,:,:],1,'r')
+           auxs[2,j,:,:] = shift(RST[j,:,:],1,'l')
+           auxs[3,j,:,:] = shift(RST[j,:,:],1,'u')
+           auxs[4,j,:,:] = shift(RST[j,:,:],1,'d')
+        RST = np.nanmean(auxs,axis=0)
+    return RST
+
+def smoother2D(mask,RST):
+    NsmoothingCycles = 2000
+    jpj,jpi = mask.shape
+    RST[mask ==0 ] = np.nan
+    auxs = np.zeros((5,jpj,jpi),np.float32)
+    for ss in range(NsmoothingCycles):
+        auxs[0,:,:] = RST
+        auxs[1,:,:] = shift(RST,1,'r')
+        auxs[2,:,:] = shift(RST,1,'l')
+        auxs[3,:,:] = shift(RST,1,'u')
+        auxs[4,:,:] = shift(RST,1,'d')
+        RST = np.nanmean(auxs,axis=0)
+    return RST
+
+args = argument()
+TheMask = Mask.from_file(args.maskfile)
+tmask_3D = TheMask._data_array.astype(float)
+csvfile = args.csvfile
+
+Surf_Values = pd.read_csv(csvfile, index_col=0)
+
+jpk, jpj, jpi = TheMask.shape
+R3l_0 = np.ones((1, jpj, jpi), np.float64) * 1.0
+R3l = np.ones((jpk, jpj, jpi), np.float64) * 1.0
+
+nav_lev=TheMask.zlevels
+
+for sub in V2.P:
+    bname = sub.name
+    if bname=='med':
+        continue
+    print(bname)
+    surf_val = Surf_Values.loc[bname].R3l
+    curr_mask = SubMask(sub, mask=TheMask)
+    R3l_0[0,curr_mask[0,:,:]] = surf_val
+
+print('smoothing...')
+R3l_s = smoother2D(tmask_3D[0,:,:], R3l_0[0,:,:])
+#R3l_s = R3l_0[0,:,:]
+
+R3l_s = CDOM_profile(R3l_s[np.newaxis,:,:], 1.0, nav_lev[:,np.newaxis,np.newaxis], 150.0)
+
+
+RSTwriter(args.rst_file, "R3l", R3l_s, TheMask)
+netcdf4.write_3d_file(R3l_s, "R3l", args.ndg_file, TheMask, compression=True)
+

IC/create_N-basins_restarts.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+. ../profile.inc
+source /g100_work/OGS23_PRACE_IT/ggalli00/VENVZ/sequence_gg.sh
+
+csvfile=$PWD/mean_R3l_from_RRS412.csv
+#my_prex_or_die "python create_R3l_nudging.py -o $csvfile"
+
+
+rst=/g100_work/OGS_test2528/V12C_QUID_neccton/PREPROC/IC/RST/RST.20220101-00:00:00.R3l.nc
+ndg=/g100_work/OGS_test2528/V12C_QUID_neccton/PREPROC/R3l/bc/R3l_yyyy0630-00:00:00.nc
+maskfile=/g100_work/OGS_devC/Benchmark/SETUP/PREPROC/MASK/meshmask.nc
+
+my_prex_or_die "python create_N-basins_restarts.py -r $rst -n $ndg -m $maskfile -c $csvfile"
+

IC/create_R3l_nudging.py

@@ -0,0 +1,177 @@
+import numpy as np
+import xarray as xr
+import pandas as pd
+import matplotlib.pyplot as plt
+from glob import glob
+#bitsea stuff
+from bitsea.basins import V2 as OGS
+from bitsea.commons.grid import RegularGrid
+import argparse
+
+def argument():
+    parser = argparse.ArgumentParser(description="""
+    creates a csv file with mean surface R3l in each subbasin,
+    with gradient estimated from satellite reflectance and
+    absolute value from our tuned R3l values in nwm and lev3.
+    these values are then fed to ANOTHER SCRIPT that extends at depth and
+    create the 3D nudging file
+
+
+    SPiani notes
+
+    for mysub in OGS.P:
+        print(mysub)
+    MM = mysub.is_inside(Lon_rrs[:,np.newaxis], Lat_rrs)
+
+    R = RegularGrid(lon=Lon_rrs, lat=Lat_rrs)
+    Area = R.area
+    """,
+    formatter_class=argparse.RawTextHelpFormatter
+    )
+
+
+    parser.add_argument('--outfile','-o ',
+                        type=str,
+                        default=None,
+                        required=True,
+                        help=''' output csv file''')
+    return parser.parse_args()
+
+
+args = argument()
+
+
+def get_submask(mysub, lat, lon, X):
+    seamask = ~np.isnan(X)
+    tmask = mysub.is_inside(lon[:,np.newaxis], lat).T
+    nanmask = tmask.astype(float)
+    nanmask[~tmask] = np.nan
+    area = RegularGrid(lon=lon, lat=lat).area
+    # mask_land
+    tmask = tmask * seamask
+    nanmask = nanmask * seamask
+    area = area * seamask * nanmask
+    return tmask, nanmask, area
+
+def saturate(Y, low_val=1.0, hig_val=1.0):
+    Y = np.maximum(Y, low_val)
+    Y = np.minimum(Y, hig_val)
+    return Y
+
+PL_dir = '/g100_work/OGS23_PRACE_IT/plazzari/COPERNICUS/ANALYZE_REFLECTANCE/'
+CA_dir = '/g100_scratch/userexternal/camadio0/Neccton_hindcast1999_2022_v18/wrkdir/POSTPROC/output/AVE_FREQ_3/RRS_DAILY/2003/'
+CS_dir = '/g100_work/OGS23_PRACE_IT/csoto/rrs_data/V11C/SAT/MONTLY/montly_20*/'
+
+#infile_nudge = '/g100_work/OGS_devC/camadio/Neccton_hindcast1999_2022_v17/R3l_yyyy0630-00:00:00.nc'
+infile_nudge = '/g100_work/OGS_devC/camadio/Neccton_hindcast1999_2022_v17/R3l_yyyy0630-00:00:00.nc'
+infile_rrs = '*_rrs_mediterranean_cmems_montly.nc'
+
+flist_rrs = sorted(glob(CS_dir+infile_rrs)) #this also gets in 2000 and 2001
+#flist_rrs = [ff for ff in flist_rrs if ]
+
+# loads mean RRS412 field
+var_rrs = 'RRS412'
+D = xr.open_mfdataset(flist_rrs, combine='nested', concat_dim='time')
+Lat_rrs = D['lat'].values[:]
+Lon_rrs = D['lon'].values[:]
+RRS = D[var_rrs].mean(dim='time').values[:]
+D.close()
+
+#RRS_range = np.nanmax(RRS) - np.nanmin(RRS)
+RRS_range = np.nanpercentile(RRS, 99) - np.nanpercentile(RRS, 1)
+
+# loads R3l nudge file, takes the surface value
+var_ndg = 'R3l'
+D = xr.open_dataset(infile_nudge)
+Lat_ndg = D['latitude'].values[:]
+Lon_ndg = D['longitude'].values[:]
+R3l = D[var_ndg].values[0,:,:]
+D.close()
+
+R3l_min = np.nanpercentile(R3l, 10) # =0.43
+R3l_max = np.nanpercentile(R3l, 90) # =0.88
+#R3l_max = 1.0
+R3l_range = R3l_max - R3l_min
+
+#X = ((RRS - np.nanmin(RRS)) * (R3l_range / RRS_range)) + np.nanmin(R3l)
+
+iRRS = 1. / RRS
+
+#hig_val = np.nanmean(iRRS[500:1100, 1100:]) # circa nwm
+#hig_val = np.nanmean(iRRS[400:800,:500])   # circa alb 
+#low_val = np.nanmean(iRRS[:400, 2500:])     # circa lev3/4
+
+tmask_alb, nanmask_alb, area_alb = get_submask(OGS.alb, Lat_rrs, Lon_rrs, RRS)
+tmask_nwm, nanmask_nwm, area_nwm = get_submask(OGS.nwm, Lat_rrs, Lon_rrs, RRS)
+tmask_lev3, nanmask_lev3, area_lev3 = get_submask(OGS.lev3, Lat_rrs, Lon_rrs, RRS)
+
+'''
+either 
+R3l_max = 1.0 & hig_val = hig_val(alb)
+or
+R3l_max = 0.88 % hig_val = hig_val(nwm)
+'''
+
+# compute regional means to decide iRRS range
+#hig_val = np.nansum(nanmask_alb * area_alb * iRRS) / np.nansum(nanmask_alb * area_alb)
+hig_val = np.nansum(nanmask_nwm * area_nwm * iRRS) / np.nansum(nanmask_nwm * area_nwm)
+low_val = np.nansum(nanmask_lev3 * area_lev3 * iRRS) / np.nansum(nanmask_lev3 * area_lev3)
+
+#iRRS = saturate(iRRS, low_val, hig_val)
+
+iRRS_range = hig_val - low_val
+#iRRS_range = np.nanmax(iRRS) - np.nanmin(iRRS)
+
+# RESCALE iRRS so that it falls into the range of R3l
+#X = ((iRRS - np.nanmin(iRRS)) * (R3l_range / iRRS_range)) + np.nanmin(R3l)
+X = ((iRRS - low_val) * (R3l_range / iRRS_range)) + np.nanmin(R3l)
+
+fig0 = plt.figure()
+ax0 = plt.axes()
+#im0 = ax0.pcolormesh(X, vmin=R3l_min, vmax=R3l_max)
+im0 = ax0.pcolormesh(X, vmin=R3l_min, vmax=1.0)
+plt.colorbar(im0)
+fig0.suptitle('R3l new')
+plt.savefig('R3l_from_iRRS.png')
+
+fig1 = plt.figure()
+ax1 = plt.axes()
+im1 = ax1.pcolormesh(R3l, vmin=R3l_min, vmax=R3l_max)
+plt.colorbar(im1)
+fig1.suptitle('R3l old')
+plt.savefig('R3l_2basin.png')
+
+fig2 = plt.figure()
+ax2 = plt.axes()
+cmin = np.nanpercentile(iRRS, 10)
+cmax = np.nanpercentile(iRRS, 90)
+im2 = ax2.pcolormesh(iRRS, vmin=cmin, vmax=cmax)
+plt.colorbar(im2)
+fig2.suptitle('iRRS')
+plt.savefig('iRRS.png')
+
+# basin mean values
+blist = [ii.name for ii in OGS.P.basin_list]
+Vals_Dict = dict.fromkeys(blist)
+for mysub in OGS.P:
+    sn = mysub.name
+    tmask, nanmask, area = get_submask(mysub, Lat_rrs, Lon_rrs, RRS)
+    sub_mean = np.nansum(nanmask * area * X) / np.nansum(nanmask * area)
+    print(f'{sn}: {sub_mean}')
+    Vals_Dict[sn] = sub_mean
+
+#manually put atl in
+Vals_Dict['atl'] = Vals_Dict['alb']
+
+#manually set aeg value because else it is too high!
+Vals_Dict['aeg'] = 0.55
+
+C = pd.DataFrame.from_dict(Vals_Dict, orient='index', columns=['R3l'])
+C.to_csv(args.outfile)
+
+
+#for mysub in OGS.P:
+#    sn = mysub.name
+#    tmask, nanmask, area = get_submask(mysub, Lat_rrs, Lon_rrs, iRRS)
+#    sub_mean = np.nansum(nanmask * area * iRRS) / np.nansum(nanmask * area)
+#    print(f'{sn}: {sub_mean}')