Changeset 1783 in lmdz_wrf

Timestamp:

Feb 28, 2018, 3:26:43 AM (7 years ago)

Author:

lfita

Message:

Adding:

• `WRFzwindMO': wind extrapolation at a given height computation using Monin-Obukhow

Location:

trunk/tools

Files:

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

trunk/tools/diag_tools.py

@@ -29 +29 @@
 # compute_WRFuasvas: Fucntion to compute geographical rotated WRF 2-meter winds
 # derivate_centered: Function to compute the centered derivate of a given field
-# def Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
+# Forcompute_cllmh: Function to compute cllmh: low/medium/hight cloud fraction following newmicro.F90 from LMDZ via Fortran subroutine
 # Forcompute_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ via a Fortran module
+# Forcompute_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F'
+# 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_zwindMO: Function to compute the wind at a given height following the Monin-Obukhov theory
 # W_diagnostic: Class to compute WRF diagnostics variables
-# Forcompute_psl_ptarget: Function to compute the sea-level pressure following target_pressure value found in `p_interp.F'
 
 # Others just providing variable values
@@ -910 +913 @@
     return zmla, zmladims, zmlavdims
 
-
 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
@@ -952 +954 @@
         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_zwindMO(ust, znt, rmol, uas, vas, sina, cosa, zval, dimns, dimvns):
+    """ Function to compute the wind at a given height following the Monin-Obukhov theory
+    Forcompute_zwind(ust, znt, rmol, uas, vas, sina, cosa, zval, dimns, dimvns)
+      [ust]: u* in similarity theory (assuming [[t],y,x]) [ms-1]
+      [znt]: thermal time-varying roughness length (assuming [[t],y,x]) [m]
+      [rmol]: inverse of Obukhov length (assuming [[t],y,x]) [m-1]
+      [uas]= x-component of unstaggered 10 m wind (assuming [[t],y,x]) [ms-1]
+      [vas]= y-component of unstaggered 10 m wind (assuming [[t],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 [uas]
+      [dimvns]= list of the name of the variables with the values of the 
+        dimensions of [uas]
+    """
+    fname = 'Forcompute_zwindMO'
+
+    vardims = dimns[:]
+    varvdims = dimvns[:]
+
+    if len(uas.shape) == 3:
+        var1= np.zeros((uas.shape[0],uas.shape[1],uas.shape[2]), dtype=np.float)
+        var2= np.zeros((uas.shape[0],uas.shape[1],uas.shape[2]), dtype=np.float)
+
+        dx = uas.shape[2]
+        dy = uas.shape[1]
+        dt = uas.shape[0]
+
+        pvar1, pvar2 = fdin.module_fordiagnostics.compute_zwindmo3d(                 \
+          ust=ust.transpose(), znt=znt[:].transpose(), rmol=rmol[:].transpose(),     \
+          uas=uas.transpose(), vas=vas.transpose(), sina=sina.transpose(),           \
+          cosa=cosa.transpose(), newz=zval, d1=dx, d2=dy, d3=dt)
+        var1 = pvar1.transpose()
+        var2 = pvar2.transpose()
+    else:
+        print errormsg
+        print '  ' + fname + ': rank', len(uas.shape), 'not ready !!'
+        print '  it only computes 3D [t,y,x] rank values'
         quit(-1)
 

trunk/tools/diagnostics.inf

@@ -35 +35 @@
 va, WRFva, U@V@SINALPHA@COSALPHA
 uavaz, WRFzwind, 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
 wds, TSwds, u@v

trunk/tools/diagnostics.py

@@ -81 +81 @@
   'WRFmrso', 'WRFp',                                                                 \
   'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
-  'WRFheightrel', 'WRFua', 'WRFva', 'WRzwind']
+  'WRFheightrel', 'WRFua', 'WRFva', 'WRFzwind', 'WRFzwindMO']
 
 methods = ['accum', 'deaccum']
@@ -1178 +1178 @@
         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]
+    elif diagn == 'WRFzwindMO':
+        var0 = ncobj.variables[depvars[0]][:]
+        var1 = ncobj.variables[depvars[1]][:]
+        var2 = ncobj.variables[depvars[2]][:]
+        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])
+
+        diagout1, diagout2, diagoutd, diagoutvd = diag.Forcompute_zwindMO(var0, var1,\
+          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)
+
     else:
         print errormsg

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_zwindMCO3D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodolog
 
 !!!
@@ -680 +681 @@
   END SUBROUTINE compute_zwind4D
 
+  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
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, d2, d3
+    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: ust, znt, rmol
+    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: uas, vas
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: sina, cosa
+    REAL(r_k), INTENT(in)                                :: newz
+    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: uznew, vznew
+ 
+! Local
+    INTEGER                                              :: i, j, it
+
+!!!!!!! Variables
+! ust: u* in similarity theory [ms-1]
+! znt: thermal time-varying roughness length [m]
+! rmol: Inverse of the Obukhov length [m-1]
+! uas: x-component 10-m wind speed [ms-1]
+! vas: y-component 10-m wind speed [ms-1]
+! sina, cosa: local sine and cosine of map rotation [1.]
+
+    fname = 'compute_zwindMO3D'
+
+    DO i=1, d1
+      DO j=1, d2
+        DO it=1, d3
+          CALL var_zwind_MOtheor(ust(i,j,it), znt(i,j,it), rmol(i,j,it), uas(i,j,it), vas(i,j,it),    &
+            sina(i,j), cosa(i,j), newz, uznew(i,j,it), vznew(i,j,it))
+        END DO
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_zwindMO3D
+
 END MODULE module_ForDiagnostics

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_MOtheor: Subroutine of wind extrapolation following Moin-Obukhov theory
 ! VirtualTemp1D: Function for 1D calculation of virtual temperaure
 ! VirtualTemperature: WRF's AFWA method to compute virtual temperature
+! stabfunc_businger: Fucntion of the stability function after Businger et al. (1971)
 
 !!!!!!! Calculations
@@ -1109 +1111 @@
   END SUBROUTINE var_zwind
 
+  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)
+
+    IMPLICIT NONE
+
+    REAL, INTENT(in)                                     :: ust, znt, rmol, u10, v10, sa, ca
+    REAL, INTENT(in)                                     :: newz
+    REAL, INTENT(out)                                    :: uznew, vznew
+
+! Local
+    REAL                                                 :: OL
+    REAL                                                 :: stability
+    REAL                                                 :: wsz, alpha
+    REAL, DIMENSION(2)                                   :: uvnewz
+
+!!!!!!! Variables
+! ust: u* in similarity theory [ms-1]
+! z0: roughness length [m]
+!!! L. Fita, CIMA. Feb. 2018
+!! NOT SURE if it should be z0 instead?
+! znt: thermal time-varying roughness length [m]
+! rmol: inverse of Obukhov length [m-1]
+! u10: x-component 10-m wind speed [ms-1]
+! v10: y-component 10-m wind speed [ms-1]
+! sa, ca: local sine and cosine of map rotation [1.]
+! 
+    fname = 'var_zwind_MOtheor'
+
+    ! Obukhov Length (using the Boussinesq approximation giving Tv from t2)
+    OL = 1/rmol
+
+    ! Wind speed at desired height
+    PRINT *,'ust:', ust, 'znt:', znt, 'OL:', OL
+
+    CALL stabfunc_businger(newz,OL,stability)
+    PRINT *,'  z/L:', newz/OL,' stabfunc:', stability, 'log:', LOG(newz/znt), 'log+stability:', LOG(newz/znt) + stability
+    PRINT *,'  ust/karman:', ust/karman
+
+    wsz = ust/karman*( LOG(newz/znt) + stability)
+    PRINT *,'  wsz:', wsz
+
+    ! Without taking into account Ekcman pumping, etc... redistributed by components unsing 10-m wind
+    !   as reference...
+    alpha = ATAN2(v10,u10)
+    uvnewz(1) = wsz*COS(alpha)
+    uvnewz(2) = wsz*SIN(alpha)
+    PRINT *,'  uvnewz:', uvnewz
+
+    ! Earth-rotation
+    uznew = uvnewz(1)*ca - uvnewz(2)*sa
+    vznew = uvnewz(1)*sa + uvnewz(2)*ca
+    PRINT *,'  uznew:', uznew, ' vznew:', vznew
+
+    RETURN
+
+  END SUBROUTINE var_zwind_MOtheor
+
+  ! L. Fita, CIMA. Feb. 2018
+  ! WRF seems to have problems with my functions, let'suse subroutine instead
+  !REAL FUNCTION stabfunc_businger(z,L)
+  SUBROUTINE stabfunc_businger(z,L,stabfunc_busingerv)
+  ! Fucntion of the stability function after Businger et al. (1971), JAS, 28(2), 181–189
+
+    IMPLICIT NONE
+
+    REAL, INTENT(in)                                     :: z,L
+    REAL, INTENT(out)                                    :: stabfunc_busingerv
+
+! Local
+    REAL                                                 :: zL, X
+
+!!!!!!! Variables
+! z: height [m]
+! L: Obukhov length [-]
+
+    fname = 'stabfunc_businger'
+
+    IF (L /= 0.) THEN
+      zL = z/L
+    ELSE
+      ! Neutral
+      zL = 0.
+    END IF
+
+    IF (zL > 0.) THEN
+    ! Stable case
+      stabfunc_busingerv = 4.7*z/L
+    ELSE IF (zL < 0.) THEN
+    ! unstable
+      X = (1. - 15.*z/L)**(0.25)
+      !stabfunc_busingerv = -2.*LOG((1.+X)/2.)-LOG((1.+X**2)/2.)+2.*ATAN(X)-piconst/2.
+      stabfunc_busingerv = LOG( ((1.+X**2)/2.)*((1.+X)/2.)**2)-2.*ATAN(X)+piconst/2.
+    ELSE
+      stabfunc_busingerv = 0.
+    END IF
+
+    RETURN
+
+!  END FUNCTION stabfunc_businger
+  END SUBROUTINE stabfunc_businger
+
 END MODULE module_ForDiagnosticsVars

trunk/tools/module_definitions.f90

@@ -57 +57 @@
 ! epislon r_d/r_v for Virtual temperature
   REAL(r_k), PARAMETER                                   :: epsilonv = 0.622
+! pi
+  REAL(r_k), PARAMETER                                   :: piconst = 3.1415926535897932384626433
+! Karman constant
+  REAL(r_k), PARAMETER                                   :: karman = 0.4
+
 
 END MODULE module_definitions