convert_pt_to_CT.py 2.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293
  1. #!/usr/bin/env python
  2. # L. Brodeau, 2015
  3. # Potential temperature to conservative temperature (TEOS 10)
  4. import sys
  5. import os
  6. import numpy as nmp
  7. from netCDF4 import Dataset
  8. from string import replace
  9. import gsw
  10. if len(sys.argv) != 5:
  11. print 'Usage: '+sys.argv[0]+' <Temperature_file_to_convert> <temperature_name> <Absolute_salinity_file> <salinity_name>'
  12. sys.exit(0)
  13. cf_temp = sys.argv[1]
  14. cv_temp = sys.argv[2]
  15. cf_sal = sys.argv[3]
  16. cv_sal = sys.argv[4]
  17. cf_out = replace(cf_temp, cv_temp, cv_temp+'_TEOS10')
  18. os.system('rm -f '+cf_out)
  19. os.system('cp '+cf_temp+' '+cf_out)
  20. print '\n'
  21. f_sal = Dataset(cf_sal) # r+ => can read and write in the file... )
  22. vcdim = f_sal.variables[cv_sal].dimensions
  23. cv_t = vcdim[0]; print ' *** record dimension is called "'+cv_t+'"'
  24. Nt = f_sal.dimensions[cv_t].size ; print ' *** There are '+str(Nt)+' time records...\n'
  25. # Inquire the shape of arrays:
  26. nb_dim = len(vcdim)
  27. print ' *** '+cf_sal+' has '+str(nb_dim)+' dimmensions!'
  28. if not nb_dim in [ 2, 3, 4 ]: print ' ERROR: unsupported number of dimmensions! =>', nb_dim ; sys.exit(0)
  29. # Opening the Netcdf output file:
  30. f_out = Dataset(cf_out, 'r+') # r+ => can read and write in the file... )
  31. print 'File ', cf_out, 'is open...\n'
  32. for jt in range(Nt):
  33. print '\n --- treating record # '+str(jt)
  34. if nb_dim==4: xsal = f_sal.variables[cv_sal][jt,:,:,:]
  35. if nb_dim==3: xsal = f_sal.variables[cv_sal][jt,:,:]
  36. if nb_dim==2: xsal = f_sal.variables[cv_sal][jt,:]
  37. # Extracting tmask at surface level:
  38. if nb_dim==4:
  39. xtemp = f_out.variables[cv_temp][jt,:,:,:]
  40. f_out.variables[cv_temp][jt,:,:,:] = gsw.CT_from_pt(xsal, xtemp)
  41. if nb_dim==3:
  42. xtemp = f_out.variables[cv_temp][jt,:,:]
  43. f_out.variables[cv_temp][jt,:,:] = gsw.CT_from_pt(xsal, xtemp)
  44. if nb_dim==2:
  45. xtemp = f_out.variables[cv_temp][jt,:]
  46. f_out.variables[cv_temp][jt,:] = gsw.CT_from_pt(xsal, xtemp)
  47. f_out.variables[cv_temp].long_name = 'Conservative Temperature (TEOS10) built from potential temperature'
  48. f_sal.close()
  49. f_out.close()
  50. print cf_out+' sucessfully created!'