Changeset 1759 in lmdz_wrf for trunk/tools


Ignore:
Timestamp:
Jan 18, 2018, 10:51:03 PM (8 years ago)
Author:
lfita
Message:

Adding:

  • CAPE, CIN, ZLFC, PLFC, LI from WRF phys/module_diag_afwa.F
Location:
trunk/tools
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/diag_tools.py

    r1758 r1759  
    300300# Diagnostics
    301301##
     302def Forcompute_cape_afwa(ta, hur, pa, zg, hgt, parcelm, dimns, dimvns):
     303    """ Function to compute the CAPE, CIN, ZLFC, PLFC, LI following WRF
     304      'phys/module_diaf_afwa.F' methodology
     305    Forcompute_cape_afwa(ta, hur, pa, hgt, zsfc, parcelm, dimns, dimvns)
     306      [ta]= air-temperature values (assuming [[t],z,y,x]) [K]
     307      [pa]= pressure values (assuming [[t],z,y,x]) [Pa]
     308      [zg]= gopotential height (assuming [[t],z,y,x]) [gpm]
     309      [hgt]= topographical height (assuming [m]
     310      [parcelm]= method of air-parcel to use
     311      [dimns]= list of the name of the dimensions of [pa]
     312      [dimvns]= list of the name of the variables with the values of the
     313        dimensions of [pa]
     314    """
     315    fname = 'Forcompute_cape_afwa'
     316
     317    psldims = dimns[:]
     318    pslvdims = dimvns[:]
     319
     320    if len(pa.shape) == 4:
     321        cape = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
     322        cin = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
     323        zlfc = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
     324        plfc = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
     325        li = np.zeros((pa.shape[0],pa.shape[2],pa.shape[3]), dtype=np.float)
     326
     327        dx = pa.shape[3]
     328        dy = pa.shape[2]
     329        dz = pa.shape[1]
     330        dt = pa.shape[0]
     331        psldims.pop(1)
     332        pslvdims.pop(1)
     333
     334        pcape,pcin,pzlfc,pplfc,pli= fdin.module_fordiagnostics.compute_cape_afwa4d(  \
     335          ta=ta[:].transpose(), hur=hur[:].transpose(), press=pa[:].transpose(),     \
     336          zg=zg[:].transpose(), hgt=hgt.transpose(), parcelmethod=parcelm, d1=dx,    \
     337          d2=dy, d3=dz, d4=dt)
     338        cape = pcape.transpose()
     339        cin = pcin.transpose()
     340        zlfc = pzlfc.transpose()
     341        plfc = pplfc.transpose()
     342        li = pli.transpose()
     343    else:
     344        print errormsg
     345        print '  ' + fname + ': rank', len(pa.shape), 'not ready !!'
     346        print '  it only computes 4D [t,z,y,x] rank values'
     347        quit(-1)
     348
     349    return cape, cin, zlfc, plfc, li, psldims, pslvdims
    302350
    303351def var_clt(cfra):

  • trunk/tools/diagnostics.inf

    r1758 r1759  
    44
    55bils, WRFbils, HFX@LH
     6cape, WRFcape_afwa, WRFt@WRFrh@WRFp@WRFgeop
    67clt, clt, CLDFRA
    78cll, cllmh, CLDFRA@WRFp

  • trunk/tools/diagnostics.py

    r1758 r1759  
    690690        ncvar.insert_variable(ncobj, 'bils', diagout, dnamesvar, dvnamesvar, newnc)
    691691
     692# WRFcape_afwa CAPE, CIN, ZLFC, PLFC, LI following WRF 'phys/module_diaf_afwa.F'
     693#   methodology as WRFt, WRFrh, WRFp, WRFgeop, HGT
     694    elif diagn == 'WRFcape_afwa':
     695        var0 = WRFt
     696        var1 = WRFrh
     697        var2 = WRFp
     698        dz = WRFgeop.shape[1]
     699        # de-staggering
     700        var3 = 0.5*(WRFgeop[:,0:dz-1,:,:]+WRFgeop[:,1:dz,:,:])*9.81
     701        var4 = ncobj.variables[depvars[4]][0,:,:]
     702
     703        dnamesvar = list(ncobj.variables['T'].dimensions)
     704        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
     705
     706        diagout = np.zeros(var0.shape, dtype=np.float)
     707        diagout1, diagout2, diagout3, diagout4, diagout5, diagoutd, diagoutvd =      \
     708          diag.Forcompute_cape_afwa(var0, var1, var2, var3, var4, 3, dnamesvar,      \
     709          dvnamesvar)
     710
     711        # Removing the nonChecking variable-dimensions from the initial list
     712        varsadd = []
     713        for nonvd in NONchkvardims:
     714            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
     715            varsadd.append(nonvd)
     716
     717        ncvar.insert_variable(ncobj, 'cape', diagout1, diagoutd, diagoutvd, newnc)
     718        ncvar.insert_variable(ncobj, 'cin', diagout2, diagoutd, diagoutvd, newnc)
     719        ncvar.insert_variable(ncobj, 'zlfc', diagout3, diagoutd, diagoutvd, newnc)
     720        ncvar.insert_variable(ncobj, 'plfc', diagout4, diagoutd, diagoutvd, newnc)
     721        ncvar.insert_variable(ncobj, 'li', diagout5, diagoutd, diagoutvd, newnc)
     722
    692723# WRFclivi WRF water vapour path WRFdens, QICE, QGRAUPEL, QHAIL
    693724    elif diagn == 'WRFclivi':
    … …  
    775806        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
    776807
    777         print 'Lluis dnamesvar:', dnamesvar
    778         print 'Lluis dvnamesvar:', dvnamesvar
    779808        diagout = np.zeros(var0.shape, dtype=np.float)
    780809        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ptarget(var0, var1, var2, \
    781810          var3, var4, 700000., dnamesvar, dvnamesvar)
    782811
    783         print 'Lluis diagoutd:', diagoutd
    784         print 'Lluis diagoutvd:', diagoutvd
    785         print 'Lluis shape:', diagout.shape
    786         # Removing the nonChecking variable-dimensions from the initial list
    787         varsadd = []
    788         for nonvd in NONchkvardims:
    789             if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
    790             varsadd.append(nonvd)
    791         print 'Lluis after:', diagoutvd
    792         print 'Lluis diagoutd:', diagoutd
    793         print 'Lluis diagoutvd:', diagoutvd
    794         print 'Lluis shape:', diagout.shape
     812        # Removing the nonChecking variable-dimensions from the initial list
     813        varsadd = []
     814        for nonvd in NONchkvardims:
     815            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
     816            varsadd.append(nonvd)
    795817
    796818        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)

  • trunk/tools/module_ForDiagnostics.f90

    r1758 r1759  
    568568  END FUNCTION VirtualTemp1D
    569569
     570! ---- BEGIN modified from module_diag_afwa.F ---- !
     571
     572  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     573  !~
     574  !~ Name:
     575  !~    Theta
     576  !~
     577  !~ Description:
     578  !~    This function calculates potential temperature as defined by
     579  !~    Poisson's equation, given temperature and pressure ( hPa ).
     580  !~
     581  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     582  FUNCTION Theta ( t, p )
     583
     584    IMPLICIT NONE
     585
     586     !~ Variable declaration
     587     !  --------------------
     588     REAL(r_k), INTENT ( IN )                            :: t
     589     REAL(r_k), INTENT ( IN )                            :: p
     590     REAL(r_k)                                           :: theta
     591
     592     ! Using WRF values
     593     !REAL :: Rd ! Dry gas constant
     594     !REAL :: Cp ! Specific heat of dry air at constant pressure
     595     !REAL :: p00 ! Standard pressure ( 1000 hPa )
     596     REAL(r_k)                                           :: Rd, p00
     597 
     598     !Rd =  287.04
     599     !Cp = 1004.67
     600     !p00 = 1000.00
     601
     602     Rd = r_d
     603     p00 = p1000mb/100.
     604
     605     !~ Poisson's equation
     606     !  ------------------
     607     theta = t * ( (p00/p)**(Rd/Cp) )
     608 
     609  END FUNCTION Theta
     610
     611  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     612  !~
     613  !~ Name:
     614  !~    Thetae
     615  !~
     616  !~ Description:
     617  !~    This function returns equivalent potential temperature using the
     618  !~    method described in Bolton 1980, Monthly Weather Review, equation 43.
     619  !~
     620  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     621  FUNCTION Thetae ( tK, p, rh, mixr )
     622
     623    IMPLICIT NONE
     624
     625     !~ Variable Declarations
     626     !  ---------------------
     627     REAL(r_k) :: tK        ! Temperature ( K )
     628     REAL(r_k) :: p         ! Pressure ( hPa )
     629     REAL(r_k) :: rh        ! Relative humidity
     630     REAL(r_k) :: mixr      ! Mixing Ratio ( kg kg^-1)
     631     REAL(r_k) :: te        ! Equivalent temperature ( K )
     632     REAL(r_k) :: thetae    ! Equivalent potential temperature
     633 
     634     ! Using WRF values
     635     !REAL, PARAMETER :: R  = 287.04         ! Universal gas constant (J/deg kg)
     636     !REAL, PARAMETER :: P0 = 1000.0         ! Standard pressure at surface (hPa)
     637     REAL(r_k)                                                :: R, p00, Lv
     638     !REAL, PARAMETER :: lv = 2.54*(10**6)   ! Latent heat of vaporization
     639                                            ! (J kg^-1)
     640     !REAL, PARAMETER :: cp = 1004.67        ! Specific heat of dry air constant
     641                                            ! at pressure (J/deg kg)
     642     REAL(r_k) :: tlc                            ! LCL temperature
     643 
     644     R = r_d
     645     p00 = p1000mb/100.
     646     lv = XLV
     647
     648     !~ Calculate the temperature of the LCL
     649     !  ------------------------------------
     650     tlc = TLCL ( tK, rh )
     651 
     652     !~ Calculate theta-e
     653     !  -----------------
     654     thetae = (tK * (p00/p)**( (R/Cp)*(1.- ( (.28E-3)*mixr*1000.) ) ) )* &
     655                 exp( (((3.376/tlc)-.00254))*&
     656                    (mixr*1000.*(1.+(.81E-3)*mixr*1000.)) )
     657 
     658  END FUNCTION Thetae
     659
     660  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     661  !~
     662  !~ Name:
     663  !~    The2T.f90
     664  !~
     665  !~ Description:
     666  !~    This function returns the temperature at any pressure level along a
     667  !~    saturation adiabat by iteratively solving for it from the parcel
     668  !~    thetae.
     669  !~
     670  !~ Dependencies:
     671  !~    function thetae.f90
     672  !~
     673  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     674  FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel )
     675
     676    IMPLICIT NONE
     677 
     678     !~ Variable Declaration
     679     !  --------------------
     680     REAL(r_k),    INTENT     ( IN ) :: thetaeK
     681     REAL(r_k),    INTENT     ( IN ) :: pres
     682     LOGICAL, INTENT ( INOUT )  :: flag
     683     REAL(r_k)                       :: tparcel
     684 
     685     REAL(r_k) :: thetaK
     686     REAL(r_k) :: tovtheta
     687     REAL(r_k) :: tcheck
     688     REAL(r_k) :: svpr, svpr2
     689     REAL(r_k) :: smixr, smixr2
     690     REAL(r_k) :: thetae_check, thetae_check2
     691     REAL(r_k) :: tguess_2, correction
     692 
     693     LOGICAL :: found
     694     INTEGER :: iter
     695 
     696     ! Using WRF values
     697     !REAL :: R     ! Dry gas constant
     698     !REAL :: Cp    ! Specific heat for dry air
     699     !REAL :: kappa ! Rd / Cp
     700     !REAL :: Lv    ! Latent heat of vaporization at 0 deg. C
     701     REAL(r_k)                                                :: R, kappa, Lv
     702
     703     R = r_d
     704     Lv = XLV
     705     !R     = 287.04
     706     !Cp    = 1004.67
     707     Kappa = R/Cp
     708     !Lv    = 2.500E+6
     709
     710     !~ Make initial guess for temperature of the parcel
     711     !  ------------------------------------------------
     712     tovtheta = (pres/100000.0)**(r/cp)
     713     tparcel  = thetaeK/exp(lv*.012/(cp*295.))*tovtheta
     714
     715     iter = 1
     716     found = .false.
     717     flag = .false.
     718
     719     DO
     720        IF ( iter > 105 ) EXIT
     721
     722        tguess_2 = tparcel + REAL ( 1 )
     723
     724        svpr   = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) )
     725        smixr  = ( 0.622*svpr ) / ( (pres/100.0)-svpr )
     726        svpr2  = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) )
     727        smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 )
     728
     729        !  ------------------------------------------------------------------ ~!
     730        !~ When this function was orinially written, the final parcel         ~!
     731        !~ temperature check was based off of the parcel temperature and      ~!
     732        !~ not the theta-e it produced.  As there are multiple temperature-   ~!
     733        !~ mixing ratio combinations that can produce a single theta-e value, ~!
     734        !~ we change the check to be based off of the resultant theta-e       ~!
     735        !~ value.  This seems to be the most accurate way of backing out      ~!
     736        !~ temperature from theta-e.                                          ~!
     737        !~                                                                    ~!
     738        !~ Rentschler, April 2010                                             ~!
     739        !  ------------------------------------------------------------------  !
     740
     741        !~ Old way...
     742        !thetaK = thetaeK / EXP (lv * smixr  /(cp*tparcel) )
     743        !tcheck = thetaK * tovtheta
     744
     745        !~ New way
     746        thetae_check  = Thetae ( tparcel,  pres/100., 100., smixr  )
     747        thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 )
     748
     749        !~ Whew doggies - that there is some accuracy...
     750        !IF ( ABS (tparcel-tcheck) < .05) THEN
     751        IF ( ABS (thetaeK-thetae_check) < .001) THEN
     752           found = .true.
     753           flag  = .true.
     754           EXIT
     755        END IF
     756
     757        !~ Old
     758        !tparcel = tparcel + (tcheck - tparcel)*.3
     759
     760        !~ New
     761        correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check )
     762        tparcel = tparcel + correction
     763
     764        iter = iter + 1
     765     END DO
     766
     767     !IF ( .not. found ) THEN
     768     !   print*, "Warning! Thetae to temperature calculation did not converge!"
     769     !   print*, "Thetae ", thetaeK, "Pressure ", pres
     770     !END IF
     771
     772  END FUNCTION The2T
     773
     774  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     775  !~
     776  !~ Name:
     777  !~    VirtualTemperature
     778  !~
     779  !~ Description:
     780  !~    This function returns virtual temperature given temperature ( K )
     781  !~    and mixing ratio.
     782  !~
     783  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     784  FUNCTION VirtualTemperature ( tK, w ) result ( Tv )
     785
     786    IMPLICIT NONE
     787
     788     !~ Variable declaration
     789     real(r_k), intent ( in ) :: tK !~ Temperature
     790     real(r_k), intent ( in ) :: w  !~ Mixing ratio ( kg kg^-1 )
     791     real(r_k)                :: Tv !~ Virtual temperature
     792
     793     Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w )
     794
     795  END FUNCTION VirtualTemperature
     796
     797  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     798  !~
     799  !~ Name:
     800  !~    SaturationMixingRatio
     801  !~
     802  !~ Description:
     803  !~    This function calculates saturation mixing ratio given the
     804  !~    temperature ( K ) and the ambient pressure ( Pa ).  Uses
     805  !~    approximation of saturation vapor pressure.
     806  !~
     807  !~ References:
     808  !~    Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10
     809  !~
     810  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     811  FUNCTION SaturationMixingRatio ( tK, p ) result ( ws )
     812
     813    IMPLICIT NONE
     814
     815    REAL(r_k), INTENT ( IN ) :: tK
     816    REAL(r_k), INTENT ( IN ) :: p
     817    REAL(r_k)                :: ws
     818
     819    REAL(r_k) :: es
     820
     821    es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) )
     822    ws = ( 0.622*es ) / ( (p/100.0)-es )
     823
     824  END FUNCTION SaturationMixingRatio
     825
     826  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     827  !~                                                                     
     828  !~ Name:                                                               
     829  !~    tlcl                                                               
     830  !~                                                                       
     831  !~ Description:                                                           
     832  !~    This function calculates the temperature of a parcel of air would have
     833  !~    if lifed dry adiabatically to it's lifting condensation level (lcl). 
     834  !~                                                                         
     835  !~ References:                                                             
     836  !~    Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22
     837  !~                                                                         
     838  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     839  FUNCTION TLCL ( tk, rh )
     840   
     841    IMPLICIT NONE
     842 
     843    REAL(r_k), INTENT ( IN ) :: tK   !~ Temperature ( K )
     844    REAL(r_k), INTENT ( IN ) :: rh   !~ Relative Humidity ( % )
     845    REAL(r_k)                :: tlcl
     846   
     847    REAL(r_k) :: denom, term1, term2
     848
     849    term1 = 1.0 / ( tK - 55.0 )
     850!! Lluis
     851!    IF ( rh > REAL (0) ) THEN
     852    IF ( rh > zeroRK ) THEN
     853      term2 = ( LOG (rh/100.0)  / 2840.0 )
     854    ELSE
     855      term2 = ( LOG (0.001/oneRK) / 2840.0 )
     856    END IF
     857    denom = term1 - term2
     858!! Lluis
     859!    tlcl = ( 1.0 / denom ) + REAL ( 55 )
     860    tlcl = ( oneRK / denom ) + 55*oneRK
     861
     862  END FUNCTION TLCL
     863
     864  FUNCTION var_cape_afwa1D(nz, tk, rhv, p, hgtv, sfc, cape, cin, zlfc, plfc, lidx, parcel) RESULT (ostat)
     865! Function to compute cape on a 1D column following implementation in phys/module_diag_afwa.F
     866
     867    IMPLICIT NONE
     868
     869    INTEGER, INTENT(in)                                  :: nz, sfc
     870    REAL(r_k), DIMENSION(nz), INTENT(in)                 :: tk, rhv, p, hgtv
     871    REAL(r_k), INTENT(out)                               :: cape, cin, zlfc, plfc, lidx
     872    INTEGER                                              :: ostat
     873    INTEGER, INTENT(in)                                  :: parcel
     874 
     875    ! Local
     876    !~ Derived profile variables
     877    !  -------------------------
     878    REAL(r_k), DIMENSION(nz)                             :: rh, hgt, ws, w, dTvK, buoy
     879    REAL(r_k)                                            :: tlclK, plcl, nbuoy, pbuoy
     880 
     881    !~ Source parcel information
     882    !  -------------------------
     883    REAL(r_k)                                            :: srctK, srcrh, srcws, srcw, srcp,          &
     884      srctheta, srcthetaeK
     885    INTEGER                                              :: srclev
     886    REAL(r_k)                                            :: spdiff
     887   
     888    !~ Parcel variables
     889    !  ----------------
     890    REAL(r_k)                                            :: ptK, ptvK, tvK, pw
     891 
     892    !~ Other utility variables
     893    !  -----------------------
     894    INTEGER                                              :: i, j, k
     895    INTEGER                                              :: lfclev
     896    INTEGER                                              :: prcl
     897    INTEGER                                              :: mlev
     898    INTEGER                                              :: lyrcnt
     899    LOGICAL                                              :: flag
     900    LOGICAL                                              :: wflag
     901    REAL(r_k)                                            :: freeze
     902    REAL(r_k)                                            :: pdiff
     903    REAL(r_k)                                            :: pm, pu, pd
     904    REAL(r_k)                                            :: lidxu
     905    REAL(r_k)                                            :: lidxd
     906 
     907    REAL(r_k), PARAMETER                                 :: Rd = r_d
     908    REAL(r_k), PARAMETER                                 :: RUNDEF = -9.999E30
     909
     910!!!!!!! Variables
     911! nz: Number of vertical levels
     912! sfc: Surface level in the profile
     913! tk: Temperature profile [K]
     914! rhv: Relative Humidity profile [1]
     915! rh: Relative Humidity profile [%]
     916! p: Pressure profile [Pa]
     917! hgtv: Height profile [m]
     918! hgt: Height profile [gpm]
     919! cape: CAPE [Jkg-1]
     920! cin: CIN [Jkg-1]
     921! zlfc: LFC Height [gpm]
     922! plfc: LFC Pressure [Pa]
     923! lidx: Lifted index
     924!   FROM: https://en.wikipedia.org/wiki/Lifted_index
     925!     lidx >= 6: Very Stable Conditions
     926!     6 > lidx > 1: Stable Conditions, Thunderstorms Not Likely
     927!     0 > lidx > -2: Slightly Unstable, Thunderstorms Possible, With Lifting Mechanism (i.e., cold front, daytime heating, ...)
     928!     -2 > lidx > -6: Unstable, Thunderstorms Likely, Some Severe With Lifting Mechanism
     929!     -6 > lidx: Very Unstable, Severe Thunderstorms Likely With Lifting Mechanism
     930! ostat: Function return status (Nonzero is bad)
     931! parcel:
     932!   Most Unstable = 1 (default)
     933!   Mean layer = 2
     934!   Surface based = 3
     935!~ Derived profile variables
     936!  -------------------------
     937! ws: Saturation mixing ratio
     938! w: Mixing ratio
     939! dTvK: Parcel / ambient Tv difference
     940! buoy: Buoyancy
     941! tlclK: LCL temperature [K]
     942! plcl: LCL pressure [Pa]
     943! nbuoy: Negative buoyancy
     944! pbuoy: Positive buoyancy
     945 
     946!~ Source parcel information
     947!  -------------------------
     948! srctK: Source parcel temperature [K]
     949! srcrh: Source parcel rh [%]
     950! srcws: Source parcel sat. mixing ratio
     951! srcw: Source parcel mixing ratio
     952! srcp: Source parcel pressure [Pa]
     953! srctheta: Source parcel theta [K]
     954! srcthetaeK: Source parcel theta-e [K]
     955! srclev: Level of the source parcel
     956! spdiff: Pressure difference
     957   
     958!~ Parcel variables
     959!  ----------------
     960! ptK: Parcel temperature [K]
     961! ptvK: Parcel virtual temperature [K]
     962! tvK: Ambient virtual temperature [K]
     963! pw: Parcel mixing ratio
     964 
     965!~ Other utility variables
     966!  -----------------------
     967! lfclev: Level of LFC
     968! prcl: Internal parcel type indicator
     969! mlev: Level for ML calculation
     970! lyrcnt: Number of layers in mean layer
     971! flag: Dummy flag
     972! wflag: Saturation flag
     973! freeze: Water loading multiplier
     974! pdiff: Pressure difference between levs
     975! pm, pu, pd: Middle, upper, lower pressures
     976! lidxu: Lifted index at upper level
     977! lidxd: Lifted index at lower level
     978
     979    fname = 'var_cape_afwa' 
     980
     981    !~ Initialize variables
     982    !  --------------------
     983    rh = rhv*100.
     984    hgt = hgtv*g
     985    ostat = 0
     986    CAPE = zeroRK
     987    CIN = zeroRK
     988    ZLFC = RUNDEF
     989    PLFC = RUNDEF
     990 
     991    !~ Look for submitted parcel definition
     992    !~ 1 = Most unstable
     993    !~ 2 = Mean layer
     994    !~ 3 = Surface based
     995    !  -------------------------------------
     996    IF ( parcel > 3 .or. parcel < 1 ) THEN
     997       prcl = 1
     998    ELSE
     999       prcl =  parcel
     1000    END IF
     1001 
     1002    !~ Initalize our parcel to be (sort of) surface based.  Because of
     1003    !~ issues we've been observing in the WRF model, specifically with
     1004    !~ excessive surface moisture values at the surface, using a true
     1005    !~ surface based parcel is resulting a more unstable environment
     1006    !~ than is actually occuring.  To address this, our surface parcel
     1007    !~ is now going to be defined as the parcel between 25-50 hPa
     1008    !~ above the surface. UPDATE - now that this routine is in WRF,
     1009    !~ going to trust surface info. GAC 20140415
     1010    !  ----------------------------------------------------------------
     1011 
     1012    !~ Compute mixing ratio values for the layer
     1013    !  -----------------------------------------
     1014    DO k = sfc, nz
     1015      ws  ( k )   = SaturationMixingRatio ( tK(k), p(k) )
     1016      w   ( k )   = ( rh(k)/100.0 ) * ws ( k )
     1017    END DO
     1018 
     1019    srclev      = sfc
     1020    srctK       = tK    ( sfc )
     1021    srcrh       = rh    ( sfc )
     1022    srcp        = p     ( sfc )
     1023    srcws       = ws    ( sfc )
     1024    srcw        = w     ( sfc )
     1025    srctheta    = Theta ( tK(sfc), p(sfc)/100.0 )
     1026   
     1027      !~ Compute the profile mixing ratio.  If the parcel is the MU parcel,
     1028      !~ define our parcel to be the most unstable parcel within the lowest
     1029      !~ 180 mb.
     1030      !  -------------------------------------------------------------------
     1031      mlev = sfc + 1
     1032      DO k = sfc + 1, nz
     1033   
     1034         !~ Identify the last layer within 100 hPa of the surface
     1035         !  -----------------------------------------------------
     1036         pdiff = ( p (sfc) - p (k) ) / REAL ( 100 )
     1037         IF ( pdiff <= REAL (100) ) mlev = k
     1038
     1039         !~ If we've made it past the lowest 180 hPa, exit the loop
     1040         !  -------------------------------------------------------
     1041         IF ( pdiff >= REAL (180) ) EXIT
     1042
     1043         IF ( prcl == 1 ) THEN
     1044            !IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
     1045            IF ( (w(k) > srcw) ) THEN
     1046               srctheta = Theta ( tK(k), p(k)/100.0 )
     1047               srcw = w ( k )
     1048               srclev  = k
     1049               srctK   = tK ( k )
     1050               srcrh   = rh ( k )
     1051               srcp    = p  ( k )
     1052            END IF
     1053         END IF
     1054   
     1055      END DO
     1056   
     1057      !~ If we want the mean layer parcel, compute the mean values in the
     1058      !~ lowest 100 hPa.
     1059      !  ----------------------------------------------------------------
     1060      lyrcnt =  mlev - sfc + 1
     1061      IF ( prcl == 2 ) THEN
     1062   
     1063         srclev   = sfc
     1064         srctK    = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt )
     1065         srcw     = SUM ( w  (sfc:mlev) ) / REAL ( lyrcnt )
     1066         srcrh    = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt )
     1067         srcp     = SUM ( p  (sfc:mlev) ) / REAL ( lyrcnt )
     1068         srctheta = Theta ( srctK, srcp/100. )
     1069   
     1070      END IF
     1071   
     1072      srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw )
     1073   
     1074      !~ Calculate temperature and pressure of the LCL
     1075      !  ---------------------------------------------
     1076      tlclK = TLCL ( tK(srclev), rh(srclev) )
     1077      plcl  = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) )
     1078   
     1079      !~ Now lift the parcel
     1080      !  -------------------
     1081   
     1082      buoy  = REAL ( 0 )
     1083      pw    = srcw
     1084      wflag = .false.
     1085      DO k  = srclev, nz
     1086         IF ( p (k) <= plcl ) THEN
     1087   
     1088            !~ The first level after we pass the LCL, we're still going to
     1089            !~ lift the parcel dry adiabatically, as we haven't added the
     1090            !~ the required code to switch between the dry adiabatic and moist
     1091            !~ adiabatic cooling.  Since the dry version results in a greater
     1092            !~ temperature loss, doing that for the first step so we don't over
     1093            !~ guesstimate the instability.
     1094            !  ----------------------------------------------------------------
     1095   
     1096            IF ( wflag ) THEN
     1097               flag  = .false.
     1098   
     1099               !~ Above the LCL, our parcel is now undergoing moist adiabatic
     1100               !~ cooling.  Because of the latent heating being undergone as
     1101               !~ the parcel rises above the LFC, must iterative solve for the
     1102               !~ parcel temperature using equivalant potential temperature,
     1103               !~ which is conserved during both dry adiabatic and
     1104               !~ pseudoadiabatic displacements.
     1105               !  --------------------------------------------------------------
     1106               ptK   = The2T ( srcthetaeK, p(k), flag )
     1107   
     1108               !~ Calculate the parcel mixing ratio, which is now changing
     1109               !~ as we condense moisture out of the parcel, and is equivalent
     1110               !~ to the saturation mixing ratio, since we are, in theory, at
     1111               !~ saturation.
     1112               !  ------------------------------------------------------------
     1113               pw = SaturationMixingRatio ( ptK, p(k) )
     1114   
     1115               !~ Now we can calculate the virtual temperature of the parcel
     1116               !~ and the surrounding environment to assess the buoyancy.
     1117               !  ----------------------------------------------------------
     1118               ptvK  = VirtualTemperature ( ptK, pw )
     1119               tvK   = VirtualTemperature ( tK (k), w (k) )
     1120   
     1121               !~ Modification to account for water loading
     1122               !  -----------------------------------------
     1123               freeze = 0.033 * ( 263.15 - pTvK )
     1124               IF ( freeze > 1.0 ) freeze = 1.0
     1125               IF ( freeze < 0.0 ) freeze = 0.0
     1126   
     1127               !~ Approximate how much of the water vapor has condensed out
     1128               !~ of the parcel at this level
     1129               !  ---------------------------------------------------------
     1130               freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7
     1131   
     1132               pTvK = pTvK - pTvK * ( srcw - pw ) + freeze
     1133               dTvK ( k ) = ptvK - tvK
     1134               buoy ( k ) = g * ( dTvK ( k ) / tvK )
     1135   
     1136            ELSE
     1137   
     1138               !~ Since the theta remains constant whilst undergoing dry
     1139               !~ adiabatic processes, can back out the parcel temperature
     1140               !~ from potential temperature below the LCL
     1141               !  --------------------------------------------------------
     1142               ptK   = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
     1143   
     1144               !~ Grab the parcel virtual temperture, can use the source
     1145               !~ mixing ratio since we are undergoing dry adiabatic cooling
     1146               !  ----------------------------------------------------------
     1147               ptvK  = VirtualTemperature ( ptK, srcw )
     1148   
     1149               !~ Virtual temperature of the environment
     1150               !  --------------------------------------
     1151               tvK   = VirtualTemperature ( tK (k), w (k) )
     1152   
     1153               !~ Buoyancy at this level
     1154               !  ----------------------
     1155               dTvK ( k ) = ptvK - tvK
     1156               buoy ( k ) = g * ( dtvK ( k ) / tvK )
     1157   
     1158               wflag = .true.
     1159   
     1160            END IF
     1161   
     1162         ELSE
     1163   
     1164            !~ Since the theta remains constant whilst undergoing dry
     1165            !~ adiabatic processes, can back out the parcel temperature
     1166            !~ from potential temperature below the LCL
     1167            !  --------------------------------------------------------
     1168            ptK   = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
     1169   
     1170            !~ Grab the parcel virtual temperture, can use the source
     1171            !~ mixing ratio since we are undergoing dry adiabatic cooling
     1172            !  ----------------------------------------------------------
     1173            ptvK  = VirtualTemperature ( ptK, srcw )
     1174   
     1175            !~ Virtual temperature of the environment
     1176            !  --------------------------------------
     1177            tvK   = VirtualTemperature ( tK (k), w (k) )
     1178   
     1179            !~ Buoyancy at this level
     1180            !  ---------------------
     1181            dTvK ( k ) = ptvK - tvK
     1182            buoy ( k ) = g * ( dtvK ( k ) / tvK )
     1183   
     1184         END IF
     1185
     1186         !~ Chirp
     1187         !  -----
     1188  !          WRITE ( *,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw*1000.,ptvK,tvK,buoy(k)
     1189   
     1190      END DO
     1191   
     1192      !~ Add up the buoyancies, find the LFC
     1193      !  -----------------------------------
     1194      flag   = .false.
     1195      lfclev = -1
     1196      nbuoy  = REAL ( 0 )
     1197      pbuoy = REAL ( 0 )
     1198      DO k = sfc + 1, nz
     1199         IF ( tK (k) < 253.15 ) EXIT
     1200         CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
     1201         CIN  = CIN  + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
     1202   
     1203         !~ If we've already passed the LFC
     1204         !  -------------------------------
     1205         IF ( flag .and. buoy (k) > REAL (0) ) THEN
     1206            pbuoy = pbuoy + buoy (k)
     1207         END IF
     1208   
     1209         !~ We are buoyant now - passed the LFC
     1210         !  -----------------------------------
     1211         IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN
     1212            flag = .true.
     1213            pbuoy = pbuoy + buoy (k)
     1214            lfclev = k
     1215         END IF
     1216   
     1217         !~ If we think we've passed the LFC, but encounter a negative layer
     1218         !~ start adding it up.
     1219         !  ----------------------------------------------------------------
     1220         IF ( flag .and. buoy (k) < REAL (0) ) THEN
     1221            nbuoy = nbuoy + buoy (k)
     1222
     1223            !~ If the accumulated negative buoyancy is greater than the
     1224            !~ positive buoyancy, then we are capped off.  Got to go higher
     1225            !~ to find the LFC. Reset positive and negative buoyancy summations
     1226            !  ----------------------------------------------------------------
     1227            IF ( ABS (nbuoy) > pbuoy ) THEN
     1228               flag   = .false.
     1229               nbuoy  = REAL ( 0 )
     1230               pbuoy  = REAL ( 0 )
     1231               lfclev = -1
     1232            END IF
     1233         END IF
     1234
     1235      END DO
     1236
     1237      !~ Calculate lifted index by interpolating difference between
     1238      !~ parcel and ambient Tv to 500mb.
     1239      !  ----------------------------------------------------------
     1240      DO k = sfc + 1, nz
     1241
     1242         pm = 50000.
     1243         pu = p ( k )
     1244         pd = p ( k - 1 )
     1245
     1246         !~ If we're already above 500mb just set lifted index to 0.
     1247         !~ --------------------------------------------------------
     1248         IF ( pd .le. pm ) THEN
     1249            lidx = zeroRK
     1250            EXIT
     1251   
     1252         ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
     1253
     1254            !~ Found trapping pressure: up, middle, down.
     1255            !~ We are doing first order interpolation. 
     1256            !  ------------------------------------------
     1257            lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
     1258            lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
     1259            lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
     1260            EXIT
     1261
     1262         ENDIF
     1263
     1264      END DO
     1265   
     1266      !~ Assuming the the LFC is at a pressure level for now
     1267      !  ---------------------------------------------------
     1268      IF ( lfclev > zeroRK ) THEN
     1269         PLFC = p   ( lfclev )
     1270         ZLFC = hgt ( lfclev )
     1271      END IF
     1272   
     1273      IF ( PLFC /= PLFC .OR. PLFC < zeroRK ) THEN
     1274         PLFC = -oneRK
     1275         ZLFC = -oneRK
     1276      END IF
     1277   
     1278      IF ( CAPE /= CAPE ) cape = zeroRK
     1279   
     1280      IF ( CIN  /= CIN  ) cin  = zeroRK
     1281
     1282      !~ Chirp
     1283      !  -----
     1284  !       WRITE ( *,* ) ' CAPE: ', cape, ' CIN:  ', cin
     1285  !       WRITE ( *,* ) ' LFC:  ', ZLFC, ' PLFC: ', PLFC
     1286  !       WRITE ( *,* ) ''
     1287  !       WRITE ( *,* ) ' Exiting buoyancy.'
     1288  !       WRITE ( *,* ) ' ==================================== '
     1289  !       WRITE ( *,* ) ''
     1290   
     1291    RETURN
     1292
     1293  END FUNCTION var_cape_afwa1D
     1294
     1295! ---- END modified from module_diag_afwa.F ---- !
     1296
     1297  SUBROUTINE compute_cape_afwa4D(ta, hur, press, zg, hgt, cape, cin, zlfc, plfc, li, parcelmethod,    &
     1298    d1, d2, d3, d4)
     1299! Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute CAPE, CIN, ZLFC, PLFC, LI
     1300
     1301    IMPLICIT NONE
     1302
     1303    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4, parcelmethod
     1304    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ta, hur, press, zg
     1305    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
     1306    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: cape, cin, zlfc, plfc, li
     1307 
     1308! Local
     1309    INTEGER                                              :: i, j, it
     1310    INTEGER                                              :: ofunc
     1311
     1312!!!!!!! Variables
     1313! ta: air temperature [K]
     1314! hur: relative humidity [%]
     1315! press: air pressure [Pa]
     1316! zg: geopotential height [gpm]
     1317! hgt: topographical height [m]
     1318! cape: Convective available potential energy [Jkg-1]
     1319! cin: Convective inhibition [Jkg-1]
     1320! zlfc: height at the Level of free convection [m]
     1321! plfc: pressure at the Level of free convection [Pa]
     1322! li: lifted index [1]
     1323! parcelmethod:
     1324!   Most Unstable = 1 (default)
     1325!   Mean layer = 2
     1326!   Surface based = 3
     1327
     1328    fname = 'compute_cape_afwa4D'
     1329
     1330    DO i=1, d1
     1331      DO j=1, d2
     1332        DO it=1, d4
     1333          ofunc = var_cape_afwa1D(d3, ta(i,j,:,it), hur(i,j,:,it), press(i,j,:,it), zg(i,j,:,it),     &
     1334            1, cape(i,j,it), cin(i,j,it), zlfc(i,j,it), plfc(i,j,it), li(i,j,it), parcelmethod)
     1335          zlfc(i,j,it) = zlfc(i,j,it)/g - hgt(i,j)
     1336        END DO
     1337      END DO
     1338    END DO
     1339
     1340    RETURN
     1341
     1342  END SUBROUTINE compute_cape_afwa4D
    5701343
    5711344END MODULE module_ForDiagnostics

  • trunk/tools/module_definitions.f90

    r1758 r1759  
    3131  REAL(r_k), PARAMETER                                   :: RCp = 0.286 ! R*Cp [-]
    3232  REAL(r_k), PARAMETER                                   :: p0ref = 100000 ! pressure reference [Pa]
     33! WRF gravity
     34  REAL(r_k), PARAMETER                                   :: g = 9.81
    3335! Ratio between molecular weights of water and dry air
    3436  REAL(r_k), PARAMETER                                   :: mol_watdry = 0.622

  • trunk/tools/variables_values.dat

    r1755 r1759  
    2222oliq, c, condensed_water_mixing_ratio, 0., 3.e-4, condensed|water|mixing|ratio, kgkg-1, BuPu
    2323OLIQ, c, condensed_water_mixing_ratio, 0., 3.e-4, condensed|water|mixing|ratio, kgkg-1, BuPu
     24cape, cape, convective_available_potential_energy, 0., 40000., convective|available|potential|energy, Jkg-1, Reds
    2425cdrag, cdrag, cdrag, 0., 0.01, drag|coefficient|for|LW|and|SH, '-', rainbow
    2526ci, ci, cloud_iced_water_mixing_ratio, 0., 0.0003, cloud|iced|water|mixing|ratio, kgkg-1, Purples
    2627iwcon, ci, cloud_iced_water_mixing_ratio, 0., 0.0003, cloud|iced|water|mixing|ratio, kgkg-1, Purples
    2728LIWCON, ci, cloud_iced_water_mixing_ratio, 0., 0.0003, cloud|iced|water|mixing|ratio, kgkg-1, Purples
     29cin, cin, convective_inhibition, 0., 40000., convective|inhibition, Jkg-1, Greens
    2830cl, cl, cloud_liquidwater_mixing_ratio, 0., 0.0003, cloud|liquid|water|mixing|ratio, kgkg-1, Blues
    2931lwcon, cl, cloud_liquidwater_mixing_ratio, 0., 0.0003, cloud|liquid|water|mixing|ratio, kgkg-1, Blues
    … …  
    236238lambda_th, lambdath, thermal_plume_vertical_velocity, -30., 30., thermal|plume|vertical|velocity, m/s, seismic
    237239LLAMBDA_TH, lambdath, thermal_plume_vertical_velocity, -30., 30., thermal|plume|vertical|velocity, m/s, seismic
    238 lev, lev, lev,  0., 100., vertical level, -, Greens
    239 iz, lev, lev,  0., 100., vertical level, -, Greens
     240lev, lev, lev,  0., 100., vertical|level, -, Greens
     241iz, lev, lev,  0., 100., vertical|level, -, Greens
     242li, li, lifted_index, -20., 20., lifted|index, -, Sesimic
    240243lmaxth, lmaxth, upper_level_thermals, 0., 100., upper|level|thermals, 1, Greens
    241244LLMAXTH, lmaxth, upper_level_thermals, 0., 100., upper|level|thermals, 1, Greens
    … …  
    252255HGT, orog, orography,  0., 3000., surface|altitude, m,terrain
    253256HGT_M, orog, orography,  0., 3000., surface|altitude, m,terrain
     257plfc, plfc, pressure_level_free_convection, 100., 101000., pressure|of|level|free|convection, Pa, BuPu
    254258pfc, pfc, pressure_free_convection, 100., 1100., pressure|free|convection, hPa, BuPu
    255259plfc, pfc, pressure_free_convection, 100., 1100., pressure|free|convection, hPa, BuPu
    … …  
    649653zg_per, zg_per, geopotential_height_per, 0., 80000., geopotential|height|perturbation, m2s-2, rainbow
    650654PH, zg_per, geopotential_height_per, 0., 80000., geopotential|height|perturbation, m2s-2, rainbow
     655zlfc, zlfc, height_level_free_convection, 100., 2000., height|of|level|free|convection, m, BuPu
    651656zmaxth, zmaxth, thermal_plume_height, 0., 4000., maximum|thermals|plume|height, m, YlOrRd
    652657zmax_th, zmaxth, thermal_plume_height, 0., 4000., maximum|thermals|plume|height, m, YlOrRd
Note: See TracChangeset for help on using the changeset viewer.