Changeset 1784 in lmdz_wrf for trunk

Timestamp:

Mar 1, 2018, 12:48:27 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `WRFzwind_log': extrapolate the wind at a given height following the 'logarithmic law' methodology

Location:

trunk/tools

Files:

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

trunk/tools/diag_tools.py

@@ -34 +34 @@
 # Forcompute_zmla_gen: Function to compute the boundary layer height following a generic method with Fortran
 # Forcompute_zwind: Function to compute the wind at a given height following the power law method
+# Forcompute_zwind_log: Function to compute the wind at a given height following the logarithmic law method
 # Forcompute_zwindMO: Function to compute the wind at a given height following the Monin-Obukhov theory
 # W_diagnostic: Class to compute WRF diagnostics variables
@@ -948 +949 @@
           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()
+    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 Forcompute_zwind_log(ua, va, z, uas, vas, sina, cosa, zval, dimns, dimvns):
+    """ Function to compute the wind at a given height following the logarithmic 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]
+      [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]
+      [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_log'
+
+    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_zwind_log4d(                \
+          ua=ua.transpose(), 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

@@ -35 +35 @@
 va, WRFva, U@V@SINALPHA@COSALPHA
 uavaz, WRFzwind, U@V@WRFz@U10@V10@SINALPHA@COSALPHA@z=100.
+uavaz, WRFzwind_log, U@V@WRFz@U10@V10@SINALPHA@COSALPHA@z=100.
 uavaz, WRFzwindMO, UST@ZNT@RMOL@U10@V10@SINALPHA@COSALPHA@z=100.
 wa, OMEGAw, vitw@pres@temp

trunk/tools/diagnostics.py

@@ -81 +81 @@
   'WRFmrso', 'WRFp',                                                                 \
   'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
-  'WRFheightrel', 'WRFua', 'WRFva', 'WRFzwind', 'WRFzwindMO']
+  'WRFheightrel', 'WRFua', 'WRFva', 'WRFzwind', 'WRFzwind_log', 'WRFzwindMO']
 
 methods = ['accum', 'deaccum']
@@ -1146 +1146 @@
         ncvar.insert_variable(ncobj, 'zmla', diagout, diagoutd, diagoutvd, newnc)
 
-# WRFzwind wind extrapolation at a given hieght computation from WRF U, V, WRFz, 
-#   U10, V10, SINALPHA, COSALPHA, z=[zval]
+# WRFzwind wind extrapolation at a given height using power law computation from WRF 
+#   U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
     elif diagn == 'WRFzwind':
         var0 = ncobj.variables[depvars[0]][:]
@@ -1178 +1178 @@
         ncvar.insert_variable(ncobj, 'vaz', diagout2, diagoutd, diagoutvd, newnc)
 
+# WRFzwind wind extrapolation at a given hieght using logarithmic law computation 
+#   from WRF U, V, WRFz, U10, V10, SINALPHA, COSALPHA, z=[zval]
+    elif diagn == 'WRFzwind_log':
+        var0 = ncobj.variables[depvars[0]][:]
+        var1 = ncobj.variables[depvars[1]][:]
+        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 = np.float(depvars[7].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_log(unvar0,  \
+          unvar1, var2, var3, var4, var5, var6, var7, 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)
+
 # WRFzwindMO wind extrapolation at a given height computation using Monin-Obukhow 
 #   theory from WRF UST, ZNT, RMOL, U10, V10, SINALPHA, COSALPHA, z=[zval]
+#   NOTE: only useful for [zval] < 80. m
     elif diagn == 'WRFzwindMO':
         var0 = ncobj.variables[depvars[0]][:]

trunk/tools/module_ForDiagnostics.f90

@@ -29 +29 @@
 ! 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
+! compute_zwind_log4D: Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
 ! compute_zwindMCO3D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodolog
 
@@ -681 +682 @@
   END SUBROUTINE compute_zwind4D
 
+  SUBROUTINE compute_zwind_log4D(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 'logarithmic law' methodology
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
+    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)              :: 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 surface [m]
+! sina, cosa: local sine and cosine of map rotation [1.]
+! zmla3D: boundary layer height from surface [m]
+
+    fname = 'compute_zwind_log4D'
+
+    DO i=1, d1
+      DO j=1, d2
+        DO it=1, d4
+          CALL var_zwind_log(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
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_zwind_log4D
+
   SUBROUTINE compute_zwindMO3D(d1, d2, d3, ust, znt, rmol, uas, vas, sina, cosa, newz, uznew, vznew)
 ! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
+!   NOTE: only usefull for newz < 80. m
 
     IMPLICIT NONE

trunk/tools/module_ForDiagnosticsVars.f90

@@ -30 +30 @@
 ! 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
+! var_zwind_log: Subroutine to extrapolate the wind at a given height following the 'logarithmic law' methodology
 ! var_zwind_MOtheor: Subroutine of wind extrapolation following Moin-Obukhov theory
 ! VirtualTemp1D: Function for 1D calculation of virtual temperaure
@@ -1111 +1112 @@
   END SUBROUTINE var_zwind
 
+  SUBROUTINE var_zwind_log(d1, u, v, z, u10, v10, sa, ca, newz, unewz, vnewz)
+! Subroutine to extrapolate the wind at a given height following the 'logarithmic 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, 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, logz0, uvnewz
+
+!!!!!!! Variables
+! u,v: vertical wind components [ms-1]
+! z: height above surface on half-mass levels [m]
+! u10,v10: 10-m wind components [ms-1]
+! 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_log'
+
+    !PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______'
+    IF (z(1) < newz ) THEN
+      DO inear = 1,d1-2
+        ! L. Fita, CIMA. Feb. 2018
+        !! Choose between extra/inter-polate. Maybe better interpolate?
+        ! Here we extrapolate from two closest lower levels
+        !zaground = z(inear+2)
+        zaground = z(inear+1)
+        !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)
+       zz(2) = z(inear+1)
+    END IF
+
+    ! Computing for each component
+    logz0 = (v2*LOG(zz(1))-v1*LOG(zz(2)))/(v2-v1)
+    
+    uvnewz = v2*(LOG(newz)-logz0)/(LOG(zz(2))-logz0)
+    ! 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_log
+
   SUBROUTINE var_zwind_MOtheor(ust, znt, rmol, u10, v10, sa, ca, newz, uznew, vznew)
   ! Subroutine of wind extrapolation following Moin-Obukhov theory R. B. Stull, 1988, 
   !   Springer (p376-383)
+  ! Only usefull for newz < 80. m
 
     IMPLICIT NONE