barentssea: script for NCEP forcing

barentssea/Input/extract_ncep_station.py

+import netCDF4
+from pylab import *
+from netcdftime import utime
+
+years='*'
+path='/work/gg0877/KST/MiMeMo/ncep'
+
+ncs = {'air':netCDF4.MFDataset('%s/air.2m.gauss.%s.nc'%(path,years)),
+       'tcdc':netCDF4.MFDataset('%s/tcdc.eatm.gauss.%s.nc'%(path,years)),
+       'pres':netCDF4.MFDataset('%s/pres.sfc.gauss.%s.nc'%(path,years)),
+       'uwnd':netCDF4.MFDataset('%s/uwnd.10m.gauss.%s.nc'%(path,years)),
+       'vwnd':netCDF4.MFDataset('%s/vwnd.10m.gauss.%s.nc'%(path,years)),
+       'shum':netCDF4.MFDataset('%s/shum.2m.gauss.%s.nc'%(path,years))}
+
+ut = utime(ncs['air'].variables['time'].units)
+dates = ut.num2date(ncs['air'].variables['time'][:])
+
+lonlat = {'barents_sea':(45.0,73.0),'greenland_shelf':(-14.,77.5), 'nansen_basin':(40.0,84.0)}
+
+if len(sys.argv)>3:
+  try:
+    stationname = sys.argv[1]
+    lonlat[stationname] = (float(sys.argv[2]),float(sys.argv[3]))
+  except:
+    print('usage: get_st_profiles.py stationname lon lat')
+    sys.exit()
+elif len(sys.argv)>1:
+  stationname=sys.argv[1]
+  if stationname not in lonlat:
+    print('known station names: %s'%str(lonlat.keys()))
+    sys.exit()
+else:
+  stationname='barents_sea'
+
+# barents sea
+i = 24
+j = 8
+
+f = open('meteo_ncepr2.dat','w')
+
+def getindex(stationname='barents_sea'):
+  import numpy as np
+  lon,lat = lonlat[stationname]
+  nceplon = ncs['air']['lon'][:]
+  nceplat = ncs['air']['lat'][:]
+  i = (np.abs(nceplon-lon)).argmin()
+  j = (np.abs(nceplat-lat)).argmin()
+  return i,j
+
+def getvar(var,i=24,j=8):
+  #if var=='tcdc' or var=='pres':
+  try:
+    return ncs[var].variables[var][:,j,i]
+  except:
+    return ncs[var].variables[var][:,0,j,i]
+
+def rh(sh,t,p):
+  es = 6.112*exp(17.67*t/(t+243.5))
+  e  = sh*p/(0.378*sh+0.622)
+  rh = max(0.0,min(1.0,e/es))
+  return rh
+
+if True:
+  i,j = getindex(stationname)
+  print('extract NCEP for station %s (i,j = %d,%d)'%(stationname,i,j))
+  air = getvar('air',i=i,j=j)-273.15
+  uwnd = getvar('uwnd',i=i,j=j)
+  vwnd = getvar('vwnd',i=i,j=j)
+  tcdc = getvar('tcdc',i=i,j=j)/100.
+  shum = getvar('shum',i=i,j=j)
+  pres = getvar('pres',i=i,j=j)/100.
+
+
+for n,date in enumerate(dates):
+  f.write('%s %6.2f %6.2f %7.2f %6.2f %7.4f %5.2f\n'%(str(date),uwnd[n],vwnd[n],pres[n],air[n],rh(shum[n],air[n],pres[n])*100.,tcdc[n]))
+
+f.close()