Changeset 339 in lmdz_wrf for trunk/tools/validation_sim.py

Timestamp:

Mar 1, 2015, 7:49:18 PM (10 years ago)

Author:

lfita

Message:

Fixing 'WRFws' and 'WRFwd' calculations

File:

• trunk/tools/validation_sim.py (modified) (9 diffs)

trunk/tools/validation_sim.py

@@ -2 +2 @@
 # L. Fita, LMD-Jussieu. February 2015
 ## e.g. sfcEneAvigon # validation_sim.py -d X@west_east@None,Y@south_north@None,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,T@time@time -k single-station -l 4.878773,43.915876,12. -o /home/lluis/DATA/obs/HyMeX/IOP15/sfcEnergyBalance_Avignon/OBSnetcdf.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v HFX@H,LH@LE,GRDFLX@G
-## e.g. AIREP # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@alti,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/AIREP/2012/10/AIREP_121018.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@t
+## e.g. AIREP # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@alti,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/AIREP/2012/10/AIREP_121018.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@t,WRFtd@td,WRFws@u,WRFwd@dd
 ## e.g. ATRCore # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@CFtime -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@altitude,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/ATRCore/V3/ATR_1Hz-HYMEXBDD-SOP1-v3_20121018_as120051.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFT@air_temperature@subc@273.15
 ## e.g. BAMED # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@CFtime -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@altitude,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/BAMED/BAMED_SOP1_B12_TOT5.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFT@tas_north
@@ -577 +577 @@
         'WRFrh': relative humidty fom WRF variables
         'WRFt': temperature from WRF variables
+        'WRFtd': dew-point temperature from WRF variables
+        'WRFws': wind speed from WRF variables
+        'WRFwd': wind direction from WRF variables
         'WRFz': height from WRF variables
     """
@@ -642 +645 @@
                 shape = ncobj.variables['P'].shape
 
-
             elif varn == 'WRFt':
 #        print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
@@ -651 +653 @@
                 dimensions = ncobj.variables['T'].dimensions
                 shape = ncobj.variables['P'].shape
+
+            elif varn == 'WRFtd':
+#        print '    ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...'
+# tacking from: http://en.wikipedia.org/wiki/Dew_point
+                p0=100000.
+                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
+
+                temp = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
+
+                qv = ncobj.variables['QVAPOR'][:]
+
+                tk = temp - 273.15
+                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
+                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
+
+                rh = qv/data2
+                
+                pa = rh * data1/100.
+                varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
+
+                dimensions = ncobj.variables['T'].dimensions
+                shape = ncobj.variables['P'].shape
+
+            elif varn == 'WRFws':
+#        print '    ' + main + ': computing wind speed from WRF as SQRT(U**2 + V**2) ...' 
+                uwind = ncobj.variables['U'][:]
+                vwind = ncobj.variables['V'][:]
+                dx = uwind.shape[3]
+                dy = vwind.shape[2]
+                 
+# de-staggering
+                ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
+                va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
+
+                varNOcheckv = np.sqrt(ua*ua + va*va)
+                dimensions = tuple(['Time','bottom_top','south_north','west_east'])
+                shape = ua.shape
+
+            elif varn == 'WRFwd':
+#        print '    ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...' 
+                uwind = ncobj.variables['U'][:]
+                vwind = ncobj.variables['V'][:]
+                dx = uwind.shape[3]
+                dy = vwind.shape[2]
+
+# de-staggering
+                ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
+                va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
+
+                theta = np.arctan2(ua, va)
+                dimensions = tuple(['Time','bottom_top','south_north','west_east'])
+                shape = ua.shape
+                varNOcheckv = 360.*(1. + theta/(2.*np.pi))
+
+                dimensions = ncobj.variables['U'].dimensions
+                shape = ncobj.variables['U'].shape
 
             elif varn == 'WRFz':
@@ -660 +718 @@
 
             else:
-                print erromsg
+                print errormsg
                 print '  ' + fname + ": variable '" + varn + "' nor ready !!"
                 quit(-1)
@@ -685 +743 @@
   '[constant] to variables values; divc,[constant]: divide by [constant] to ' +      \
   'variables values'
-varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'WRFt', 'WRFz']
+varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'WRFt', 'WRFtd', 'WRFws',        \
+  'WRFwd', 'WRFz']
 varNOcheckinf = "'WRFdens': air density from WRF variables; 'WRFght': geopotential"+ \
   " height from WRF variables; 'WRFp': pressure from WRF variables; 'WRFrh': " +     \
- "relative humidty fom WRF variables; 'WRFt': temperature from WRF variables; " +    \
- "'WRFz': height from WRF variables"
+  "relative humidty fom WRF variables; 'WRFt': temperature from WRF variables; " +   \
+  "'WRFtd': dew-point temperature from WRF variables; 'WRFws': wind speed from " +   \
+  "WRF variables; 'WRFwd': wind speed direction from WRF variables; 'WRFz': " +     \
+  "height from WRF variables"
 
 dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
@@ -1092 +1153 @@
               dims['T'][0]+':'+str(coindtvalues[it][0])
             slicevar, dimslice = slice_variable(ovsim, slicev)
-            simobsvalues.append([ slicevar, ovobs[coindtvalues[it,1]]])
             slicev = dims['X'][0] + ':' + str(trajpos[2,it]-Ngrid) + '@' +           \
               str(trajpos[2,it]+Ngrid) + '|' + dims['Y'][0] + ':' +                  \
@@ -1162 +1222 @@
                 print simobsvalues[varsimobs][:][it]
 
+    arrayvals = np.array(simobsvalues)
+    if len(valvars[ivar]) > 2:
+        const=np.float(valvars[ivar][3])
+        if valvars[ivar][2] == 'sumc':
+            simobsSvalues = simobsSvalues + const
+            arrayvals[:,0] = arrayvals[:,0] + const
+        elif valvars[ivar][2] == 'subc':
+            simobsSvalues = simobsSvalues - const
+            arrayvals[:,0] = arrayvals[:,0] - const
+        elif valvars[ivar][2] == 'mulc':
+            simobsSvalues = simobsSvalues * const
+            arrayvals[:,0] = arrayvals[:,0] * const
+        elif valvars[ivar][2] == 'divc':
+            simobsSvalues = simobsSvalues / const
+            arrayvals[:,0] = arrayvals[:,0] / const
+        else:
+            print errormsg
+            print '  ' + fname + ": operation '" + valvars[ivar][2] + "' not ready!!"
+            quit(-1)
+
 # Statistics around values
     aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
@@ -1178 +1258 @@
       valvars[ivar][1]
     basicvardef(newvar, varsimobs, descvar, ovobs.getncattr('units'))
-    newvar[:] = np.array(simobsvalues)
+    newvar[:] = arrayvals
 
 # Around values
     if dims.has_key('Z'):
-        newdim = onewnc.createDimension('zaround',Ngrid*2+1)
+        if not onewnc.dimensions.has_key('zaround'):
+            newdim = onewnc.createDimension('zaround',Ngrid*2+1)
+
         newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f',            \
           ('time','zaround','yaround','xaround'), fill_value=fillValueF)