Changeset 1773 in lmdz_wrf

Timestamp:

Feb 12, 2018, 6:09:39 PM (7 years ago)

Author:

lfita

Message:

Adding `zmla_generic': Computation of boundary layer height following a generic method from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim. applied also in Garcia-Diez et al., 2013, QJRMS

Location:

trunk/tools

Files:

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

trunk/tools/diag_tools.py

@@ -349 +349 @@
     return cape, cin, zlfc, plfc, li, psldims, pslvdims
 
+
+
 def var_clt(cfra):
     """ Function to compute the total cloud fraction following 'newmicro.F90' from 
@@ -867 +869 @@
 
     return psl, psldims, pslvdims
+
+def Forcompute_zmla_gen(theta, qratio, zpl, hgt, dimns, dimvns):
+    """ Function to compute the boundary layer height following a generic method 
+      with Fortran
+    Forcompute_zmla_gen(theta, qratio, zpl, hgt, zmla, dimns, dimvns)
+      [theta]= potential air-temperature values (assuming [[t],z,y,x]) [K]
+      [qratio]= water mixing ratio (assuming [[t],z,y,x]) [kgkg-1]
+      [zpl]= height from sea level (assuming [[t],z,y,x]) [m]
+      [hgt]= topographical height (assuming [m]
+      [zmla]= boundary layer height [m]
+      [dimns]= list of the name of the dimensions of [theta]
+      [dimvns]= list of the name of the variables with the values of the 
+        dimensions of [theta]
+    """
+    fname = 'Forcompute_zmla_gen'
+
+    zmladims = dimns[:]
+    zmlavdims = dimvns[:]
+
+    if len(theta.shape) == 4:
+        zmla= np.zeros((theta.shape[0],theta.shape[2],theta.shape[3]), dtype=np.float)
+
+        dx = theta.shape[3]
+        dy = theta.shape[2]
+        dz = theta.shape[1]
+        dt = theta.shape[0]
+        zmladims.pop(1)
+        zmlavdims.pop(1)
+
+        pzmla= fdin.module_fordiagnostics.compute_zmla_generic4d(                    \
+          tpot=theta[:].transpose(), qratio=qratio[:].transpose(),                   \
+          z=zpl[:].transpose(), hgt=hgt.transpose(), d1=dx, d2=dy, d3=dz, d4=dt)
+        zmla = pzmla.transpose()
+    else:
+        print errormsg
+        print '  ' + fname + ': rank', len(theta.shape), 'not ready !!'
+        print '  it only computes 4D [t,z,y,x] rank values'
+        quit(-1)
+
+    return zmla, zmladims, zmlavdims
+
 
 def compute_OMEGAw(omega, p, t, dimns, dimvns):

trunk/tools/diagnostics.inf

@@ -44 +44 @@
 wssturb, turbulence, wss
 zg, WRFght, PH@PHB
+zmla, WRFzmlagen, T@QVAPOR@WRFgeop@HGT
 

trunk/tools/diagnostics.py

@@ -1104 +1104 @@
         ncvar.insert_variable(ncobj, 'zhgtrel', diagout, dnames, diagoutvd, newnc)
 
+# WRFzmla_gen generic boundary layer hieght computation from WRF theta, QVAPOR, WRFgeop, HGT, 
+    elif diagn == 'WRFzmlagen':
+        var0 = ncobj.variables[depvars[0]][:]+300.
+        var1 = ncobj.variables[depvars[1]][:]
+        dimz = var0.shape[1]
+        var2 = WRFgeop[:,1:dimz+1,:,:]/9.8
+        var3 = ncobj.variables[depvars[3]][0,:,:]
+
+        diagout, diagoutd, diagoutvd = diag.Forcompute_zmla_gen(var0,var1,var2,var3, \
+          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, 'zmla', diagout, diagoutd, diagoutvd, newnc)
+
     else:
         print errormsg

trunk/tools/module_ForDiagnostics.f90

@@ -27 +27 @@
 ! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates
 ! 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
+
 !!!
 ! Calculations
@@ -597 +599 @@
   END SUBROUTINE compute_cape_afwa4D
 
+  SUBROUTINE compute_zmla_generic4D(tpot, qratio, z, hgt, zmla3D, d1, d2, d3, d4)
+! Subroutine to compute pbl-height following a generic method
+!    from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
+!    applied also in Garcia-Diez et al., 2013, QJRMS
+!   where 
+!     "The technique identifies the ML height as a threshold increase of potential temperature from 
+!       its minimum value within the boundary layer."
+!   here applied similarly to Garcia-Diez et al. where 
+!      zmla = "...first level where potential temperature exceeds the minimum potential temperature
+!        reached in the mixed layer by more than 1.5 K"
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: tpot, qratio, z
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
+    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: zmla3D
+ 
+! 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]
+! zmla3D: boundary layer height from surface [m]
+
+    fname = 'compute_zmla_generic4D'
+
+    DO i=1, d1
+      DO j=1, d2
+        DO it=1, d4
+          CALL var_zmla_generic(d3, qratio(i,j,:,it), tpot(i,j,:,it), z(i,j,:,it), hgt(i,j),          &
+            zmla3D(i,j,it))
+        END DO
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_zmla_generic4D
+
+
 END MODULE module_ForDiagnostics

trunk/tools/module_ForDiagnosticsVars.f90

@@ -11 +11 @@
 
   IMPLICIT NONE
-
-  REAL(r_k), PARAMETER                                   :: ZEPSEC=1.0D-12
-! Low limit pressure for medium clouds [Pa]
-  REAL(r_k), PARAMETER                                   :: prmhc = 44000.d0
-! Low limit pressure for High clouds [Pa]
-  REAL(r_k), PARAMETER                                   :: prmlc = 68000.d0
 
   CONTAINS
@@ -34 +28 @@
 ! var_cllmh: low, medium, high-cloud [0,1]
 ! var_clt: total cloudiness [0,1]
+! var_zmla_generic: Subroutine to compute pbl-height following a generic method
 ! VirtualTemp1D: Function for 1D calculation of virtual temperaure
 ! VirtualTemperature: WRF's AFWA method to compute virtual temperature
@@ -967 +962 @@
 ! ---- END modified from module_diag_afwa.F ---- !
 
+  SUBROUTINE var_zmla_generic(dz, qv, tpot, z, topo, zmla)
+!  Subroutine to compute pbl-height following a generic method
+!    from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
+!    applied also in Garcia-Diez et al., 2013, QJRMS
+!   where 
+!     "The technique identifies the ML height as a threshold increase of potential temperature from 
+!       its minimum value within the boundary layer."
+!   here applied similarly to Garcia-Diez et al. where 
+!      zmla = "...first level where potential temperature exceeds the minimum potential temperature
+!        reached in the mixed layer by more than 1.5 K"
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dz
+    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: qv, tpot, z
+    REAL(r_k), INTENT(in)                                :: topo
+    REAL(r_k), INTENT(out)                               :: zmla
+
+! Local
+    INTEGER                                              :: i
+    INTEGER                                              :: mldlev, bllev
+    REAL(r_k)                                            :: dqvar, tpotmin, refdt
+
+!!!!!!! Variables
+! qv: water vapour mixing ratio
+! tpot: potential temperature [K]
+! z: height above sea level [m]
+! topo: topographic height [m]
+! zmla: boundary layer height [m]
+
+    fname = 'var_zmla_generic'
+
+    ! Pecentage of difference of mixing ratio used to determine Mixed layer depth
+    dqvar = 0.1
+
+    ! MLD = Mixed layer with no substantial variation of mixing ratio /\qv < 10% ?
+    !PRINT *,'  Mixed layer mixing ratios qv[1] lev qv[lev] dqvar% _______'
+    DO mldlev = 2, dz
+      IF (ABS(qv(mldlev)-qv(1))/qv(1) > dqvar ) EXIT
+    !  PRINT *,qv(1), mldlev, qv(mldlev), ABS(qv(mldlev)-qv(1))/qv(1)
+    END DO
+
+    ! Looking for the minimum potential temperature within the MLD [tpotmin = min(tpot)_0^MLD]
+    tpotmin = MINVAL(tpot(1:mldlev))
+
+    ! Change in temperature to determine boundary layer height
+    refdt = 1.5
+
+    ! Determine the first level where tpot > tpotmin + 1.5 K
+    !PRINT *,'  Mixed layer tpotmin lev tpotmin[lev] dtpot _______'
+    DO bllev = 1, dz
+      IF (ABS(tpot(bllev)-tpotmin) > refdt ) EXIT
+    !  PRINT *,tpotmin, bllev, tpot(bllev), ABS(tpot(bllev)-tpotmin)
+    END DO
+    
+    !PRINT *,'   height end MLD:', z(mldlev)
+    !PRINT *,'   pbl height:', z(bllev)
+
+    zmla = z(bllev) - topo
+
+    RETURN
+
+  END SUBROUTINE var_zmla_generic
 
 END MODULE module_ForDiagnosticsVars

trunk/tools/module_definitions.f90

@@ -25 +25 @@
 
 ! Scientific constants
+!!
+
+  REAL(r_k), PARAMETER                                   :: ZEPSEC=1.0D-12
+! Low limit pressure for medium clouds [Pa]
+  REAL(r_k), PARAMETER                                   :: prmhc = 44000.d0
+! Low limit pressure for High clouds [Pa]
+  REAL(r_k), PARAMETER                                   :: prmlc = 68000.d0
+
   REAL(r_k), PARAMETER                                   :: grav = 9.81 ! gravity [ms-1]
   REAL(r_k), PARAMETER                                   :: Cpmol = 29.07 ! Molar gas heat capacity at constant pressure [Jmol-1K-1]