Changeset 1776 in lmdz_wrf for trunk

Timestamp:

Feb 15, 2018, 7:18:41 PM (7 years ago)

Author:

lfita

Message:

Adding `zwind': wind at a given height following power-law methodology

Location:

trunk/tools

Files:

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

trunk/tools/diag_tools.py

@@ -910 +910 @@
     return zmla, zmladims, zmlavdims
 
+
+def Forcompute_zwind(ua, va, zsl, uas, vas, hgt, 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]
+      [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.]
+      [zval]= desired height for winds [m]
+      [dimns]= list of the name of the dimensions of [ua]
+      [dimvns]= list of the name of the variables with the values of the 
+        dimensions of [ua]
+    """
+    fname = 'Forcompute_zwind'
+
+    vardims = dimns[:]
+    varvdims = dimvns[:]
+
+    if len(ua.shape) == 4:
+        var1= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)
+        var2= np.zeros((ua.shape[0],ua.shape[2],ua.shape[3]), dtype=np.float)
+
+        dx = ua.shape[3]
+        dy = ua.shape[2]
+        dz = ua.shape[1]
+        dt = ua.shape[0]
+        vardims.pop(1)
+        varvdims.pop(1)
+
+        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)
+        var1 = pvar1.transpose()
+        var2 = pvar2.transpose()
+    else:
+        print errormsg
+        print '  ' + fname + ': rank', len(ua.shape), 'not ready !!'
+        print '  it only computes 4D [t,z,y,x] rank values'
+        quit(-1)
+
+    return var1, var2, vardims, varvdims
 
 def compute_OMEGAw(omega, p, t, dimns, dimvns):

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.
 wa, OMEGAw, vitw@pres@temp
 wds, TSwds, u@v
@@ -45 +46 @@
 zg, WRFght, PH@PHB
 zmla, WRFzmlagen, T@QVAPOR@WRFgeop@HGT
-

trunk/tools/diagnostics.py

@@ -81 +81 @@
   'WRFmrso', 'WRFp',                                                                 \
   'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
-  'WRFheightrel', 'WRFua', 'WRFva']
+  'WRFheightrel', 'WRFua', 'WRFva', 'WRzwind']
 
 methods = ['accum', 'deaccum']
@@ -344 +344 @@
               gen.searchInlist(NONcheckingvars, depv) and                            \
               not gen.searchInlist(methods, depv) and not depvars[0] == 'deaccum'    \
-              and not depvars[0] == 'accum':
+              and not depvars[0] == 'accum' and not depv[0:2] == 'z=':
                 print errormsg
                 print '  ' + main + ": file '" + opts.ncfile +                       \
@@ -1123 +1123 @@
         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]
+    elif diagn == 'WRFzwind':
+        var0 = ncobj.variables[depvars[0]][:]
+        var1 = ncobj.variables[depvars[1]][:]
+        dimz = var0.shape[1]
+        var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
+        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])
+
+        # un-staggering 3D winds
+        unstgdims = [var0.shape[0], var0.shape[1], var0.shape[2], var0.shape[3]-1]
+        va = np.zeros(tuple(unstgdims), dtype=np.float)
+        unvar0 = np.zeros(tuple(unstgdims), dtype=np.float)
+        unvar1 = np.zeros(tuple(unstgdims), dtype=np.float)
+        unvar0 = 0.5*(var0[:,:,:,0:var0.shape[3]-1] + var0[:,:,:,1:var0.shape[3]])
+        unvar1 = 0.5*(var1[:,:,0:var1.shape[2]-1,:] + var1[:,:,1:var1.shape[2],:])
+
+        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwind(unvar0,      \
+          unvar1, var2, var3, var4, var5, var6, var7, var8, dnames, dvnames)
+
+        # Removing the nonChecking variable-dimensions from the initial list
+        varsadd = []
+        for nonvd in NONchkvardims:
+            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
+            varsadd.append(nonvd)
+
+        ncvar.insert_variable(ncobj, 'uaz', diagout1, diagoutd, diagoutvd, newnc)
+        ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
+
     else:
         print errormsg

trunk/tools/module_ForDiagnostics.f90

@@ -28 +28 @@
 ! compute_zint4D: Subroutine to vertically integrate a 4D variable in any vertical coordinates
 ! compute_zmla_generic4D: Subroutine to compute pbl-height following a generic method
+! compute_zwind4D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
 
 !!!
@@ -642 +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 to compute extrapolate the wind at a given height following the 'power law' methodology
+
+    IMPLICIT NONE
+
+    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,d4), INTENT(in)           :: uas, vas
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt, sina, cosa
+    REAL(r_k), INTENT(in)                                :: zextrap
+    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: uaz, vaz
+ 
+! Local
+    INTEGER                                              :: i, j, it
+
+!!!!!!! Variables
+! tpot: potential air temperature [K]
+! qratio: water vapour mixing ratio [kgkg-1]
+! z: height above sea level [m]
+! hgt: terrain height [m]
+! sina, cosa: local sine and cosine of map rotation [1.]
+! zmla3D: boundary layer height from surface [m]
+
+    fname = 'compute_zwind4D'
+
+    DO i=1, d1
+      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))
+        END DO
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_zwind4D
 
 END MODULE module_ForDiagnostics

trunk/tools/module_ForDiagnosticsVars.f90

@@ -29 +29 @@
 ! var_clt: total cloudiness [0,1]
 ! var_zmla_generic: Subroutine to compute pbl-height following a generic method
+! var_zwind: Subroutine to extrapolate the wind at a given height following the 'power law' methodology
 ! VirtualTemp1D: Function for 1D calculation of virtual temperaure
 ! VirtualTemperature: WRF's AFWA method to compute virtual temperature
@@ -1026 +1027 @@
   END SUBROUTINE var_zmla_generic
 
+  SUBROUTINE var_zwind(d1, u, v, z, u10, v10, sa, ca, topo, newz, unewz, vnewz)
+! Subroutine to extrapolate the wind at a given height following the 'power law' methodology
+!    wss[newz] = wss[z1]*(newz/z1)**alpha
+!    alpha = (ln(wss[z2])-ln(wss[z1]))/(ln(z2)-ln(z1))
+! AFTER: Phd Thesis: 
+!   Benedicte Jourdier. Ressource eolienne en France metropolitaine : methodes d’evaluation du 
+!   potentiel, variabilite et tendances. Climatologie. Ecole Doctorale Polytechnique, 2015. French
+!
+    IMPLICIT NONE
+
+    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(out)                               :: unewz, vnewz
+
+! Local
+    INTEGER                                              :: inear
+    REAL(r_k)                                            :: zaground
+    REAL(r_k), DIMENSION(2)                              :: v1, v2, zz, alpha, uvnewz
+
+!!!!!!! Variables
+! u,v: vertical wind components [ms-1]
+! z: height above surface [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
+! unewz,vnewz: Wind compoonents at the given height [ms-1]
+
+    fname = 'var_zwind'
+
+    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)
+        IF ( zaground >= newz) EXIT
+      END DO
+    ELSE
+      PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz'
+      inear = d1 - 2
+    END IF
+
+    IF (inear == d1-2) THEN
+    ! No vertical pair of levels is below newz, using 10m wind as first value and the first level as the second
+       v1(1) = u10
+       v1(2) = v10
+       v2(1) = u(1)
+       v2(2) = v(1)
+       zz(1) = 10.
+       zz(2) = z(1)
+    ELSE
+       v1(1) = u(inear)
+       v1(2) = v(inear)
+       v2(1) = u(inear+1)
+       v2(2) = v(inear+1)
+       zz(1) = z(inear) - topo
+       zz(2) = z(inear+1) - topo
+    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)
+    
+    uvnewz = v1*(newz/zz(1))**alpha
+    ! Earth-rotation
+    unewz = uvnewz(1)*ca - uvnewz(2)*sa
+    vnewz = uvnewz(1)*sa + uvnewz(2)*ca
+
+    PRINT *,'  result vz:', uvnewz
+
+    !STOP
+
+    RETURN
+
+  END SUBROUTINE var_zwind
+
 END MODULE module_ForDiagnosticsVars

trunk/tools/variables_values.dat

@@ -462 +462 @@
 U10, uas, eastward_wind, -30., 30., eastward|10m|wind, ms-1, seismic
 Vent zonal 10m, uas, eastward_wind, -30., 30., eastward|10m|wind, ms-1, seismic
+uaz, uaz, eastward_wind, -30., 30., eastward|wind|at|z|value, ms-1, seismic
 va, va, northward_wind, -30., 30., northward|wind, ms-1, seismic
 vitv, va, northward_wind, -30., 30., northward|wind, ms-1, seismic
@@ -472 +473 @@
 V10, vas, northward_wind, -30., 30., northward|10m|wind, ms-1, seismic
 Vent meridien 10m, vas, northward_wind, -30., 30., northward|10m|wind, ms-1, seismic
+vaz, vaz, northward_wind, -30., 30., northward|wind|at|z|value, ms-1, seismic
 veget, veget, vegetation_type, 1, 13, type|of|plant|functional|type, 1, Reds
 wakedeltaq, wakedeltaq, wake_delta_vapor, -0.003, 0.003, wake|delta|mixing|ratio, kgkg-1s-1, seismic