Changeset 1777 in lmdz_wrf for trunk

Timestamp:

Feb 16, 2018, 3:17:02 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `WRFz': as height above surface unstaggering WRFgeop
• Fixing zwind computation

Location:

trunk/tools

Files:

• diag_tools.py (modified) (2 diffs)
• diagnostics.inf (modified) (1 diff)
• diagnostics.py (modified) (9 diffs)
• module_ForDiagnostics.f90 (modified) (4 diffs)
• module_ForDiagnosticsVars.f90 (modified) (6 diffs)
• variables_values.dat (modified) (1 diff)

trunk/tools/diag_tools.py

@@ -911 +911 @@
 
 
-def Forcompute_zwind(ua, va, zsl, uas, vas, hgt, sina, cosa, zval, dimns, dimvns):
+def Forcompute_zwind(ua, va, z, uas, vas, sina, cosa, zval, dimns, dimvns):
     """ Function to compute the wind at a given height following the power law method
     Forcompute_zwind(ua, va, zsl, uas, vas, hgt, sina, cosa, zval, dimns, dimvns)
       [ua]= x-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
       [va]= y-component of unstaggered 3D wind (assuming [[t],z,y,x]) [ms-1]
-      [zsl]= height above sea level [m]
+      [z]= height above surface [m]
       [uas]= x-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
       [vas]= y-component of unstaggered 10 m wind (assuming [[t],z,y,x]) [ms-1]
-      [hgt]= topographical height (assuming [m]
       [sina]= local sine of map rotation [1.]
       [cosa]= local cosine of map rotation [1.]
@@ -944 +943 @@
 
         pvar1, pvar2= fdin.module_fordiagnostics.compute_zwind4d(ua=ua.transpose(),  \
-          va=va[:].transpose(), zsl=zsl[:].transpose(), uas=uas.transpose(),         \
-          vas=vas.transpose(), hgt=hgt.transpose(), sina=sina.transpose(),           \
-          cosa=cosa.transpose(), zextrap=zval, d1=dx, d2=dy, d3=dz, d4=dt)
+          va=va[:].transpose(), z=z[:].transpose(), uas=uas.transpose(),             \
+          vas=vas.transpose(), sina=sina.transpose(), cosa=cosa.transpose(),         \
+          zextrap=zval, d1=dx, d2=dy, d3=dz, d4=dt)
         var1 = pvar1.transpose()
         var2 = pvar2.transpose()

trunk/tools/diagnostics.inf

@@ -34 +34 @@
 ua, WRFua, U@V@SINALPHA@COSALPHA
 va, WRFva, U@V@SINALPHA@COSALPHA
-uavaz, WRFzwind, U@V@WRFgeop@U10@V10@HGT@SINALPHA@COSALPHA@z=100.
+uavaz, WRFzwind, U@V@WRFz@U10@V10@SINALPHA@COSALPHA@z=100.
 wa, OMEGAw, vitw@pres@temp
 wds, TSwds, u@v

trunk/tools/diagnostics.py

@@ -6 +6 @@
 #      ciclad: f2py --f90flags="-fPIC" --f90exec=/usr/bin/gfortran -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -L/usr/lib64/ -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -m module_ForDiagnostics -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90 >& run_f2py.log
 
-## e.g. # diagnostics.py -d 'Time@time,bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -v 'clt|CLDFRA,cllmh|CLDFRA@WRFp,RAINTOT|RAINC@RAINNC@XTIME' -f WRF_LMDZ/NPv31/wrfout_d01_1980-03-01_00:00:00
+## e.g. # diagnostics.py -d 'Time@WRFtime,bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -v 'clt|CLDFRA,cllmh|CLDFRA@WRFp,RAINTOT|RAINC@RAINNC@XTIME' -f WRF_LMDZ/NPv31/wrfout_d01_1980-03-01_00:00:00
 ## e.g. # diagnostics.py -f /home/lluis/PY/diagnostics.inf -d variable_combo -v WRFprc
 
@@ -90 +90 @@
   'WRFp', 'WRFtd',                                                                   \
   'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs', 'WRFrhs', 'WRFrvors', \
-  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight']
+  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz']
 
 NONchkvardims = ['WRFtime']
@@ -124 +124 @@
 WRFpos_compute = False
 WRFtime_compute = False
+WRFz_compute = False
 
 # File creation
@@ -147 +148 @@
     if dnv == 'WRFpos': WRFpos_compute = True
     if dnv == 'WRFtime': WRFtime_compute = True
+    if dnv == 'WRFz':WRFz_compute = True
 
 # diagnostics to compute
@@ -163 +165 @@
         if depvars == 'WRFpos': WRFpos_compute = True
         if depvars == 'WRFtime': WRFtime_compute = True
+        if depvars == 'WRFz': WRFz_compute = True
     else:
         depvars = diags[idiag].split('|')[1].split('@')
@@ -173 +176 @@
         if gen.searchInlist(depvars, 'WRFpos'): WRFpos_compute = True
         if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
+        if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True
 
 # Dictionary with the new computed variables to be able to add them
@@ -325 +329 @@
     dictcompvars['WRFtime'] = {'name': 'time', 'standard_name': 'time',              \
       'long_name': 'time', 'units': tunits, 'calendar': 'gregorian'}
+
+if WRFz_compute:
+    print '  ' + main + ': Retrieving z: height above surface value from WRF as ' +  \
+      'unstagger(PH + PHB)/9.8-hgt'
+    dimv = ncobj.variables['PH'].shape
+    WRFzg = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/9.8
+
+    unzgd = (dimv[0], dimv[1]-1, dimv[2], dimv[3])
+    unzg = np.zeros(unzgd, dtype=np.float)
+    unzg = 0.5*(WRFzg[:,0:dimv[1]-1,:,:] + WRFzg[:,1:dimv[1],:,:])
+
+    WRFz = np.zeros(unzgd, dtype=np.float)
+    for iz in range(dimv[1]-1):
+        WRFz[:,iz,:,:] = unzg[:,iz,:,:] - ncobj.variables['HGT'][:]
+
+    # Attributes of the variable
+    Vvals = gen.variables_values('WRFz')
+    dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
+      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
 
 ### ## #
@@ -1123 +1146 @@
         ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
 
-# WRFzwind wind extrapolation at a given hieght computation from WRF U, V, WRFgeop, 
-#   U10, V10, HGT, SINALPHA, COSALPHA, z=[zval]
+# WRFzwind wind extrapolation at a given hieght computation from WRF U, V, WRFz, 
+#   U10, V10, SINALPHA, COSALPHA, z=[zval]
     elif diagn == 'WRFzwind':
         var0 = ncobj.variables[depvars[0]][:]
         var1 = ncobj.variables[depvars[1]][:]
-        dimz = var0.shape[1]
-        var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
+        var2 = WRFz
         var3 = ncobj.variables[depvars[3]][:]
         var4 = ncobj.variables[depvars[4]][:]
         var5 = ncobj.variables[depvars[5]][0,:,:]
         var6 = ncobj.variables[depvars[6]][0,:,:]
-        var7 = ncobj.variables[depvars[7]][0,:,:]
-        var8 = np.float(depvars[8].split('=')[1])
+        var7 = np.float(depvars[7].split('=')[1])
 
         # un-staggering 3D winds
@@ -1146 +1167 @@
 
         diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0,      \
-          unvar1, var2, var3, var4, var5, var6, var7, var8, dnames, dvnames)
+          unvar1, var2, var3, var4, var5, var6, var7, dnames, dvnames)
 
         # Removing the nonChecking variable-dimensions from the initial list

trunk/tools/module_ForDiagnostics.f90

@@ -643 +643 @@
   END SUBROUTINE compute_zmla_generic4D
 
-  SUBROUTINE compute_zwind4D(ua, va, zsl, uas, vas, hgt, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
+  SUBROUTINE compute_zwind4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
 ! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
 
@@ -649 +649 @@
 
     INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
-    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ua, va, zsl
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ua, va, z
     REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: uas, vas
-    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt, sina, cosa
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: sina, cosa
     REAL(r_k), INTENT(in)                                :: zextrap
     REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: uaz, vaz
@@ -661 +661 @@
 ! tpot: potential air temperature [K]
 ! qratio: water vapour mixing ratio [kgkg-1]
-! z: height above sea level [m]
-! hgt: terrain height [m]
+! z: height above surface [m]
 ! sina, cosa: local sine and cosine of map rotation [1.]
 ! zmla3D: boundary layer height from surface [m]
@@ -671 +670 @@
       DO j=1, d2
         DO it=1, d4
-          CALL var_zwind(d3, ua(i,j,:,it), va(i,j,:,it), zsl(i,j,:,it), uas(i,j,it), vas(i,j,it),    &
-            hgt(i,j), sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
+          CALL var_zwind(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it),       &
+            sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
         END DO
       END DO

trunk/tools/module_ForDiagnosticsVars.f90

@@ -1027 +1027 @@
   END SUBROUTINE var_zmla_generic
 
-  SUBROUTINE var_zwind(d1, u, v, z, u10, v10, sa, ca, topo, newz, unewz, vnewz)
+  SUBROUTINE var_zwind(d1, u, v, z, u10, v10, sa, ca, newz, unewz, vnewz)
 ! Subroutine to extrapolate the wind at a given height following the 'power law' methodology
 !    wss[newz] = wss[z1]*(newz/z1)**alpha
@@ -1039 +1039 @@
     INTEGER, INTENT(in)                                  :: d1
     REAL(r_k), DIMENSION(d1), INTENT(in)                 :: u,v,z
-    REAL(r_k), INTENT(in)                                :: u10, v10, topo, sa, ca, newz
+    REAL(r_k), INTENT(in)                                :: u10, v10, sa, ca, newz
     REAL(r_k), INTENT(out)                               :: unewz, vnewz
 
@@ -1049 +1049 @@
 !!!!!!! Variables
 ! u,v: vertical wind components [ms-1]
-! z: height above surface [m]
+! z: height above surface on half-mass levels [m]
 ! u10,v10: 10-m wind components [ms-1]
-! topo: topographical height [m]
 ! sa, ca: local sine and cosine of map rotation [1.]
 ! newz: desired height above grpund of extrapolation
@@ -1058 +1057 @@
     fname = 'var_zwind'
 
-    PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______'
+    !PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______'
     IF (z(1) < newz ) THEN
       DO inear = 1,d1-2
         zaground = z(inear+2)
-        PRINT *, inear, z(inear), newz, z(inear+1), z(inear+2)
+        !PRINT *, inear, z(inear), newz, z(inear+1), z(inear+2)
         IF ( zaground >= newz) EXIT
       END DO
     ELSE
-      PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz'
+      !PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz'
       inear = d1 - 2
     END IF
@@ -1083 +1082 @@
        v2(1) = u(inear+1)
        v2(2) = v(inear+1)
-       zz(1) = z(inear) - topo
-       zz(2) = z(inear+1) - topo
+       zz(1) = z(inear)
+       zz(2) = z(inear+1)
     END IF
 
     ! Computing for each component
     alpha = (LOG(ABS(v2))-LOG(ABS(v1)))/(LOG(zz(2))-LOG(zz(1)))
-    PRINT *,' Computing with v1:', v1, ' ms-1 v2:', v2, ' ms-1'
-    PRINT *,' z1:', zz(1), 'm z2:', zz(2), ' m'
-    PRINT *,' alhpa u:', alpha(1), ' alpha 2:', alpha(2)
+    !PRINT *,' Computing with v1:', v1, ' ms-1 v2:', v2, ' ms-1'
+    !PRINT *,' z1:', zz(1), 'm z2:', zz(2), ' m'
+    !PRINT *,' alhpa u:', alpha(1), ' alpha 2:', alpha(2)
     
     uvnewz = v1*(newz/zz(1))**alpha
@@ -1098 +1097 @@
     vnewz = uvnewz(1)*sa + uvnewz(2)*ca
 
-    PRINT *,'  result vz:', uvnewz
+    !PRINT *,'  result vz:', uvnewz
 
     !STOP

trunk/tools/variables_values.dat

@@ -643 +643 @@
 Y, y, y,  0., 100., y, -, Blues
 z, z, z,  0., 100., z, -, Greens
+WRFz, z, height, 0., 80000., height|above|surface, m, rainbow
 Z, z, z,  0., 100., z, -, Greens
 zhgt, zhgt, height,  0., 50000., sea|level|altitude, m, rainbow