diff --git a/example/calcul_pressure-lagrangien-2-eulerien-patch-correctif.py b/example/calcul_pressure-lagrangien-2-eulerien-patch-correctif.py new file mode 100644 index 0000000000000000000000000000000000000000..a7d953753876e85ac0669bb069917ca596861397 --- /dev/null +++ b/example/calcul_pressure-lagrangien-2-eulerien-patch-correctif.py @@ -0,0 +1,186 @@ +#!/usr/bin/python + +import sys +import numpy as np +import netCDF4 as nc +from numpy import arange, dtype +import math + +#******************************************************************************* +def main(args): + + """USE + {0} mask.nc input.nc output.nc + +DESCRIPTION + Computes pressure from density and e3t (that changes with SSH). + +""" + if len(args)<6 or args[1]=='-h' or args[1]=='--help': + print( main.func_doc.format(args[0]) ) + return 0 + + + +#=============================================================================== +# input +#=============================================================================== + + file_pressure = args[1] + fichierpression = nc.Dataset(file_pressure) + + pbc_a=fichierpression.variables['Pbc_a'][:,:,:] + pbc_G=fichierpression.variables['Pbc_G'][:,:,:] + depth=fichierpression.variables['dept'][:] + lon=fichierpression.variables['lont'][:] + lat=fichierpression.variables['latt'][:] + + file_ssh = args[2] + fichierssh = nc.Dataset(file_ssh) + + ssh_a=fichierssh.variables['ssh_a'][:,:] + ssh_G=fichierssh.variables['ssh_G'][:,:] + + jj=pbc_a.shape[1] + ii=pbc_a.shape[2] + kk=pbc_a.shape[0] + + file_mean = args[3] + fichierrhop = nc.Dataset(file_mean) + rhop_m=fichierrhop.variables['rhop_mean'][0,:,:] + + file_mask = args[5] + fichiermask = nc.Dataset(file_mask) + gdepw=fichiermask.variables['gdepw'][0,:,:,:] + gdept_0_all=fichiermask.variables['gdept_0'][0,74,:,:] + gdept_0=fichiermask.variables['gdept_0'][0,0,:,:] + tmask = fichiermask.variables['tmask'][0, :, :, :] + + + + +#=============================================================================== +# output +#=============================================================================== + + p_file=args[4] + print('Creating '+p_file) + + # this will overwrite the file + ncfile = nc.Dataset(p_file,'w', format='NETCDF3_CLASSIC') + setattr(ncfile,'history',' '.join(args)) + ncfile.createDimension('x',ii) + ncfile.createDimension('y',jj) + ncfile.createDimension('dept',kk) + + + data = ncfile.createVariable('lont',dtype('float64').char,('y','x')) + data[:] = lon + setattr(data, 'units',"degrees_north") + setattr(data, 'standard_name',"longitude_at_T_location") + + data = ncfile.createVariable('latt',dtype('float64').char,('y','x')) + data[:] = lat + setattr(data, 'units',"degrees_east") + setattr(data, 'standard_name',"latitude_at_T_location") + + data = ncfile.createVariable('dept',dtype('float64').char,('dept')) + data[:] = depth + setattr(data, 'units',"m") + setattr(data, 'positive',"down") + setattr(data, 'standard_name',"depth_at_T_location") + setattr(data, 'axis',"Z") + + p_corr_a = ncfile.createVariable('P_corr_A',dtype('float32').char,('dept','y','x')) + setattr(p_corr_a, 'standard_name',"Amplitude Of Pressure Total") + setattr(p_corr_a, 'units',"m-1.kg.s-2") + setattr(p_corr_a, 'coordinates', "dept latt lont") + + p_corr_G = ncfile.createVariable('P_corr_G',dtype('float32').char,('dept','y','x')) + setattr(p_corr_G, 'standard_name',"Phase Of Pressure Total") + setattr(p_corr_G, 'units',"m-1.kg.s-2") + setattr(p_corr_G, 'coordinates', "dept latt lont") + +#=============================================================================== +# process +#=============================================================================== + + print( 'Dimensions: kk={} jj={} ii={} '.format(kk,jj,ii) ) + + p_zr = np.zeros(shape=(kk,jj,ii)) + p_zi = np.zeros(shape=(kk,jj,ii)) + + ssh_zr = np.zeros(shape=(jj,ii)) + ssh_zi = np.zeros(shape=(jj,ii)) + + grav=9.81 + + + for k in range (0,75): + for j in range (0,484): + for i in range (0,388): + + if math.isnan(pbc_a[k, j, i])==True: + p_zr[k, j, i] == np.NaN; + p_zi[k, j, i] == np.NaN; + else: + p_zr[k, j, i]=pbc_a[k, j, i]*math.cos(math.radians(pbc_G[k, j, i])) + p_zi[k, j, i]=pbc_a[k, j, i]*math.sin(math.radians(pbc_G[k, j, i])) + + for j in range (0,484): + for i in range (0,388): + + if math.isnan(ssh_a[j, i])==True: + ssh_zr[j, i] == np.NaN; + ssh_zi[j, i] == np.NaN; + else: + ssh_zr[j, i]=ssh_a[j, i]*math.cos(math.radians(ssh_G[j, i])) + ssh_zi[j, i]=ssh_a[j, i]*math.sin(math.radians(ssh_G[j, i])) + + p_corr_zr = np.zeros(shape=(kk,jj,ii)) + p_corr_zi = np.zeros(shape=(kk,jj,ii)) + zz = np.zeros(shape=(jj,ii)) + + for k in range (0,75): + for j in range(0, 484): + for i in range(0, 388): + + zz[j, i] = (0.5 * gdept_0[j,i]) + gdepw[k, j, i] + p_corr_zr[k, j, i] = p_zr[k, j, i] + grav * rhop_m[j, i] * ssh_zr[j, i] * ((gdept_0_all[j,i] - zz[j,i]) / gdept_0_all[j,i]) + p_corr_zi[k, j, i] = p_zi[k, j, i] + grav * rhop_m[j, i] * ssh_zi[j, i] * ((gdept_0_all[j,i] - zz[j,i]) / gdept_0_all[j,i]) + + + p_a = np.zeros(shape=(kk,jj,ii)) + p_G = np.zeros(shape=(kk,jj,ii)) + + + for k in range (0,75): + for j in range (0,484): + for i in range (0,388): + p_a[k, j, i] = (p_corr_zr[k, j, i] * p_corr_zr[k, j, i] + p_corr_zi[k, j, i] * p_corr_zi[k, j, i])**(0.5) + p_G[k, j, i] = math.degrees(math.atan2(p_corr_zi[k, j, i], p_corr_zr[k, j, i])) + + p_a[tmask == 0] = nc.default_fillvals['f4'] + p_G[tmask == 0] = nc.default_fillvals['f4'] + p_corr_a[:,:,:] = p_a + p_corr_G[:,:,:] = p_G + + + + +#=============================================================================== +# clean-up +#=============================================================================== + + fichierpression.close() + fichierssh.close() + fichierrhop.close() + fichiermask.close() + ncfile.close() + + +#******************************************************************************* +if __name__ == '__main__': + #file:///usr/share/doc/packages/python/html/python-2.7-docs-html/library/sys.html#sys.exit + sys.exit(main(sys.argv)) + diff --git a/example/calcul_pressure-lagrangien.py b/example/calcul_pressure-lagrangien.py new file mode 100644 index 0000000000000000000000000000000000000000..5f01e81bcdfff26f390f1b32dea5d84a109549a6 --- /dev/null +++ b/example/calcul_pressure-lagrangien.py @@ -0,0 +1,138 @@ +#!/usr/bin/python + +import sys +import numpy as np +import netCDF4 as nc +from numpy import arange, dtype + +#******************************************************************************* +def main(args): + + """USE + {0} mask.nc input.nc output.nc + +DESCRIPTION + Computes pressure from density and e3t (that changes with SSH). + +""" + if len(args)<4 or args[1]=='-h' or args[1]=='--help': + print( main.func_doc.format(args[0]) ) + return 0 + +# =============================================================================== +# mask file +# =============================================================================== + + fichiermask=args[1] + + mask_file=nc.Dataset(fichiermask) + depth=mask_file.variables['gdept_1d'][:] + tmask=mask_file.variables['tmask'][0,:,:,:] + mask_file.close() + +#=============================================================================== +# input +#=============================================================================== + + fichierjulien=args[2] + lefichier = nc.Dataset(fichierjulien) + + e3tv=lefichier.variables['e3t'] + ji=e3tv.shape[3] + jj=e3tv.shape[2] + jk=e3tv.shape[1] + jl=e3tv.shape[0] + + e3wv=lefichier.variables['e3w'] + rhopv=lefichier.variables['rhop'] + itv=lefichier.variables['time_counter'] + + lon=lefichier.variables['nav_lon_grid_T'][:] + lat=lefichier.variables['nav_lat_grid_T'][:] + +#=============================================================================== +# output +#=============================================================================== + + p_file=args[3] + print('Creating '+p_file) + + # this will overwrite the file + ncfile = nc.Dataset(p_file,'w', format='NETCDF3_CLASSIC') + setattr(ncfile,'history',' '.join(args)) + ncfile.createDimension('x',ji) + ncfile.createDimension('y',jj) + ncfile.createDimension('dept',jk) + ncfile.createDimension('time') # unlimited + + data = ncfile.createVariable('lont',dtype('float64').char,('y','x')) + data[:] = lon + setattr(data, 'units',"degrees_north") + setattr(data, 'standard_name',"longitude_at_T_location") + + data = ncfile.createVariable('latt',dtype('float64').char,('y','x')) + data[:] = lat + setattr(data, 'units',"degrees_east") + setattr(data, 'standard_name',"latitude_at_T_location") + + data = ncfile.createVariable('dept',dtype('float64').char,('dept')) + data[:] = depth + setattr(data, 'units',"m") + setattr(data, 'positive',"down") + setattr(data, 'standard_name',"depth_at_T_location") + setattr(data, 'axis',"Z") + + otv = ncfile.createVariable('time',dtype('float64').char,('time')) + setattr(otv, 'units',itv.units) + setattr(otv, 'standard_name',"time") + setattr(otv, 'axis',"T") + + pv = ncfile.createVariable('Pbc',dtype('float32').char,('time','dept','y','x')) + setattr(pv, 'standard_name',"Amplitude Of Pressure Total") + setattr(pv, 'units',"m-1.kg.s-2") + setattr(pv, 'coordinates', "time dept latt lont") + +#=============================================================================== +# process +#=============================================================================== + + print( 'Dimensions: jl={} jk={} jj={} ji={}'.format(jl,jk,jj,ji) ) + + p_bc=np.zeros(shape=(jk,jj,ji)) + + grav=9.81 + + for l in range(0,jl): + + time=itv[l] + print( 'Time frame {}/{}'.format(l,jl) ) + + e3t=e3tv[l, 0, :, :] + e3w=e3wv[l, :, :, :] + rhop=rhopv[l, :, :, :] + + for k in range(0,jk): + if k==0: + p0=0.5*grav*e3t*rhop[k,:,:] + else: + p0+=0.5*grav*e3w[k,:,:]*(rhop[k,:,:]+rhop[k-1,:,:]) + + p0*=tmask[k,:,:] + p_bc[k,:,:]=p0 + + p_bc[tmask==0]=nc.default_fillvals['f4'] + pv[l,:,:,:] = p_bc + otv[l] = time + +#=============================================================================== +# clean-up +#=============================================================================== + + lefichier.close() + ncfile.close() + + +#******************************************************************************* +if __name__ == '__main__': + #file:///usr/share/doc/packages/python/html/python-2.7-docs-html/library/sys.html#sys.exit + sys.exit(main(sys.argv)) diff --git a/example/calcul_pressure_eulerien.py b/example/calcul_pressure_eulerien.py new file mode 100644 index 0000000000000000000000000000000000000000..4c39d629fb103f4a3942b459f53fa12b82c6b023 --- /dev/null +++ b/example/calcul_pressure_eulerien.py @@ -0,0 +1,143 @@ +#!/usr/bin/python + +import sys +import numpy as np +import netCDF4 as nc +from numpy import arange, dtype + +#******************************************************************************* +def main(args): + + """USE + {0} mask.nc input.nc output.nc + +DESCRIPTION + Computes pressure from density and e3t (that changes with SSH). + +""" + if len(args)<4 or args[1]=='-h' or args[1]=='--help': + print( main.func_doc.format(args[0]) ) + return 0 + +# =============================================================================== +# mask file +# =============================================================================== + + fichiermask=args[1] + + mask_file=nc.Dataset(fichiermask) + depth=mask_file.variables['gdept_1d'][:] + tmask=mask_file.variables['tmask'][0,:,:,:] + e3t = mask_file.variables['e3t'][:,:,:,:] + e3w = mask_file.variables['e3w'][:,:,:,:] + mask_file.close() + +#=============================================================================== +# input +#=============================================================================== + + fichierjulien=args[2] + lefichier = nc.Dataset(fichierjulien) + + rhopv = lefichier.variables['rhop'] + ji=rhopv.shape[3] + jj=rhopv.shape[2] + jk=rhopv.shape[1] + jl=rhopv.shape[0] + + ssh1 = lefichier.variables['ssh'] + + + + itv=lefichier.variables['time_counter'] + + lon=lefichier.variables['nav_lon_grid_T'][:] + lat=lefichier.variables['nav_lat_grid_T'][:] + +#=============================================================================== +# output +#=============================================================================== + + p_file=args[3] + print('Creating '+p_file) + + # this will overwrite the file + ncfile = nc.Dataset(p_file,'w', format='NETCDF3_CLASSIC') + setattr(ncfile,'history',' '.join(args)) + ncfile.createDimension('x',ji) + ncfile.createDimension('y',jj) + ncfile.createDimension('dept',jk) + ncfile.createDimension('time') # unlimited + + data = ncfile.createVariable('lont',dtype('float64').char,('y','x')) + data[:] = lon + setattr(data, 'units',"degrees_north") + setattr(data, 'standard_name',"longitude_at_T_location") + + data = ncfile.createVariable('latt',dtype('float64').char,('y','x')) + data[:] = lat + setattr(data, 'units',"degrees_east") + setattr(data, 'standard_name',"latitude_at_T_location") + + data = ncfile.createVariable('dept',dtype('float64').char,('dept')) + data[:] = depth + setattr(data, 'units',"m") + setattr(data, 'positive',"down") + setattr(data, 'standard_name',"depth_at_T_location") + setattr(data, 'axis',"Z") + + otv = ncfile.createVariable('time',dtype('float64').char,('time')) + setattr(otv, 'units',itv.units) + setattr(otv, 'standard_name',"time") + setattr(otv, 'axis',"T") + + pv = ncfile.createVariable('Ptot',dtype('float32').char,('time','dept','y','x')) + setattr(pv, 'standard_name',"Amplitude Of Pressure Total") + setattr(pv, 'units',"m-1.kg.s-2") + setattr(pv, 'coordinates', "time dept latt lont") + +#=============================================================================== +# process +#=============================================================================== + + print( 'Dimensions: jl={} jk={} jj={} ji={}'.format(jl,jk,jj,ji) ) + + p_tot=np.zeros(shape=(jk,jj,ji)) + + grav=9.81 + + for l in range(0,jl): + + time=itv[l] + print( 'Time frame {}/{}'.format(l,jl) ) + + e3t1=e3t[0, 0, :, :] + e3w1=e3w[0, :, :, :] + rhop=rhopv[l, :, :, :] + ssh=ssh1[l, :, :] + + for k in range(0,jk): + if k==0: + p0=0.5*grav*e3t1*rhop[k,:,:]+(rhop[k,:,:]*grav*ssh[:,:]) + else: + p0+=0.5*grav*e3w1[k,:,:]*(rhop[k,:,:]+rhop[k-1,:,:]) + + p0*=tmask[k,:,:] + p_tot[k,:,:]=p0 + + p_tot[tmask==0]=nc.default_fillvals['f4'] + pv[l,:,:,:] = p_tot + otv[l] = time + +#=============================================================================== +# clean-up +#=============================================================================== + + lefichier.close() + ncfile.close() + + +#******************************************************************************* +if __name__ == '__main__': + #file:///usr/share/doc/packages/python/html/python-2.7-docs-html/library/sys.html#sys.exit + sys.exit(main(sys.argv)) diff --git a/example/diff_eulerien-lagrangien.py b/example/diff_eulerien-lagrangien.py new file mode 100644 index 0000000000000000000000000000000000000000..febcc5002b231d45873f4cbc812251cb8d94be61 --- /dev/null +++ b/example/diff_eulerien-lagrangien.py @@ -0,0 +1,113 @@ +#!/usr/bin/python + +import sys +import numpy as np +import netCDF4 as nc +from numpy import arange, dtype +import math + +#******************************************************************************* +def main(args): + + """USE + {0} mask.nc input.nc output.nc + +DESCRIPTION + Computes pressure from density and e3t (that changes with SSH). + +""" + if len(args)<4 or args[1]=='-h' or args[1]=='--help': + print( main.func_doc.format(args[0]) ) + return 0 + + + +#=============================================================================== +# input +#=============================================================================== + + file_pressure_eulerien = args[1] + fichierpression1 = nc.Dataset(file_pressure_eulerien) + + ptot_a=fichierpression1.variables['Ptot_a'][:,:,:] + ptot_G=fichierpression1.variables['Ptot_G'][:,:,:] + depth=fichierpression1.variables['dept'][:] + lon=fichierpression1.variables['lont'][:] + lat=fichierpression1.variables['latt'][:] + + jj=ptot_a.shape[1] + ii=ptot_a.shape[2] + kk=ptot_a.shape[0] + + + file_pressure_lagpatch = args[2] + fichierpression2 = nc.Dataset(file_pressure_lagpatch) + P_corr_A=fichierpression2.variables['P_corr_A'][:,:,:] + P_corr_G=fichierpression2.variables['P_corr_G'][:,:,:] + + + +#=============================================================================== +# output +#=============================================================================== + + p_file=args[3] + print('Creating '+p_file) + + # this will overwrite the file + ncfile = nc.Dataset(p_file,'w', format='NETCDF3_CLASSIC') + setattr(ncfile,'history',' '.join(args)) + ncfile.createDimension('x',ii) + ncfile.createDimension('y',jj) + ncfile.createDimension('dept',kk) + + + data = ncfile.createVariable('lont',dtype('float64').char,('y','x')) + data[:] = lon + setattr(data, 'units',"degrees_north") + setattr(data, 'standard_name',"longitude_at_T_location") + + data = ncfile.createVariable('latt',dtype('float64').char,('y','x')) + data[:] = lat + setattr(data, 'units',"degrees_east") + setattr(data, 'standard_name',"latitude_at_T_location") + + data = ncfile.createVariable('dept',dtype('float64').char,('dept')) + data[:] = depth + setattr(data, 'units',"m") + setattr(data, 'positive',"down") + setattr(data, 'standard_name',"depth_at_T_location") + setattr(data, 'axis',"Z") + + p_diff_a = ncfile.createVariable('P_corr_A',dtype('float32').char,('dept','y','x')) + setattr(p_diff_a, 'standard_name',"Amplitude Of Pressure Total") + setattr(p_diff_a, 'units',"m-1.kg.s-2") + setattr(p_diff_a, 'coordinates', "dept latt lont") + + p_diff_G = ncfile.createVariable('P_corr_G',dtype('float32').char,('dept','y','x')) + setattr(p_diff_G, 'standard_name',"Phase Of Pressure Total") + setattr(p_diff_G, 'units',"m-1.kg.s-2") + setattr(p_diff_G, 'coordinates', "dept latt lont") + +#=============================================================================== +# process +#=============================================================================== + + p_diff_a[:,:,:] = ptot_a[:,:,:] - P_corr_A[:,:,:] + p_diff_G[:,:,:] = ptot_G[:,:,:] - P_corr_G[:,:,:] + + +#=============================================================================== +# clean-up +#=============================================================================== + + fichierpression1.close() + fichierpression2.close() + ncfile.close() + + +#******************************************************************************* +if __name__ == '__main__': + #file:///usr/share/doc/packages/python/html/python-2.7-docs-html/library/sys.html#sys.exit + sys.exit(main(sys.argv)) +