From a09408af0dbc0c12fbb5e4e50755c4b3b1ff7f78 Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Wed, 19 Nov 2025 12:16:50 +0100 Subject: [PATCH 01/15] Importing from https://gitlab.hpc.cineca.it/OGS/preproc/-/tree/free-surface-atm --- FORCINGS/degradation/create_meshmask_nc.py | 311 ++++++++++ .../degradation/create_meshmask_nc_MITgcm.py | 317 ++++++++++ .../degradation/create_meshmask_nc_nemo.py | 321 ++++++++++ .../degradation/create_meshmask_reduced.py | 26 + FORCINGS/degradation/forcing_reducer.py | 293 +++++++++ FORCINGS/degradation/forcingswriter.py | 180 ++++++ FORCINGS/degradation/main.py | 28 + FORCINGS/degradation/main_MITgcm.py | 16 + FORCINGS/degradation/reducer.py | 579 ++++++++++++++++++ FORCINGS/degradation/restart_interpolator.py | 113 ++++ FORCINGS/degradation/river_reindexer.py | 71 +++ IC/Ext_without_dres.py | 42 +- 12 files changed, 2282 insertions(+), 15 deletions(-) create mode 100644 FORCINGS/degradation/create_meshmask_nc.py create mode 100644 FORCINGS/degradation/create_meshmask_nc_MITgcm.py create mode 100644 FORCINGS/degradation/create_meshmask_nc_nemo.py create mode 100644 FORCINGS/degradation/create_meshmask_reduced.py create mode 100644 FORCINGS/degradation/forcing_reducer.py create mode 100644 FORCINGS/degradation/forcingswriter.py create mode 100644 FORCINGS/degradation/main.py create mode 100644 FORCINGS/degradation/main_MITgcm.py create mode 100644 FORCINGS/degradation/reducer.py create mode 100644 FORCINGS/degradation/restart_interpolator.py create mode 100644 FORCINGS/degradation/river_reindexer.py diff --git a/FORCINGS/degradation/create_meshmask_nc.py b/FORCINGS/degradation/create_meshmask_nc.py new file mode 100644 index 0000000..6bccc82 --- /dev/null +++ b/FORCINGS/degradation/create_meshmask_nc.py @@ -0,0 +1,311 @@ +# this script requires netcdf 4 +# to load them it is needed the following virtual environment: +# source /gpfs/work/IscrC_MYMEDBIO/COPERNICUS/sequence.sh +# LOAD PACKAGES + + +import numpy as np + +#import scipy.io.netcdf as NC +import netCDF4 as NC + +# Script to create meshmask file + +class OrigMask(): + + def __init__(self,inputfile): + NCin=NC.Dataset(inputfile,"r") + self.Lon = (NCin.variables['glamt'][0,0,:]).copy() + self.Lat = (NCin.variables['gphit'][0,:,0]).copy() + self.nav_lev = (NCin.variables['nav_lev'][:]).copy().astype(np.float); + self.nc_handler = NCin + + + + def get_wp(self): + ''' + Prints some information about the number of waterpooints in longitude + The aim is to provide some information about where to cut this mesh to get the ogstm meshmask. + ''' + NCin=self.nc_handler + jpi=NCin.dimensions['x'].size + jpj=NCin.dimensions['y'].size + tmask= (NCin.variables['tmask'][0,:,:,:]).copy().astype(np.int64) #shape is (121, 380, 1307) + # Biscay correction + for j in range(jpj): + for i in range(jpi): + if(self.Lon[i]<0.0 and self.Lat[j]>42): + tmask[:,j,i]=0. + waterpoints_longitude = tmask.sum(axis=0).sum(axis=0) + + Gibraltar_lon = -5.75 + imed = self.Lon>Gibraltar_lon + med_waterpoints = waterpoints_longitude[imed].sum() + print "mediterranean waterpoints:", med_waterpoints + print "atlantic waterpoints:" + for lon in np.arange(-9,-8,1./16): + iatl = (self.Lon < Gibraltar_lon) & (self.Lon> lon) + atl_waterpoints = waterpoints_longitude[iatl].sum() + print lon, atl_waterpoints + + def getCutLocation(self,lon): + ind =np.argmin(np.abs(self.Lon-lon)) + return ind + + + +time = 1; +z_a = 1; +y_a = 1; +x_a = 1; + + + + +def create_meshmask_nc(OrigMaskobj,outfile,lon_cut,depth_cut,biscay_land = True,free_surface=False): + + NCin=OrigMaskobj.nc_handler + + jpi=NCin.dimensions['x'].size-lon_cut + jpj=NCin.dimensions['y'].size + jpk=NCin.dimensions['z'].size-depth_cut #116 + +# double coastp(y, x) ; + coastp = np.zeros((jpj,jpi),np.double) +# double fmask(time, z, y, x) ; + fmask = np.ones((time,jpk,jpj,jpi),np.double); + fmask[0,:,:,:]=(NCin.variables['fmask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# fmask[0,:,1, 0] = 0.; +# fmask[0,:,0, 1] = 0.; +# fmask[0,:,-1,0] = 0.; +# fmask[0,:,0,-1] = 0.; + +# double gdept(time, z, y_a, x_a) ; + gdept = np.zeros((time,jpk,y_a,x_a),np.double); + gdept[0,0:jpk,0,0] = (NCin.variables['gdept_1d'][0,:jpk]).copy().astype(np.double); + +# double gdepw(time, z, y_a, x_a) ; + gdepw = np.zeros((time,jpk,y_a,x_a),np.double); + gdepw[0,0:jpk,0,0] = (NCin.variables['gdepw_1d'][0,:jpk]).copy().astype(np.double); + + +# double glamt(time, z_a, y, x) ; + glamt = np.ones((time,z_a,jpj,jpi),np.double); + glamt[0,0,:,:]=(NCin.variables['glamt'][0,:,lon_cut:]).copy().astype(np.double); + +# double glamu(time, z_a, y, x) ; + glamu = np.ones((time,z_a,jpj,jpi),np.double); + glamu[0,0,:,:]=(NCin.variables['glamu'][0,:,lon_cut:]).copy().astype(np.double); + +# double glamv(time, z_a, y, x) ; + glamv = np.ones((time,z_a,jpj,jpi),np.double); + glamv[0,0,:,:]=(NCin.variables['glamv'][0,:,lon_cut:]).copy().astype(np.double); + +# double glamf(time, z_a, y, x) ; + glamf = np.ones((time,z_a,jpj,jpi),np.double); + glamf[0,0,:,:]=(NCin.variables['glamf'][0,:,lon_cut:]).copy().astype(np.double); + +# double gphit(time, z_a, y, x) ; + gphit = np.ones((time,z_a,jpj,jpi),np.double); + gphit[0,0,:,:]=(NCin.variables['gphit'][0,:,lon_cut:]).copy().astype(np.double); + +# double gphiu(time, z_a, y, x) ; + gphiu = np.ones((time,z_a,jpj,jpi),np.double); + gphiu[0,0,:,:]=(NCin.variables['gphiu'][0,:,lon_cut:]).copy().astype(np.double); + +# double gphiv(time, z_a, y, x) ; + gphiv = np.ones((time,z_a,jpj,jpi),np.double); + gphiv[0,0,:,:]=(NCin.variables['gphiv'][0,:,lon_cut:]).copy().astype(np.double); + +# double gphif(time, z_a, y, x) ; + gphif = np.ones((time,z_a,jpj,jpi),np.double); + gphif[0,0,:,:]=(NCin.variables['gphif'][0,:,lon_cut:]).copy().astype(np.double); + +# double ff(time, z_a, y, x) ; + ff = np.ones((time,z_a,jpj,jpi),np.double); + ff[0,0,:,:]=(NCin.variables['ff'][0,:,lon_cut:]).copy().astype(np.double); + +# float nav_lat(y, x) ; + nav_lat = np.ones((jpj,jpi),np.float); + nav_lat[:,:] = (NCin.variables['nav_lat'][:,lon_cut:]).copy().astype(np.float); + +# float nav_lev(z) ; + nav_lev = np.ones((jpk,),np.float); + nav_lev[:] = (NCin.variables['nav_lev'][:jpk]).copy().astype(np.float); + +# float nav_lon(y, x) ; + nav_lon = np.ones((jpj,jpi),np.float); + nav_lon[:,:] = (NCin.variables['nav_lon'][:,lon_cut:]).copy().astype(np.float); + + +# double e1f(time, z_a, y, x) ; + e1f = np.ones((time,z_a,jpj,jpi),np.double); + e1f[0,0,:,:] = (NCin.variables['e1f'][0,:,lon_cut:]).copy().astype(np.double) + +# double e1t(time, z_a, y, x) ; + e1t = np.ones((time,z_a,jpj,jpi),np.double); + e1t[0,0,:,:] = (NCin.variables['e1t'][0,:,lon_cut:]).copy().astype(np.double) + +# double e1u(time, z_a, y, x) ; + e1u = np.ones((time,z_a,jpj,jpi),np.double); + e1u[0,0,:,:] = (NCin.variables['e1u'][0,:,lon_cut:]).copy().astype(np.double) + +# double e1v(time, z_a, y, x) ; + e1v = np.ones((time,z_a,jpj,jpi),np.double); + e1v[0,0,:,:] = (NCin.variables['e1v'][0,:,lon_cut:]).copy().astype(np.double) + + +# double e2f(time, z_a, y, x) ; + e2f = np.ones((time,z_a,jpj,jpi),np.double); + e2f[0,0,:,:] = (NCin.variables['e2f'][0,:,lon_cut:]).copy().astype(np.double) + + +# double e2t(time, z_a, y, x) ; + e2t = np.ones((time,z_a,jpj,jpi),np.double) + e2t[0,0,:,:] = (NCin.variables['e2t'][0,:,lon_cut:]).copy().astype(np.double) + +# double e2u(time, z_a, y, x) ; + e2u = np.ones((time,z_a,jpj,jpi),np.double); + e2u[0,0,:,:] = (NCin.variables['e2u'][0,:,lon_cut:]).copy().astype(np.double) + +# double e2v(time, z_a, y, x) ; + e2v = np.ones((time,z_a,jpj,jpi),np.double); + e2v[0,0,:,:] = (NCin.variables['e2v'][0,:,lon_cut:]).copy().astype(np.double) + + if free_surface: + +# double e3t(time, z, y, x) ; + e3t = np.ones((time,jpk,jpj,jpi),np.double) + e3t[0,:,:,:] = (NCin.variables['e3t_0'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# double e3u(time, z, y, x) ; + e3u = np.ones((time,jpk,jpj,jpi),np.double); + e3u[0,:,:,:] = (NCin.variables['e3u_0'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# double e3v(time, z, y, x) ; + e3v = np.ones((time,jpk,jpj,jpi),np.double); + e3v[0,:,:,:] = (NCin.variables['e3v_0'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# double e3w(time, z, y, x) ; + e3w = np.ones((time,jpk,jpj,jpi),np.double); + e3w[0,:,:,:] = (NCin.variables['e3w_0'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + else: +# double e3t(time, z, y_a, x_a) ; + e3t = np.ones((time,jpk,y_a,x_a),np.double); + e3t[0,:,0,0] = (NCin.variables['e3t_0'][0,:jpk]).copy().astype(np.float); + +# double e3w(time, z, y_a, x_a) ; + e3w = np.zeros((time,jpk,y_a,x_a),np.double); + e3w[0,:,0,0] = (NCin.variables['e3w_0'][0,:jpk]).copy().astype(np.float); + +# double tmask(time, z, y, x) ; + tmask = np.ones((time,jpk,jpj,jpi),np.double); + tmask[0,:,:,:] = (NCin.variables['tmask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + tmask[0,:,0, :] =0.; + tmask[0,:,:, 0] =0.; + tmask[0,:,jpj-1,:] =0.; + tmask[0,:,:,jpi-1] =0.; + tmask[0,jpk-1,:,:] =0.; + + if biscay_land: + Lon=glamt[0,0,0,:]; + Lat=gphit[0,0,:,0]; + for j in range(jpj): + for i in range(jpi): + if(Lon[i]<0.0 and Lat[j]>42 ): tmask[0,:,j,i]=0. + if(Lon[i]<-6. and Lat[j]>37.25): tmask[0,:,j,i]=0. + +# double umask(time, z, y, x) ; + umask = np.ones((time,jpk,jpj,jpi),np.double); + umask[0,:,:,:] = (NCin.variables['umask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + umask[0,:,0, :] =0.; + umask[0,:,:, 0] =0.; + umask[0,:,jpj-1,:] =0.; + umask[0,:,:,jpi-1] =0.; +# umask[0,:,:,-2] =0.; To be checked !! + umask[0,jpk-1,:,:] =0.; + + if biscay_land: + Lon=glamu[0,0,0,:]; + Lat=gphiu[0,0,:,0]; + for j in range(jpj): + for i in range(jpi): + if(Lon[i]<0.0 and Lat[j]>42 ): umask[0,:,j,i]=0. + if(Lon[i]<-6. and Lat[j]>37.25): umask[0,:,j,i]=0. +# double vmask(time, z, y, x) ; + + vmask = np.ones((time,jpk,jpj,jpi),np.double); + vmask[0,:,:,:] = (NCin.variables['vmask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + vmask[0,:,:, 0] =0.; + vmask[0,:,0, :] =0.; + vmask[0,:,jpj-1,:] =0.; + vmask[0,:,:,jpi-1] =0.; +# vmask[0,:,-2,:] =0.; To be checked !! + vmask[0,jpk-1,:,:] =0.; + + if biscay_land: + Lon=glamv[0,0,0,:]; + Lat=gphiv[0,0,:,0]; + for j in range(jpj): + for i in range(jpi): + if(Lon[i]<0.0 and Lat[j]>42 ): vmask[0,:,j,i]=0. + if(Lon[i]<-6. and Lat[j]>37.25): vmask[0,:,j,i]=0. + + ############################################################## + # write meshmask netcdf file ! + ############################################################## + ncOUT=NC.Dataset(outfile,"w"); + + ncOUT.createDimension('x',jpi); + ncOUT.createDimension('y',jpj); + ncOUT.createDimension('z',jpk); + ncOUT.createDimension('time',time) + + ncOUT.createDimension('x_a',x_a); + ncOUT.createDimension('y_a',y_a); + ncOUT.createDimension('z_a',z_a); + + ncvar = ncOUT.createVariable('coastp','d',('y','x')) ; ncvar[:] = coastp; + ncvar = ncOUT.createVariable('e1f' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e1f ; + ncvar = ncOUT.createVariable('e1t' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e1t ; + ncvar = ncOUT.createVariable('e1u' ,'d',('time','z_a', 'y', 'x') ); ncvar[:] = e1u ; + ncvar = ncOUT.createVariable('e1v' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e1v ; + ncvar = ncOUT.createVariable('e2f' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e2f ; + ncvar = ncOUT.createVariable('e2t' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e2t ; + ncvar = ncOUT.createVariable('e2u' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e2u ; + ncvar = ncOUT.createVariable('e2v' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e2v ; + if free_surface: + ncvar = ncOUT.createVariable('e3t_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3t ; + ncvar = ncOUT.createVariable('e3u_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3u ; + ncvar = ncOUT.createVariable('e3v_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3v ; + ncvar = ncOUT.createVariable('e3w_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3w ; + else: + ncvar = ncOUT.createVariable('e3t' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = e3t ; + ncvar = ncOUT.createVariable('e3w' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = e3w ; + ncvar = ncOUT.createVariable('ff' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = ff ; + ncvar = ncOUT.createVariable('fmask' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = fmask ; + ncvar = ncOUT.createVariable('gdept' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = gdept ; + ncvar = ncOUT.createVariable('gdepw' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = gdepw ; + ncvar = ncOUT.createVariable('glamf' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamf ; + ncvar = ncOUT.createVariable('glamt' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamt ; + ncvar = ncOUT.createVariable('glamu' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamu ; + ncvar = ncOUT.createVariable('glamv' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamv ; + ncvar = ncOUT.createVariable('gphif' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphif ; + ncvar = ncOUT.createVariable('gphit' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphit ; + ncvar = ncOUT.createVariable('gphiu' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphiu ; + ncvar = ncOUT.createVariable('gphiv' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphiv ; + ncvar = ncOUT.createVariable('nav_lat','f',('y','x')) ; ncvar[:] = nav_lat; + ncvar = ncOUT.createVariable('nav_lev' ,'f',('z',)) ; ncvar[:] = nav_lev; + ncvar = ncOUT.createVariable('nav_lon','f',('y','x')) ; ncvar[:] = nav_lon; +# float time(time) ; +# short time_steps(time) ; + ncvar = ncOUT.createVariable('tmask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = tmask + ncvar = ncOUT.createVariable('umask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = umask + ncvar = ncOUT.createVariable('vmask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = vmask + ncOUT.close() + + + + + diff --git a/FORCINGS/degradation/create_meshmask_nc_MITgcm.py b/FORCINGS/degradation/create_meshmask_nc_MITgcm.py new file mode 100644 index 0000000..82b129b --- /dev/null +++ b/FORCINGS/degradation/create_meshmask_nc_MITgcm.py @@ -0,0 +1,317 @@ +# this script requires netcdf 4 +# to load them it is needed the following virtual environment: +# source /gpfs/work/IscrC_MYMEDBIO/COPERNICUS/sequence.sh +# LOAD PACKAGES + + +import numpy as np + +#import scipy.io.netcdf as NC +import netCDF4 as NC + +# Script to create meshmask file + +class OrigMask(): + + def __init__(self,inputfile): + NCin=NC.Dataset(inputfile,"r") + self.Lon = (NCin.variables['XC'][0,:]).copy() + self.Lat = (NCin.variables['YC'][:,0]).copy() + self.nav_lev = - (NCin.variables['Z'][:]).copy().astype(np.float); + self.nc_handler = NCin + + def get_wp(self): + ''' + Prints some information about the number of waterpooints in longitude + The aim is to provide some information about where to cut this mesh to get the ogstm meshmask. + ''' + NCin=self.nc_handler + jpi=NCin.dimensions['X'].size + jpj=NCin.dimensions['Y'].size + jpk=NCin.dimensions['Z'].size + tmask =np.zeros((jpk,jpj,jpi),dtype=np.int64) + HFacC = (NCin.variables['HFacC'][:,:,:]).copy() #shape is (121, 380, 1307) + idx_wet = HFacC > 0. + tmask[idx_wet] = 1 + + waterpoints_longitude = tmask.sum(axis=0).sum(axis=0) + + Gibraltar_lon = -5.75 + imed = self.Lon>Gibraltar_lon + med_waterpoints = waterpoints_longitude[imed].sum() + print "mediterranean waterpoints:", med_waterpoints + print "atlantic waterpoints:" + for lon in np.arange(-9,-8,1./16): + iatl = (self.Lon < Gibraltar_lon) & (self.Lon> lon) + atl_waterpoints = waterpoints_longitude[iatl].sum() + print lon, atl_waterpoints + + def getCutLocation(self,lon): + ind =np.argmin(np.abs(self.Lon-lon)) + return ind + + + +time = 1; +z_a = 1; +y_a = 1; +x_a = 1; + + + + +def create_meshmask_nc(OrigMaskobj,outfile,free_surface=False): + + NCin=OrigMaskobj.nc_handler + + jpi=NCin.dimensions['X'].size + jpj=NCin.dimensions['Y'].size + jpk=NCin.dimensions['Z'].size + time = 1 + +# double fmask(time, z, y, x) ; + fmask = np.zeros((time,jpk,jpj,jpi),np.double); + print("fmask not used, assigned to 0.") +# fmask[0,:,:,:]=(NCin.variables['fmask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# fmask[0,:,1, 0] = 0.; +# fmask[0,:,0, 1] = 0.; +# fmask[0,:,-1,0] = 0.; +# fmask[0,:,0,-1] = 0.; + +# double gdept(time, z, y_a, x_a) ; + gdept = np.zeros((time,jpk,y_a,x_a),np.double); + gdept[0,0:jpk,0,0] = - (NCin.variables['Z'][:jpk]).copy().astype(np.double); + +# double gdepw(time, z, y_a, x_a) ; + gdepw = np.zeros((time,jpk,y_a,x_a),np.double); + gdepw[0,0:jpk,0,0] = - (NCin.variables['Zp1'][0:jpk]).copy().astype(np.double); + + +# double glamt(time, z_a, y, x) ; + glamt = np.ones((time,z_a,jpj,jpi),np.double); + glamt[0,0,:,:]=(NCin.variables['XC'][:,:]).copy().astype(np.double); + +# double glamu(time, z_a, y, x) ; + glamu = np.ones((time,z_a,jpj,jpi),np.double); + glamu[0,0,:,:]=(NCin.variables['XG'][0:jpj,1:jpi+1]).copy().astype(np.double); + +# double glamv(time, z_a, y, x) ; + glamv = np.ones((time,z_a,jpj,jpi),np.double); + glamv[0,0,:,:]=(NCin.variables['XC'][:,:]).copy().astype(np.double); + +# double glamf(time, z_a, y, x) ; + glamf = np.ones((time,z_a,jpj,jpi),np.double); + glamf[0,0,:,:]=(NCin.variables['XG'][1:jpj+1,1:jpi+1]).copy().astype(np.double); + + +# double gphit(time, z_a, y, x) ; + gphit = np.ones((time,z_a,jpj,jpi),np.double); + gphit[0,0,:,:]=(NCin.variables['YC'][:,:]).copy().astype(np.double); + +# double gphiu(time, z_a, y, x) ; + gphiu = np.ones((time,z_a,jpj,jpi),np.double); + gphiu[0,0,:,:]=(NCin.variables['YC'][:,:]).copy().astype(np.double); + +# double gphiv(time, z_a, y, x) ; + gphiv = np.ones((time,z_a,jpj,jpi),np.double); + gphiv[0,0,:,:]=(NCin.variables['YG'][1:jpj+1,0:jpi]).copy().astype(np.double); + +# double gphif(time, z_a, y, x) ; + gphif = np.ones((time,z_a,jpj,jpi),np.double); + gphif[0,0,:,:]=(NCin.variables['YG'][1:jpj+1,1:jpi+1]).copy().astype(np.double); + +# double ff(time, z_a, y, x) ; + ff = np.zeros((time,z_a,jpj,jpi),np.double); +# ff[0,0,:,:]=(NCin.variables['ff'][0,:,lon_cut:]).copy().astype(np.double); + + +# float nav_lat(y, x) ; + nav_lat = np.ones((jpj,jpi),np.float); + nav_lat[:,:] = (NCin.variables['YC'][:,:]).copy().astype(np.float); + +# float nav_lev(z) ; + nav_lev = np.ones((jpk,),np.float); + nav_lev[:] = -(NCin.variables['Z'][:]).copy().astype(np.float); + +# float nav_lon(y, x) ; + nav_lon = np.ones((jpj,jpi),np.float); + nav_lon[:,:] = (NCin.variables['XC'][:,:]).copy().astype(np.float); + + +# double e1f(time, z_a, y, x) ; + e1f = np.ones((time,z_a,jpj,jpi),np.double); + e1f[0,0,:,:] = (NCin.variables['dxV'][1:jpj+1,1:jpi+1]).copy().astype(np.double) + +# double e1t(time, z_a, y, x) ; + e1t = np.ones((time,z_a,jpj,jpi),np.double); + e1t[0,0,:,:] = (NCin.variables['dxF'][:,:]).copy().astype(np.double) + +# double e1u(time, z_a, y, x) ; + e1u = np.ones((time,z_a,jpj,jpi),np.double); + e1u[0,0,:,:] = (NCin.variables['dxC'][:,1:jpi+1]).copy().astype(np.double) + +# double e1v(time, z_a, y, x) ; + e1v = np.ones((time,z_a,jpj,jpi),np.double); + e1v[0,0,:,:] = (NCin.variables['dxG'][1:jpj+1,:]).copy().astype(np.double) + + +# double e2f(time, z_a, y, x) ; + e2f = np.ones((time,z_a,jpj,jpi),np.double); + e2f[0,0,:,:] = (NCin.variables['dyU'][1:jpj+1,1:jpi+1]).copy().astype(np.double) + + +# double e2t(time, z_a, y, x) ; + e2t = np.ones((time,z_a,jpj,jpi),np.double) + e2t[0,0,:,:] = (NCin.variables['dyF'][:,:]).copy().astype(np.double) + +# double e2u(time, z_a, y, x) ; + e2u = np.ones((time,z_a,jpj,jpi),np.double); + e2u[0,0,:,:] = (NCin.variables['dyG'][:,1:jpi+1]).copy().astype(np.double) + +# double e2v(time, z_a, y, x) ; + e2v = np.ones((time,z_a,jpj,jpi),np.double); + e2v[0,0,:,:] = (NCin.variables['dyC'][1:jpj+1,:]).copy().astype(np.double) + + e3t_1d=np.zeros((jpk)) + Zp1 = - (NCin.variables['Zp1'][:]).copy().astype(np.double) + e3t_1d=np.diff(Zp1) + + e3w_1d=np.zeros((jpk)) + Z = - (NCin.variables['Z'][:]).copy().astype(np.double) + + e3w_1d[0]=Z[0]*2. + for k in range(1,jpk): + e3w_1d[k] = Z[k] -Z[k-1] + + + if free_surface: + +# double e3t(time, z, y, x) ; + e3t = np.ones((time,jpk,jpj,jpi),np.double) + HFacC = (NCin.variables['HFacC'][:,:,:]).copy().astype(np.double) + for jj in range(jpj): + for ji in range(jpi): + e3t[0,:,jj,ji] = e3t_1d[:] * HFacC[:,jj,ji] + +# double e3u(time, z, y, x) ; + e3u = np.ones((time,jpk,jpj,jpi),np.double); + HFacW = (NCin.variables['HFacW'][:,:,1:jpi+1]).copy().astype(np.double) + for jj in range(jpj): + for ji in range(jpi): + e3u[0,:,jj,ji] = e3t_1d[:] * HFacW[:,jj,ji] + +# double e3v(time, z, y, x) ; + e3v = np.ones((time,jpk,jpj,jpi),np.double); + HFacS = (NCin.variables['HFacS'][:,1:jpj+1,:]).copy().astype(np.double) + for jj in range(jpj): + for ji in range(jpi): + e3v[0,:,jj,ji] = e3t_1d[:] * HFacS[:,jj,ji] + +# double e3w(time, z, y, x) ; + e3w = np.ones((time,jpk,jpj,jpi),np.double); + for jj in range(jpj): + for ji in range(jpi): + e3w[0,:,jj,ji] = e3w_1d[:] * HFacC[:,jj,ji] + else: +# double e3t(time, z, y_a, x_a) ; + e3t = np.ones((time,jpk,y_a,x_a),np.double); + e3t[0,:,0,0] = e3t_1d[:] + +# double e3w(time, z, y_a, x_a) ; + e3w = np.zeros((time,jpk,y_a,x_a),np.double); + e3w[0,:,0,0] = e3w_1d[:] + +# double tmask(time, z, y, x) ; + + tmask =np.zeros((time,jpk,jpj,jpi),dtype=np.double) + HFacC = (NCin.variables['HFacC'][:,:,:]).copy() #shape is (121, 380, 1307) + idx_wet = HFacC > 0. + tmask[0,idx_wet] = 1. + + tmask[0,:,0, :] =0.; + tmask[0,:,:, 0] =0.; + tmask[0,:,jpj-1,:] =0.; + tmask[0,:,:,jpi-1] =0.; + tmask[0,jpk-1,:,:] =0.; + + +# double umask(time, z, y, x) ; + umask = np.zeros((time,jpk,jpj,jpi),np.double); + HFacW = (NCin.variables['HFacW'][:,:,1:jpi+1]).copy() #shape is (121, 380, 1307) + idx_wet = HFacW > 0. + umask[0,idx_wet] = 1. + + umask[0,:,0, :] =0.; + umask[0,:,:, 0] =0.; + umask[0,:,jpj-1,:] =0.; + umask[0,:,:,jpi-1] =0.; +# umask[0,:,:,-2] =0.; To be checked !! + umask[0,jpk-1,:,:] =0.; + +# double vmask(time, z, y, x) ; + + vmask = np.ones((time,jpk,jpj,jpi),np.double); + HFacS = (NCin.variables['HFacS'][:,1:jpj+1,:]).copy() #shape is (121, 380, 1307) + idx_wet = HFacS > 0. + vmask[0,idx_wet] = 1. + vmask[0,:,:, 0] =0.; + vmask[0,:,0, :] =0.; + vmask[0,:,jpj-1,:] =0.; + vmask[0,:,:,jpi-1] =0.; +# vmask[0,:,-2,:] =0.; To be checked !! + vmask[0,jpk-1,:,:] =0.; + + ############################################################## + # write meshmask netcdf file ! + ############################################################## + ncOUT=NC.Dataset(outfile,"w"); + + ncOUT.createDimension('x',jpi); + ncOUT.createDimension('y',jpj); + ncOUT.createDimension('z',jpk); + ncOUT.createDimension('time',time) + + ncOUT.createDimension('x_a',x_a); + ncOUT.createDimension('y_a',y_a); + ncOUT.createDimension('z_a',z_a); + +# ncvar = ncOUT.createVariable('coastp','d',('y','x')) ; ncvar[:] = coastp; + ncvar = ncOUT.createVariable('e1f' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e1f ; + ncvar = ncOUT.createVariable('e1t' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e1t ; + ncvar = ncOUT.createVariable('e1u' ,'d',('time','z_a', 'y', 'x') ); ncvar[:] = e1u ; + ncvar = ncOUT.createVariable('e1v' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e1v ; + ncvar = ncOUT.createVariable('e2f' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e2f ; + ncvar = ncOUT.createVariable('e2t' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e2t ; + ncvar = ncOUT.createVariable('e2u' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e2u ; + ncvar = ncOUT.createVariable('e2v' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e2v ; + if free_surface: + ncvar = ncOUT.createVariable('e3t_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3t ; + ncvar = ncOUT.createVariable('e3u_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3u ; + ncvar = ncOUT.createVariable('e3v_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3v ; + ncvar = ncOUT.createVariable('e3w_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3w ; + else: + ncvar = ncOUT.createVariable('e3t' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = e3t ; + ncvar = ncOUT.createVariable('e3w' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = e3w ; + ncvar = ncOUT.createVariable('ff' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = ff ; + ncvar = ncOUT.createVariable('fmask' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = fmask ; + ncvar = ncOUT.createVariable('gdept' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = gdept ; + ncvar = ncOUT.createVariable('gdepw' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = gdepw ; + ncvar = ncOUT.createVariable('glamf' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamf ; + ncvar = ncOUT.createVariable('glamt' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamt ; + ncvar = ncOUT.createVariable('glamu' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamu ; + ncvar = ncOUT.createVariable('glamv' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamv ; + ncvar = ncOUT.createVariable('gphif' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphif ; + ncvar = ncOUT.createVariable('gphit' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphit ; + ncvar = ncOUT.createVariable('gphiu' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphiu ; + ncvar = ncOUT.createVariable('gphiv' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphiv ; + ncvar = ncOUT.createVariable('nav_lat','f',('y','x')) ; ncvar[:] = nav_lat; + ncvar = ncOUT.createVariable('nav_lev' ,'f',('z',)) ; ncvar[:] = nav_lev; + ncvar = ncOUT.createVariable('nav_lon','f',('y','x')) ; ncvar[:] = nav_lon; +# float time(time) ; +# short time_steps(time) ; + ncvar = ncOUT.createVariable('tmask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = tmask + ncvar = ncOUT.createVariable('umask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = umask + ncvar = ncOUT.createVariable('vmask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = vmask + ncOUT.close() + diff --git a/FORCINGS/degradation/create_meshmask_nc_nemo.py b/FORCINGS/degradation/create_meshmask_nc_nemo.py new file mode 100644 index 0000000..5503487 --- /dev/null +++ b/FORCINGS/degradation/create_meshmask_nc_nemo.py @@ -0,0 +1,321 @@ +# For nemo3.4 +# this script requires netcdf 4 +# to load them it is needed the following virtual environment: +# source /gpfs/work/IscrC_MYMEDBIO/COPERNICUS/sequence.sh +# LOAD PACKAGES + + +import numpy as np + +#import scipy.io.netcdf as NC +import netCDF4 as NC + +# Script to create meshmask file + +class OrigMask(): + + def __init__(self,inputfile): + NCin=NC.Dataset(inputfile,"r") + self.Lon = (NCin.variables['glamt'][0,0,:]).copy() + self.Lat = (NCin.variables['gphit'][0,:,0]).copy() + self.nav_lev = (NCin.variables['nav_lev'][:]).copy().astype(np.float); + self.nc_handler = NCin + + + + def get_wp(self): + ''' + Prints some information about the number of waterpooints in longitude + The aim is to provide some information about where to cut this mesh to get the ogstm meshmask. + ''' + NCin=self.nc_handler + jpi=NCin.dimensions['x'].size + jpj=NCin.dimensions['y'].size + tmask= (NCin.variables['tmask'][0,:,:,:]).copy().astype(np.int64) #shape is (121, 380, 1307) + # Biscay correction + for j in range(jpj): + for i in range(jpi): + if(self.Lon[i]<0.0 and self.Lat[j]>42): + tmask[:,j,i]=0. + waterpoints_longitude = tmask.sum(axis=0).sum(axis=0) + + Gibraltar_lon = -5.75 + imed = self.Lon>Gibraltar_lon + med_waterpoints = waterpoints_longitude[imed].sum() + print "mediterranean waterpoints:", med_waterpoints + print "atlantic waterpoints:" + for lon in np.arange(-9,-8,1./16): + iatl = (self.Lon < Gibraltar_lon) & (self.Lon> lon) + atl_waterpoints = waterpoints_longitude[iatl].sum() + print lon, atl_waterpoints + + def getCutLocation(self,lon): + ind =np.argmin(np.abs(self.Lon-lon)) + return ind + + + +time = 1; +z_a = 1; +y_a = 1; +x_a = 1; + + + + +def create_meshmask_nc(OrigMaskobj,outfile,lon_cut,depth_cut,biscay_land = True,free_surface=False): + + NCin=OrigMaskobj.nc_handler + + jpi=NCin.dimensions['x'].size-lon_cut + jpj=NCin.dimensions['y'].size + jpk=NCin.dimensions['z'].size-depth_cut #116 + +# double coastp(y, x) ; + coastp = np.zeros((jpj,jpi),np.double) +# double fmask(time, z, y, x) ; + fmask = np.ones((time,jpk,jpj,jpi),np.double); + fmask[0,:,:,:]=(NCin.variables['fmask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# fmask[0,:,1, 0] = 0.; +# fmask[0,:,0, 1] = 0.; +# fmask[0,:,-1,0] = 0.; +# fmask[0,:,0,-1] = 0.; + +# double gdept(time, z, y_a, x_a) ; + gdept = np.zeros((time,jpk,y_a,x_a),np.double); + gdept[0,0:jpk,0,0] = (NCin.variables['gdept_0'][0,:jpk]).copy().astype(np.double); + +# double gdepw(time, z, y_a, x_a) ; + gdepw = np.zeros((time,jpk,y_a,x_a),np.double); + gdepw[0,0:jpk,0,0] = (NCin.variables['gdepw_0'][0,:jpk]).copy().astype(np.double); + + +# double glamt(time, z_a, y, x) ; + glamt = np.ones((time,z_a,jpj,jpi),np.double); + glamt[0,0,:,:]=(NCin.variables['glamt'][0,:,lon_cut:]).copy().astype(np.double); + +# double glamu(time, z_a, y, x) ; + glamu = np.ones((time,z_a,jpj,jpi),np.double); + glamu[0,0,:,:]=(NCin.variables['glamu'][0,:,lon_cut:]).copy().astype(np.double); + +# double glamv(time, z_a, y, x) ; + glamv = np.ones((time,z_a,jpj,jpi),np.double); + glamv[0,0,:,:]=(NCin.variables['glamv'][0,:,lon_cut:]).copy().astype(np.double); + +# double glamf(time, z_a, y, x) ; + glamf = np.ones((time,z_a,jpj,jpi),np.double); + glamf[0,0,:,:]=(NCin.variables['glamf'][0,:,lon_cut:]).copy().astype(np.double); + +# double gphit(time, z_a, y, x) ; + gphit = np.ones((time,z_a,jpj,jpi),np.double); + gphit[0,0,:,:]=(NCin.variables['gphit'][0,:,lon_cut:]).copy().astype(np.double); + +# double gphiu(time, z_a, y, x) ; + gphiu = np.ones((time,z_a,jpj,jpi),np.double); + gphiu[0,0,:,:]=(NCin.variables['gphiu'][0,:,lon_cut:]).copy().astype(np.double); + +# double gphiv(time, z_a, y, x) ; + gphiv = np.ones((time,z_a,jpj,jpi),np.double); + gphiv[0,0,:,:]=(NCin.variables['gphiv'][0,:,lon_cut:]).copy().astype(np.double); + +# double gphif(time, z_a, y, x) ; + gphif = np.ones((time,z_a,jpj,jpi),np.double); + gphif[0,0,:,:]=(NCin.variables['gphif'][0,:,lon_cut:]).copy().astype(np.double); + +# double ff(time, z_a, y, x) ; + ff = np.ones((time,z_a,jpj,jpi),np.double); + ff[0,0,:,:]=(NCin.variables['ff'][0,:,lon_cut:]).copy().astype(np.double); + +# float nav_lat(y, x) ; + nav_lat = np.ones((jpj,jpi),np.float); + nav_lat[:,:] = (NCin.variables['nav_lat'][:,lon_cut:]).copy().astype(np.float); + +# float nav_lev(z) ; + nav_lev = np.ones((jpk,),np.float); + nav_lev[:] = (NCin.variables['nav_lev'][:jpk]).copy().astype(np.float); + +# float nav_lon(y, x) ; + nav_lon = np.ones((jpj,jpi),np.float); + nav_lon[:,:] = (NCin.variables['nav_lon'][:,lon_cut:]).copy().astype(np.float); + + +# double e1f(time, z_a, y, x) ; + e1f = np.ones((time,z_a,jpj,jpi),np.double); + e1f[0,0,:,:] = (NCin.variables['e1f'][0,:,lon_cut:]).copy().astype(np.double) + +# double e1t(time, z_a, y, x) ; + e1t = np.ones((time,z_a,jpj,jpi),np.double); + e1t[0,0,:,:] = (NCin.variables['e1t'][0,:,lon_cut:]).copy().astype(np.double) + +# double e1u(time, z_a, y, x) ; + e1u = np.ones((time,z_a,jpj,jpi),np.double); + e1u[0,0,:,:] = (NCin.variables['e1u'][0,:,lon_cut:]).copy().astype(np.double) + +# double e1v(time, z_a, y, x) ; + e1v = np.ones((time,z_a,jpj,jpi),np.double); + e1v[0,0,:,:] = (NCin.variables['e1v'][0,:,lon_cut:]).copy().astype(np.double) + + +# double e2f(time, z_a, y, x) ; + e2f = np.ones((time,z_a,jpj,jpi),np.double); + e2f[0,0,:,:] = (NCin.variables['e2f'][0,:,lon_cut:]).copy().astype(np.double) + + +# double e2t(time, z_a, y, x) ; + e2t = np.ones((time,z_a,jpj,jpi),np.double) + e2t[0,0,:,:] = (NCin.variables['e2t'][0,:,lon_cut:]).copy().astype(np.double) + +# double e2u(time, z_a, y, x) ; + e2u = np.ones((time,z_a,jpj,jpi),np.double); + e2u[0,0,:,:] = (NCin.variables['e2u'][0,:,lon_cut:]).copy().astype(np.double) + +# double e2v(time, z_a, y, x) ; + e2v = np.ones((time,z_a,jpj,jpi),np.double); + e2v[0,0,:,:] = (NCin.variables['e2v'][0,:,lon_cut:]).copy().astype(np.double) + + if free_surface: +# double e3t_0(time, z, y_a, x_a) ; + e3t_0 = np.ones((time,jpk,y_a,x_a),np.double); + e3t_0[0,:,0,0] = (NCin.variables['e3t_0'][0,:jpk]).copy().astype(np.double); + +# double e3w_0(time, z, y_a, x_a) ; + e3w_0 = np.zeros((time,jpk,y_a,x_a),np.double); + e3w_0[0,:,0,0] = (NCin.variables['e3w_0'][0,:jpk]).copy().astype(np.double); + +# double e3t(time, z, y, x) ; + e3t = np.ones((time,jpk,jpj,jpi),np.double) + e3t[0,:,:,:] = (NCin.variables['e3t'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# double e3u(time, z, y, x) ; + e3u = np.ones((time,jpk,jpj,jpi),np.double); + e3u[0,:,:,:] = (NCin.variables['e3u'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# double e3v(time, z, y, x) ; + e3v = np.ones((time,jpk,jpj,jpi),np.double); + e3v[0,:,:,:] = (NCin.variables['e3v'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + +# double e3w(time, z, y, x) ; + e3w = np.ones((time,jpk,jpj,jpi),np.double); + e3w[0,:,:,:] = (NCin.variables['e3w'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + else: +# double e3t(time, z, y_a, x_a) ; + e3t = np.ones((time,jpk,y_a,x_a),np.double); + e3t[0,:,0,0] = (NCin.variables['e3t_0'][0,:jpk]).copy().astype(np.float); + +# double e3w(time, z, y_a, x_a) ; + e3w = np.zeros((time,jpk,y_a,x_a),np.double); + e3w[0,:,0,0] = (NCin.variables['e3w_0'][0,:jpk]).copy().astype(np.float); + +# double tmask(time, z, y, x) ; + tmask = np.ones((time,jpk,jpj,jpi),np.double); + tmask[0,:,:,:] = (NCin.variables['tmask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + tmask[0,:,0, :] =0.; + tmask[0,:,:, 0] =0.; + tmask[0,:,jpj-1,:] =0.; + tmask[0,:,:,jpi-1] =0.; + tmask[0,jpk-1,:,:] =0.; + + if biscay_land: + Lon=glamt[0,0,0,:]; + Lat=gphit[0,0,:,0]; + for j in range(jpj): + for i in range(jpi): + if(Lon[i]<0.0 and Lat[j]>42 ): tmask[0,:,j,i]=0. + if(Lon[i]<-6. and Lat[j]>37.25): tmask[0,:,j,i]=0. + +# double umask(time, z, y, x) ; + umask = np.ones((time,jpk,jpj,jpi),np.double); + umask[0,:,:,:] = (NCin.variables['umask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + #umask[0,:,0, :] =0.; + #umask[0,:,:, 0] =0.; + #umask[0,:,jpj-1,:] =0.; + #umask[0,:,:,jpi-1] =0.; +# umask[0,:,:,-2] =0.; To be checked !! + umask[0,jpk-1,:,:] =0.; + + if biscay_land: + Lon=glamu[0,0,0,:]; + Lat=gphiu[0,0,:,0]; + for j in range(jpj): + for i in range(jpi): + if(Lon[i]<0.0 and Lat[j]>42 ): umask[0,:,j,i]=0. + if(Lon[i]<-6. and Lat[j]>37.25): umask[0,:,j,i]=0. +# double vmask(time, z, y, x) ; + + vmask = np.ones((time,jpk,jpj,jpi),np.double); + vmask[0,:,:,:] = (NCin.variables['vmask'][0,:jpk,:,lon_cut:]).copy().astype(np.double) + #vmask[0,:,:, 0] =0.; + #vmask[0,:,0, :] =0.; + #vmask[0,:,jpj-1,:] =0.; + #vmask[0,:,:,jpi-1] =0.; +# vmask[0,:,-2,:] =0.; To be checked !! + vmask[0,jpk-1,:,:] =0.; + + if biscay_land: + Lon=glamv[0,0,0,:]; + Lat=gphiv[0,0,:,0]; + for j in range(jpj): + for i in range(jpi): + if(Lon[i]<0.0 and Lat[j]>42 ): vmask[0,:,j,i]=0. + if(Lon[i]<-6. and Lat[j]>37.25): vmask[0,:,j,i]=0. + + ############################################################## + # write meshmask netcdf file ! + ############################################################## + ncOUT=NC.Dataset(outfile,"w"); + + ncOUT.createDimension('x',jpi); + ncOUT.createDimension('y',jpj); + ncOUT.createDimension('z',jpk); + ncOUT.createDimension('time',time) + + ncOUT.createDimension('x_a',x_a); + ncOUT.createDimension('y_a',y_a); + ncOUT.createDimension('z_a',z_a); + + ncvar = ncOUT.createVariable('coastp','d',('y','x')) ; ncvar[:] = coastp; + ncvar = ncOUT.createVariable('e1f' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e1f ; + ncvar = ncOUT.createVariable('e1t' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e1t ; + ncvar = ncOUT.createVariable('e1u' ,'d',('time','z_a', 'y', 'x') ); ncvar[:] = e1u ; + ncvar = ncOUT.createVariable('e1v' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e1v ; + ncvar = ncOUT.createVariable('e2f' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e2f ; + ncvar = ncOUT.createVariable('e2t' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[:] = e2t ; + ncvar = ncOUT.createVariable('e2u' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e2u ; + ncvar = ncOUT.createVariable('e2v' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = e2v ; + if free_surface: +# ncvar = ncOUT.createVariable('e3t_0' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = e3t_0 ; +# ncvar = ncOUT.createVariable('e3w_0' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = e3w_0 ; + ncvar = ncOUT.createVariable('e3t_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3t ; + ncvar = ncOUT.createVariable('e3u_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3u ; + ncvar = ncOUT.createVariable('e3v_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3v ; + ncvar = ncOUT.createVariable('e3w_0' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = e3w ; + else: + ncvar = ncOUT.createVariable('e3t' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = e3t ; + ncvar = ncOUT.createVariable('e3w' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = e3w ; + ncvar = ncOUT.createVariable('ff' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = ff ; + ncvar = ncOUT.createVariable('fmask' ,'d',('time','z', 'y', 'x')) ; ncvar[:] = fmask ; + ncvar = ncOUT.createVariable('gdept' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = gdept ; + ncvar = ncOUT.createVariable('gdepw' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[:] = gdepw ; + ncvar = ncOUT.createVariable('glamf' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamf ; + ncvar = ncOUT.createVariable('glamt' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamt ; + ncvar = ncOUT.createVariable('glamu' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamu ; + ncvar = ncOUT.createVariable('glamv' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = glamv ; + ncvar = ncOUT.createVariable('gphif' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphif ; + ncvar = ncOUT.createVariable('gphit' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphit ; + ncvar = ncOUT.createVariable('gphiu' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphiu ; + ncvar = ncOUT.createVariable('gphiv' ,'d',('time','z_a', 'y', 'x')) ; ncvar[:] = gphiv ; + ncvar = ncOUT.createVariable('nav_lat','f',('y','x')) ; ncvar[:] = nav_lat; + ncvar = ncOUT.createVariable('nav_lev' ,'f',('z',)) ; ncvar[:] = nav_lev; + ncvar = ncOUT.createVariable('nav_lon','f',('y','x')) ; ncvar[:] = nav_lon; +# float time(time) ; +# short time_steps(time) ; + ncvar = ncOUT.createVariable('tmask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = tmask + ncvar = ncOUT.createVariable('umask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = umask + ncvar = ncOUT.createVariable('vmask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = vmask + ncOUT.close() + + + + + diff --git a/FORCINGS/degradation/create_meshmask_reduced.py b/FORCINGS/degradation/create_meshmask_reduced.py new file mode 100644 index 0000000..23a3425 --- /dev/null +++ b/FORCINGS/degradation/create_meshmask_reduced.py @@ -0,0 +1,26 @@ +import reducer +import numpy as np +isFreeSurface=True + +inputmesh="meshmask.nc" + +R =reducer.reducer(isFreeSurface) +R.read_fine(inputmesh) +passo=4 +Mred,I,J = R.reduce(passo) +Mexp,I,J = R.expand(Mred,I,J) + + +outputmesh="meshmask_472.nc" +R.dumpfile(Mexp, outputmesh) + +np.savetxt('South_West_I_indexes.txt',I,'%d') # of mesh expanded +np.savetxt('South_West_J_indexes.txt',J,'%d') + +#Mexp= reducer.mesh("meshmask_472.nc") +M_med, cut = R.cut_med(Mexp, lon_cut=-8.875, depth_cut=2,biscay_land=True) +R.dumpfile(M_med, "meshmask_470.nc") + +# For BC - input for river_reindexer +np.savetxt('ogstm_South_West_I_indexes.txt',I[cut:],'%d') # of M_med mesh +np.savetxt('ogstm_South_West_J_indexes.txt',J,'%d') diff --git a/FORCINGS/degradation/forcing_reducer.py b/FORCINGS/degradation/forcing_reducer.py new file mode 100644 index 0000000..e362872 --- /dev/null +++ b/FORCINGS/degradation/forcing_reducer.py @@ -0,0 +1,293 @@ +import argparse + +def argument(): + parser = argparse.ArgumentParser(description = ''' + Degradates forcings INGV files in another resolution defined for ogstm model. + New resolution can be 2,3,4,... coarser than the original one. + In order to apply the reduction, 2 files are needed: + ogstm_South_West_I_indexes.txt + ogstm_South_West_J_indexes.txt + ''', formatter_class=argparse.RawTextHelpFormatter) + + + parser.add_argument( '--inputdir','-i', + type = str, + required = True, + help = 'Path of INGV forcings directory') + parser.add_argument( '--outputdir','-o', + type = str, + required = True, + help = 'Path of OGSTM forcings directory') + parser.add_argument( '--ingvmask','-M', + type = str, + required = True, + help = 'Path of the INGV mask') + + parser.add_argument( '--maskfile', '-m', + type = str, + default = None, + required = True, + help = ''' Path of maskfile''') + parser.add_argument( '--starttime','-st', + type = str, + required = False, + default = '19500101', + help = 'start date in yyyymmdd format') + parser.add_argument( '--endtime','-et', + type = str, + required = False, + default = '21000101', + help = 'end date in yyyymmdd format') + + + + return parser.parse_args() + +args = argument() + +from reducer import mesh +import numpy as np +from commons.mask import Mask +from commons.dataextractor import DataExtractor +import forcingswriter +from commons.Timelist import TimeInterval, TimeList +from commons.utils import addsep +import os + +try: + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nranks =comm.size + isParallel = True +except: + rank = 0 + nranks = 1 + isParallel = False + + + +def weighted_mean(values,weights): + return (values*weights).sum()/weights.sum() + +def interp2d(M2dfine, Maskfine,MaskCoarse,I,J): + ''' + Arguments: + * M2dfine * 3d numpy array + * Maskfine * Mask object + * MaskCoarse * Mask object + * I * numpy array 1d of integers + * J * idem + ''' + step = I[3]-I[2] + _, jpj, jpi = MaskCoarse.shape + + M2dCoarse=np.ones((jpj,jpi),np.float32)*1.e+20 + for ji in range(jpi): + for jj in range(jpj): + if MaskCoarse.mask[0,jj,ji]: + i=I[ji] + j=J[jj] + mask = Maskfine.mask[0,j:j+step,i:i+step] + values= M2dfine[j:j+step,i:i+step] + area = Maskfine.area[j:j+step,i:i+step] + M2dCoarse[jj,ji] = weighted_mean(values[mask], area[mask]) + return M2dCoarse + + +def interp3d(M3dfine, Maskfine,MaskCoarse,I,J): + ''' + Degradates a 3d field defined on cell center, using area weighted mean. + + Arguments: + * M3dfine * 3d numpy array + * Maskfine * Mask object + * MaskCoarse * Mask object + * I * numpy array 1d of integers, coarse_lon[k] = fine_lon[I[k]] + * J * idem + + Returns: + * M3dCoarse * 3d numpy array + ''' + step = I[3]-I[2] + jpk, jpj, jpi = MaskCoarse.shape + + M3dCoarse=np.zeros((jpk,jpj,jpi),np.float32) + for jk in range(jpk): + for ji in range(jpi): + for jj in range(jpj): + if MaskCoarse.mask[jk,jj,ji]: + i=I[ji] + j=J[jj] + mask = Maskfine.mask[jk,j:j+step,i:i+step] + values= M3dfine[jk,j:j+step,i:i+step] + area = Maskfine.area[ j:j+step,i:i+step] + M3dCoarse[jk,jj,ji] = weighted_mean(values[mask], area[mask]) + return M3dCoarse + +def interp3d_uv(Ufine,Vfine, Maskfine,MaskCoarse,I,J): + ''' + Generates degradated u,v fields by conserving divercence. + + Arguments: + * Ufine * 3d numpy array + * Vfine * 3d numpy array + * Maskfine * mesh object + * MaskCoarse * mesh object + * I * numpy array 1d of integers, coarse_lon[k] = fine_lon[I[k]] + * J * idem + + Returns: + * Ucoarse * 3d numpy array on coarse mask + * Vcoarse * + ''' + + Ufine[Maskfine.umask==0] = 0.0 + Vfine[Maskfine.vmask==0] = 0.0 + step = I[3]-I[2] + jpi = MaskCoarse.jpi + jpj = MaskCoarse.jpj + jpk = MaskCoarse.jpk + UCoarse=np.zeros((jpk,jpj,jpi),np.float32) + VCoarse=np.zeros((jpk,jpj,jpi),np.float32) + + for jk in range(jpk): + for ji in range(jpi): + for jj in range(jpj): + if (MaskCoarse.umask[jk,jj,ji]==1): + i=I[ji] + j=J[jj] + u = Ufine[jk,j:j+step,i+step-1] + e2u= Maskfine.e2u[ j:j+step,i+step-1] + e3u= Maskfine.e3u_0[jk,j:j+step,i+step-1] + flux = (u*e2u*e3u).sum() + norm = MaskCoarse.e2u[jj,ji]*MaskCoarse.e3u_0[jk,jj,ji] + UCoarse[jk,jj,ji] = flux/norm + + for jk in range(jpk): + for ji in range(jpi): + for jj in range(jpj): + if (MaskCoarse.vmask[jk,jj,ji]==1): + i=I[ji] + j=J[jj] + v = Vfine[jk,j+step-1,i:i+step] + e1v= Maskfine.e1v[ j+step-1,i:i+step] + e3v= Maskfine.e3v_0[jk,j+step-1,i:i+step] + flux = (v*e1v*e3v).sum() + norm = MaskCoarse.e1v[jj,ji]*MaskCoarse.e3v_0[jk,jj,ji] + VCoarse[jk,jj,ji] = flux/norm + return UCoarse,VCoarse + + +def interp2d_tau(tauxFine,tauyFine, Maskfine,MaskCoarse,I,J): + ''' + Generates degradated taux,tauy fields by conserving stress. + + Arguments: + * tauxFine * 2d numpy array + * tauyFine * idem + * Maskfine * mesh object + * MaskCoarse * mesh object + * I * numpy array 1d of integers, coarse_lon[k] = fine_lon[I[k]] + * J * idem + + Returns: + * tauxCoarse * 2d numpy array on coarse mask + * tauxCoarse * idem + ''' + + tauxFine[Maskfine.umask[0,:,:]==0] = 0.0 + tauyFine[Maskfine.vmask[0,:,:]==0] = 0.0 + step = I[3]-I[2] + jpi = MaskCoarse.jpi + jpj = MaskCoarse.jpj + tauxCoarse=np.zeros((jpj,jpi),np.float32) + tauyCoarse=np.zeros((jpj,jpi),np.float32) + + for ji in range(jpi): + for jj in range(jpj): + if (MaskCoarse.umask[0,jj,ji]==1): + i=I[ji] + j=J[jj] + u = tauxFine[j:j+step,i+step-1] + e2u= Maskfine.e2u[j:j+step,i+step-1] + force = (u*e2u).sum() + tauxCoarse[jj,ji] = force/MaskCoarse.e2u[jj,ji] + + for ji in range(jpi): + for jj in range(jpj): + if (MaskCoarse.vmask[0,jj,ji]==1): + i=I[ji] + j=J[jj] + v = tauyFine[j+step-1,i:i+step] + e1v= Maskfine.e1v[j+step-1,i:i+step] + force = (v*e1v).sum() + tauyCoarse[jj,ji] = force/MaskCoarse.e1v[jj,ji] + return tauxCoarse,tauyCoarse + + +INPUTDIR = addsep(args.inputdir) +OUTDIR = addsep(args.outputdir) + +Mfine =mesh(filename=args.ingvmask,isFreeSurface=True) +Mcoarse=mesh(filename=args.maskfile,isFreeSurface=True) +MaskINGV = Mask(args.ingvmask) +MaskCoarse = Mask(args.maskfile) + + + +I = np.loadtxt('South_West_I_indexes.txt',np.int32) +J = np.loadtxt('South_West_J_indexes.txt',np.int32) + + +TI = TimeInterval(args.starttime,args.endtime,"%Y%m%d") +TL = TimeList.fromfilenames(TI, INPUTDIR, "*U.nc", prefix="", dateformat="%Y%m%d") + + + +for filename in TL.filelist[rank::nranks]: + basename = os.path.basename(filename).rsplit("U.nc")[0] + + infileU = INPUTDIR + basename + "U.nc" + infileV = INPUTDIR + basename + "V.nc" + infileW = INPUTDIR + basename + "W.nc" + infileT = INPUTDIR + basename + "T.nc" + # same names + outfileU = OUTDIR + basename + "U.nc" + outfileV = OUTDIR + basename + "V.nc" + outfileW = OUTDIR + basename + "W.nc" + outfileT = OUTDIR + basename + "T.nc" + + print outfileU + + U = DataExtractor(MaskINGV,infileU,'vozocrtx').values + V = DataExtractor(MaskINGV,infileV,'vomecrty').values + TAUX = DataExtractor(MaskINGV,infileU,'sozotaux', dimvar=2).values + TAUY = DataExtractor(MaskINGV,infileV,'sometauy', dimvar=2).values + u,v = interp3d_uv(U, V, Mfine, Mcoarse, I, J) + taux, tauy = interp2d_tau(TAUX, TAUY, Mfine, Mcoarse, I, J) + forcingswriter.writefileU(outfileU, MaskCoarse, u, taux) + forcingswriter.writefileV(outfileV, MaskCoarse, v, tauy) + + + + T = DataExtractor(MaskINGV,infileT,'votemper').values + S = DataExtractor(MaskINGV,infileT,'vosaline').values + E = DataExtractor(MaskINGV,infileT,'sossheig', dimvar=2).values + SW= DataExtractor(MaskINGV,infileT,'soshfldo', dimvar=2).values + + t = interp3d(T, MaskINGV, MaskCoarse, I, J) + s = interp3d(S, MaskINGV, MaskCoarse, I, J) + e = interp2d(E, MaskINGV, MaskCoarse, I, J) + sw= interp2d(SW,MaskINGV, MaskCoarse, I, J) + forcingswriter.writefileT(outfileT, MaskCoarse, t, s, e,sw) + + + #W = DataExtractor(MaskINGV,infileW,'vovecrtz').values + K = DataExtractor(MaskINGV,infileW,'votkeavt').values + #w = interp3d(W, MaskINGV, MaskCoarse, I, J) + k = interp3d(K, MaskINGV, MaskCoarse, I, J) + + forcingswriter.writefileW(outfileW, MaskCoarse, None, k) + + diff --git a/FORCINGS/degradation/forcingswriter.py b/FORCINGS/degradation/forcingswriter.py new file mode 100644 index 0000000..4ec4b53 --- /dev/null +++ b/FORCINGS/degradation/forcingswriter.py @@ -0,0 +1,180 @@ +import netCDF4 as NC +from commons.mask import Mask +#Maskfile_OGS="meshmask_6125.nc" +#TheMask = Mask(Maskfile_OGS) + +def writefileU(filename, mask, U, taux): + + jpk,jpj,jpi = mask.shape + + + ncOUT = NC.Dataset(filename,"w",format="NETCDF4") + ncOUT.createDimension('x',jpi) + ncOUT.createDimension('y',jpj) + ncOUT.createDimension('depthu',jpk) + ncOUT.createDimension('time_counter',1) + + ncvar = ncOUT.createVariable('nav_lon','f',('y','x')) + ncvar[:]=mask.xlevels + ncvar = ncOUT.createVariable('nav_lat','f',('y','x')) + ncvar[:]=mask.ylevels + ncvar = ncOUT.createVariable('depthu','f',('depthu',)) + ncvar[:]=mask.zlevels + + ncvar=ncOUT.createVariable('vozocrtx', 'f' ,('time_counter', 'depthu', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=U + setattr(ncvar,'standard_name',"sea_water_x_velocity") + setattr(ncvar,'long_name',"sea_water_x_velocity") + setattr(ncvar,'units',"m/s" ) + setattr(ncvar,'online_operation',"instant") + setattr(ncvar,'interval_operation',"1 d") + setattr(ncvar,'interval_write',"1 d") + setattr(ncvar,'cell_methods',"time: point") + setattr(ncvar,'coordinates',"time_instant depthu nav_lon nav_lat") + + ncvar = ncOUT.createVariable('sozotaux', 'f', ('time_counter', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=taux + setattr(ncvar,'standard_name',"surface_downward_x_stress") + setattr(ncvar,'long_name',"surface_downward_x_stress") + setattr(ncvar,'units',"N/m2" ) + ncOUT.close() + +def writefileV(filename, mask, V, tauy): + + jpk,jpj,jpi = mask.shape + + + ncOUT = NC.Dataset(filename,"w",format="NETCDF4") + ncOUT.createDimension('x',jpi) + ncOUT.createDimension('y',jpj) + ncOUT.createDimension('depthv',jpk) + ncOUT.createDimension('time_counter',1) + + ncvar = ncOUT.createVariable('nav_lon','f',('y','x')) + ncvar[:]=mask.xlevels + ncvar = ncOUT.createVariable('nav_lat','f',('y','x')) + ncvar[:]=mask.ylevels + ncvar = ncOUT.createVariable('depthv','f',('depthv',)) + ncvar[:]=mask.zlevels + + ncvar=ncOUT.createVariable('vomecrty', 'f' ,('time_counter', 'depthv', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=V + setattr(ncvar,'standard_name',"sea_water_y_velocity") + setattr(ncvar,'long_name',"sea_water_y_velocity") + setattr(ncvar,'units',"m/s" ) + setattr(ncvar,'online_operation',"instant") + setattr(ncvar,'interval_operation',"1 d") + setattr(ncvar,'interval_write',"1 d") + setattr(ncvar,'cell_methods',"time: point") + setattr(ncvar,'coordinates',"time_instant depthv nav_lon nav_lat") + + ncvar = ncOUT.createVariable('sometauy', 'f', ('time_counter', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=tauy + setattr(ncvar,'standard_name',"surface_downward_y_stress") + setattr(ncvar,'long_name',"surface_downward_y_stress") + setattr(ncvar,'units',"N/m2" ) + ncOUT.close() + +def writefileW(filename, mask, W, K): + + jpk,jpj,jpi = mask.shape + + + ncOUT = NC.Dataset(filename,"w",format="NETCDF4") + ncOUT.createDimension('x',jpi) + ncOUT.createDimension('y',jpj) + ncOUT.createDimension('depthw',jpk) + ncOUT.createDimension('time_counter',1) + + ncvar = ncOUT.createVariable('nav_lon','f',('y','x')) + ncvar[:]=mask.xlevels + ncvar = ncOUT.createVariable('nav_lat','f',('y','x')) + ncvar[:]=mask.ylevels + ncvar = ncOUT.createVariable('depthw','f',('depthw',)) + ncvar[:]=mask.zlevels + if W is not None: + ncvar=ncOUT.createVariable('vovecrtz', 'f' ,('time_counter', 'depthw', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=W + setattr(ncvar,'standard_name',"sea_water_z_velocity") + setattr(ncvar,'long_name',"sea_water_z_velocity") + setattr(ncvar,'units',"m/s" ) + setattr(ncvar,'online_operation',"instant") + setattr(ncvar,'interval_operation',"1 d") + setattr(ncvar,'interval_write',"1 d") + setattr(ncvar,'cell_methods',"time: point") + setattr(ncvar,'coordinates', "time_instant depthu nav_lon nav_lat") + + + ncvar=ncOUT.createVariable('votkeavt', 'f' ,('time_counter', 'depthw', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=K + + setattr(ncvar,'standard_name',"sea_water_z_velocity") + setattr(ncvar,'long_name',"sea_water_z_velocity") + setattr(ncvar,'units',"m/s" ) + setattr(ncvar,'online_operation',"instant") + setattr(ncvar,'interval_operation',"1 d") + setattr(ncvar,'interval_write',"1 d") + setattr(ncvar,'cell_methods',"time: point") + + setattr(ncvar,'coordinates',"time_instant depthu nav_lon nav_lat") + + + ncOUT.close() + +def writefileT(filename, mask, T, S, ETA, SW_Rad): + + jpk,jpj,jpi = mask.shape + + + ncOUT = NC.Dataset(filename,"w",format="NETCDF4") + ncOUT.createDimension('x',jpi) + ncOUT.createDimension('y',jpj) + ncOUT.createDimension('deptht',jpk) + ncOUT.createDimension('time_counter',1) + + ncvar = ncOUT.createVariable('nav_lon','f',('y','x')) + ncvar[:]=mask.xlevels + ncvar = ncOUT.createVariable('nav_lat','f',('y','x')) + ncvar[:]=mask.ylevels + ncvar = ncOUT.createVariable('deptht','f',('deptht',)) + ncvar[:]=mask.zlevels + + ncvar=ncOUT.createVariable('vosaline', 'f' ,('time_counter', 'deptht', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=S + setattr(ncvar,'standard_name',"sea_water_z_velocity") + setattr(ncvar,'long_name',"sea_water_z_velocity") + setattr(ncvar,'units',"m/s" ) + setattr(ncvar,'online_operation',"instant") + setattr(ncvar,'interval_operation',"1 d") + setattr(ncvar,'interval_write',"1 d") + setattr(ncvar,'cell_methods',"time: point") + setattr(ncvar,'coordinates', "time_instant depthu nav_lon nav_lat") + + + ncvar=ncOUT.createVariable('votemper', 'f' ,('time_counter', 'deptht', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=T + + setattr(ncvar,'standard_name',"sea_water_z_velocity") + setattr(ncvar,'long_name',"sea_water_z_velocity") + setattr(ncvar,'units',"m/s" ) + setattr(ncvar,'online_operation',"instant") + setattr(ncvar,'interval_operation',"1 d") + setattr(ncvar,'interval_write',"1 d") + setattr(ncvar,'cell_methods',"time: point") + setattr(ncvar,'coordinates',"time_instant depthu nav_lon nav_lat") + + ncvar = ncOUT.createVariable('sossheig', 'f', ('time_counter', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=ETA + setattr(ncvar,'standard_name',"sea_surface_height_above_geoid") + setattr(ncvar,'long_name',"sea_surface_height_above_geoid") + setattr(ncvar,'units',"m" ) + + ncvar = ncOUT.createVariable('soshfldo', 'f', ('time_counter', 'y', 'x'),fill_value=1.0e+20) + ncvar[:]=SW_Rad + setattr(ncvar,'standard_name',"Shortwave Radiation") + setattr(ncvar,'long_name',"Shortwave Radiation") + setattr(ncvar,'units',"W/m2" ) + + ncOUT.close() + + diff --git a/FORCINGS/degradation/main.py b/FORCINGS/degradation/main.py new file mode 100644 index 0000000..20d2cb1 --- /dev/null +++ b/FORCINGS/degradation/main.py @@ -0,0 +1,28 @@ +# SECTION FOR REDUCTION OF 1/16 meshmask +# nemo3.4 class reader +import create_meshmask_nc_nemo as c_mask + +infile = '/gss/gss_work/DRES_OGS_BiGe/gbolzon/masks/V4/mesh_mask_V1INGV.nc' + +outfile = 'meshmask_INGVfor_ogstm.nc' + +M = c_mask.OrigMask(infile) + +lon_cut = 0 +depth_cut=0 + +free_surface=True +Biscay_land = False +c_mask.create_meshmask_nc(M,outfile,lon_cut,depth_cut,Biscay_land,free_surface) + +outfile = 'meshmask.nc' + +lon_cut = 149 # or M.getCutLocation(-8.875) +depth_cut=2 + +Biscay_land = True +c_mask.create_meshmask_nc(M,outfile,lon_cut,depth_cut,Biscay_land,free_surface) + + +M.nc_handler.close() +print('DONE!') diff --git a/FORCINGS/degradation/main_MITgcm.py b/FORCINGS/degradation/main_MITgcm.py new file mode 100644 index 0000000..94b6046 --- /dev/null +++ b/FORCINGS/degradation/main_MITgcm.py @@ -0,0 +1,16 @@ +# SECTION FOR REDUCTION OF 1/16 meshmask +# nemo3.4 class reader +import create_meshmask_nc_MITgcm as c_mask + +infile = '/home/plazzari/data/ENEA/grid.t001.nc' + +outfile = 'meshmask_MITgcm.nc' + +M = c_mask.OrigMask(infile) + +free_surface=True +c_mask.create_meshmask_nc(M,outfile,free_surface) + + +M.nc_handler.close() +print('DONE!') diff --git a/FORCINGS/degradation/reducer.py b/FORCINGS/degradation/reducer.py new file mode 100644 index 0000000..0d58199 --- /dev/null +++ b/FORCINGS/degradation/reducer.py @@ -0,0 +1,579 @@ +import numpy as np +import netCDF4 + + +class mesh(): + def __init__(self,filename=None, isFreeSurface=True): + if filename is not None: + NCin = netCDF4.Dataset(filename,'r') + self.jpi=NCin.dimensions['x'].size + self.jpj=NCin.dimensions['y'].size + self.jpk =NCin.dimensions['z'].size + + self.glamt = np.array(NCin.variables['glamt'][0,0,:,:]) + self.glamu = np.array(NCin.variables['glamu'][0,0,:,:]) + self.glamv = np.array(NCin.variables['glamv'][0,0,:,:]) + self.gphit = np.array(NCin.variables['gphit'][0,0,:,:]) + self.gphiu = np.array(NCin.variables['gphiu'][0,0,:,:]) + self.gphiv = np.array(NCin.variables['gphiv'][0,0,:,:]) + + + self.tmask = np.array(NCin.variables['tmask'][0,:,:,:]).astype(np.bool) + self.umask = np.array(NCin.variables['umask'][0,:,:,:]).astype(np.bool) + self.vmask = np.array(NCin.variables['vmask'][0,:,:,:]).astype(np.bool) + + if isFreeSurface: + self.e3t_0 = np.array(NCin.variables['e3t_0'][0,:,:,:])#_0 a 1/24 + self.e3u_0 = np.array(NCin.variables['e3u_0'][0,:,:,:]) + self.e3v_0 = np.array(NCin.variables['e3v_0'][0,:,:,:]) + self.e3w_0 = np.array(NCin.variables['e3w_0'][0,:,:,:]) + + self.e1t = np.array(NCin.variables['e1t'][0,0,:,:]) + self.e1u = np.array(NCin.variables['e1u'][0,0,:,:]) + self.e1v = np.array(NCin.variables['e1v'][0,0,:,:]) + self.e2t = np.array(NCin.variables['e2t'][0,0,:,:]) + self.e2u = np.array(NCin.variables['e2u'][0,0,:,:]) + self.e2v = np.array(NCin.variables['e2v'][0,0,:,:]) + self.gdept = np.array(NCin.variables['gdept'][0,:,0,0]) + self.gdepw = np.array(NCin.variables['gdepw'][0,:,0,0]) + self.nav_lev = np.array(NCin.variables['nav_lev']) + NCin.close() + else: + self.jpi=None + self.jpj=None + self.jpk=None + + self.glamt = None + self.glamu = None + self.glamv = None + self.gphit = None + self.gphiu = None + self.gphiv = None + + + self.gmask = None + self.gmask = None + self.gmask = None + self.e3t_0 = None + if isFreeSurface: + self.e3u_0 = None + self.e3v_0 = None + self.e3w_0 = None + self.e1t = None + self.e1u = None + self.e1v = None + self.e2t = None + self.e2u = None + self.e2v = None + self.gdept = None + self.gdepw = None + self.nav_lev = None + self.finecells = None + + +class reducer(): + def __init__(self, isFreeSurface=True): + self.time = 1 + self.z_a = 1 + self.x_a = 1 + self.y_a = 1 + self.isFreeSurface = isFreeSurface + self.jpi = None + self.jpj = None + self.jpk = None + + def read_fine(self,filename): + self.Mfine = mesh(filename=filename,isFreeSurface=self.isFreeSurface) + + + + + def reducing_by_boxes_2d_sum(self, M2d, I, J, direction=1): + ''' Valid for e2u, e2v etc''' + jpi = len(I) + jpj = len(J) + step = self.step + m2d = np.zeros((jpj,jpi),np.double) + for ji in range(jpi): + for jj in range(jpj): + j=J[jj] + i=I[ji] + if direction ==1: + m2d[jj,ji] = np.sum(M2d[j+step-1,i:i+step]) + if direction ==2: + m2d[jj,ji] = np.sum(M2d[j:j+step,i+step-1]) + return m2d + + def reducing_by_boxes_3d_mean(self, M3d, W3d, I, J, direction='u'): + ''' Valid for e3u, e3v ''' + jpi = len(I) + jpj = len(J) + jpk = self.jpk + step = self.step + m3d = np.zeros((jpk,jpj,jpi),np.double) + + for jk in range(jpk): + for ji in range(jpi): + for jj in range(jpj): + j=J[jj] + i=I[ji] + if direction =="u": + e = M3d[jk,j:j+step,i+step-1] + w = W3d[j:j+step,i+step-1] + ew = (e*w).sum() + m3d[jk,jj,ji] = ew/w.sum() + if direction =="v": + e = M3d[jk,j+step-1,i:i+step] + w = W3d[j+step-1,i:i+step] + ew = (e*w).sum() + m3d[jk,jj,ji] = ew/w.sum() + return m3d + + + def reducing_by_boxes_double(self, M3d, I, J): + ''' + valid for e3u_0, e3v_0 + mean over box + meglio fare solo i punti d'acqua? + ''' + jpi = self.jpi + jpj = self.jpj + jpk = self.jpk + step = self.step + m3d = np.zeros((jpk,jpj,jpi),np.double) + for ji in range(jpi): + for jj in range(jpj): + j=J[jj] + i=I[ji] + for jk in range(jpk): + m3d[jk,jj,ji] = np.mean(M3d[jk,j:j+step,i:i+step]) + return m3d + + + + def reducing_by_boxes_bool(self,Bool3d, I, J, direction='u'): + ''' + Reduction of mask arrays + ''' + jpi = self.jpi + jpj = self.jpj + jpk = self.jpk + step = self.step + bool3d = np.zeros((jpk,jpj,jpi),np.bool) + count3d= np.zeros((jpk,jpj,jpi),np.int32) + for ji in range(jpi): + for jj in range(jpj): + j=J[jj] + i=I[ji] + if direction == 'u': + for jk in range(jpk): + bool3d[jk,jj,ji] = np.any(Bool3d[jk,j:j+step,i+step-1]) + if direction == 'v': + for jk in range(jpk): + bool3d[jk,jj,ji] = np.any(Bool3d[jk,j+step-1,i:i+step]) + if direction == 't': + for jk in range(jpk): + bool3d[jk,jj,ji] = np.any(Bool3d[jk,j:j+step,i:i+step]) + count3d[jk,jj,ji]= np.sum(Bool3d[jk,j:j+step,i:i+step]) + return bool3d, count3d + +# def reducing_by_boxes_2d_sum(self, M2d, I, J, direction="u"): +# ''' Valid for e2u, e2v etc''' +# +# m2d = np.zeros((jpj,jpi),np.double) +# for ji in range(jpi): +# for jj in range(jpj): +# j=J[jj] +# i=I[ji] +# if direction =="u": +# m2d[jj,ji] = np.sum(M2d[j:j+step,i+step-1]) +# if direction =="v": +# m2d[jj,ji] = np.sum(M2d[j+step-1,i:i+step]) +# return m2d + def reducing_by_boxes_2d_mean(self, M2d, I, J, direction="lam"): + ''' + Valid for glamu, glamt + lam = longitude index + phy = latitude index + U box + + * V box * * * T box: * * * + + + * + + + * * * + + + * + + + * * * + ''' + step = self.step + m2d = np.zeros((self.jpj,self.jpi),np.double) + for ji in range(self.jpi): + for jj in range(self.jpj): + j=J[jj] + i=I[ji] + if direction =="lam": + m2d[jj,ji] = np.mean(M2d[j:j+step, i+step-1]) + if direction =="phi": + m2d[jj,ji] = np.mean(M2d[j+step-1 ,i:i+step]) +# if direction =="t": +# m2d[jj,ji] = np.mean(M2d[j:j+step, i:i+step]) + + return m2d + + def reduce(self, step): + ''' + Generates a mask object reduced by step + Arguments: + * step * integer + + The coarse mesh matches exactly the fine one at East and North. + Returns: + * m * mymesh object + * I * numpy array of integers, such as coarse_lon[k] = fine_lon[I[k]] + I[k] is the left index of the box. + * J * idem, for latitude. + ''' + + print "Start Mesh reduction" + m = mesh(filename=None,isFreeSurface=self.isFreeSurface) + self.step = step + + m.jpi, i_start = divmod(self.Mfine.jpi,step) + m.jpj, j_start = divmod(self.Mfine.jpj,step) + m.jpk = self.Mfine.jpk + + self.jpi = m.jpi + self.jpj = m.jpj + self.jpk = m.jpk + + I_glo = np.arange(self.Mfine.jpi) + J_glo = np.arange(self.Mfine.jpj) + I = I_glo[i_start::step] + J = J_glo[j_start::step] + + + m.tmask, m.finecells = self.reducing_by_boxes_bool(self.Mfine.tmask, I, J, 't') + m.umask, _ = self.reducing_by_boxes_bool(self.Mfine.umask, I, J, 'u') + m.vmask, _ = self.reducing_by_boxes_bool(self.Mfine.vmask, I, J, 'v') + + + if self.isFreeSurface: + print "Start e3 section" + m.e3t_0 = self.reducing_by_boxes_double(self.Mfine.e3t_0, I, J) + m.e3u_0 = self.reducing_by_boxes_3d_mean(self.Mfine.e3u_0, self.Mfine.e2u, I, J, 'u') +# m.e3u_0 = self.reducing_by_boxes_double(self.Mfine.e3u_0, I, J) + m.e3v_0 = self.reducing_by_boxes_3d_mean(self.Mfine.e3v_0, self.Mfine.e1v, I, J, 'v') +# m.e3v_0 = self.reducing_by_boxes_double(self.Mfine.e3v_0, I, J) + m.e3w_0 = self.reducing_by_boxes_double(self.Mfine.e3w_0, I, J) + + print "Start e1 e2 section" + m.e1t = self.reducing_by_boxes_2d_sum(self.Mfine.e1t, I, J, 1) + m.e1u = self.reducing_by_boxes_2d_sum(self.Mfine.e1u, I, J, 1) + m.e1v = self.reducing_by_boxes_2d_sum(self.Mfine.e1v, I, J, 1) + m.e2t = self.reducing_by_boxes_2d_sum(self.Mfine.e2t, I, J, 2) + m.e2u = self.reducing_by_boxes_2d_sum(self.Mfine.e2u, I, J, 2) + m.e2v = self.reducing_by_boxes_2d_sum(self.Mfine.e2v, I, J, 2) + + print "Start phi lam section" + m.glamt = self.reducing_by_boxes_2d_mean(self.Mfine.glamt, I, J,'lam') + m.glamu = self.reducing_by_boxes_2d_mean(self.Mfine.glamu, I, J,'lam') + m.glamv = self.reducing_by_boxes_2d_mean(self.Mfine.glamv, I, J,'lam') + m.gphit = self.reducing_by_boxes_2d_mean(self.Mfine.gphit, I, J,'phi') + m.gphiu = self.reducing_by_boxes_2d_mean(self.Mfine.gphiu, I, J,'phi') + m.gphiv = self.reducing_by_boxes_2d_mean(self.Mfine.gphiv, I, J,'phi') + + m.gdept = self.Mfine.gdept + m.gdepw = self.Mfine.gdepw + m.nav_lev = self.Mfine.nav_lev + return m,I,J + + def nearest_extend3D(self,M3d): + ''' + Argument + * M3d * a numpy array of dimension (jpk, jpj,jpi) + + Returns: + * out * a numpy array of dimension (jpk,jpj+2,jpi+2), + with frame of one cell all around, in longitude and latitude + ''' + jpk, jpj, jpi = M3d.shape + out = np.ones((jpk,jpj+2,jpi+2), dtype=M3d.dtype) + out[:,1:-1, 1:-1] = M3d + out[:, 1:-1 ,0] = M3d[:,:, 0] + out[:, 1:-1,-1] = M3d[:,:,-1] + out[:, 0 ,1:-1] = M3d[:, 0,:] + out[:,-1 ,1:-1] = M3d[:,-1,:] + return out + + def nearest_extend1D(self,array): + nP=len(array)+2 + out = np.zeros((nP),dtype=array.dtype) + out[0] = array[ 0] + out[1:-1] = array + out[-1] = array[-1] + return out + + def nearest_extend2D(self,M2d): + ''' + Argument: + * M2d * a numpy array having dimension (jpj,jpi) + + Returns: + * out * a numpy array having dimension (jpj+2, jpi+2) + with a frame of one cell all around, obtainded by + nearest extrapolation + ''' + + jpj, jpi = M2d.shape + out = np.zeros((jpj+2,jpi+2), dtype=M2d.dtype) + out[1:-1, 1:-1] = M2d + out[: ,0] = self.nearest_extend1D(M2d[:, 0]) + out[:,-1] = self.nearest_extend1D(M2d[:,-1]) + out[0, :] = self.nearest_extend1D(M2d[ 0,:]) + out[-1,:] = self.nearest_extend1D(M2d[-1,:]) + return out + def extrap_2d(self,M2d): + ''' + Argument: + * M2d * a numpy array having dimension (jpj,jpi) + + Returns: + * out * a numpy array having dimension (jpj+2, jpi+2) + with a frame of one cell all around, obtainded by + linear extrapolation + The idea is: + + F(n+1) = f(n) + ( f(n) - (f(n-1)) ) = 2*f(n) - f(n-1) + F(-1) = f(0) - ( f(1) - f(0) ) = 2*f(0) - f(1) + ''' + jpj, jpi = M2d.shape + out = np.zeros((jpj+2,jpi+2), dtype=M2d.dtype) + out[1:-1, 1:-1] = M2d + out[:, 0] = self.nearest_extend1D( 2*M2d[:, 0] - M2d[:, 1] ) + out[:,-1] = self.nearest_extend1D( 2*M2d[:,-1] - M2d[:,-2] ) + out[0, :] = self.nearest_extend1D( 2*M2d[0, :] - M2d[1, :] ) + out[-1,:] = self.nearest_extend1D( 2*M2d[-1,:] - M2d[-2,:] ) + return out + + def expand(self,m,I,J): + ''' + Argument: + * m * a mymesh object, having dimensions jpk,jpj,jpi + * I * numpy array[jpi] of indexes + * J * numpy array[jpi] of indexes + Returns: + + * M * a mymesh object, having dimensions jpk, jpj+2, jpi+2 + obtained by generating a frame of one cell all around, + in longitude and latitude. + * I * numpy array of indexes + ''' + M = mesh(filename=None, isFreeSurface=self.isFreeSurface) + M.jpi = m.jpi+2 + M.jpj = m.jpj+2 + M.jpk = m.jpk + + I_exp = np.ones((M.jpi),np.int32)*(-999) + J_exp = np.ones((M.jpj),np.int32)*(-999) + I_exp[1:-1]=I + J_exp[1:-1]=J + + # here pure frame of zeros + M.tmask = np.zeros((M.jpk,M.jpj,M.jpi),np.bool) + M.umask = np.zeros((M.jpk,M.jpj,M.jpi),np.bool) + M.vmask = np.zeros((M.jpk,M.jpj,M.jpi),np.bool) + M.tmask[:,1:-1, 1:-1] = m.tmask + M.umask[:,1:-1, 1:-1] = m.umask + M.vmask[:,1:-1, 1:-1] = m.vmask + M.finecells = np.zeros((M.jpk,M.jpj,M.jpi),np.int32) + M.finecells[:,1:-1, 1:-1] = m.finecells + #------------------------------------ + if self.isFreeSurface: + M.e3t_0 = self.nearest_extend3D(m.e3t_0) + M.e3u_0 = self.nearest_extend3D(m.e3u_0) + M.e3v_0 = self.nearest_extend3D(m.e3v_0) + M.e3w_0 = self.nearest_extend3D(m.e3w_0) + M.e1t = self.nearest_extend2D(m.e1t) + M.e1u = self.nearest_extend2D(m.e1u) + M.e1v = self.nearest_extend2D(m.e1v) + M.e2t = self.nearest_extend2D(m.e2t) + M.e2u = self.nearest_extend2D(m.e2u) + M.e2v = self.nearest_extend2D(m.e2v) + + + M.glamu = self.extrap_2d(m.glamu) + M.glamv = self.extrap_2d(m.glamv) + M.glamt = self.extrap_2d(m.glamt) + M.gphiu = self.extrap_2d(m.gphiu) + M.gphiv = self.extrap_2d(m.gphiv) + M.gphit = self.extrap_2d(m.gphit) + + M.gdept = m.gdept + M.gdepw = m.gdepw + M.nav_lev = m.nav_lev + + return M,I_exp,J_exp + + + def cut_med(self,m,lon_cut=-8.875,depth_cut=0,biscay_land=True): + ''' + Cuts the atlantic buffer and unuseful bottom of the Med Sea + + Arguments : + * m * a mesh object + * lon_cut * float indicating the first longitude of returned mesh + * depth_cut * integer, number of layers of bottom which are land points + * biscay_land * logical flag; if True, Biscay Gulf will be set as land + + Returns : + * M * mesh object, ready for ogstm + * lon_ind_cut * integer + ''' + lon_ind_cut =np.argmin(np.abs(m.glamt[0,:]-lon_cut)) + + M = mesh(filename=None, isFreeSurface=self.isFreeSurface) + + M.jpi = m.jpi - lon_ind_cut + M.jpj = m.jpj + M.jpk = m.jpk - depth_cut + + + + M.tmask = m.tmask[:M.jpk,:,lon_ind_cut:] + M.umask = m.umask[:M.jpk,:,lon_ind_cut:] + M.vmask = m.vmask[:M.jpk,:,lon_ind_cut:] + + M.tmask[:,:,0] = False + M.vmask[:,:,0] = False + M.tmask[:,:,0] = False + M.finecells = m.finecells[:M.jpk,:,lon_ind_cut:] + if self.isFreeSurface: + M.e3t_0 = m.e3t_0[:M.jpk,:,lon_ind_cut:] + M.e3u_0 = m.e3u_0[:M.jpk,:,lon_ind_cut:] + M.e3v_0 = m.e3v_0[:M.jpk,:,lon_ind_cut:] + M.e3w_0 = m.e3w_0[:M.jpk,:,lon_ind_cut:] + M.e1t = m.e1t[:,lon_ind_cut:] + M.e1u = m.e1u[:,lon_ind_cut:] + M.e1v = m.e1v[:,lon_ind_cut:] + M.e2t = m.e2t[:,lon_ind_cut:] + M.e2u = m.e2u[:,lon_ind_cut:] + M.e2v = m.e2v[:,lon_ind_cut:] + + + M.glamu = m.glamu[:,lon_ind_cut:] + M.glamv = m.glamv[:,lon_ind_cut:] + M.glamt = m.glamt[:,lon_ind_cut:] + M.gphiu = m.gphiu[:,lon_ind_cut:] + M.gphiv = m.gphiv[:,lon_ind_cut:] + M.gphit = m.gphit[:,lon_ind_cut:] + M.gdept = m.gdept[:M.jpk] + M.gdepw = m.gdepw[:M.jpk] + M.nav_lev = m.nav_lev[:M.jpk] + + if biscay_land: + bool1 = (M.glamt < 0) & (M.gphit > 42) + bool2 = (M.glamt < -6.) & (M.gphit > 37.25) + land = bool1 | bool2 + + for k in range(M.jpk): + tmask = M.tmask[k,:,:] + umask = M.umask[k,:,:] + vmask = M.vmask[k,:,:] + tmask[land] = False + umask[land] = False + vmask[land] = False + return M, lon_ind_cut + + + def dumpfile(self, m, outfile): + ncOUT=netCDF4.Dataset(outfile,"w"); + + ncOUT.createDimension('x',m.jpi); + ncOUT.createDimension('y',m.jpj); + ncOUT.createDimension('z',m.jpk); + ncOUT.createDimension('time',self.time) + + ncOUT.createDimension('x_a',self.x_a); + ncOUT.createDimension('y_a',self.y_a); + ncOUT.createDimension('z_a',self.z_a); + + + ncvar = ncOUT.createVariable('e1t' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[0,0,:,:] = m.e1t ; + ncvar = ncOUT.createVariable('e1u' ,'d',('time','z_a', 'y', 'x') ); ncvar[0,0,:,:] = m.e1u ; + ncvar = ncOUT.createVariable('e1v' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.e1v ; + + ncvar = ncOUT.createVariable('e2t' ,'d',('time','z_a', 'y', 'x') ) ; ncvar[0,0,:,:] = m.e2t ; + ncvar = ncOUT.createVariable('e2u' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.e2u ; + ncvar = ncOUT.createVariable('e2v' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.e2v ; + + if self.isFreeSurface: + ncvar = ncOUT.createVariable('e3t_0' ,'d',('time','z', 'y', 'x')) ; ncvar[0,:,:,:] = m.e3t_0 + ncvar = ncOUT.createVariable('e3u_0' ,'d',('time','z', 'y', 'x')) ; ncvar[0,:,:,:] = m.e3u_0 + ncvar = ncOUT.createVariable('e3v_0' ,'d',('time','z', 'y', 'x')) ; ncvar[0,:,:,:] = m.e3v_0 + ncvar = ncOUT.createVariable('e3w_0' ,'d',('time','z', 'y', 'x')) ; ncvar[0,:,:,:] = m.e3w_0 + + + + ncvar = ncOUT.createVariable('gdept' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[0,:,0,0] = m.gdept + ncvar = ncOUT.createVariable('gdepw' ,'d',('time','z', 'y_a', 'x_a')) ; ncvar[0,:,0,0] = m.gdepw + + ncvar = ncOUT.createVariable('glamt' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.glamt + ncvar = ncOUT.createVariable('glamu' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.glamu + ncvar = ncOUT.createVariable('glamv' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.glamv + + ncvar = ncOUT.createVariable('gphit' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.gphit + ncvar = ncOUT.createVariable('gphiu' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.gphiu + ncvar = ncOUT.createVariable('gphiv' ,'d',('time','z_a', 'y', 'x')) ; ncvar[0,0,:,:] = m.gphiv + ncvar = ncOUT.createVariable('nav_lat','f',('y','x')) ; ncvar[:] = m.gphit + ncvar = ncOUT.createVariable('nav_lev' ,'f',('z',)) ; ncvar[:] = m.nav_lev + ncvar = ncOUT.createVariable('nav_lon','f',('y','x')) ; ncvar[:] = m.glamt + ncvar = ncOUT.createVariable('tmask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = m.tmask + ncvar = ncOUT.createVariable('umask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = m.umask + ncvar = ncOUT.createVariable('vmask' ,'d',('time','z', 'y', 'x') ) ; ncvar[:] = m.vmask + ncvar = ncOUT.createVariable('finecells' ,'i',('time','z', 'y', 'x') ); ncvar[:] = m.finecells + ncOUT.close() + + def nearest_2d_cutting_along_slices(self, M2d,I,J): + ''' + Performs a nearest 2d interp over two regular grids, based on precalculated indexes. + Returns a reduced 2d array, having dimension (len(J),len(I)) + + Arguments: + * M2d * 2d numpy array + * I * integer numpy array + M2d longitude[I[k]] is the nearest of M2d_reduced longitude[k] + * J * integer numpy array, same meaning. + + + Returns: + * M2d_reduced * 2d numpy array + + ''' + jpi = len(I) + jpj = len(J) + M2d_reduced = np.zeros((jpj,jpi), dtype=M2d.dtype) + for ji in range(jpi): + for jj in range(jpj): + M2d_reduced[jj,ji] = M2d[J[jj],I[ji]] + return M2d_reduced + +if __name__=="__main__": + import numpy as np + isFreeSurface=True + inputmesh="meshmask.nc" + import pickle + + # just to save time for tests + #fid=open("temp.dat") + #LIST=pickle.load(fid) + #[R, Mred, I, J ] = LIST + #fid.close() + + R =reducer(isFreeSurface) + R.read_fine(inputmesh) + passo=4 + Mred,I,J = R.reduce(passo) + Mexp,I,J = R.expand(Mred,I,J) + outputmesh="meshmask_6125.nc" + R.dumpfile(Mexp, outputmesh) + np.savetxt('South_West_I_indexes.txt',I,'%d') # of mesh not expanded + np.savetxt('South_West_J_indexes.txt',J,'%d') + + + + + diff --git a/FORCINGS/degradation/restart_interpolator.py b/FORCINGS/degradation/restart_interpolator.py new file mode 100644 index 0000000..f20d0a7 --- /dev/null +++ b/FORCINGS/degradation/restart_interpolator.py @@ -0,0 +1,113 @@ +import numpy as np +from commons.mask import Mask +from commons.dataextractor import DataExtractor +import netCDF4 +import glob,os + +def weighted_mean(values,weights): + return (values*weights).sum()/weights.sum() +def interp3d(M3dfine, Maskfine,MaskCoarse,I,J): + ''' + Degradates a 3d field defined on cell center, using area weighted mean. + + Arguments: + * M3dfine * 3d numpy array + * Maskfine * Mask object + * MaskCoarse * Mask object + * I * numpy array 1d of integers, coarse_lon[k] = fine_lon[I[k]] + * J * idem + + Returns: + * M3dCoarse * 3d numpy array + ''' + step = I[3]-I[2] + jpk, jpj, jpi = MaskCoarse.shape + + M3dCoarse=np.zeros((jpk,jpj,jpi),np.float64) + for jk in range(jpk): + for ji in range(jpi): + for jj in range(jpj): + if MaskCoarse.mask[jk,jj,ji]: + i=I[ji] + j=J[jj] + mask = Maskfine.mask[jk,j:j+step,i:i+step] + values= M3dfine[jk,j:j+step,i:i+step] + area = Maskfine.area[ j:j+step,i:i+step] + M3dCoarse[jk,jj,ji] = weighted_mean(values[mask], area[mask]) + return M3dCoarse + + +def RSTwriter(outfile, var,rst, TheMask): + rst[~TheMask.mask] = 1.e+20 + jpk, jpj, jpi = TheMask.shape + ncOUT=netCDF4.Dataset(outfile,"w", format="NETCDF4") + ncOUT.createDimension('x',jpi); + ncOUT.createDimension('y',jpj); + ncOUT.createDimension('z',jpk); + ncOUT.createDimension('time',1) + + TRN = 'TRN' + var; + ncvar = ncOUT.createVariable('nav_lon' ,'d',('y','x') ); ncvar[:] = TheMask.xlevels + ncvar = ncOUT.createVariable('nav_lat' ,'d',('y','x') ); ncvar[:] = TheMask.ylevels + ncvar = ncOUT.createVariable('nav_lev' ,'d',('z') ); ncvar[:] = TheMask.zlevels + ncvar = ncOUT.createVariable('time' ,'d',('time',) ); ncvar = 1.; + ncvar = ncOUT.createVariable(TRN ,'d',('time','z','y','x') ); ncvar[:] = rst; + + + setattr(ncOUT.variables[TRN] ,'missing_value',1.e+20 ); + setattr(ncOUT.variables['time'],'Units' ,'seconds since 1582-10-15 00:00:00'); + setattr(ncOUT ,'TimeString' ,'20010101-00:00:00'); + ncOUT.close() + + + + +I = np.loadtxt('ogstm_South_West_I_indexes.txt', np.int32) # of M_med mesh +J = np.loadtxt('ogstm_South_West_J_indexes.txt', np.int32) + +delta = 871 - 722 +I = I-delta +I[I<0] = 0 + + + +MaskCoarse = Mask("/marconi_work/OGS_dev_0/DEGRADATION_4_70/PREPROC/preproc_meshgen_forcings/mesh_gen/meshmask_470.nc") +Mask16 = Mask("/marconi_scratch/userexternal/gbolzon0/2017_RA16/wrkdir/MODEL/meshmask.nc") +# correction of spain point, setting tmask=1 where 1/4 mask is 1 +orig_tmask=Mask16.mask.copy() +J_corr=[149, 149, 149, 150, 150, 151, 151, 152, 153] +I_corr=[137, 138, 139, 138, 139, 138, 139, 139, 139] + +ncorr=len(J_corr) +for i in range(ncorr): Mask16.mask[:5, J_corr[i], I_corr[i] ] = True + + +INPUTDIR="/marconi_scratch/userexternal/gbolzon0/REA_KD490/INIT/" +OUTPUTDIR="/marconi_work/OGS_dev_0/DEGRADATION_4_70/PREPROC/preproc_meshgen_forcings/mesh_gen/INIT_INTERP/" + +RST_LIST=glob.glob(INPUTDIR +"RST*nc") + + +for filename in RST_LIST: + basename=os.path.basename(filename) + outfile = OUTPUTDIR + basename + var = basename.rsplit(".")[2] + print outfile + + RST = DataExtractor(Mask16, filename, "TRN" + var).values + values=RST[orig_tmask] + print "ORIG: " , values.min(), values.max() +# correction of spain point, setting values where tmask was 0 ####### + nearest_values= RST[:5,149,140] + for i in range(ncorr): + RST[:5, J_corr[i], I_corr[i] ] = nearest_values +###################################################################### + + rst = interp3d(RST, Mask16, MaskCoarse, I, J) + values = rst[MaskCoarse.mask] + if np.isnan(values).any(): + raise ValueError + print "INTERPOLATED: " , values.min(), values.max() + + + RSTwriter(outfile, var, rst, MaskCoarse) diff --git a/FORCINGS/degradation/river_reindexer.py b/FORCINGS/degradation/river_reindexer.py new file mode 100644 index 0000000..d16daa6 --- /dev/null +++ b/FORCINGS/degradation/river_reindexer.py @@ -0,0 +1,71 @@ +# This script prints as standart output two columns of indexes +# that have to be putted in rivers xls file for reduced mesh + +import numpy as np +class indexer(): + def __init__(self,I,J): + self.I = I + self.J = J + self.jpi = len(I) + self.jpj = len(J) + def switch_lon_index(self,ifine): + for k in range(self.jpi): + if self.I[k]>ifine: + coarse_ind = k + break + return coarse_ind-1 + def switch_lat_index(self,jfine): + for k in range(self.jpj): + if self.J[k]>jfine: + coarse_ind = k + break + return coarse_ind-1 + + + + + +I_ind = np.loadtxt('ogstm_South_West_I_indexes.txt',np.int32) +J_ind = np.loadtxt('ogstm_South_West_J_indexes.txt',np.int32) +Indexer = indexer(I_ind,J_ind) + + +# Prese le colonne I,J del file river_8_on_meshgrid_v2.xlsx e messe in nei files I_river_16_from_xls.txt J_river_16_from_xls.txt +I_river_16 = np.loadtxt("I_river_16_from_xls.txt",np.int32) -1 + 222 # because plazzari removed 222 +J_river_16 = np.loadtxt("J_river_16_from_xls.txt",np.int32) -1 + + +riverPoints = len(I_river_16) +I_red = np.zeros((riverPoints,),np.int32) +J_red = np.zeros((riverPoints,),np.int32) + +for iR in range(riverPoints): + i = Indexer.switch_lon_index(I_river_16[iR]) + j = Indexer.switch_lat_index(J_river_16[iR]) + print "%d\t%d" %(i+1,j+1) + I_red[iR] = i + J_red[iR] = j + + +from commons.mask import Mask +TheMask = Mask("meshmask_469.nc") + +from layer_integral.mapplot import mapplot +tmask = TheMask.mask_at_level(0) +fig, ax = mapplot({'data':tmask, 'clim':[0,1]}) + +ax.plot(I_red, J_red,'w.') +fig.show() + + + + + + + + + + + + + diff --git a/IC/Ext_without_dres.py b/IC/Ext_without_dres.py index 6b22050..bc55174 100644 --- a/IC/Ext_without_dres.py +++ b/IC/Ext_without_dres.py @@ -4,6 +4,7 @@ from bitsea.commons.mask import Mask from bitsea.commons.dataextractor import DataExtractor from IC import RSTwriter, smoother +from dataclasses import dataclass date_out="20220701-00:00:00" date_in="20230701-00:00:00" @@ -12,66 +13,77 @@ MDS_dir=Path("/g100_work/OGS_devC/RA_24/PREPROC/RESTARTS/MDS") -Int_dir=Path("/g100_work/OGS_devC/RA_24/PREPROC/RESTARTS/Interim_2023") +Interim_dir=Path("/g100_work/OGS_devC/RA_24/PREPROC/RESTARTS/Interim_2023") # generated with /g100_work/OGS_prodC/OPA/Interim-dev/archive/202306/MODEL/RST.20230701.tar OUTDIR=Path("/g100_work/OGS_devC/RA_24/PREPROC/RESTARTS/generated") -DICT = {"chl":["P1l","P2l","P3l","P4l"], "phyc":["P1c","P2c","P3c","P4c"]} +@dataclass +class mds_var: + name: str + varlist: list + conversion : float -for var_mds in DICT.keys(): - VARLIST=DICT[var_mds] +MDS_VARLIST=[ + mds_var("chl",["P1l","P2l","P3l","P4l"],1.0), + mds_var("phyc",["P1c","P2c","P3c","P4c"],12.0) + ] + + +for var_mds in MDS_VARLIST: + VARLIST=var_mds.varlist Int=np.zeros(TheMask.shape,np.float64) for var in VARLIST: - Int += DataExtractor(TheMask,Int_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values - MDS = netcdf4.readfile(MDS_dir/ f"{date_in_mds}_d-OGS--PFTC-MedBFM3-MED-b20240216_re-sv05.00.nc", var_mds)[0,:] - ratio = MDS / Int[:,:,80:] + Int += DataExtractor(TheMask,Interim_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values + MDS = netcdf4.readfile(MDS_dir/ f"{date_in_mds}_d-OGS--PFTC-MedBFM3-MED-b20240216_re-sv05.00.nc", var_mds.name)[0,:] + ratio = MDS * var_mds.conversion / Int[:,:,80:] for var in VARLIST: - Int = DataExtractor(TheMask,Int_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values + Int = DataExtractor(TheMask,Interim_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values M = Int.copy() M[:,:,80:] = Int[:,:,80:]*ratio RSTwriter(OUTDIR / f"RST.{date_out}.{var}.nc",var,M,TheMask) - +import sys +sys.exit() var="O2o" -Int =DataExtractor(TheMask,Int_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values +Int =DataExtractor(TheMask,Interim_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values MDS = netcdf4.readfile(MDS_dir/ f"{date_in_mds}_d-OGS--BIOL-MedBFM3-MED-b20240216_re-sv05.00.nc", "o2") M=Int.copy() M[:,:,80:]=MDS RSTwriter(OUTDIR / f"RST.{date_out}.{var}.nc",var,M,TheMask) var="N3n" -Int =DataExtractor(TheMask,Int_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values +Int =DataExtractor(TheMask,Interim_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values MDS = netcdf4.readfile(MDS_dir/ f"{date_in_mds}_d-OGS--NUTR-MedBFM3-MED-b20240216_re-sv05.00.nc", "no3") M=Int.copy() M[:,:,80:]=MDS RSTwriter(OUTDIR / f"RST.{date_out}.{var}.nc",var,M,TheMask) var="N1p" -Int =DataExtractor(TheMask,Int_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values +Int =DataExtractor(TheMask,Interim_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values MDS = netcdf4.readfile(MDS_dir/ f"{date_in_mds}_d-OGS--NUTR-MedBFM3-MED-b20240216_re-sv05.00.nc", "po4") M=Int.copy() M[:,:,80:]=MDS RSTwriter(OUTDIR / f"RST.{date_out}.{var}.nc",var,M,TheMask) var="N4n" -Int =DataExtractor(TheMask,Int_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values +Int =DataExtractor(TheMask,Interim_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values MDS = netcdf4.readfile(MDS_dir/ f"{date_in_mds}_d-OGS--NUTR-MedBFM3-MED-b20240216_re-sv05.00.nc", "nh4") M=Int.copy() M[:,:,80:]=MDS RSTwriter(OUTDIR / f"RST.{date_out}.{var}.nc",var,M,TheMask) var="O3c" -Int =DataExtractor(TheMask,Int_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values +Int =DataExtractor(TheMask,Interim_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values MDS = netcdf4.readfile(MDS_dir/ f"{date_in_mds}_d-OGS--CARB-MedBFM3-MED-b20240216_re-sv05.00.nc", "dissic")*12000 M=Int.copy() M[:,:,80:]=MDS RSTwriter(OUTDIR / f"RST.{date_out}.{var}.nc",var,M,TheMask) var="O3h" -Int =DataExtractor(TheMask,Int_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values +Int =DataExtractor(TheMask,Interim_dir/f"RST.{date_in}.{var}.nc", f"TRN{var}").values MDS = netcdf4.readfile(MDS_dir/ f"{date_in_mds}_d-OGS--CARB-MedBFM3-MED-b20240216_re-sv05.00.nc", "talk")*1000 M=Int.copy() M[:,:,80:]=MDS From bf927e948ecc0e4a10aa6e3fe57dcf3887e59fce Mon Sep 17 00:00:00 2001 From: Giovanni Galli Date: Wed, 10 Dec 2025 10:53:10 +0100 Subject: [PATCH 02/15] new scripts: degrade_mesh.py and degrade_forcings.py do what they say in the title, degrade horizontal resolution of physics forcings and meshmask by an integer number. --- FORCINGS/degradation/degrade_forcings.py | 490 +++++++++++++++++++++++ FORCINGS/degradation/degrade_mesh.py | 467 +++++++++++++++++++++ 2 files changed, 957 insertions(+) create mode 100644 FORCINGS/degradation/degrade_forcings.py create mode 100644 FORCINGS/degradation/degrade_mesh.py diff --git a/FORCINGS/degradation/degrade_forcings.py b/FORCINGS/degradation/degrade_forcings.py new file mode 100644 index 0000000..5eb0fce --- /dev/null +++ b/FORCINGS/degradation/degrade_forcings.py @@ -0,0 +1,490 @@ +import numpy as np +import xarray as xr +from glob import glob +from mpi4py import MPI +from argparse import ArgumentParser +import yaml +from itertools import chain +# my stuff +import degrade_mesh as dm + +''' +degrades resolution of OGSTM physics forcings +on T, U, V, W grids +needs both original and degraded resolution meshmask.nc +(run degrade_mesh.py before) +input parameters (like directories etc.) go into a .yaml file +usage: python degrade_forcings.py -y degrade_forcings.yaml + +yaml file example: +maskfile: '/path/to/meshmask_in.nc' #original mesh +maskfile_d: '/path/to/coarse/meshmask_out.nc' #degraded mesh +ffdir: '/path/to/wrkdir/MODEL/FORCINGS/' +outdir: '/path/to/output/' +ndeg: 6 #nr of cells to join in i, j +y0: 2000 #start year +yE: 2020 #end year + +(takes about 3min per year of daily files on 1 node (256GB) and 16 cores, +if it doesn't die by OOM, which sometimes does and sometimes not. +If it does go with 8 cores instead) +''' + +def argument(): + parser = ArgumentParser() + parser.add_argument('--yamlfile','-y', + type=str, + required=True, + help='file with all parameters') + return parser.parse_args() + +def load_parameters(): + yamlfile=argument().yamlfile + # yamlfile = 'degrade_forcings.yaml' + with open(yamlfile) as f: + Params = yaml.load(f, Loader=yaml.Loader) + return(Params) + +def load_mesh_light(maskfile, ndeg=1): + ''' + load meshmask, expand fields + just those needed to degrade forcings + ''' + M = xr.open_dataset(maskfile) + M1 = {} + # + M1["e1t"] = dm.xpnd_wrap(M["e1t"], 'interp', ndeg) + M1["e1u"] = dm.xpnd_wrap(M["e1u"], 'interp', ndeg) + M1["e1v"] = dm.xpnd_wrap(M["e1v"], 'interp', ndeg) + M1["e2t"] = dm.xpnd_wrap(M["e2t"], 'edge', ndeg) + M1["e2u"] = dm.xpnd_wrap(M["e2u"], 'edge', ndeg) + M1["e2v"] = dm.xpnd_wrap(M["e2v"], 'edge', ndeg) + M1["e3t_0"] = dm.xpnd_wrap(M["e3t_0"], 'edge', ndeg) + M1["e3u_0"] = dm.xpnd_wrap(M["e3u_0"], 'edge', ndeg) + M1["e3v_0"] = dm.xpnd_wrap(M["e3v_0"], 'edge', ndeg) + M1["glamt"] = dm.xpnd_wrap(M["glamt"], 'interp', ndeg) + M1["glamu"] = dm.xpnd_wrap(M["glamu"], 'interp', ndeg) + M1["glamv"] = dm.xpnd_wrap(M["glamv"], 'interp', ndeg) + M1["gphit"] = dm.xpnd_wrap(M["gphit"], 'interp', ndeg) + M1["gphiu"] = dm.xpnd_wrap(M["gphiu"], 'interp', ndeg) + M1["gphiv"] = dm.xpnd_wrap(M["gphiv"], 'interp', ndeg) + M1["tmask"] = dm.xpnd_wrap(M["tmask"], 'edge', ndeg) #there is some error here! + M1["umask"] = dm.xpnd_wrap(M["umask"], 'edge', ndeg) #there is some error here! + M1["vmask"] = dm.xpnd_wrap(M["vmask"], 'edge', ndeg) #there is some error here! + # Bonus variable for degrading forcings + M1['h_column_t'] = M1['e3t_0'].sum(dim='z') + # + M.close() + M1 = xr.Dataset(M1) + return M1 + + +def load_coords_degraded(maskfile_d): + ''' + loads lat, lon coordinates on T, U, V grids + from already degraded mask + to avoid having to degrade each time the same arrays + need to run degrade_mesh.py first! + ''' + Md = xr.open_dataset(maskfile_d) + C = Md[["glamt", "glamu", "glamv", "gphit", "gphiu", "gphiv"]] + C['glamw'] = C['glamt'] + C['gphiw'] = C['gphit'] + return C + +def load_tfile(infile, ndeg=1): + ''' + loads t-grid forcing file, pads a rim all around it + to make j, i multiples of ndeg, + populates rim with edge values + (ugly, but it's the atlantic buffer that gets populated, + the model won't use those values anyway) + ''' + F = xr.open_dataset(infile) + F1 = {} + # most of these ogstm does not need (commented) + F1['nav_lat'] = dm.xpnd_wrap(F['nav_lat'], 'interp', ndeg) + F1['nav_lon'] = dm.xpnd_wrap(F['nav_lon'], 'interp', ndeg) + F1['votemper'] = dm.xpnd_wrap(F['votemper'], 'edge', ndeg) + F1['vosaline'] = dm.xpnd_wrap(F['vosaline'], 'edge', ndeg) + F1['sossheig'] = dm.xpnd_wrap(F['sossheig'], 'edge', ndeg) + #F1['sossh_ib'] = dm.xpnd_wrap(F['sossh_ib'], 'edge', ndeg) + #F1['sowaflup'] = dm.xpnd_wrap(F['sowaflup'], 'edge', ndeg) + #F1['soevapor'] = dm.xpnd_wrap(F['soevapor'], 'edge', ndeg) + #F1['soprecip'] = dm.xpnd_wrap(F['soprecip'], 'edge', ndeg) + F1['sorunoff'] = dm.xpnd_wrap(F['sorunoff'], 'edge', ndeg) + F1['soshfldo'] = dm.xpnd_wrap(F['soshfldo'], 'edge', ndeg) + #F1['sohefldo'] = dm.xpnd_wrap(F['sohefldo'], 'edge', ndeg) + #F1['solofldo'] = dm.xpnd_wrap(F['solofldo'], 'edge', ndeg) + #F1['sosefldo'] = dm.xpnd_wrap(F['sosefldo'], 'edge', ndeg) + #F1['solafldo'] = dm.xpnd_wrap(F['solafldo'], 'edge', ndeg) + #F1['somxl010'] = dm.xpnd_wrap(F['somxl010'], 'edge', ndeg) + # overwrite coordinates, else it complains + for vv in F1.keys(): + F1[vv].coords['nav_lat'].values[:] = F1['nav_lat'].values[:] + F1[vv].coords['nav_lon'].values[:] = F1['nav_lon'].values[:] + # + F1['deptht'] = F['deptht'] + F1['deptht_bounds'] = F['deptht_bounds'] + F1['time_instant'] = F['time_instant'] + F1['time_instant_bounds'] = F['time_instant_bounds'] + F1['time_counter'] = F['time_counter'] + F1['time_counter_bounds'] = F['time_counter_bounds'] + F1['time_centered'] = F['time_centered'] + F1['time_centered_bounds'] = F['time_centered_bounds'] + # + F1 = xr.Dataset(F1) + F.close() + return F1 + +def load_ufile(infile, ndeg=1): + ''' + as in load_tfile() + ''' + F = xr.open_dataset(infile) + F1 = {} + # + F1['nav_lat'] = dm.xpnd_wrap(F['nav_lat'], 'interp', ndeg) + F1['nav_lon'] = dm.xpnd_wrap(F['nav_lon'], 'interp', ndeg) + F1['vozocrtx'] = dm.xpnd_wrap(F['vozocrtx'], 'edge', ndeg) + F1['sozotaux'] = dm.xpnd_wrap(F['sozotaux'], 'edge', ndeg) + vlist = list(F1.keys()) + # overwrite coordinates, else it complains + for vv in vlist: + F1[vv].coords['nav_lat'].values[:] = F1['nav_lat'].values[:] + F1[vv].coords['nav_lon'].values[:] = F1['nav_lon'].values[:] + F1['depthu'] = F['depthu'] + F1['depthu_bounds'] = F['depthu_bounds'] + F1['time_instant'] = F['time_instant'] + F1['time_instant_bounds'] = F['time_instant_bounds'] + F1['time_counter'] = F['time_counter'] + F1['time_counter_bounds'] = F['time_counter_bounds'] + F1['time_centered'] = F['time_centered'] + F1['time_centered_bounds'] = F['time_centered_bounds'] + # + F1 = xr.Dataset(F1) + F.close() + return F1 + +def load_vfile(infile, ndeg=1): + ''' + as in load_tfile(), load_ufile + ''' + F = xr.open_dataset(infile) + F1 = {} + # + F1['nav_lat'] = dm.xpnd_wrap(F['nav_lat'], 'interp', ndeg) + F1['nav_lon'] = dm.xpnd_wrap(F['nav_lon'], 'interp', ndeg) + F1['vomecrty'] = dm.xpnd_wrap(F['vomecrty'], 'edge', ndeg) + F1['sometauy'] = dm.xpnd_wrap(F['sometauy'], 'edge', ndeg) + vlist = list(F1.keys()) + # overwrite coordinates, else it complains + for vv in vlist: + F1[vv].coords['nav_lat'].values[:] = F1['nav_lat'].values[:] + F1[vv].coords['nav_lon'].values[:] = F1['nav_lon'].values[:] + F1['depthv'] = F['depthv'] + F1['depthv_bounds'] = F['depthv_bounds'] + F1['time_instant'] = F['time_instant'] + F1['time_instant_bounds'] = F['time_instant_bounds'] + F1['time_counter'] = F['time_counter'] + F1['time_counter_bounds'] = F['time_counter_bounds'] + F1['time_centered'] = F['time_centered'] + F1['time_centered_bounds'] = F['time_centered_bounds'] + # + F1 = xr.Dataset(F1) + F.close() + return F1 + +def load_wfile(infile, ndeg=1): + ''' + as in load_tfile(), load_ufile(), load vfile() + ''' + F = xr.open_dataset(infile) + F1 = {} + # + F1['nav_lat'] = dm.xpnd_wrap(F['nav_lat'], 'interp', ndeg) + F1['nav_lon'] = dm.xpnd_wrap(F['nav_lon'], 'interp', ndeg) + F1['vovecrtz'] = dm.xpnd_wrap(F['vovecrtz'], 'edge', ndeg) + F1['votkeavt'] = dm.xpnd_wrap(F['votkeavt'], 'edge', ndeg) + # overwrite coordinates, else it complains + for vv in F1.keys(): + F1[vv].coords['nav_lat'].values[:] = F1['nav_lat'].values[:] + F1[vv].coords['nav_lon'].values[:] = F1['nav_lon'].values[:] + F1['depthw'] = F['depthw'] + F1['depthw_bounds'] = F['depthw_bounds'] + #F1['time_instant'] = F['time_instant'] + #F1['time_instant_bounds'] = F['time_instant_bounds'] + F1['time_counter'] = F['time_counter'] + F1['time_counter_bounds'] = F['time_counter_bounds'] + F1['time_centered'] = F['time_centered'] + F1['time_centered_bounds'] = F['time_centered_bounds'] + # + F1 = xr.Dataset(F1) + F.close() + return F1 + + +def load_fffile(infile, tuv, ndeg=1): + if tuv == 'T': + F1 = load_tfile(infile, ndeg) + elif tuv == 'U': + F1 = load_ufile(infile, ndeg) + elif tuv == 'V': + F1 = load_vfile(infile, ndeg) + elif tuv == 'W': + F1 = load_wfile(infile, ndeg) + return F1 + +def get_mask_fields(M): + ''' + gets the arrays needed to compute cell surface areas and volumes + those from the mask that do not change + so it doesn't have to retrieve them each time + ''' + h_column_t = M['h_column_t'].values[:] + tmask = M['tmask'].values[:] + umask = M['umask'].values[:] + vmask = M['vmask'].values[:] + e3t_0 = M['e3t_0'].values[:] + e3u_0 = M['e3u_0'].values[:] + e3v_0 = M['e3v_0'].values[:] + e1u = M['e1u'].values[:] + e2u = M['e2u'].values[:] + e1v = M['e1v'].values[:] + e2v = M['e2v'].values[:] + e1t = M['e1t'].values[:] + e2t = M['e2t'].values[:] + At = e1t * e2t #e1v * e2u + Aw = np.repeat(At, repeats=tmask.shape[1], axis=1) + return h_column_t, tmask, umask, vmask, e3t_0, e3u_0, e3v_0, e1u, e2u, e1v, e2v, e1t, e2t, At, Aw + +def get_nanmask(mask): + nanmask = np.ones(mask.shape) + nanmask[mask==0.0] = np.nan + return nanmask + +def get_weights(M, T): + ''' + given the meshmask and the current t-grid file (with free surface) + computes cell surface areas on U, V grids and cell volume on T grid + and surface area on T / W grid + algorithm from ogstm/src/IO/forcing_phys.f90 + (NB, this one pure numpy because xarray is not good at broadcasting) + ''' + # THIS GOES OUTSIDE THE FUNCTION, SO I DON'T DO IT EACH TIME I CALL IT + h_column_t, tmask, umask, vmask, e3t_0, e3u_0, e3v_0, e1u, e2u, e1v, e2v, e1t, e2t, At, Aw = get_mask_fields(M) + nan_tmask = get_nanmask(tmask) + nan_umask = get_nanmask(umask) + nan_vmask = get_nanmask(vmask) + At = At * nan_tmask[0,0,:,:] #used for degrading T-grid, needs NaN (weight only valid values) + # THIS GOES OUTSIDE THE FUNCTION, SO I DON'T DO IT EACH TIME I CALL IT + ssh = tmask * T['sossheig'].values[:] + correction_e3t = 1.0 + (ssh / h_column_t) + e3t = tmask * e3t_0 * correction_e3t + # + e1u_x_e2u = (e1u * e2u) + e1v_x_e2v = (e1v * e2v) + e1t_x_e2t = (e1t * e2t) + diff_e3t = e3t - e3t_0 + # + s0 = e1t_x_e2t * diff_e3t + s1 = np.zeros(s0.shape) + s2 = np.zeros(s0.shape) + s1[:,:,:,1:] = e1t_x_e2t[:,:,:,1:] * diff_e3t[:,:,:,1:] #THIS INDEXING SHOULD BE CORRECT + s2[:,:,1:,:] = e1t_x_e2t[:,:,1:,:] * diff_e3t[:,:,1:,:] #THIS INDEXING SHOULD BE CORRECT + # + e3u = e3u_0 + (0.5 * (umask / e1u_x_e2u) * (s0 + s1)) + e3v = e3v_0 + (0.5 * (vmask / e1v_x_e2v) * (s0 + s2)) + # + B = {} + B['e1v'] = xr.DataArray(e1v, name='e1v' , dims=('time','z_a','y','x')) + B['e2u'] = xr.DataArray(e2u, name='e2u' , dims=('time','z_a','y','x')) + B['At'] = xr.DataArray(At, name='At' , dims=('time','z_a','y','x')) + B['Aw'] = xr.DataArray(Aw, name='Aw' , dims=('time','z','y','x')) + B['Au'] = xr.DataArray((e2u * e3u), name='Au' , dims=('time','z','y','x')) + B['Av'] = xr.DataArray((e1v * e3v), name='Av' , dims=('time','z','y','x')) + B['V'] = xr.DataArray((nan_tmask * e1v * e2u * e3t), name='V' , dims=('time','z','y','x')) #for degrading T-grid, don't weight land values + B = xr.Dataset(B) + return B + +def degrade_wrap(X, degr_op, B=1.0, ndeg=1): + print('---1---') + X = dm.reshape_blocks(X, ndeg) + print('---2---') + X = degr_op(X, ndeg, B) + print('...') + return X + +def init_ds(D, C, tuv): + ''' + adds to degraded dataset the variables that are always the same: + nav_lat, nav_lon, depth, + and those that do not need regridding, time + ''' + Dd = {} + Dd[f'nav_lon'] = C[f'glam{tuv}'] + Dd[f'nav_lat'] = C[f'gphi{tuv}'] + Dd[f'depth{tuv}'] = D[f'depth{tuv}'] + Dd[f'depth{tuv}_bounds'] = D[f'depth{tuv}_bounds'] + #Dd['time_instant'] = D['time_instant'] + Dd['time_counter'] = D['time_counter'] + Dd['time_centered'] = D['time_centered'] + #Dd['time_instant_bounds'] = D['time_instant_bounds'] + Dd['time_counter_bounds'] = D['time_counter_bounds'] + Dd['time_centered_bounds'] = D['time_centered_bounds'] + if tuv != 'w': + # for some reason these are not in W-grid files + Dd['time_instant'] = D['time_instant'] + Dd['time_instant_bounds'] = D['time_instant_bounds'] + return Dd + +def degrade_V(V, B, C, ndeg=1): + ''' + degrades resolution of V-grid file + vomecrty + sometauy + ''' + B = B.rename({'time':'time_counter', 'z':'depthv'}) + #Vd = {} + Vd = init_ds(V, C, 'v') + # + Vd['vomecrty'] = degrade_wrap(V['vomecrty'], dm.vawmean_jstep, B['Av'], ndeg) + Vd['sometauy'] = degrade_wrap(V['sometauy'], dm.e1vwmean_jstep, B['e1v'], ndeg) + # + Vd = xr.Dataset(Vd) + return Vd + +def degrade_U(U, B, C, ndeg=1): + ''' + degrades resolution of U-grid file + vozocrtx + sozotaux + ''' + B = B.rename({'time':'time_counter', 'z':'depthu'}) + #Ud = {} + Ud = init_ds(U, C, 'u') + # + Ud['vozocrtx'] = degrade_wrap(U['vozocrtx'], dm.uawmean_istep, B['Au'], ndeg) + Ud['sozotaux'] = degrade_wrap(U['sozotaux'], dm.e2uwmean_istep, B['e2u'], ndeg) + # + Ud = xr.Dataset(Ud) + return Ud + +def degrade_W(W, B, C, ndeg=1): + ''' + degrades resolution of W-grid file + vovecrtz + votkeavt + ''' + B = B.rename({'time':'time_counter', 'z':'depthw'}) + #Wd = {} + Wd = init_ds(W, C, 'w') + # + # use dm.vwmean(), but actually uses 3D Aw as weight + Wd['vovecrtz'] = degrade_wrap(W['vovecrtz'], dm.vwmean, B['Aw'], ndeg) + Wd['votkeavt'] = degrade_wrap(W['votkeavt'], dm.vwmean, B['Aw'], ndeg) + # + Wd = xr.Dataset(Wd) + return Wd + + +def degrade_T(T, B, C, ndeg=1): + ''' + degrades resolution of T-grid file + ''' + B = B.rename({'time':'time_counter', 'z':'deptht'}) + Td = init_ds(T, C, 't') + #Td = {} + # + Td['votemper'] = degrade_wrap(T['votemper'], dm.vwmean, B['V'], ndeg) + Td['vosaline'] = degrade_wrap(T['vosaline'], dm.vwmean, B['V'], ndeg) + Td['sossheig'] = degrade_wrap(T['sossheig'], dm.awmean, B['At'], ndeg) + Td['soshfldo'] = degrade_wrap(T['soshfldo'], dm.awmean, B['At'], ndeg) + Td['sorunoff'] = degrade_wrap(T['sorunoff'], dm.awmean, B['At'], ndeg) + # + Td = xr.Dataset(Td) + return Td + +def degrade(F, tuv, B, C, outdir, outfile, ndeg=1): + if tuv == 'T': + Fd = degrade_T(F, B, C, ndeg) + elif tuv == 'U': + Fd = degrade_U(F, B, C, ndeg) + elif tuv == 'V': + Fd = degrade_V(F, B, C, ndeg) + elif tuv == 'W': + Fd = degrade_W(F, B, C, ndeg) + outfile = outdir+outfile + dm.dump_mesh(Fd, outfile) + return Fd + +def get_flist(tuvw, Params): + ffdir = Params['ffdir'] + y0 = Params['y0'] + yE = Params['yE'] + flist = [glob(f'{ffdir}/{YYYY}/??/{tuvw}*.nc') for YYYY in range(y0, yE+1)] + flist = sorted(list(chain.from_iterable(flist))) + return flist + +if __name__=='__main__': + # + try: + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nranks = comm.size + except: + comm = None + rank = 0 + nranks = 1 + print (rank,nranks) + # + Params = load_parameters() + # + maskfile = Params['maskfile'] + maskfile_d = Params['maskfile_d'] + ffdir = Params['ffdir'] + outdir = Params['outdir'] + y0 = Params['y0'] + yE = Params['yE'] + # + ndeg = Params['ndeg'] + # + flistt = get_flist('T', Params) + flistu = get_flist('U', Params) + flistv = get_flist('V', Params) + flistw = get_flist('W', Params) + # + # load mesh and expand lat, lon to make them multiples of ndeg + if rank==0: + print('loading mask') + M = load_mesh_light(maskfile, ndeg) #NB, here Au, Av, V, are calculated with e3tuv_0 + C = load_coords_degraded(maskfile_d) + else: + M = None + C = None + # + if nranks > 1: + M = comm.bcast(M, root=0) + C = comm.bcast(C, root=0) + # + # LOOP ON FILES + print('loading T, U, V files') + ziter = list(zip(flistt, flistu, flistv, flistw)) + for tfile, ufile, vfile, wfile in ziter[rank::nranks]: + V = load_fffile(vfile, 'V', ndeg) + U = load_fffile(ufile, 'U', ndeg) + W = load_fffile(wfile, 'W', ndeg) + T = load_fffile(tfile, 'T', ndeg) + B = get_weights(M, T) + # degrade mesh + print('degrado') + outft = tfile.split('/')[-1] + outfu = ufile.split('/')[-1] + outfv = vfile.split('/')[-1] + outfw = wfile.split('/')[-1] + degrade(T, 'T', B, C, outdir, outft, ndeg) + degrade(U, 'U', B, C, outdir, outfu, ndeg) + degrade(V, 'V', B, C, outdir, outfv, ndeg) + degrade(W, 'W', B, C, outdir, outfw, ndeg) + + diff --git a/FORCINGS/degradation/degrade_mesh.py b/FORCINGS/degradation/degrade_mesh.py new file mode 100644 index 0000000..c4634b7 --- /dev/null +++ b/FORCINGS/degradation/degrade_mesh.py @@ -0,0 +1,467 @@ +import numpy as np +import xarray as xr +from argparse import ArgumentParser +import yaml + +''' +generates a reduced horizontal resolution mesh from an original mesh (with Atlantic buffer) +vertical resolution is unchanged. +resolution is degraded by merging cells in (ndeg x ndeg) blocks. +functining: +0. load original mesh, +1. expand lat, lon to make the new size a multiple of ndeg, extrapolate values properly +2. compute cell surface areas and volume (for weighted means) +3. degrade resolution of original variables +4. dump new meshmask.nc +usage: python degrade_mesh.py -y degrade_mesh.yaml + +yaml file example: +maskfile: '/path/to/meshmask_in.nc' +outfile: '/path/to/meshmask_out.nc' +outfile_med: '/path/to/meshmask_out_med.nc' +ndeg: 6 #nr of cells to join in i, j +thresh: 1 #water point per block to assign to waterpoint in destination grid + +''' + +def argument(): + parser = ArgumentParser() + parser.add_argument('--yamlfile','-y', + type=str, + required=True, + help='file with all parameters') + return parser.parse_args() + +def load_parameters(): + yamlfile=argument().yamlfile + # yamlfile = 'degrade_mesh.yaml' + with open(yamlfile) as f: + Params = yaml.load(f, Loader=yaml.Loader) + return(Params) + + +# LOAD MESHMASK + +def adjust_dims(X, n=4): + ''' + add singleton dimensions in order to have it 4D + (t, z, y, x) + ''' + cdim = 0 + nd = X.ndim + while X.ndim < n: + new_dim = f'j{cdim}' + X = X.expand_dims({new_dim: 1}) + cdim = cdim + 1 + return X, nd +# +def deadjust_dims(X, nd): + ''' + remove some singleton dimensions + in order to get back to original dimensions + ''' + while X.ndim > nd and X.shape[0] == 1: + X = X.squeeze(axis=0) + return X + +def fill_rim_lin_2D(X): + ''' + fills an outer rim of NaN by linear extrapolation of internal values + in x and y direction (works also with constant values) + for glamt, gphit, etc. + ''' + t, w, j, i = X.shape + # detect rim + X2d = X.squeeze().values + isnan = np.isnan(X2d) + south_rim = np.argmax((~isnan).any(axis=1)) # first valid row + north_rim = np.argmax((~isnan[::-1]).any(axis=1)) # last valid row from top + north_rim = j - north_rim # last valid row + west_rim = 0 # first valid column (Atl open boundary) + east_rim = np.argmin((~isnan[::-1]).any(axis=0)) # last valid colum + # linear interpolation + filled = X2d.copy() + dy = filled[north_rim-1,:] - filled[north_rim-2,:] + # fill south + filled[:south_rim, :] = filled[south_rim,:] - np.outer(np.arange(south_rim)[::-1]+1, dy) + # fill north + filled[north_rim:, :] = filled[north_rim-1,:] + np.outer(np.arange(j-north_rim)+1, dy) + # fill west + dx = filled[:,1] - filled[:,0] + filled[:, :west_rim] = filled[:, west_rim][:,None] - np.outer(dx, np.arange(west_rim)[::-1]+1) + # fill east + filled[:, east_rim:] = filled[:, east_rim-1][:,None] + np.outer(dx, np.arange(i-east_rim)) + # out + X[...,:,:] = filled + return X + +def expand_array(X, pad_op, ndeg=1): + ''' + pad array with linear extrapolation (pad_op='interp') + or edge vaslues (pad_op='edge') + to make j, i multiples of ndeg. + also add a rim of size ndeg, so I'm sure to have + land points at the N, S, E boundaries + when degrading resolution + ''' + t, w, j, i = X.shape + dj = (ndeg - j % ndeg) % ndeg + di = (ndeg - i % ndeg) % ndeg + pw = {'y':(dj+ndeg, ndeg), 'x':(0, di+ndeg)} + if pad_op=='edge': + X = X.pad(pad_width=pw, mode='edge') + elif pad_op=='interp': + X = X.pad(pad_width=pw, mode='constant', constant_values=np.nan) + X = fill_rim_lin_2D(X) + return X + +def xpnd_wrap(X, pad_op, ndeg=1): + X, nd = adjust_dims(X) + X = expand_array(X, pad_op, ndeg) + X = deadjust_dims(X, nd) + return X + +def load_mesh(maskfile, ndeg=1): + ''' + load meshmask, expand fields + ''' + M = xr.open_dataset(maskfile) + M1 = {} + # + M1["coastp"] = xpnd_wrap(M["coastp"], 'edge', ndeg) + M1["e1f"] = xpnd_wrap(M["e1f"], 'interp', ndeg) + M1["e1t"] = xpnd_wrap(M["e1t"], 'interp', ndeg) + M1["e1u"] = xpnd_wrap(M["e1u"], 'interp', ndeg) + M1["e1v"] = xpnd_wrap(M["e1v"], 'interp', ndeg) + M1["e2f"] = xpnd_wrap(M["e2f"], 'edge', ndeg) + M1["e2t"] = xpnd_wrap(M["e2t"], 'edge', ndeg) + M1["e2u"] = xpnd_wrap(M["e2u"], 'edge', ndeg) + M1["e2v"] = xpnd_wrap(M["e2v"], 'edge', ndeg) + M1["e3t_0"] = xpnd_wrap(M["e3t_0"], 'edge', ndeg) + M1["e3u_0"] = xpnd_wrap(M["e3u_0"], 'edge', ndeg) + M1["e3v_0"] = xpnd_wrap(M["e3v_0"], 'edge', ndeg) + M1["e3w_0"] = xpnd_wrap(M["e3w_0"], 'edge', ndeg) + M1["ff"] = xpnd_wrap(M["ff"], 'interp', ndeg) + M1["gdept"] = M["gdept"] + M1["gdepw"] = M["gdepw"] + M1["glamf"] = xpnd_wrap(M["glamf"], 'interp', ndeg) + M1["glamt"] = xpnd_wrap(M["glamt"], 'interp', ndeg) + M1["glamu"] = xpnd_wrap(M["glamu"], 'interp', ndeg) + M1["glamv"] = xpnd_wrap(M["glamv"], 'interp', ndeg) + M1["gphif"] = xpnd_wrap(M["gphif"], 'interp', ndeg) + M1["gphit"] = xpnd_wrap(M["gphit"], 'interp', ndeg) + M1["gphiu"] = xpnd_wrap(M["gphiu"], 'interp', ndeg) + M1["gphiv"] = xpnd_wrap(M["gphiv"], 'interp', ndeg) + M1["nav_lat"] = xpnd_wrap(M["nav_lat"], 'interp', ndeg) + M1["nav_lon"] = xpnd_wrap(M["nav_lon"], 'interp', ndeg) + M1["nav_lev"] = M["nav_lev"] + M1["fmask"] = xpnd_wrap(M["fmask"], 'edge', ndeg) #there is some error here! + M1["tmask"] = xpnd_wrap(M["tmask"], 'edge', ndeg) #there is some error here! + M1["umask"] = xpnd_wrap(M["umask"], 'edge', ndeg) #there is some error here! + M1["vmask"] = xpnd_wrap(M["vmask"], 'edge', ndeg) #there is some error here! + # Compute cell surface areas and volume + M1['At'] = xr.DataArray(data = M1['e1t'].values[:] * M1['e2t'].values[:], + dims = ('time','z_a','y','x'), + name = 'At') + M1['Au'] = xr.DataArray(data = M1['e2u'].values[:] * M1['e3u_0'].values[:], + dims = ('time','z','y','x'), + name = 'Au') + M1['Av'] = xr.DataArray(data = M1['e2v'].values[:] * M1['e3v_0'].values[:], + dims = ('time','z','y','x'), + name = 'Av') + M1['V'] = xr.DataArray(data = M1['e1v'].values[:] * M1['e2u'].values[:] * M1['e3t_0'].values[:], + dims = ('time','z','y','x'), + name = 'V') + # Bonus variable for degrading forcings + M1['h_column_t'] = (M1['tmask'] * M1['e3t_0']).sum(dim='z') + # + M.close() + M1 = xr.Dataset(M1) + return M1 + +# DEGRADE RESOLUTION + +def reshape_blocks(X, ndeg=1): + ''' + reshape in blocks of shape ndeg along j, i + xarray's convoluted way of doing np.reshape, + assume I'm always reshaping dims called y, x + chunking is necessary to avoid OOM + (still cuz xarray's way of reshaping is not super smart) + ''' + cs = ndeg * 5 + X = X.chunk({'y':cs, 'x':cs}) + X = X.coarsen(y=ndeg, x=ndeg).construct(y=('y', 'y_b'), x=('x', 'x_b')) + return X + +# REGRIDDING FUNCTIONS +def isum_jmean(X, ndeg=1, W=1.0): + ''' + sum along longitude (i) & mean along latitude (j) + no need to weight by area! for e1t + ''' + X = X.sum(dim='x_b', skipna=True) #blocks sum along lon (i) + X = X.mean(dim='y_b', skipna=True) #block mean along lat (j) + return X + +def imean_jstep(X, ndeg=1, W=1.0): + ''' + mean along i (lon) + & one element each ndeg along j (lat) + glamt + ''' + X = X.isel({'y_b':-1}) + X = X.mean(dim='x_b', skipna=True) + return X + +def jmean_istep(X, ndeg=1, W=1.0): + ''' + mean along j (lat) + & one element each ndeg along i (lon) + gphit + ''' + X = X.isel({'x_b':-1}) + X = X.mean(dim='y_b', skipna=True) + return X + +def isum_jmean(X, ndeg=1, W=1.0): + ''' + sum along longitude (i) & mean along latitude (j) + no need to weight by area! + e1t + ''' + X = X.sum(dim='x_b', skipna=True) #blocks sum along lon (i) + X = X.mean(dim='y_b', skipna=True) #block mean along lat (j) + return X + +def jsum_imean(X, ndeg=1, W=1.0): + ''' + sum along latitude (j) & mean along longitude (i) + no need to weight by area! + e2t + ''' + X = X.sum(dim='y_b', skipna=True) #blocks sum along lat (j) + X = X.mean(dim='x_b', skipna=True) #block mean along lon (i) + return X + +def isum_jstep(X, ndeg=1, W=1.0): + ''' + sum along longitude (i) & one point each ndeg along latitude (j) + ''' + X = X.sum(dim='x_b', skipna=True) + X = X.isel({'y_b': -1}) + return X + +def jsum_istep(X, ndeg=1, W=1.0): + ''' + sum along latitude (j) & one point each ndeg along longitude (i) + ''' + X = X.sum(dim='y_b', skipna=True) + X = X.isel({'x_b': -1}) + return X + +def awmean(X, ndeg=1, W=1.0): + ''' + mean, weighted by surface area (t-grid) + ''' + At = reshape_blocks(W, ndeg) + X = (At * X).sum(dim=('y_b', 'x_b'), skipna=True) / At.sum(dim=('y_b', 'x_b'), skipna=True) + return X + +def vwmean(X, ndeg=1, W=1.0): + ''' + mean, weighted by surface area (t-grid) + ''' + V = reshape_blocks(W, ndeg) + X = (V * X).sum(dim=('y_b', 'x_b'), skipna=True) / V.sum(dim=('y_b', 'x_b'), skipna=True) + return X + +def uawmean_istep(X, ndeg=1, W=1.0): + ''' + mean, weighted by lateral u-grid area along lat (j) + & one element each ndeg along lon (i) + ''' + Au = reshape_blocks(W, ndeg) + Au = Au.isel({'x_b':-1}) + X = X.isel({'x_b':-1}) + X = (Au * X).sum(dim='y_b', skipna=True) / Au.sum(dim='y_b', skipna=True) + return X + +def vawmean_jstep(X, ndeg=1, W=1.0): + ''' + mean, weighted by lateral u-grid area along lat (j) + & one element each ndeg along lon (i) + ''' + Av = reshape_blocks(W, ndeg) + Av = Av.isel({'y_b':-1}) + X = X.isel({'y_b':-1}) + X = (Av * X).sum(dim='x_b', skipna=True) / Av.sum(dim='x_b', skipna=True) + return X + +def e1vwmean_jstep(X, ndeg=1, W=1.0): + ''' + mean, weighted by e1v, v-grid length along lon (i) + & one element each ndeg along lat (j) + sometauy (from phy forcings) + ''' + e1v = reshape_blocks(W, ndeg) + e1v = e1v.isel({'y_b':-1}) + X = X.isel({'y_b':-1}) + X = (e1v * X).sum(dim='x_b', skipna=True) / e1v.sum(dim='x_b', skipna=True) + return X + +def e2uwmean_istep(X, ndeg=1, W=1.0): + ''' + mean, weighted by e1u, u-grid length along lat (j) + & one element each ndeg along lon (i) + sozotaux (from phy forcings) + ''' + e2u = reshape_blocks(W, ndeg) + e2u = e2u.isel({'x_b':-1}) + X = X.isel({'x_b':-1}) + X = (e2u * X).sum(dim='y_b', skipna=True) / e2u.sum(dim='y_b', skipna=True) + return X + +def waterpt_thresh(X, thresh=1 ,ndeg=1): + ''' + 1 or 0 if enough water points (or volume?) + the higher thresh, the least water points needed + to classify a cell as water + thresh=ndeg**2: all 36 points must be water + thresh=1 (36) just one point needed + tmask + ''' + X = reshape_blocks(X, ndeg) + thr = ndeg**2 - thresh + 1 + X = X.sum(dim=('x_b', 'y_b')) + X = (X >= (ndeg**2 // thr)).astype(int) + return X + +def from_tmask(tmask): + ''' + computes umask, vmask, fmask from tmask + after Madec 2015 NEMO Manual, pp. 67 + NB Madec's fmask definition is not correct! + tmask.shape: (1, 141, 380, 1307) + ''' + # CHECK I'M PADDING AT THE RIGHT END OF THE ARRAY!!! + tmask_ipad = tmask.shift(x=1, fill_value=0.0) + tmask_jpad = tmask.shift(y=1, fill_value=0.0) + tmask_ijpad = tmask.shift(x=1, y=1, fill_value=0.0) + # + umask = tmask * tmask_ipad + vmask = tmask * tmask_jpad + # ERROR: THIS fmask DOES NOT CORRESPOND WITH WHAT'S IN meshmask.nc + fmask = tmask * tmask_ipad * tmask_jpad * tmask_ijpad + # + return umask, vmask, fmask + +def noop(X): + ''' + no operation, for variables that need no regridding + e.g. nav_lev + ''' + return X + +# END REGRIDDING OPERATIONS + +def degr_wrap(X, degr_op, ndeg=1, W=1.0): + X, nd = adjust_dims(X) + X = reshape_blocks(X, ndeg) + X = degr_op(X, ndeg, W) + X = deadjust_dims(X, nd) + return X + +def degrade_mesh(M1, thresh=1, ndeg=1): + ''' + generates new reduced mesh for each variable + ''' + M2 = {} + # + M2["coastp"] = degr_wrap(M1["coastp"], awmean, ndeg, W=M1['At']) + M2["e1f"] = degr_wrap(M1["e1f"], isum_jstep, ndeg, W=None) + M2["e1t"] = degr_wrap(M1["e1t"], isum_jmean, ndeg, W=None) + M2["e1u"] = degr_wrap(M1["e1u"], isum_jmean, ndeg, W=None) + M2["e1v"] = degr_wrap(M1["e1v"], isum_jstep, ndeg, W=None) + M2["e2f"] = degr_wrap(M1["e2f"], jsum_istep, ndeg, W=None) + M2["e2t"] = degr_wrap(M1["e2t"], jsum_istep, ndeg, W=None) + M2["e2u"] = degr_wrap(M1["e2u"], jsum_istep, ndeg, W=None) + M2["e2v"] = degr_wrap(M1["e2v"], jsum_istep, ndeg, W=None) + M2["e3t_0"] = degr_wrap(M1["e3t_0"], vwmean, ndeg, W=M1['V']) + M2["e3u_0"] = degr_wrap(M1["e3u_0"], uawmean_istep, ndeg, W=M1['Au']) + M2["e3v_0"] = degr_wrap(M1["e3v_0"], vawmean_jstep, ndeg, W=M1['Av']) + M2["e3w_0"] = degr_wrap(M1["e3w_0"], vwmean, ndeg, W=M1['V']) + M2["ff"] = degr_wrap(M1["ff"], jmean_istep, ndeg, W=None) + M2["gdept"] = noop(M1["gdept"]) + M2["gdepw"] = noop(M1["gdepw"]) + M2["glamf"] = degr_wrap(M1["glamf"], jmean_istep, ndeg, W=None) + M2["glamt"] = degr_wrap(M1["glamt"], imean_jstep, ndeg, W=None) + M2["glamu"] = degr_wrap(M1["glamu"], jmean_istep, ndeg, W=None) + M2["glamv"] = degr_wrap(M1["glamv"], imean_jstep, ndeg, W=None) + M2["gphif"] = degr_wrap(M1["gphif"], imean_jstep, ndeg, W=None) + M2["gphit"] = degr_wrap(M1["gphit"], jmean_istep, ndeg, W=None) + M2["gphiu"] = degr_wrap(M1["gphiu"], jmean_istep, ndeg, W=None) + M2["gphiv"] = degr_wrap(M1["gphiv"], imean_jstep, ndeg, W=None) + M2["nav_lat"] = degr_wrap(M1["nav_lat"], jmean_istep, ndeg, W=None) + M2["nav_lon"] = degr_wrap(M1["nav_lon"], imean_jstep, ndeg, W=None) + M2["nav_lev"] = noop(M1["nav_lev"]) + M2["tmask"] = waterpt_thresh(M1["tmask"], thresh ,ndeg) + umask, vmask, fmask = from_tmask(M2["tmask"]) + M2["umask"] = umask + M2["vmask"] = vmask + M2["fmask"] = fmask + # + M2 = xr.Dataset(M2) + return M2 + +def cut_med(M2, lon_cut=-8.875, depth_cut=4153.0, biscay_land=True): + ''' + cut out atlantic buffer, + mask bay of biscay + ''' + lon1D = M2['glamt'].isel(time=0, z_a=0, y=0) + z1D = M2['nav_lev'] + # + x_sel = np.arange(len(lon1D))[lon1D >= lon_cut] + z_sel = np.arange(len(z1D))[z1D <= depth_cut] + # + MMed = M2.isel(z=z_sel, x=x_sel) + if biscay_land: + bool1 = (MMed.glamt < 0) & (MMed.gphit > 42) + bool2 = (MMed.glamt < -6.) & (MMed.gphit > 37.25) + land = (bool1 | bool2).squeeze() + MMed['tmask'] = MMed['tmask'] * ~land + MMed['umask'] = MMed['umask'] * ~land + MMed['vmask'] = MMed['vmask'] * ~land + MMed['fmask'] = MMed['fmask'] * ~land + return MMed + +def dump_mesh(M2, outfile): + print(f'saving: {outfile}') + M2.to_netcdf(outfile, mode='w', unlimited_dims='time', format='NETCDF4_CLASSIC') + return + +if __name__=='__main__': + Params = load_parameters() + maskfile = Params['maskfile'] + outfile = Params['outfile'] + outfile_med = Params['outfile_med'] + # by how much cells to degrade (nr of cells to join in i, j) + ndeg = Params['ndeg'] + # threshold of points to assign waterpoint to degraded mesh + thresh = Params['thresh'] #THE ONLY OPTION THAT ENSURES NO RIVERS END UP ON LAND + # load meshmask and expand dimensions + print('---1---') + M1 = load_mesh(maskfile, ndeg) + # + print('---2---') + M2 = degrade_mesh(M1, thresh, ndeg) + # + print('---3---') + dump_mesh(M2, outfile) + # + print('---4---') + MMed = cut_med(M2) + dump_mesh(MMed, outfile_med) + + From efd89430eec20bd9bcf5521c61e64e3c00692a5e Mon Sep 17 00:00:00 2001 From: Giovanni Galli Date: Fri, 12 Dec 2025 11:28:25 +0100 Subject: [PATCH 03/15] IC/degrade_restart.py degrades the resolution of restart files by an integer value. --- IC/degrade_restart.py | 140 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 140 insertions(+) create mode 100644 IC/degrade_restart.py diff --git a/IC/degrade_restart.py b/IC/degrade_restart.py new file mode 100644 index 0000000..6a7dd67 --- /dev/null +++ b/IC/degrade_restart.py @@ -0,0 +1,140 @@ +import numpy as np +import xarray as xr +from argparse import ArgumentParser +import yaml +from glob import glob +# +import degrade_mesh as dm +import degrade_forcings as df + +''' +degrades the horizontal resolution of BFM restarts by an integer value +joins boxes of ndeg x ndeg cells, does volume weighted mean. +needs both original and degraded-resolution meshmask +(run degrade mesh.py before) + +usage: +python degrade_mesh.py -y degrade_mesh.yaml + +yaml file example: +mesh_in: '/path/to/meshmask_original.nc' +mesh_out: '/path/to/meshmask_degraded.nc' +indir: '/path/to/restarts/rst.nc' #also with wildcards '*' +outdir: '/path/to/output_dir/' +infile: 'RST.*-00:00:00.__VNAME__.nc' #__VNAME__ gets replaced and is needed to read the variable +ndeg: 6 #size of boxes for resolution degradation +''' + +def argument(): + parser = ArgumentParser() + parser.add_argument('--yamlfile','-y', + type=str, + required=True, + help='file with all parameters') + return parser.parse_args() + +def load_parameters(): + yamlfile=argument().yamlfile + # yamlfile = 'degrade_restart.yaml' + with open(yamlfile) as f: + Params = yaml.load(f, Loader=yaml.Loader) + return(Params) + +def get_Volume(M): + ''' + gets cell volume for weighted mean + ''' + V0 = M['e1t'].values[:] * M['e2t'].values[:] * M['e3t_0'] + nanmask = M['tmask'].values[:] + nanmask[nanmask==0.0] = np.nan + V0 = V0 * nanmask #else you see border effects at the coast + return V0 + +def load_rst(infile, vname, ndeg=1): + ''' + loads the data array from a variable in the restart + ''' + D = xr.open_dataset(infile) + D1 = {} + # + D1[vname] = dm.xpnd_wrap(D[vname], 'edge', ndeg) + # + D1['nav_lat'] = dm.xpnd_wrap(D['nav_lat'], 'interp', ndeg) + D1['nav_lon'] = dm.xpnd_wrap(D['nav_lon'], 'interp', ndeg) + #D1[vname].coords['nav_lat'].values[:] = D1['nav_lat'].values[:] + #D1[vname].coords['nav_lon'].values[:] = D1['nav_lon'].values[:] + D1['time'] = D['time'] + D1['nav_lev'] = D['nav_lev'] + # + D1 = xr.Dataset(D1) + D1 = D1.assign_attrs(D.attrs) + D.close() + return D1 + +def init_rst(C): + ''' + lat, lon are the same for all files + ''' + R = {} + R['nav_lon'] = C['glamt'] + R['nav_lat'] = C['gphit'] + return R + +def degrade_bgc(DI, V0, C, vname, ndeg=1): + ''' + degrades resolution of BFM restart file + ''' + R = init_rst(C) + R[vname] = df.degrade_wrap(DI[vname], dm.vwmean, V0, ndeg) + R['nav_lev'] = DI['nav_lev'] + R = xr.Dataset(R) + R = R.assign_attrs(DI.attrs) + return R + +def get_flist(Params): + indir = Params['indir'] + infile = Params['infile'] + flist = glob(indir+infile.replace('__VNAME__', '*')) + itrbl = [(ff.split('.')[-2], ff) for ff in flist] + return itrbl + +if __name__=='__main__': + try: + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nranks = comm.size + except: + comm = None + rank = 0 + nranks = 1 + # + Params = load_parameters() + mesh_in = Params['mesh_in'] + mesh_out = Params['mesh_out'] + indir = Params['indir'] + infile = Params['infile'] + outdir = Params['outdir'] + ndeg = Params['ndeg'] + # + if rank==0: + M = df.load_mesh_light(mesh_in, ndeg) + C = df.load_coords_degraded(mesh_out) + V0 = get_Volume(M) + else: + M = None + C = None + V0 = None + if nranks > 1: + M = comm.bcast(M, root=0) + C = comm.bcast(C, root=0) + V0 = comm.bcast(V0, root=0) + # HERE LOOP ON FILES AND VARIABLES + itrbl = get_flist(Params) + for vname, fname in itrbl[rank::nranks]: + vname = 'TRN'+vname + print(f'degrado: {vname}') + # + DI = load_rst(fname, vname, ndeg) + Dd = degrade_bgc(DI, V0, C, vname, ndeg) + outfile = outdir+fname.split('/')[-1] + dm.dump_mesh(Dd, outfile) From b2ed7af77083a11723f8fb82366247ee93580ad5 Mon Sep 17 00:00:00 2001 From: Giovanni Galli Date: Tue, 16 Dec 2025 16:51:51 +0100 Subject: [PATCH 04/15] degrade_mesh.py reconstruct gdept, gdepw from e3t_0, e3w_0 --- FORCINGS/degradation/degrade_mesh.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/FORCINGS/degradation/degrade_mesh.py b/FORCINGS/degradation/degrade_mesh.py index c4634b7..8aecbb5 100644 --- a/FORCINGS/degradation/degrade_mesh.py +++ b/FORCINGS/degradation/degrade_mesh.py @@ -121,6 +121,15 @@ def xpnd_wrap(X, pad_op, ndeg=1): X = deadjust_dims(X, nd) return X +def e3_2_gdep(M, tuvw): + ''' + reconstructs gdept from e3t_0 + cause no gdept in the original Nemo mesh + ''' + e3 = M[f'e3{tuvw}_0'] + gdep = e3.cumsum(dim="z") - (e3 / 2.0) + return gdep + def load_mesh(maskfile, ndeg=1): ''' load meshmask, expand fields @@ -142,8 +151,8 @@ def load_mesh(maskfile, ndeg=1): M1["e3v_0"] = xpnd_wrap(M["e3v_0"], 'edge', ndeg) M1["e3w_0"] = xpnd_wrap(M["e3w_0"], 'edge', ndeg) M1["ff"] = xpnd_wrap(M["ff"], 'interp', ndeg) - M1["gdept"] = M["gdept"] - M1["gdepw"] = M["gdepw"] + M1["gdept"] = xpnd_wrap(e3_2_gdep(M, 't'), 'edge', ndeg) #M["gdept"] + M1["gdepw"] = xpnd_wrap(e3_2_gdep(M, 'w'), 'edge', ndeg) #M["gdepw"] M1["glamf"] = xpnd_wrap(M["glamf"], 'interp', ndeg) M1["glamt"] = xpnd_wrap(M["glamt"], 'interp', ndeg) M1["glamu"] = xpnd_wrap(M["glamu"], 'interp', ndeg) @@ -392,8 +401,10 @@ def degrade_mesh(M1, thresh=1, ndeg=1): M2["e3v_0"] = degr_wrap(M1["e3v_0"], vawmean_jstep, ndeg, W=M1['Av']) M2["e3w_0"] = degr_wrap(M1["e3w_0"], vwmean, ndeg, W=M1['V']) M2["ff"] = degr_wrap(M1["ff"], jmean_istep, ndeg, W=None) - M2["gdept"] = noop(M1["gdept"]) - M2["gdepw"] = noop(M1["gdepw"]) + #M2["gdept"] = noop(M1["gdept"]) + #M2["gdepw"] = noop(M1["gdepw"]) + M2["gdept"] = degr_wrap(M1["gdept"], vwmean, ndeg, W=M1['V']) + M2["gdepw"] = degr_wrap(M1["gdepw"], vwmean, ndeg, W=M1['V']) M2["glamf"] = degr_wrap(M1["glamf"], jmean_istep, ndeg, W=None) M2["glamt"] = degr_wrap(M1["glamt"], imean_jstep, ndeg, W=None) M2["glamu"] = degr_wrap(M1["glamu"], jmean_istep, ndeg, W=None) From 82ae5e94d67bfd1bf22c83c08b66aac3e263e583 Mon Sep 17 00:00:00 2001 From: Giovanni Galli Date: Wed, 17 Dec 2025 14:30:31 +0100 Subject: [PATCH 05/15] added somxl010 to T-grid forcing files --- FORCINGS/degradation/degrade_forcings.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/FORCINGS/degradation/degrade_forcings.py b/FORCINGS/degradation/degrade_forcings.py index 5eb0fce..0bdd3a0 100644 --- a/FORCINGS/degradation/degrade_forcings.py +++ b/FORCINGS/degradation/degrade_forcings.py @@ -118,7 +118,7 @@ def load_tfile(infile, ndeg=1): #F1['solofldo'] = dm.xpnd_wrap(F['solofldo'], 'edge', ndeg) #F1['sosefldo'] = dm.xpnd_wrap(F['sosefldo'], 'edge', ndeg) #F1['solafldo'] = dm.xpnd_wrap(F['solafldo'], 'edge', ndeg) - #F1['somxl010'] = dm.xpnd_wrap(F['somxl010'], 'edge', ndeg) + F1['somxl010'] = dm.xpnd_wrap(F['somxl010'], 'edge', ndeg) # overwrite coordinates, else it complains for vv in F1.keys(): F1[vv].coords['nav_lat'].values[:] = F1['nav_lat'].values[:] @@ -401,6 +401,7 @@ def degrade_T(T, B, C, ndeg=1): Td['sossheig'] = degrade_wrap(T['sossheig'], dm.awmean, B['At'], ndeg) Td['soshfldo'] = degrade_wrap(T['soshfldo'], dm.awmean, B['At'], ndeg) Td['sorunoff'] = degrade_wrap(T['sorunoff'], dm.awmean, B['At'], ndeg) + Td['somxl010'] = degrade_wrap(T['somxl010'], dm.awmean, B['At'], ndeg) # Td = xr.Dataset(Td) return Td From 28045888e06925d466874a166dea8f397c8b8f1e Mon Sep 17 00:00:00 2001 From: Giovanni Galli Date: Thu, 18 Dec 2025 10:51:49 +0100 Subject: [PATCH 06/15] degrade_forcings.py saves files under /YYYY/MM/[TUVW]YYYYMM*.nc --- FORCINGS/degradation/degrade_forcings.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/FORCINGS/degradation/degrade_forcings.py b/FORCINGS/degradation/degrade_forcings.py index 0bdd3a0..15c7e40 100644 --- a/FORCINGS/degradation/degrade_forcings.py +++ b/FORCINGS/degradation/degrade_forcings.py @@ -1,6 +1,7 @@ import numpy as np import xarray as xr from glob import glob +import os from mpi4py import MPI from argparse import ArgumentParser import yaml @@ -406,6 +407,18 @@ def degrade_T(T, B, C, ndeg=1): Td = xr.Dataset(Td) return Td +def make_outdir(outdir, outfile): + #/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/qDEG_SETUP/FORCINGS/ + #T20020308-12:00:00.nc + yyyy = outfile[1:5] + mm = outfile[5:7] + outdir = f'{outdir}/{yyyy}/{mm}/' + if not os.path.exists(outdir): + os.makedirs(outdir) + else: + pass + return outdir + def degrade(F, tuv, B, C, outdir, outfile, ndeg=1): if tuv == 'T': Fd = degrade_T(F, B, C, ndeg) @@ -415,6 +428,7 @@ def degrade(F, tuv, B, C, outdir, outfile, ndeg=1): Fd = degrade_V(F, B, C, ndeg) elif tuv == 'W': Fd = degrade_W(F, B, C, ndeg) + outdir = make_outdir(outdir, outfile) outfile = outdir+outfile dm.dump_mesh(Fd, outfile) return Fd From 12aa1072e8971f160f4553c5367bf0d61d0f156f Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Thu, 18 Dec 2025 18:33:52 +0100 Subject: [PATCH 07/15] improving code readability --- FORCINGS/degradation/degrade_mesh.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/FORCINGS/degradation/degrade_mesh.py b/FORCINGS/degradation/degrade_mesh.py index 8aecbb5..268db0b 100644 --- a/FORCINGS/degradation/degrade_mesh.py +++ b/FORCINGS/degradation/degrade_mesh.py @@ -201,6 +201,8 @@ def reshape_blocks(X, ndeg=1): cs = ndeg * 5 X = X.chunk({'y':cs, 'x':cs}) X = X.coarsen(y=ndeg, x=ndeg).construct(y=('y', 'y_b'), x=('x', 'x_b')) + + # X has the dimension y_b, x_b of size ndeg return X # REGRIDDING FUNCTIONS @@ -279,6 +281,7 @@ def awmean(X, ndeg=1, W=1.0): def vwmean(X, ndeg=1, W=1.0): ''' + W is a dataarray with cell volumes 3D mean, weighted by surface area (t-grid) ''' V = reshape_blocks(W, ndeg) @@ -331,8 +334,10 @@ def e2uwmean_istep(X, ndeg=1, W=1.0): X = (e2u * X).sum(dim='y_b', skipna=True) / e2u.sum(dim='y_b', skipna=True) return X -def waterpt_thresh(X, thresh=1 ,ndeg=1): +def waterpt_thresh(X24, thresh=1 ,ndeg=1): ''' + Args: + X: tmask 1 or 0 if enough water points (or volume?) the higher thresh, the least water points needed to classify a cell as water @@ -340,11 +345,11 @@ def waterpt_thresh(X, thresh=1 ,ndeg=1): thresh=1 (36) just one point needed tmask ''' - X = reshape_blocks(X, ndeg) + X24 = reshape_blocks(X24, ndeg) thr = ndeg**2 - thresh + 1 - X = X.sum(dim=('x_b', 'y_b')) - X = (X >= (ndeg**2 // thr)).astype(int) - return X + X_degr = X24.sum(dim=('x_b', 'y_b')) + X_degr= (X_degr >= (ndeg**2 // thr)).astype(int) + return X_degr def from_tmask(tmask): ''' @@ -374,12 +379,12 @@ def noop(X): # END REGRIDDING OPERATIONS -def degr_wrap(X, degr_op, ndeg=1, W=1.0): - X, nd = adjust_dims(X) - X = reshape_blocks(X, ndeg) - X = degr_op(X, ndeg, W) - X = deadjust_dims(X, nd) - return X +def degr_wrap(X24, degr_op, ndeg=1, W=1.0): + X24, nd = adjust_dims(X24) # make 4D + X24_dims = reshape_blocks(X24, ndeg) + X_degr = degr_op(X24_dims, ndeg, W) + X_degr = deadjust_dims(X_degr, nd) # back to original dims + return X_degr def degrade_mesh(M1, thresh=1, ndeg=1): ''' From 26c6a7ea1d00082cee3dab0ab9362a707eabd3c4 Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Fri, 19 Dec 2025 18:26:25 +0100 Subject: [PATCH 08/15] Adding mesh_degradation.py as functional main, in order to move there argparse, yaml and main calls, Taking degrade_mesh as library, with the __main__ section to test them Improving readability and docstrings --- FORCINGS/degradation/degrade_mesh.py | 216 +++++++++++++---------- FORCINGS/degradation/mesh_degradation.py | 63 +++++++ 2 files changed, 184 insertions(+), 95 deletions(-) create mode 100644 FORCINGS/degradation/mesh_degradation.py diff --git a/FORCINGS/degradation/degrade_mesh.py b/FORCINGS/degradation/degrade_mesh.py index 268db0b..e8de622 100644 --- a/FORCINGS/degradation/degrade_mesh.py +++ b/FORCINGS/degradation/degrade_mesh.py @@ -1,59 +1,27 @@ import numpy as np import xarray as xr -from argparse import ArgumentParser -import yaml - -''' -generates a reduced horizontal resolution mesh from an original mesh (with Atlantic buffer) -vertical resolution is unchanged. -resolution is degraded by merging cells in (ndeg x ndeg) blocks. -functining: -0. load original mesh, -1. expand lat, lon to make the new size a multiple of ndeg, extrapolate values properly -2. compute cell surface areas and volume (for weighted means) -3. degrade resolution of original variables -4. dump new meshmask.nc -usage: python degrade_mesh.py -y degrade_mesh.yaml - -yaml file example: -maskfile: '/path/to/meshmask_in.nc' -outfile: '/path/to/meshmask_out.nc' -outfile_med: '/path/to/meshmask_out_med.nc' -ndeg: 6 #nr of cells to join in i, j -thresh: 1 #water point per block to assign to waterpoint in destination grid - -''' - -def argument(): - parser = ArgumentParser() - parser.add_argument('--yamlfile','-y', - type=str, - required=True, - help='file with all parameters') - return parser.parse_args() - -def load_parameters(): - yamlfile=argument().yamlfile - # yamlfile = 'degrade_mesh.yaml' - with open(yamlfile) as f: - Params = yaml.load(f, Loader=yaml.Loader) - return(Params) - - -# LOAD MESHMASK - -def adjust_dims(X, n=4): - ''' - add singleton dimensions in order to have it 4D - (t, z, y, x) + + +def adjust_dims(X: xr.DataArray, n: int = 4): + ''' + add singleton dimensions in front of existing ones in order to have n dimensions + + Args: + X: xarray DataArray, with dims < = n + n: desired number of dimensions (default 4) + + Returns: + X_adjusted: 4D xarray DataArray with dims (t, z, y, x) + orig_dim: int, original number of dimensions of X + ''' cdim = 0 - nd = X.ndim + orig_dim = X.ndim while X.ndim < n: new_dim = f'j{cdim}' X = X.expand_dims({new_dim: 1}) cdim = cdim + 1 - return X, nd + return X, orig_dim # def deadjust_dims(X, nd): ''' @@ -95,31 +63,58 @@ def fill_rim_lin_2D(X): X[...,:,:] = filled return X -def expand_array(X, pad_op, ndeg=1): +def expand_array(X: xr.DataArray, pad_op:str="edge", ndeg:int=1): ''' + Args: + X: xarray DataArray with dims (t, z, y, x) + pad_op: string, 'edge' or 'interp'. + ndeg: int, number of cells to join in i, j + + Uses xarray.DataArray.pad pad array with linear extrapolation (pad_op='interp') or edge vaslues (pad_op='edge') - to make j, i multiples of ndeg. + to have the sizes y, x multiples of ndeg. + also add a rim of size ndeg, so I'm sure to have land points at the N, S, E boundaries when degrading resolution - ''' - t, w, j, i = X.shape - dj = (ndeg - j % ndeg) % ndeg - di = (ndeg - i % ndeg) % ndeg - pw = {'y':(dj+ndeg, ndeg), 'x':(0, di+ndeg)} + Returns: + X_padded: xarray DataArray with dims (t, z, y_padded, x_padded) + ''' + jpt, jpk, jpj, jpi = X.shape + # compute padding sizes, integers 0 <= d <= ndeg -1 + dj = (ndeg - jpj % ndeg) % ndeg + di = (ndeg - jpi % ndeg) % ndeg + + padding_south = dj + ndeg # applying delta + the future border of one cell + padding_north = ndeg # just the border + padding_east = 0 # nothing on Atlantic buffer + padding_west = di + ndeg # applying delta + the future border of one cell + + # generate pad_width dict + pw = {'y':(padding_south, padding_north), 'x':(padding_east, padding_west)} + if pad_op=='edge': - X = X.pad(pad_width=pw, mode='edge') + X_padded = X.pad(pad_width = pw, mode='edge') elif pad_op=='interp': - X = X.pad(pad_width=pw, mode='constant', constant_values=np.nan) - X = fill_rim_lin_2D(X) - return X + X_padded = X.pad(pad_width = pw, mode='constant', constant_values=np.nan) + X_padded = fill_rim_lin_2D(X_padded) + return X_padded -def xpnd_wrap(X, pad_op, ndeg=1): - X, nd = adjust_dims(X) - X = expand_array(X, pad_op, ndeg) - X = deadjust_dims(X, nd) - return X +def xpnd_wrap(X: xr.DataArray, pad_op:str="edge", ndeg:int=1): + ''' + Args: + X: xarray DataArray with dims (... y, x) or less + pad_op: string, 'edge' or 'interp'. + ndeg: int, number of cells to join in i, j + + Returns: + X_expanded: xarray DataArray with dims (... y_padded, x_padded) + ''' + X_4d, orig_dims = adjust_dims(X) + X_4d_expand = expand_array(X_4d, pad_op, ndeg) + X_expand = deadjust_dims(X_4d_expand, orig_dims) + return X_expand def e3_2_gdep(M, tuvw): ''' @@ -193,17 +188,27 @@ def load_mesh(maskfile, ndeg=1): def reshape_blocks(X, ndeg=1): ''' reshape in blocks of shape ndeg along j, i - xarray's convoluted way of doing np.reshape, + xarray's convoluted way of doing np.reshape, esticatsi! assume I'm always reshaping dims called y, x chunking is necessary to avoid OOM (still cuz xarray's way of reshaping is not super smart) + + Args: + X: xarray DataArray with dims (... y, x) + ndeg: int, number of cells to join in y, x directions + x and y sizes must be multiples of ndeg + Returns: + X_reshaped: xarray DataArray with dims + (..., y, y_b, x, x_b) with dims sizes + (..., y/ndeg, ndeg, x/ndeg, ndeg) + So, X_reshaped[... 0,:,0,:] is the first block of size ndeg x ndeg + ''' cs = ndeg * 5 X = X.chunk({'y':cs, 'x':cs}) - X = X.coarsen(y=ndeg, x=ndeg).construct(y=('y', 'y_b'), x=('x', 'x_b')) + X_reshaped = X.coarsen(y=ndeg, x=ndeg).construct(y=('y', 'y_b'), x=('x', 'x_b')) - # X has the dimension y_b, x_b of size ndeg - return X + return X_reshaped # REGRIDDING FUNCTIONS def isum_jmean(X, ndeg=1, W=1.0): @@ -379,10 +384,37 @@ def noop(X): # END REGRIDDING OPERATIONS -def degr_wrap(X24, degr_op, ndeg=1, W=1.0): - X24, nd = adjust_dims(X24) # make 4D - X24_dims = reshape_blocks(X24, ndeg) - X_degr = degr_op(X24_dims, ndeg, W) +def degr_wrap(X24, degr_op:callable, ndeg=1, W=1.0): + ''' + Wrapper to degrade resolution of 2D/3D/4D arrays. + + Parameters + ---------- + X24 : xarray.DataArray + Input array with dimensions (..., y, x). + degr_op : Callable + Regridding function that performs the degradation operation. + Expected signature: degr_op(X24_dims, ndeg, W) -> xarray.DataArray + ndeg : int, optional + Number of cells to join in i, j directions. Default is 1. + W : xarray.DataArray, optional + Weights for weighted means computation. Default is 1.0. + + Returns + ------- + xarray.DataArray + Degraded array with the same dimensions as the input X24. + + Notes + ----- + The function internally transforms the input to 4D format, applies + the degradation operation on reshaped blocks, and then restores + the original dimensionality. + + ''' + X24_4D, nd = adjust_dims(X24) # make 4D + X24_6D = reshape_blocks(X24_4D, ndeg) + X_degr = degr_op(X24_6D, ndeg, W) X_degr = deadjust_dims(X_degr, nd) # back to original dims return X_degr @@ -458,26 +490,20 @@ def dump_mesh(M2, outfile): return if __name__=='__main__': - Params = load_parameters() - maskfile = Params['maskfile'] - outfile = Params['outfile'] - outfile_med = Params['outfile_med'] - # by how much cells to degrade (nr of cells to join in i, j) - ndeg = Params['ndeg'] - # threshold of points to assign waterpoint to degraded mesh - thresh = Params['thresh'] #THE ONLY OPTION THAT ENSURES NO RIVERS END UP ON LAND - # load meshmask and expand dimensions - print('---1---') - M1 = load_mesh(maskfile, ndeg) - # - print('---2---') - M2 = degrade_mesh(M1, thresh, ndeg) - # - print('---3---') - dump_mesh(M2, outfile) - # - print('---4---') - MMed = cut_med(M2) - dump_mesh(MMed, outfile_med) - - + data = np.array([[1, 2], + [3, 4], + [5, 6]]) # shape (3, 2) + arr= xr.DataArray(data, dims=('x','y'), coords={'x':[0,1,2], 'y':[10,20]}) + arr.pad({'y':(2,2),'x':(2,2)},mode='edge') + X, nd = adjust_dims(arr) + print(X) + + X_exp = xpnd_wrap(arr, pad_op='edge', ndeg=3) # su array piccoli fa una cosa strana + data = np.ones((20,20), dtype=np.float32) + arr= xr.DataArray(data, dims=('x','y'), coords={'x':np.arange(20), 'y':np.arange(20,40)}) + print(xpnd_wrap(arr, pad_op='edge', ndeg=4).sizes) + + X = xpnd_wrap(arr, pad_op='edge', ndeg=3) + X_4D, nd = adjust_dims(X) + X24_6D = reshape_blocks(X_4D, ndeg=3) + X_degr = jsum_istep(X24_6D, ndeg=3, W=None) diff --git a/FORCINGS/degradation/mesh_degradation.py b/FORCINGS/degradation/mesh_degradation.py new file mode 100644 index 0000000..f938647 --- /dev/null +++ b/FORCINGS/degradation/mesh_degradation.py @@ -0,0 +1,63 @@ +import numpy as np +import xarray as xr +from argparse import ArgumentParser +import yaml +from degrade_mesh import load_mesh, degrade_mesh, dump_mesh, cut_med + +''' +generates a reduced horizontal resolution mesh from an original mesh (with Atlantic buffer) +vertical resolution is unchanged. +resolution is degraded by merging cells in (ndeg x ndeg) blocks. +functining: +0. load original mesh, +1. expand lat, lon to make the new size a multiple of ndeg, extrapolate values properly +2. compute cell surface areas and volume (for weighted means) +3. degrade resolution of original variables +4. dump new meshmask.nc +usage: python degrade_mesh.py -y degrade_mesh.yaml + +yaml file example: +maskfile: '/path/to/meshmask_in.nc' +outfile: '/path/to/meshmask_out.nc' +outfile_med: '/path/to/meshmask_out_med.nc' +ndeg: 6 #nr of cells to join in i, j +thresh: 1 #water point per block to assign to waterpoint in destination grid + +''' + +def argument(): + parser = ArgumentParser() + parser.add_argument('--yamlfile','-y', + type=str, + required=True, + help='file with all parameters') + return parser.parse_args() + +def load_parameters(): + yamlfile=argument().yamlfile + # yamlfile = 'degrade_mesh.yaml' + with open(yamlfile) as f: + Params = yaml.load(f, Loader=yaml.Loader) + return(Params) + +Params = load_parameters() +maskfile = Params['maskfile'] +outfile = Params['outfile'] +outfile_med = Params['outfile_med'] +# by how much cells to degrade (nr of cells to join in i, j) +ndeg = Params['ndeg'] +# threshold of points to assign waterpoint to degraded mesh +thresh = Params['thresh'] #THE ONLY OPTION THAT ENSURES NO RIVERS END UP ON LAND +# load meshmask and expand dimensions +print('---1---') +M1 = load_mesh(maskfile, ndeg) +# +print('---2---') +M2 = degrade_mesh(M1, thresh, ndeg) +# +print('---3---') +dump_mesh(M2, outfile) +# +print('---4---') +MMed = cut_med(M2) +dump_mesh(MMed, outfile_med) \ No newline at end of file From 51a38422f8fb73c4f77d9989bea126deb234b763 Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Fri, 19 Dec 2025 18:29:29 +0100 Subject: [PATCH 09/15] adding test script --- FORCINGS/degradation/test.py | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 FORCINGS/degradation/test.py diff --git a/FORCINGS/degradation/test.py b/FORCINGS/degradation/test.py new file mode 100644 index 0000000..fb4f3d6 --- /dev/null +++ b/FORCINGS/degradation/test.py @@ -0,0 +1,7 @@ +import xarray as xr +from degrade_mesh import xpnd_wrap, degr_wrap +from degrade_mesh import jsum_istep +M = xr.open_dataset(maskfile) +ndeg=6 +a = xpnd_wrap(M["e2t"], 'edge', ndeg) # in load_mesh +e2t = degr_wrap(a, jsum_istep, ndeg, W=None) From f4e5b5aa8b842c6e4f067a3d7c24034638c2f545 Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Mon, 22 Dec 2025 11:51:33 +0100 Subject: [PATCH 10/15] Moving old scripts in previous/ dir --- FORCINGS/degradation/{ => previous}/create_meshmask_nc.py | 0 FORCINGS/degradation/{ => previous}/create_meshmask_nc_MITgcm.py | 0 FORCINGS/degradation/{ => previous}/create_meshmask_nc_nemo.py | 0 FORCINGS/degradation/{ => previous}/create_meshmask_reduced.py | 0 FORCINGS/degradation/{ => previous}/forcing_reducer.py | 0 FORCINGS/degradation/{ => previous}/forcingswriter.py | 0 FORCINGS/degradation/{ => previous}/main.py | 0 FORCINGS/degradation/{ => previous}/main_MITgcm.py | 0 FORCINGS/degradation/{ => previous}/reducer.py | 0 FORCINGS/degradation/{ => previous}/restart_interpolator.py | 0 FORCINGS/degradation/{ => previous}/river_reindexer.py | 0 11 files changed, 0 insertions(+), 0 deletions(-) rename FORCINGS/degradation/{ => previous}/create_meshmask_nc.py (100%) rename FORCINGS/degradation/{ => previous}/create_meshmask_nc_MITgcm.py (100%) rename FORCINGS/degradation/{ => previous}/create_meshmask_nc_nemo.py (100%) rename FORCINGS/degradation/{ => previous}/create_meshmask_reduced.py (100%) rename FORCINGS/degradation/{ => previous}/forcing_reducer.py (100%) rename FORCINGS/degradation/{ => previous}/forcingswriter.py (100%) rename FORCINGS/degradation/{ => previous}/main.py (100%) rename FORCINGS/degradation/{ => previous}/main_MITgcm.py (100%) rename FORCINGS/degradation/{ => previous}/reducer.py (100%) rename FORCINGS/degradation/{ => previous}/restart_interpolator.py (100%) rename FORCINGS/degradation/{ => previous}/river_reindexer.py (100%) diff --git a/FORCINGS/degradation/create_meshmask_nc.py b/FORCINGS/degradation/previous/create_meshmask_nc.py similarity index 100% rename from FORCINGS/degradation/create_meshmask_nc.py rename to FORCINGS/degradation/previous/create_meshmask_nc.py diff --git a/FORCINGS/degradation/create_meshmask_nc_MITgcm.py b/FORCINGS/degradation/previous/create_meshmask_nc_MITgcm.py similarity index 100% rename from FORCINGS/degradation/create_meshmask_nc_MITgcm.py rename to FORCINGS/degradation/previous/create_meshmask_nc_MITgcm.py diff --git a/FORCINGS/degradation/create_meshmask_nc_nemo.py b/FORCINGS/degradation/previous/create_meshmask_nc_nemo.py similarity index 100% rename from FORCINGS/degradation/create_meshmask_nc_nemo.py rename to FORCINGS/degradation/previous/create_meshmask_nc_nemo.py diff --git a/FORCINGS/degradation/create_meshmask_reduced.py b/FORCINGS/degradation/previous/create_meshmask_reduced.py similarity index 100% rename from FORCINGS/degradation/create_meshmask_reduced.py rename to FORCINGS/degradation/previous/create_meshmask_reduced.py diff --git a/FORCINGS/degradation/forcing_reducer.py b/FORCINGS/degradation/previous/forcing_reducer.py similarity index 100% rename from FORCINGS/degradation/forcing_reducer.py rename to FORCINGS/degradation/previous/forcing_reducer.py diff --git a/FORCINGS/degradation/forcingswriter.py b/FORCINGS/degradation/previous/forcingswriter.py similarity index 100% rename from FORCINGS/degradation/forcingswriter.py rename to FORCINGS/degradation/previous/forcingswriter.py diff --git a/FORCINGS/degradation/main.py b/FORCINGS/degradation/previous/main.py similarity index 100% rename from FORCINGS/degradation/main.py rename to FORCINGS/degradation/previous/main.py diff --git a/FORCINGS/degradation/main_MITgcm.py b/FORCINGS/degradation/previous/main_MITgcm.py similarity index 100% rename from FORCINGS/degradation/main_MITgcm.py rename to FORCINGS/degradation/previous/main_MITgcm.py diff --git a/FORCINGS/degradation/reducer.py b/FORCINGS/degradation/previous/reducer.py similarity index 100% rename from FORCINGS/degradation/reducer.py rename to FORCINGS/degradation/previous/reducer.py diff --git a/FORCINGS/degradation/restart_interpolator.py b/FORCINGS/degradation/previous/restart_interpolator.py similarity index 100% rename from FORCINGS/degradation/restart_interpolator.py rename to FORCINGS/degradation/previous/restart_interpolator.py diff --git a/FORCINGS/degradation/river_reindexer.py b/FORCINGS/degradation/previous/river_reindexer.py similarity index 100% rename from FORCINGS/degradation/river_reindexer.py rename to FORCINGS/degradation/previous/river_reindexer.py From 700d5580b9a22050529cf7c7f87cf9f144a2d15f Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Mon, 22 Dec 2025 12:38:51 +0100 Subject: [PATCH 11/15] moving all regridding functions (working on blocks) in regridding.py --- FORCINGS/degradation/degrade_mesh.py | 107 ++++++--------------------- FORCINGS/degradation/regridding.py | 65 ++++++++++++++++ 2 files changed, 88 insertions(+), 84 deletions(-) create mode 100644 FORCINGS/degradation/regridding.py diff --git a/FORCINGS/degradation/degrade_mesh.py b/FORCINGS/degradation/degrade_mesh.py index e8de622..09e2954 100644 --- a/FORCINGS/degradation/degrade_mesh.py +++ b/FORCINGS/degradation/degrade_mesh.py @@ -1,5 +1,6 @@ import numpy as np import xarray as xr +import regridding as rg def adjust_dims(X: xr.DataArray, n: int = 4): @@ -210,71 +211,7 @@ def reshape_blocks(X, ndeg=1): return X_reshaped -# REGRIDDING FUNCTIONS -def isum_jmean(X, ndeg=1, W=1.0): - ''' - sum along longitude (i) & mean along latitude (j) - no need to weight by area! for e1t - ''' - X = X.sum(dim='x_b', skipna=True) #blocks sum along lon (i) - X = X.mean(dim='y_b', skipna=True) #block mean along lat (j) - return X - -def imean_jstep(X, ndeg=1, W=1.0): - ''' - mean along i (lon) - & one element each ndeg along j (lat) - glamt - ''' - X = X.isel({'y_b':-1}) - X = X.mean(dim='x_b', skipna=True) - return X - -def jmean_istep(X, ndeg=1, W=1.0): - ''' - mean along j (lat) - & one element each ndeg along i (lon) - gphit - ''' - X = X.isel({'x_b':-1}) - X = X.mean(dim='y_b', skipna=True) - return X - -def isum_jmean(X, ndeg=1, W=1.0): - ''' - sum along longitude (i) & mean along latitude (j) - no need to weight by area! - e1t - ''' - X = X.sum(dim='x_b', skipna=True) #blocks sum along lon (i) - X = X.mean(dim='y_b', skipna=True) #block mean along lat (j) - return X -def jsum_imean(X, ndeg=1, W=1.0): - ''' - sum along latitude (j) & mean along longitude (i) - no need to weight by area! - e2t - ''' - X = X.sum(dim='y_b', skipna=True) #blocks sum along lat (j) - X = X.mean(dim='x_b', skipna=True) #block mean along lon (i) - return X - -def isum_jstep(X, ndeg=1, W=1.0): - ''' - sum along longitude (i) & one point each ndeg along latitude (j) - ''' - X = X.sum(dim='x_b', skipna=True) - X = X.isel({'y_b': -1}) - return X - -def jsum_istep(X, ndeg=1, W=1.0): - ''' - sum along latitude (j) & one point each ndeg along longitude (i) - ''' - X = X.sum(dim='y_b', skipna=True) - X = X.isel({'x_b': -1}) - return X def awmean(X, ndeg=1, W=1.0): ''' @@ -425,33 +362,33 @@ def degrade_mesh(M1, thresh=1, ndeg=1): M2 = {} # M2["coastp"] = degr_wrap(M1["coastp"], awmean, ndeg, W=M1['At']) - M2["e1f"] = degr_wrap(M1["e1f"], isum_jstep, ndeg, W=None) - M2["e1t"] = degr_wrap(M1["e1t"], isum_jmean, ndeg, W=None) - M2["e1u"] = degr_wrap(M1["e1u"], isum_jmean, ndeg, W=None) - M2["e1v"] = degr_wrap(M1["e1v"], isum_jstep, ndeg, W=None) - M2["e2f"] = degr_wrap(M1["e2f"], jsum_istep, ndeg, W=None) - M2["e2t"] = degr_wrap(M1["e2t"], jsum_istep, ndeg, W=None) - M2["e2u"] = degr_wrap(M1["e2u"], jsum_istep, ndeg, W=None) - M2["e2v"] = degr_wrap(M1["e2v"], jsum_istep, ndeg, W=None) + M2["e1f"] = degr_wrap(M1["e1f"], rg.isum_jstep, ndeg, W=None) + M2["e1t"] = degr_wrap(M1["e1t"], rg.isum_jmean, ndeg, W=None) + M2["e1u"] = degr_wrap(M1["e1u"], rg.isum_jmean, ndeg, W=None) + M2["e1v"] = degr_wrap(M1["e1v"], rg.isum_jstep, ndeg, W=None) + M2["e2f"] = degr_wrap(M1["e2f"], rg.jsum_istep, ndeg, W=None) + M2["e2t"] = degr_wrap(M1["e2t"], rg.jsum_istep, ndeg, W=None) + M2["e2u"] = degr_wrap(M1["e2u"], rg.jsum_istep, ndeg, W=None) + M2["e2v"] = degr_wrap(M1["e2v"], rg.jsum_istep, ndeg, W=None) M2["e3t_0"] = degr_wrap(M1["e3t_0"], vwmean, ndeg, W=M1['V']) M2["e3u_0"] = degr_wrap(M1["e3u_0"], uawmean_istep, ndeg, W=M1['Au']) M2["e3v_0"] = degr_wrap(M1["e3v_0"], vawmean_jstep, ndeg, W=M1['Av']) M2["e3w_0"] = degr_wrap(M1["e3w_0"], vwmean, ndeg, W=M1['V']) - M2["ff"] = degr_wrap(M1["ff"], jmean_istep, ndeg, W=None) + M2["ff"] = degr_wrap(M1["ff"], rg.jmean_istep, ndeg, W=None) #M2["gdept"] = noop(M1["gdept"]) #M2["gdepw"] = noop(M1["gdepw"]) M2["gdept"] = degr_wrap(M1["gdept"], vwmean, ndeg, W=M1['V']) M2["gdepw"] = degr_wrap(M1["gdepw"], vwmean, ndeg, W=M1['V']) - M2["glamf"] = degr_wrap(M1["glamf"], jmean_istep, ndeg, W=None) - M2["glamt"] = degr_wrap(M1["glamt"], imean_jstep, ndeg, W=None) - M2["glamu"] = degr_wrap(M1["glamu"], jmean_istep, ndeg, W=None) - M2["glamv"] = degr_wrap(M1["glamv"], imean_jstep, ndeg, W=None) - M2["gphif"] = degr_wrap(M1["gphif"], imean_jstep, ndeg, W=None) - M2["gphit"] = degr_wrap(M1["gphit"], jmean_istep, ndeg, W=None) - M2["gphiu"] = degr_wrap(M1["gphiu"], jmean_istep, ndeg, W=None) - M2["gphiv"] = degr_wrap(M1["gphiv"], imean_jstep, ndeg, W=None) - M2["nav_lat"] = degr_wrap(M1["nav_lat"], jmean_istep, ndeg, W=None) - M2["nav_lon"] = degr_wrap(M1["nav_lon"], imean_jstep, ndeg, W=None) + M2["glamf"] = degr_wrap(M1["glamf"], rg.jmean_istep, ndeg, W=None) + M2["glamt"] = degr_wrap(M1["glamt"], rg.imean_jstep, ndeg, W=None) + M2["glamu"] = degr_wrap(M1["glamu"], rg.jmean_istep, ndeg, W=None) + M2["glamv"] = degr_wrap(M1["glamv"], rg.imean_jstep, ndeg, W=None) + M2["gphif"] = degr_wrap(M1["gphif"], rg.imean_jstep, ndeg, W=None) + M2["gphit"] = degr_wrap(M1["gphit"], rg.jmean_istep, ndeg, W=None) + M2["gphiu"] = degr_wrap(M1["gphiu"], rg.jmean_istep, ndeg, W=None) + M2["gphiv"] = degr_wrap(M1["gphiv"], rg.imean_jstep, ndeg, W=None) + M2["nav_lat"] = degr_wrap(M1["nav_lat"], rg.jmean_istep, ndeg, W=None) + M2["nav_lon"] = degr_wrap(M1["nav_lon"], rg.imean_jstep, ndeg, W=None) M2["nav_lev"] = noop(M1["nav_lev"]) M2["tmask"] = waterpt_thresh(M1["tmask"], thresh ,ndeg) umask, vmask, fmask = from_tmask(M2["tmask"]) @@ -506,4 +443,6 @@ def dump_mesh(M2, outfile): X = xpnd_wrap(arr, pad_op='edge', ndeg=3) X_4D, nd = adjust_dims(X) X24_6D = reshape_blocks(X_4D, ndeg=3) - X_degr = jsum_istep(X24_6D, ndeg=3, W=None) + X_degr = rg.jsum_istep(X24_6D, ndeg=3, W=None) + + diff --git a/FORCINGS/degradation/regridding.py b/FORCINGS/degradation/regridding.py new file mode 100644 index 0000000..2542607 --- /dev/null +++ b/FORCINGS/degradation/regridding.py @@ -0,0 +1,65 @@ +# REGRIDDING FUNCTIONS +def isum_jmean(X, ndeg=1, W=1.0): + ''' + sum along longitude (i) & mean along latitude (j) + no need to weight by area! for e1t + ''' + X = X.sum(dim='x_b', skipna=True) #blocks sum along lon (i) + X = X.mean(dim='y_b', skipna=True) #block mean along lat (j) + return X + +def imean_jstep(X, ndeg=1, W=1.0): + ''' + mean along i (lon) + & one element each ndeg along j (lat) + glamt + ''' + X = X.isel({'y_b':-1}) + X = X.mean(dim='x_b', skipna=True) + return X + +def jmean_istep(X, ndeg=1, W=1.0): + ''' + mean along j (lat) + & one element each ndeg along i (lon) + gphit + ''' + X = X.isel({'x_b':-1}) + X = X.mean(dim='y_b', skipna=True) + return X + +def isum_jmean(X, ndeg=1, W=1.0): + ''' + sum along longitude (i) & mean along latitude (j) + no need to weight by area! + e1t + ''' + X = X.sum(dim='x_b', skipna=True) #blocks sum along lon (i) + X = X.mean(dim='y_b', skipna=True) #block mean along lat (j) + return X + +def jsum_imean(X, ndeg=1, W=1.0): + ''' + sum along latitude (j) & mean along longitude (i) + no need to weight by area! + e2t + ''' + X = X.sum(dim='y_b', skipna=True) #blocks sum along lat (j) + X = X.mean(dim='x_b', skipna=True) #block mean along lon (i) + return X + +def isum_jstep(X, ndeg=1, W=1.0): + ''' + sum along longitude (i) & one point each ndeg along latitude (j) + ''' + X = X.sum(dim='x_b', skipna=True) + X = X.isel({'y_b': -1}) + return X + +def jsum_istep(X, ndeg=1, W=1.0): + ''' + sum along latitude (j) & one point each ndeg along longitude (i) + ''' + X = X.sum(dim='y_b', skipna=True) + X = X.isel({'x_b': -1}) + return X From 270e7c641651961f0af2e53988f81f1669a512fd Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Mon, 22 Dec 2025 14:37:54 +0100 Subject: [PATCH 12/15] regridding.py : start improving docstring and __main__ section --- FORCINGS/degradation/regridding.py | 38 ++++++++++++++++++++++++++---- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/FORCINGS/degradation/regridding.py b/FORCINGS/degradation/regridding.py index 2542607..25407dc 100644 --- a/FORCINGS/degradation/regridding.py +++ b/FORCINGS/degradation/regridding.py @@ -1,12 +1,29 @@ +import xarray as xr + # REGRIDDING FUNCTIONS -def isum_jmean(X, ndeg=1, W=1.0): +def isum_jmean(X: xr.DataArray, ndeg=1, W=1.0): ''' sum along longitude (i) & mean along latitude (j) no need to weight by area! for e1t + + The algorithm on a single block works as follows: + [ a11 a12 a13 --> sum these ones + a21 a22 a23 + a31 a32 a33 ] + + For e1t, we do prefer to repeat the sum for each j_b, and then average along j_b + + Arguments: + X: xarray DataArray with dims + (..., y, y_b, x, x_b) with dims sizes + + Returns: + Xr: xarray DataArray with dims (..., y, x ) + ''' - X = X.sum(dim='x_b', skipna=True) #blocks sum along lon (i) - X = X.mean(dim='y_b', skipna=True) #block mean along lat (j) - return X + X_yb = X.sum(dim='x_b', skipna=True) #blocks sum along lon (i), returns (..., y, y_b, x) + Xr = X_yb.mean(dim='y_b', skipna=True) #block mean along lat (j) + return Xr def imean_jstep(X, ndeg=1, W=1.0): ''' @@ -63,3 +80,16 @@ def jsum_istep(X, ndeg=1, W=1.0): X = X.sum(dim='y_b', skipna=True) X = X.isel({'x_b': -1}) return X + +if __name__ == "__main__": + import numpy as np + + data = np.ones((20,4,20,4), dtype=np.float32) + arr= xr.DataArray(data, dims=('y','y_b','x','x_b'), coords={'x':np.arange(20), 'y':np.arange(20)}) + a=isum_jmean(arr, ndeg=4, W=None) + + data= np.arange(20).reshape((5,4)).astype(np.float32) + arr= xr.DataArray(data, dims=('y_b','x_b')) + isum_jmean(arr, ndeg=4, W=None) + + From 787bd1c702e70b13c486ebd406d66813a790a4fe Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Mon, 22 Dec 2025 18:45:50 +0100 Subject: [PATCH 13/15] umask and vmask calculated from fine mesh, looking at the cells on right and upper side of the coarse mesh Removing from_tmask() Removing unused fmask, ff --- FORCINGS/degradation/degrade_mesh.py | 60 +++++++++------------------- 1 file changed, 18 insertions(+), 42 deletions(-) diff --git a/FORCINGS/degradation/degrade_mesh.py b/FORCINGS/degradation/degrade_mesh.py index 09e2954..c69de11 100644 --- a/FORCINGS/degradation/degrade_mesh.py +++ b/FORCINGS/degradation/degrade_mesh.py @@ -160,7 +160,6 @@ def load_mesh(maskfile, ndeg=1): M1["nav_lat"] = xpnd_wrap(M["nav_lat"], 'interp', ndeg) M1["nav_lon"] = xpnd_wrap(M["nav_lon"], 'interp', ndeg) M1["nav_lev"] = M["nav_lev"] - M1["fmask"] = xpnd_wrap(M["fmask"], 'edge', ndeg) #there is some error here! M1["tmask"] = xpnd_wrap(M["tmask"], 'edge', ndeg) #there is some error here! M1["umask"] = xpnd_wrap(M["umask"], 'edge', ndeg) #there is some error here! M1["vmask"] = xpnd_wrap(M["vmask"], 'edge', ndeg) #there is some error here! @@ -293,35 +292,7 @@ def waterpt_thresh(X24, thresh=1 ,ndeg=1): X_degr= (X_degr >= (ndeg**2 // thr)).astype(int) return X_degr -def from_tmask(tmask): - ''' - computes umask, vmask, fmask from tmask - after Madec 2015 NEMO Manual, pp. 67 - NB Madec's fmask definition is not correct! - tmask.shape: (1, 141, 380, 1307) - ''' - # CHECK I'M PADDING AT THE RIGHT END OF THE ARRAY!!! - tmask_ipad = tmask.shift(x=1, fill_value=0.0) - tmask_jpad = tmask.shift(y=1, fill_value=0.0) - tmask_ijpad = tmask.shift(x=1, y=1, fill_value=0.0) - # - umask = tmask * tmask_ipad - vmask = tmask * tmask_jpad - # ERROR: THIS fmask DOES NOT CORRESPOND WITH WHAT'S IN meshmask.nc - fmask = tmask * tmask_ipad * tmask_jpad * tmask_ijpad - # - return umask, vmask, fmask - -def noop(X): - ''' - no operation, for variables that need no regridding - e.g. nav_lev - ''' - return X - -# END REGRIDDING OPERATIONS - -def degr_wrap(X24, degr_op:callable, ndeg=1, W=1.0): +def degr_wrap(X24:xr.DataArray, degr_op:callable, ndeg=1, W=1.0): ''' Wrapper to degrade resolution of 2D/3D/4D arrays. @@ -355,13 +326,23 @@ def degr_wrap(X24, degr_op:callable, ndeg=1, W=1.0): X_degr = deadjust_dims(X_degr, nd) # back to original dims return X_degr -def degrade_mesh(M1, thresh=1, ndeg=1): +def degrade_mesh(M1:xr.DataArray, thresh:int=1, ndeg:int=1): ''' generates new reduced mesh for each variable + + Arguments: + M1 : xr.DataArray where for example e2t field is + M1["e2t"] = xpnd_wrap(M["e2t"], 'edge', ndeg) + so, it has dimensions (... y,x ) multipes of ndeg + thresh : integer, number of water points needed to classify a cell as water + ndeg : integer, number of cells to join in i, j directions + + Returns: + M2: xr.DataArray + with dimensions y/ndeg, x/ndeg ''' M2 = {} # - M2["coastp"] = degr_wrap(M1["coastp"], awmean, ndeg, W=M1['At']) M2["e1f"] = degr_wrap(M1["e1f"], rg.isum_jstep, ndeg, W=None) M2["e1t"] = degr_wrap(M1["e1t"], rg.isum_jmean, ndeg, W=None) M2["e1u"] = degr_wrap(M1["e1u"], rg.isum_jmean, ndeg, W=None) @@ -374,9 +355,7 @@ def degrade_mesh(M1, thresh=1, ndeg=1): M2["e3u_0"] = degr_wrap(M1["e3u_0"], uawmean_istep, ndeg, W=M1['Au']) M2["e3v_0"] = degr_wrap(M1["e3v_0"], vawmean_jstep, ndeg, W=M1['Av']) M2["e3w_0"] = degr_wrap(M1["e3w_0"], vwmean, ndeg, W=M1['V']) - M2["ff"] = degr_wrap(M1["ff"], rg.jmean_istep, ndeg, W=None) - #M2["gdept"] = noop(M1["gdept"]) - #M2["gdepw"] = noop(M1["gdepw"]) + M2["gdept"] = degr_wrap(M1["gdept"], vwmean, ndeg, W=M1['V']) M2["gdepw"] = degr_wrap(M1["gdepw"], vwmean, ndeg, W=M1['V']) M2["glamf"] = degr_wrap(M1["glamf"], rg.jmean_istep, ndeg, W=None) @@ -389,13 +368,11 @@ def degrade_mesh(M1, thresh=1, ndeg=1): M2["gphiv"] = degr_wrap(M1["gphiv"], rg.imean_jstep, ndeg, W=None) M2["nav_lat"] = degr_wrap(M1["nav_lat"], rg.jmean_istep, ndeg, W=None) M2["nav_lon"] = degr_wrap(M1["nav_lon"], rg.imean_jstep, ndeg, W=None) - M2["nav_lev"] = noop(M1["nav_lev"]) + M2["nav_lev"] = M1["nav_lev"] M2["tmask"] = waterpt_thresh(M1["tmask"], thresh ,ndeg) - umask, vmask, fmask = from_tmask(M2["tmask"]) - M2["umask"] = umask - M2["vmask"] = vmask - M2["fmask"] = fmask - # + M2['umask'] = degr_wrap(M1['umask'], rg.jsum_istep, ndeg, W=None) + M2['vmask'] = degr_wrap(M1['vmask'], rg.isum_jstep, ndeg, W=None) + M2 = xr.Dataset(M2) return M2 @@ -418,7 +395,6 @@ def cut_med(M2, lon_cut=-8.875, depth_cut=4153.0, biscay_land=True): MMed['tmask'] = MMed['tmask'] * ~land MMed['umask'] = MMed['umask'] * ~land MMed['vmask'] = MMed['vmask'] * ~land - MMed['fmask'] = MMed['fmask'] * ~land return MMed def dump_mesh(M2, outfile): From f6c5679eba192b0d925af4b85ce5a76c88f3d4f2 Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Mon, 22 Dec 2025 18:50:54 +0100 Subject: [PATCH 14/15] putting load_mesh() close to degrade_mesh() for readability --- FORCINGS/degradation/degrade_mesh.py | 114 ++++++++++++++------------- 1 file changed, 58 insertions(+), 56 deletions(-) diff --git a/FORCINGS/degradation/degrade_mesh.py b/FORCINGS/degradation/degrade_mesh.py index c69de11..27c2fcd 100644 --- a/FORCINGS/degradation/degrade_mesh.py +++ b/FORCINGS/degradation/degrade_mesh.py @@ -126,62 +126,6 @@ def e3_2_gdep(M, tuvw): gdep = e3.cumsum(dim="z") - (e3 / 2.0) return gdep -def load_mesh(maskfile, ndeg=1): - ''' - load meshmask, expand fields - ''' - M = xr.open_dataset(maskfile) - M1 = {} - # - M1["coastp"] = xpnd_wrap(M["coastp"], 'edge', ndeg) - M1["e1f"] = xpnd_wrap(M["e1f"], 'interp', ndeg) - M1["e1t"] = xpnd_wrap(M["e1t"], 'interp', ndeg) - M1["e1u"] = xpnd_wrap(M["e1u"], 'interp', ndeg) - M1["e1v"] = xpnd_wrap(M["e1v"], 'interp', ndeg) - M1["e2f"] = xpnd_wrap(M["e2f"], 'edge', ndeg) - M1["e2t"] = xpnd_wrap(M["e2t"], 'edge', ndeg) - M1["e2u"] = xpnd_wrap(M["e2u"], 'edge', ndeg) - M1["e2v"] = xpnd_wrap(M["e2v"], 'edge', ndeg) - M1["e3t_0"] = xpnd_wrap(M["e3t_0"], 'edge', ndeg) - M1["e3u_0"] = xpnd_wrap(M["e3u_0"], 'edge', ndeg) - M1["e3v_0"] = xpnd_wrap(M["e3v_0"], 'edge', ndeg) - M1["e3w_0"] = xpnd_wrap(M["e3w_0"], 'edge', ndeg) - M1["ff"] = xpnd_wrap(M["ff"], 'interp', ndeg) - M1["gdept"] = xpnd_wrap(e3_2_gdep(M, 't'), 'edge', ndeg) #M["gdept"] - M1["gdepw"] = xpnd_wrap(e3_2_gdep(M, 'w'), 'edge', ndeg) #M["gdepw"] - M1["glamf"] = xpnd_wrap(M["glamf"], 'interp', ndeg) - M1["glamt"] = xpnd_wrap(M["glamt"], 'interp', ndeg) - M1["glamu"] = xpnd_wrap(M["glamu"], 'interp', ndeg) - M1["glamv"] = xpnd_wrap(M["glamv"], 'interp', ndeg) - M1["gphif"] = xpnd_wrap(M["gphif"], 'interp', ndeg) - M1["gphit"] = xpnd_wrap(M["gphit"], 'interp', ndeg) - M1["gphiu"] = xpnd_wrap(M["gphiu"], 'interp', ndeg) - M1["gphiv"] = xpnd_wrap(M["gphiv"], 'interp', ndeg) - M1["nav_lat"] = xpnd_wrap(M["nav_lat"], 'interp', ndeg) - M1["nav_lon"] = xpnd_wrap(M["nav_lon"], 'interp', ndeg) - M1["nav_lev"] = M["nav_lev"] - M1["tmask"] = xpnd_wrap(M["tmask"], 'edge', ndeg) #there is some error here! - M1["umask"] = xpnd_wrap(M["umask"], 'edge', ndeg) #there is some error here! - M1["vmask"] = xpnd_wrap(M["vmask"], 'edge', ndeg) #there is some error here! - # Compute cell surface areas and volume - M1['At'] = xr.DataArray(data = M1['e1t'].values[:] * M1['e2t'].values[:], - dims = ('time','z_a','y','x'), - name = 'At') - M1['Au'] = xr.DataArray(data = M1['e2u'].values[:] * M1['e3u_0'].values[:], - dims = ('time','z','y','x'), - name = 'Au') - M1['Av'] = xr.DataArray(data = M1['e2v'].values[:] * M1['e3v_0'].values[:], - dims = ('time','z','y','x'), - name = 'Av') - M1['V'] = xr.DataArray(data = M1['e1v'].values[:] * M1['e2u'].values[:] * M1['e3t_0'].values[:], - dims = ('time','z','y','x'), - name = 'V') - # Bonus variable for degrading forcings - M1['h_column_t'] = (M1['tmask'] * M1['e3t_0']).sum(dim='z') - # - M.close() - M1 = xr.Dataset(M1) - return M1 # DEGRADE RESOLUTION @@ -326,6 +270,64 @@ def degr_wrap(X24:xr.DataArray, degr_op:callable, ndeg=1, W=1.0): X_degr = deadjust_dims(X_degr, nd) # back to original dims return X_degr +def load_mesh(maskfile, ndeg=1): + ''' + load meshmask, expand fields + ''' + M = xr.open_dataset(maskfile) + M1 = {} + # + M1["coastp"] = xpnd_wrap(M["coastp"], 'edge', ndeg) + M1["e1f"] = xpnd_wrap(M["e1f"], 'interp', ndeg) + M1["e1t"] = xpnd_wrap(M["e1t"], 'interp', ndeg) + M1["e1u"] = xpnd_wrap(M["e1u"], 'interp', ndeg) + M1["e1v"] = xpnd_wrap(M["e1v"], 'interp', ndeg) + M1["e2f"] = xpnd_wrap(M["e2f"], 'edge', ndeg) + M1["e2t"] = xpnd_wrap(M["e2t"], 'edge', ndeg) + M1["e2u"] = xpnd_wrap(M["e2u"], 'edge', ndeg) + M1["e2v"] = xpnd_wrap(M["e2v"], 'edge', ndeg) + M1["e3t_0"] = xpnd_wrap(M["e3t_0"], 'edge', ndeg) + M1["e3u_0"] = xpnd_wrap(M["e3u_0"], 'edge', ndeg) + M1["e3v_0"] = xpnd_wrap(M["e3v_0"], 'edge', ndeg) + M1["e3w_0"] = xpnd_wrap(M["e3w_0"], 'edge', ndeg) + M1["ff"] = xpnd_wrap(M["ff"], 'interp', ndeg) + M1["gdept"] = xpnd_wrap(e3_2_gdep(M, 't'), 'edge', ndeg) #M["gdept"] + M1["gdepw"] = xpnd_wrap(e3_2_gdep(M, 'w'), 'edge', ndeg) #M["gdepw"] + M1["glamf"] = xpnd_wrap(M["glamf"], 'interp', ndeg) + M1["glamt"] = xpnd_wrap(M["glamt"], 'interp', ndeg) + M1["glamu"] = xpnd_wrap(M["glamu"], 'interp', ndeg) + M1["glamv"] = xpnd_wrap(M["glamv"], 'interp', ndeg) + M1["gphif"] = xpnd_wrap(M["gphif"], 'interp', ndeg) + M1["gphit"] = xpnd_wrap(M["gphit"], 'interp', ndeg) + M1["gphiu"] = xpnd_wrap(M["gphiu"], 'interp', ndeg) + M1["gphiv"] = xpnd_wrap(M["gphiv"], 'interp', ndeg) + M1["nav_lat"] = xpnd_wrap(M["nav_lat"], 'interp', ndeg) + M1["nav_lon"] = xpnd_wrap(M["nav_lon"], 'interp', ndeg) + M1["nav_lev"] = M["nav_lev"] + M1["tmask"] = xpnd_wrap(M["tmask"], 'edge', ndeg) #there is some error here! + M1["umask"] = xpnd_wrap(M["umask"], 'edge', ndeg) #there is some error here! + M1["vmask"] = xpnd_wrap(M["vmask"], 'edge', ndeg) #there is some error here! + # Compute cell surface areas and volume + M1['At'] = xr.DataArray(data = M1['e1t'].values[:] * M1['e2t'].values[:], + dims = ('time','z_a','y','x'), + name = 'At') + M1['Au'] = xr.DataArray(data = M1['e2u'].values[:] * M1['e3u_0'].values[:], + dims = ('time','z','y','x'), + name = 'Au') + M1['Av'] = xr.DataArray(data = M1['e2v'].values[:] * M1['e3v_0'].values[:], + dims = ('time','z','y','x'), + name = 'Av') + M1['V'] = xr.DataArray(data = M1['e1v'].values[:] * M1['e2u'].values[:] * M1['e3t_0'].values[:], + dims = ('time','z','y','x'), + name = 'V') + # Bonus variable for degrading forcings + M1['h_column_t'] = (M1['tmask'] * M1['e3t_0']).sum(dim='z') + # + M.close() + M1 = xr.Dataset(M1) + return M1 + + def degrade_mesh(M1:xr.DataArray, thresh:int=1, ndeg:int=1): ''' generates new reduced mesh for each variable From 53e29bd209ffeb5c79154ba2891bb40894a0f952 Mon Sep 17 00:00:00 2001 From: Giorgio Bolzon Date: Mon, 22 Dec 2025 19:13:35 +0100 Subject: [PATCH 15/15] Adding degrade_restart.py and yaml files in Giovanni Galli's version --- FORCINGS/degradation/degrade_forcings.yaml | 7 ++ FORCINGS/degradation/degrade_mesh.yaml | 8 ++ FORCINGS/degradation/degrade_restart.py | 140 +++++++++++++++++++++ FORCINGS/degradation/degrade_restart.yaml | 6 + 4 files changed, 161 insertions(+) create mode 100644 FORCINGS/degradation/degrade_forcings.yaml create mode 100644 FORCINGS/degradation/degrade_mesh.yaml create mode 100644 FORCINGS/degradation/degrade_restart.py create mode 100644 FORCINGS/degradation/degrade_restart.yaml diff --git a/FORCINGS/degradation/degrade_forcings.yaml b/FORCINGS/degradation/degrade_forcings.yaml new file mode 100644 index 0000000..f56b78e --- /dev/null +++ b/FORCINGS/degradation/degrade_forcings.yaml @@ -0,0 +1,7 @@ +maskfile: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/MASK24/meshmask_141.nc' +maskfile_d: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/qDEG_SETUP/MASKS/meshmask_025_z141.nc' +ffdir: '/leonardo_work/OGS_devC_0/Neccton_hindcast1999_2022_v22/wrkdir/MODEL/FORCINGS' +outdir: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/qDEG_SETUP/FORCINGS/' +ndeg: 6 +y0: 1999 +yE: 2001 diff --git a/FORCINGS/degradation/degrade_mesh.yaml b/FORCINGS/degradation/degrade_mesh.yaml new file mode 100644 index 0000000..3698192 --- /dev/null +++ b/FORCINGS/degradation/degrade_mesh.yaml @@ -0,0 +1,8 @@ +maskfile: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/MASK24/meshmask_141.nc' +outfile: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/qDEG_SETUP/MASKS/meshmask_025_z141.nc' +outfile_med: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/qDEG_SETUP/MASKS/meshmask_025_z125.nc' +# by how much cells to degrade (nr of cells to join in i, j) +ndeg: 6 +# threshold of points to assign waterpoint to degraded mesh +thresh: 1 #THE ONLY OPTION THAT ENSURES NO RIVERS END UP ON LAND + diff --git a/FORCINGS/degradation/degrade_restart.py b/FORCINGS/degradation/degrade_restart.py new file mode 100644 index 0000000..6a7dd67 --- /dev/null +++ b/FORCINGS/degradation/degrade_restart.py @@ -0,0 +1,140 @@ +import numpy as np +import xarray as xr +from argparse import ArgumentParser +import yaml +from glob import glob +# +import degrade_mesh as dm +import degrade_forcings as df + +''' +degrades the horizontal resolution of BFM restarts by an integer value +joins boxes of ndeg x ndeg cells, does volume weighted mean. +needs both original and degraded-resolution meshmask +(run degrade mesh.py before) + +usage: +python degrade_mesh.py -y degrade_mesh.yaml + +yaml file example: +mesh_in: '/path/to/meshmask_original.nc' +mesh_out: '/path/to/meshmask_degraded.nc' +indir: '/path/to/restarts/rst.nc' #also with wildcards '*' +outdir: '/path/to/output_dir/' +infile: 'RST.*-00:00:00.__VNAME__.nc' #__VNAME__ gets replaced and is needed to read the variable +ndeg: 6 #size of boxes for resolution degradation +''' + +def argument(): + parser = ArgumentParser() + parser.add_argument('--yamlfile','-y', + type=str, + required=True, + help='file with all parameters') + return parser.parse_args() + +def load_parameters(): + yamlfile=argument().yamlfile + # yamlfile = 'degrade_restart.yaml' + with open(yamlfile) as f: + Params = yaml.load(f, Loader=yaml.Loader) + return(Params) + +def get_Volume(M): + ''' + gets cell volume for weighted mean + ''' + V0 = M['e1t'].values[:] * M['e2t'].values[:] * M['e3t_0'] + nanmask = M['tmask'].values[:] + nanmask[nanmask==0.0] = np.nan + V0 = V0 * nanmask #else you see border effects at the coast + return V0 + +def load_rst(infile, vname, ndeg=1): + ''' + loads the data array from a variable in the restart + ''' + D = xr.open_dataset(infile) + D1 = {} + # + D1[vname] = dm.xpnd_wrap(D[vname], 'edge', ndeg) + # + D1['nav_lat'] = dm.xpnd_wrap(D['nav_lat'], 'interp', ndeg) + D1['nav_lon'] = dm.xpnd_wrap(D['nav_lon'], 'interp', ndeg) + #D1[vname].coords['nav_lat'].values[:] = D1['nav_lat'].values[:] + #D1[vname].coords['nav_lon'].values[:] = D1['nav_lon'].values[:] + D1['time'] = D['time'] + D1['nav_lev'] = D['nav_lev'] + # + D1 = xr.Dataset(D1) + D1 = D1.assign_attrs(D.attrs) + D.close() + return D1 + +def init_rst(C): + ''' + lat, lon are the same for all files + ''' + R = {} + R['nav_lon'] = C['glamt'] + R['nav_lat'] = C['gphit'] + return R + +def degrade_bgc(DI, V0, C, vname, ndeg=1): + ''' + degrades resolution of BFM restart file + ''' + R = init_rst(C) + R[vname] = df.degrade_wrap(DI[vname], dm.vwmean, V0, ndeg) + R['nav_lev'] = DI['nav_lev'] + R = xr.Dataset(R) + R = R.assign_attrs(DI.attrs) + return R + +def get_flist(Params): + indir = Params['indir'] + infile = Params['infile'] + flist = glob(indir+infile.replace('__VNAME__', '*')) + itrbl = [(ff.split('.')[-2], ff) for ff in flist] + return itrbl + +if __name__=='__main__': + try: + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nranks = comm.size + except: + comm = None + rank = 0 + nranks = 1 + # + Params = load_parameters() + mesh_in = Params['mesh_in'] + mesh_out = Params['mesh_out'] + indir = Params['indir'] + infile = Params['infile'] + outdir = Params['outdir'] + ndeg = Params['ndeg'] + # + if rank==0: + M = df.load_mesh_light(mesh_in, ndeg) + C = df.load_coords_degraded(mesh_out) + V0 = get_Volume(M) + else: + M = None + C = None + V0 = None + if nranks > 1: + M = comm.bcast(M, root=0) + C = comm.bcast(C, root=0) + V0 = comm.bcast(V0, root=0) + # HERE LOOP ON FILES AND VARIABLES + itrbl = get_flist(Params) + for vname, fname in itrbl[rank::nranks]: + vname = 'TRN'+vname + print(f'degrado: {vname}') + # + DI = load_rst(fname, vname, ndeg) + Dd = degrade_bgc(DI, V0, C, vname, ndeg) + outfile = outdir+fname.split('/')[-1] + dm.dump_mesh(Dd, outfile) diff --git a/FORCINGS/degradation/degrade_restart.yaml b/FORCINGS/degradation/degrade_restart.yaml new file mode 100644 index 0000000..e32d6b0 --- /dev/null +++ b/FORCINGS/degradation/degrade_restart.yaml @@ -0,0 +1,6 @@ +mesh_in: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/MASK24/meshmask.nc' +mesh_out: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/qDEG_SETUP/MASKS/meshmask_025_z125.nc' +indir: '/leonardo_work/OGS_devC_0/RESTARTS_NECCTON_1999/' +outdir: '/leonardo_work/OGS23_PRACE_IT_0/ggalli00/OGSTM-BFM/qDEG_SETUP/RESTARTS/' +infile: 'RST.*-00:00:00.__VNAME__.nc' +ndeg: 6