Changeset 1769 in lmdz_wrf


Ignore:
Timestamp:
Feb 8, 2018, 9:10:45 PM (7 years ago)
Author:
lfita
Message:

Re-ordering subroutines/functions to their apropriate module

Location:
trunk/tools
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/module_ForDiagnostics.f90

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

  • trunk/tools/module_ForDiagnosticsVars.f90

    r1608 r1769  
    2121
    2222!!!!!!! Variables
     23! compute_psl_ptarget4d2: Compute sea level pressure using a target pressure. Similar to the Benjamin
     24!   and Miller (1990). Method found in p_interp.F90
     25! compute_tv4d: 4D calculation of virtual temperaure
     26! SaturationMixingRatio: WRF's AFWA method to compute the saturation mixing ratio
     27! The2T: WRF's AFWA method to compute the temperature at any pressure level along a saturation adiabat
     28!   by iteratively solving for it from the parcel thetae.
     29! Theta: WRF's AFWA method to compute potential temperature
     30! Thetae: WRF's AFWA method to compute equivalent potential temperature
     31! TLCL: WRF's AFWA method to compute the temperature of a parcel of air would have if lifed dry
     32!   adiabatically to it's lifting condensation level (lcl)
     33! var_cape_afwa1D: WRF's AFWA method to compute cape, cin, fclp, fclz and li
    2334! var_cllmh: low, medium, high-cloud [0,1]
    2435! var_clt: total cloudiness [0,1]
     36! VirtualTemp1D: Function for 1D calculation of virtual temperaure
     37! VirtualTemperature: WRF's AFWA method to compute virtual temperature
    2538
    2639!!!!!!! Calculations
    … …  
    115128  END FUNCTION var_clt
    116129
     130
     131  SUBROUTINE compute_psl_ptarget4d2(press, ps, hgt, ta, qv, ptarget, psl, d1, d2, d3, d4)
     132    ! Subroutine to compute sea level pressure using a target pressure. Similar to the Benjamin
     133    !   and Miller (1990). Method found in p_interp.F90
     134
     135    IMPLICIT NONE
     136
     137    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
     138    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: press, ta, qv
     139    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: ps
     140    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
     141    REAL(r_k), INTENT(in)                                :: ptarget
     142    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: psl
     143
     144! Local
     145    INTEGER                                              :: i, j, it
     146    INTEGER                                              :: kin
     147    INTEGER                                              :: kupper
     148    REAL(r_k)                                            :: dpmin, dp, tbotextrap,   &
     149      tvbotextrap, virtual
     150    ! Exponential related to standard atmosphere lapse rate r_d*gammav/g
     151    REAL(r_k), PARAMETER                                 :: expon=r_d*gammav/grav
     152
     153!!!!!!! Variables
     154! press: Atmospheric pressure [Pa]
     155! ps: surface pressure [Pa]
     156! hgt: surface height
     157! ta: temperature [K]
     158! qv: water vapor mixing ratio
     159! dz: number of vertical levels
     160! psl: sea-level pressure
     161
     162    fname = 'compute_psl_ptarget4d2'
     163
     164    ! Minimal distance between pressures [Pa]
     165    dpmin=1.e4
     166    psl=0.
     167
     168    DO i=1,d1
     169      DO j=1,d2
     170        IF (hgt(i,j) /= 0.) THEN
     171          DO it=1,d4
     172
     173            ! target pressure to be used for the extrapolation [Pa] (defined in namelist.input)
     174            !   ptarget = 70000. default value
     175
     176            ! We are below both the ground and the lowest data level.
     177
     178            !      First, find the model level that is closest to a "target" pressure
     179            !        level, where the "target" pressure is delta-p less that the local
     180            !        value of a horizontally smoothed surface pressure field.  We use
     181            !        delta-p = 150 hPa here. A standard lapse rate temperature profile
     182            !        passing through the temperature at this model level will be used
     183            !        to define the temperature profile below ground.  This is similar
     184            !        to the Benjamin and Miller (1990) method, using 
     185            !        700 hPa everywhere for the "target" pressure.
     186
     187            kupper = 0
     188            loop_kIN: DO kin=d3,1,-1
     189              kupper = kin
     190              dp=abs( press(i,j,kin,it) - ptarget )
     191              IF (dp .GT. dpmin) EXIT loop_kIN
     192              dpmin=min(dpmin,dp)
     193            ENDDO loop_kIN
     194
     195            tbotextrap=ta(i,j,kupper,it)*(ps(i,j,it)/ptarget)**expon
     196            tvbotextrap=virtualTemp1D(tbotextrap,qv(i,j,kupper,it))
     197
     198            psl(i,j,it) = ps(i,j,it)*((tvbotextrap+gammav*hgt(i,j))/tvbotextrap)**(1/expon)
     199          END DO
     200        ELSE
     201          psl(i,j,:) = ps(i,j,:)
     202        END IF
     203      END DO
     204    END DO
     205
     206    RETURN
     207
     208  END SUBROUTINE compute_psl_ptarget4d2
     209
     210  SUBROUTINE compute_tv4d(ta,qv,tv,d1,d2,d3,d4)
     211! 4D calculation of virtual temperaure
     212
     213    IMPLICIT NONE
     214
     215    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
     216    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ta, qv
     217    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out)       :: tv
     218
     219! Variables
     220! ta: temperature [K]
     221! qv: mixing ratio [kgkg-1]
     222! tv: virtual temperature
     223
     224    tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
     225
     226  END SUBROUTINE compute_tv4d
     227
     228  FUNCTION VirtualTemp1D (ta,qv) result (tv)
     229! 1D calculation of virtual temperaure
     230
     231    IMPLICIT NONE
     232
     233    REAL(r_k), INTENT(in)                                :: ta, qv
     234    REAL(r_k)                                            :: tv
     235
     236! Variables
     237! ta: temperature [K]
     238! qv: mixing ratio [kgkg-1]
     239
     240    tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
     241
     242  END FUNCTION VirtualTemp1D
     243
     244! ---- BEGIN modified from module_diag_afwa.F ---- !
     245
     246  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     247  !~
     248  !~ Name:
     249  !~    Theta
     250  !~
     251  !~ Description:
     252  !~    This function calculates potential temperature as defined by
     253  !~    Poisson's equation, given temperature and pressure ( hPa ).
     254  !~
     255  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     256  FUNCTION Theta ( t, p )
     257
     258    IMPLICIT NONE
     259
     260     !~ Variable declaration
     261     !  --------------------
     262     REAL(r_k), INTENT ( IN )                            :: t
     263     REAL(r_k), INTENT ( IN )                            :: p
     264     REAL(r_k)                                           :: theta
     265
     266     ! Using WRF values
     267     !REAL :: Rd ! Dry gas constant
     268     !REAL :: Cp ! Specific heat of dry air at constant pressure
     269     !REAL :: p00 ! Standard pressure ( 1000 hPa )
     270     REAL(r_k)                                           :: Rd, p00
     271 
     272     !Rd =  287.04
     273     !Cp = 1004.67
     274     !p00 = 1000.00
     275
     276     Rd = r_d
     277     p00 = p1000mb/100.
     278
     279     !~ Poisson's equation
     280     !  ------------------
     281     theta = t * ( (p00/p)**(Rd/Cp) )
     282 
     283  END FUNCTION Theta
     284
     285  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     286  !~
     287  !~ Name:
     288  !~    Thetae
     289  !~
     290  !~ Description:
     291  !~    This function returns equivalent potential temperature using the
     292  !~    method described in Bolton 1980, Monthly Weather Review, equation 43.
     293  !~
     294  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     295  FUNCTION Thetae ( tK, p, rh, mixr )
     296
     297    IMPLICIT NONE
     298
     299     !~ Variable Declarations
     300     !  ---------------------
     301     REAL(r_k) :: tK        ! Temperature ( K )
     302     REAL(r_k) :: p         ! Pressure ( hPa )
     303     REAL(r_k) :: rh        ! Relative humidity
     304     REAL(r_k) :: mixr      ! Mixing Ratio ( kg kg^-1)
     305     REAL(r_k) :: te        ! Equivalent temperature ( K )
     306     REAL(r_k) :: thetae    ! Equivalent potential temperature
     307 
     308     ! Using WRF values
     309     !REAL, PARAMETER :: R  = 287.04         ! Universal gas constant (J/deg kg)
     310     !REAL, PARAMETER :: P0 = 1000.0         ! Standard pressure at surface (hPa)
     311     REAL(r_k)                                                :: R, p00, Lv
     312     !REAL, PARAMETER :: lv = 2.54*(10**6)   ! Latent heat of vaporization
     313                                            ! (J kg^-1)
     314     !REAL, PARAMETER :: cp = 1004.67        ! Specific heat of dry air constant
     315                                            ! at pressure (J/deg kg)
     316     REAL(r_k) :: tlc                            ! LCL temperature
     317 
     318     R = r_d
     319     p00 = p1000mb/100.
     320     lv = XLV
     321
     322     !~ Calculate the temperature of the LCL
     323     !  ------------------------------------
     324     tlc = TLCL ( tK, rh )
     325 
     326     !~ Calculate theta-e
     327     !  -----------------
     328     thetae = (tK * (p00/p)**( (R/Cp)*(1.- ( (.28E-3)*mixr*1000.) ) ) )* &
     329                 exp( (((3.376/tlc)-.00254))*&
     330                    (mixr*1000.*(1.+(.81E-3)*mixr*1000.)) )
     331 
     332  END FUNCTION Thetae
     333
     334  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     335  !~
     336  !~ Name:
     337  !~    The2T.f90
     338  !~
     339  !~ Description:
     340  !~    This function returns the temperature at any pressure level along a
     341  !~    saturation adiabat by iteratively solving for it from the parcel
     342  !~    thetae.
     343  !~
     344  !~ Dependencies:
     345  !~    function thetae.f90
     346  !~
     347  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     348  FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel )
     349
     350    IMPLICIT NONE
     351 
     352     !~ Variable Declaration
     353     !  --------------------
     354     REAL(r_k),    INTENT     ( IN ) :: thetaeK
     355     REAL(r_k),    INTENT     ( IN ) :: pres
     356     LOGICAL, INTENT ( INOUT )  :: flag
     357     REAL(r_k)                       :: tparcel
     358 
     359     REAL(r_k) :: thetaK
     360     REAL(r_k) :: tovtheta
     361     REAL(r_k) :: tcheck
     362     REAL(r_k) :: svpr, svpr2
     363     REAL(r_k) :: smixr, smixr2
     364     REAL(r_k) :: thetae_check, thetae_check2
     365     REAL(r_k) :: tguess_2, correction
     366 
     367     LOGICAL :: found
     368     INTEGER :: iter
     369 
     370     ! Using WRF values
     371     !REAL :: R     ! Dry gas constant
     372     !REAL :: Cp    ! Specific heat for dry air
     373     !REAL :: kappa ! Rd / Cp
     374     !REAL :: Lv    ! Latent heat of vaporization at 0 deg. C
     375     REAL(r_k)                                                :: R, kappa, Lv
     376
     377     R = r_d
     378     Lv = XLV
     379     !R     = 287.04
     380     !Cp    = 1004.67
     381     Kappa = R/Cp
     382     !Lv    = 2.500E+6
     383
     384     !~ Make initial guess for temperature of the parcel
     385     !  ------------------------------------------------
     386     tovtheta = (pres/100000.0)**(r/cp)
     387     tparcel  = thetaeK/exp(lv*.012/(cp*295.))*tovtheta
     388
     389     iter = 1
     390     found = .false.
     391     flag = .false.
     392
     393     DO
     394        IF ( iter > 105 ) EXIT
     395
     396        tguess_2 = tparcel + REAL ( 1 )
     397
     398        svpr   = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) )
     399        smixr  = ( 0.622*svpr ) / ( (pres/100.0)-svpr )
     400        svpr2  = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) )
     401        smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 )
     402
     403        !  ------------------------------------------------------------------ ~!
     404        !~ When this function was orinially written, the final parcel         ~!
     405        !~ temperature check was based off of the parcel temperature and      ~!
     406        !~ not the theta-e it produced.  As there are multiple temperature-   ~!
     407        !~ mixing ratio combinations that can produce a single theta-e value, ~!
     408        !~ we change the check to be based off of the resultant theta-e       ~!
     409        !~ value.  This seems to be the most accurate way of backing out      ~!
     410        !~ temperature from theta-e.                                          ~!
     411        !~                                                                    ~!
     412        !~ Rentschler, April 2010                                             ~!
     413        !  ------------------------------------------------------------------  !
     414
     415        !~ Old way...
     416        !thetaK = thetaeK / EXP (lv * smixr  /(cp*tparcel) )
     417        !tcheck = thetaK * tovtheta
     418
     419        !~ New way
     420        thetae_check  = Thetae ( tparcel,  pres/100., 100., smixr  )
     421        thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 )
     422
     423        !~ Whew doggies - that there is some accuracy...
     424        !IF ( ABS (tparcel-tcheck) < .05) THEN
     425        IF ( ABS (thetaeK-thetae_check) < .001) THEN
     426           found = .true.
     427           flag  = .true.
     428           EXIT
     429        END IF
     430
     431        !~ Old
     432        !tparcel = tparcel + (tcheck - tparcel)*.3
     433
     434        !~ New
     435        correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check )
     436        tparcel = tparcel + correction
     437
     438        iter = iter + 1
     439     END DO
     440
     441     !IF ( .not. found ) THEN
     442     !   print*, "Warning! Thetae to temperature calculation did not converge!"
     443     !   print*, "Thetae ", thetaeK, "Pressure ", pres
     444     !END IF
     445
     446  END FUNCTION The2T
     447
     448  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     449  !~
     450  !~ Name:
     451  !~    VirtualTemperature
     452  !~
     453  !~ Description:
     454  !~    This function returns virtual temperature given temperature ( K )
     455  !~    and mixing ratio.
     456  !~
     457  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     458  FUNCTION VirtualTemperature ( tK, w ) result ( Tv )
     459
     460    IMPLICIT NONE
     461
     462     !~ Variable declaration
     463     real(r_k), intent ( in ) :: tK !~ Temperature
     464     real(r_k), intent ( in ) :: w  !~ Mixing ratio ( kg kg^-1 )
     465     real(r_k)                :: Tv !~ Virtual temperature
     466
     467     Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w )
     468
     469  END FUNCTION VirtualTemperature
     470
     471  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     472  !~
     473  !~ Name:
     474  !~    SaturationMixingRatio
     475  !~
     476  !~ Description:
     477  !~    This function calculates saturation mixing ratio given the
     478  !~    temperature ( K ) and the ambient pressure ( Pa ).  Uses
     479  !~    approximation of saturation vapor pressure.
     480  !~
     481  !~ References:
     482  !~    Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10
     483  !~
     484  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     485  FUNCTION SaturationMixingRatio ( tK, p ) result ( ws )
     486
     487    IMPLICIT NONE
     488
     489    REAL(r_k), INTENT ( IN ) :: tK
     490    REAL(r_k), INTENT ( IN ) :: p
     491    REAL(r_k)                :: ws
     492
     493    REAL(r_k) :: es
     494
     495    es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) )
     496    ws = ( 0.622*es ) / ( (p/100.0)-es )
     497
     498  END FUNCTION SaturationMixingRatio
     499
     500  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     501  !~                                                                     
     502  !~ Name:                                                               
     503  !~    tlcl                                                               
     504  !~                                                                       
     505  !~ Description:                                                           
     506  !~    This function calculates the temperature of a parcel of air would have
     507  !~    if lifed dry adiabatically to it's lifting condensation level (lcl). 
     508  !~                                                                         
     509  !~ References:                                                             
     510  !~    Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22
     511  !~                                                                         
     512  !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
     513  FUNCTION TLCL ( tk, rh )
     514   
     515    IMPLICIT NONE
     516 
     517    REAL(r_k), INTENT ( IN ) :: tK   !~ Temperature ( K )
     518    REAL(r_k), INTENT ( IN ) :: rh   !~ Relative Humidity ( % )
     519    REAL(r_k)                :: tlcl
     520   
     521    REAL(r_k) :: denom, term1, term2
     522
     523    term1 = 1.0 / ( tK - 55.0 )
     524!! Lluis
     525!    IF ( rh > REAL (0) ) THEN
     526    IF ( rh > zeroRK ) THEN
     527      term2 = ( LOG (rh/100.0)  / 2840.0 )
     528    ELSE
     529      term2 = ( LOG (0.001/oneRK) / 2840.0 )
     530    END IF
     531    denom = term1 - term2
     532!! Lluis
     533!    tlcl = ( 1.0 / denom ) + REAL ( 55 )
     534    tlcl = ( oneRK / denom ) + 55*oneRK
     535
     536  END FUNCTION TLCL
     537
     538  FUNCTION var_cape_afwa1D(nz, tk, rhv, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, parcel) RESULT (ostat)
     539! Function to compute cape on a 1D column following implementation in phys/module_diag_afwa.F
     540
     541    IMPLICIT NONE
     542
     543    INTEGER, INTENT(in)                                  :: nz, sfc
     544    REAL(r_k), DIMENSION(nz), INTENT(in)                 :: tk, rhv, p, hgt
     545    REAL(r_k), INTENT(out)                               :: cape, cin, zlfc, plfc, lidx
     546    INTEGER                                              :: ostat
     547    INTEGER, INTENT(in)                                  :: parcel
     548 
     549    ! Local
     550    !~ Derived profile variables
     551    !  -------------------------
     552    REAL(r_k), DIMENSION(nz)                             :: rh, ws, w, dTvK, buoy
     553    REAL(r_k)                                            :: tlclK, plcl, nbuoy, pbuoy
     554 
     555    !~ Source parcel information
     556    !  -------------------------
     557    REAL(r_k)                                            :: srctK, srcrh, srcws, srcw, srcp,          &
     558      srctheta, srcthetaeK
     559    INTEGER                                              :: srclev
     560    REAL(r_k)                                            :: spdiff
     561   
     562    !~ Parcel variables
     563    !  ----------------
     564    REAL(r_k)                                            :: ptK, ptvK, tvK, pw
     565 
     566    !~ Other utility variables
     567    !  -----------------------
     568    INTEGER                                              :: i, j, k
     569    INTEGER                                              :: lfclev
     570    INTEGER                                              :: prcl
     571    INTEGER                                              :: mlev
     572    INTEGER                                              :: lyrcnt
     573    LOGICAL                                              :: flag
     574    LOGICAL                                              :: wflag
     575    REAL(r_k)                                            :: freeze
     576    REAL(r_k)                                            :: pdiff
     577    REAL(r_k)                                            :: pm, pu, pd
     578    REAL(r_k)                                            :: lidxu
     579    REAL(r_k)                                            :: lidxd
     580 
     581    REAL(r_k), PARAMETER                                 :: Rd = r_d
     582    REAL(r_k), PARAMETER                                 :: RUNDEF = -9.999E30
     583
     584!!!!!!! Variables
     585! nz: Number of vertical levels
     586! sfc: Surface level in the profile
     587! tk: Temperature profile [K]
     588! rhv: Relative Humidity profile [1]
     589! rh: Relative Humidity profile [%]
     590! p: Pressure profile [Pa]
     591! hgt: Geopotential height profile [gpm]
     592! cape: CAPE [Jkg-1]
     593! cin: CIN [Jkg-1]
     594! zlfc: LFC Height [gpm]
     595! plfc: LFC Pressure [Pa]
     596! lidx: Lifted index
     597!   FROM: https://en.wikipedia.org/wiki/Lifted_index
     598!     lidx >= 6: Very Stable Conditions
     599!     6 > lidx > 1: Stable Conditions, Thunderstorms Not Likely
     600!     0 > lidx > -2: Slightly Unstable, Thunderstorms Possible, With Lifting Mechanism (i.e., cold front, daytime heating, ...)
     601!     -2 > lidx > -6: Unstable, Thunderstorms Likely, Some Severe With Lifting Mechanism
     602!     -6 > lidx: Very Unstable, Severe Thunderstorms Likely With Lifting Mechanism
     603! ostat: Function return status (Nonzero is bad)
     604! parcel:
     605!   Most Unstable = 1 (default)
     606!   Mean layer = 2
     607!   Surface based = 3
     608!~ Derived profile variables
     609!  -------------------------
     610! ws: Saturation mixing ratio
     611! w: Mixing ratio
     612! dTvK: Parcel / ambient Tv difference
     613! buoy: Buoyancy
     614! tlclK: LCL temperature [K]
     615! plcl: LCL pressure [Pa]
     616! nbuoy: Negative buoyancy
     617! pbuoy: Positive buoyancy
     618 
     619!~ Source parcel information
     620!  -------------------------
     621! srctK: Source parcel temperature [K]
     622! srcrh: Source parcel rh [%]
     623! srcws: Source parcel sat. mixing ratio
     624! srcw: Source parcel mixing ratio
     625! srcp: Source parcel pressure [Pa]
     626! srctheta: Source parcel theta [K]
     627! srcthetaeK: Source parcel theta-e [K]
     628! srclev: Level of the source parcel
     629! spdiff: Pressure difference
     630   
     631!~ Parcel variables
     632!  ----------------
     633! ptK: Parcel temperature [K]
     634! ptvK: Parcel virtual temperature [K]
     635! tvK: Ambient virtual temperature [K]
     636! pw: Parcel mixing ratio
     637 
     638!~ Other utility variables
     639!  -----------------------
     640! lfclev: Level of LFC
     641! prcl: Internal parcel type indicator
     642! mlev: Level for ML calculation
     643! lyrcnt: Number of layers in mean layer
     644! flag: Dummy flag
     645! wflag: Saturation flag
     646! freeze: Water loading multiplier
     647! pdiff: Pressure difference between levs
     648! pm, pu, pd: Middle, upper, lower pressures
     649! lidxu: Lifted index at upper level
     650! lidxd: Lifted index at lower level
     651
     652    fname = 'var_cape_afwa' 
     653
     654    !~ Initialize variables
     655    !  --------------------
     656    rh = rhv*100.
     657    ostat = 0
     658    CAPE = zeroRK
     659    CIN = zeroRK
     660    ZLFC = RUNDEF
     661    PLFC = RUNDEF
     662 
     663    !~ Look for submitted parcel definition
     664    !~ 1 = Most unstable
     665    !~ 2 = Mean layer
     666    !~ 3 = Surface based
     667    !  -------------------------------------
     668    IF ( parcel > 3 .or. parcel < 1 ) THEN
     669       prcl = 1
     670    ELSE
     671       prcl =  parcel
     672    END IF
     673 
     674    !~ Initalize our parcel to be (sort of) surface based.  Because of
     675    !~ issues we've been observing in the WRF model, specifically with
     676    !~ excessive surface moisture values at the surface, using a true
     677    !~ surface based parcel is resulting a more unstable environment
     678    !~ than is actually occuring.  To address this, our surface parcel
     679    !~ is now going to be defined as the parcel between 25-50 hPa
     680    !~ above the surface. UPDATE - now that this routine is in WRF,
     681    !~ going to trust surface info. GAC 20140415
     682    !  ----------------------------------------------------------------
     683 
     684    !~ Compute mixing ratio values for the layer
     685    !  -----------------------------------------
     686    DO k = sfc, nz
     687      ws  ( k )   = SaturationMixingRatio ( tK(k), p(k) )
     688      w   ( k )   = ( rh(k)/100.0 ) * ws ( k )
     689    END DO
     690 
     691    srclev      = sfc
     692    srctK       = tK    ( sfc )
     693    srcrh       = rh    ( sfc )
     694    srcp        = p     ( sfc )
     695    srcws       = ws    ( sfc )
     696    srcw        = w     ( sfc )
     697    srctheta    = Theta ( tK(sfc), p(sfc)/100.0 )
     698   
     699      !~ Compute the profile mixing ratio.  If the parcel is the MU parcel,
     700      !~ define our parcel to be the most unstable parcel within the lowest
     701      !~ 180 mb.
     702      !  -------------------------------------------------------------------
     703      mlev = sfc + 1
     704      DO k = sfc + 1, nz
     705   
     706         !~ Identify the last layer within 100 hPa of the surface
     707         !  -----------------------------------------------------
     708         pdiff = ( p (sfc) - p (k) ) / REAL ( 100 )
     709         IF ( pdiff <= REAL (100) ) mlev = k
     710
     711         !~ If we've made it past the lowest 180 hPa, exit the loop
     712         !  -------------------------------------------------------
     713         IF ( pdiff >= REAL (180) ) EXIT
     714
     715         IF ( prcl == 1 ) THEN
     716            !IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
     717            IF ( (w(k) > srcw) ) THEN
     718               srctheta = Theta ( tK(k), p(k)/100.0 )
     719               srcw = w ( k )
     720               srclev  = k
     721               srctK   = tK ( k )
     722               srcrh   = rh ( k )
     723               srcp    = p  ( k )
     724            END IF
     725         END IF
     726   
     727      END DO
     728   
     729      !~ If we want the mean layer parcel, compute the mean values in the
     730      !~ lowest 100 hPa.
     731      !  ----------------------------------------------------------------
     732      lyrcnt =  mlev - sfc + 1
     733      IF ( prcl == 2 ) THEN
     734   
     735         srclev   = sfc
     736         srctK    = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt )
     737         srcw     = SUM ( w  (sfc:mlev) ) / REAL ( lyrcnt )
     738         srcrh    = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt )
     739         srcp     = SUM ( p  (sfc:mlev) ) / REAL ( lyrcnt )
     740         srctheta = Theta ( srctK, srcp/100. )
     741   
     742      END IF
     743   
     744      srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw )
     745   
     746      !~ Calculate temperature and pressure of the LCL
     747      !  ---------------------------------------------
     748      tlclK = TLCL ( tK(srclev), rh(srclev) )
     749      plcl  = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) )
     750   
     751      !~ Now lift the parcel
     752      !  -------------------
     753   
     754      buoy  = REAL ( 0 )
     755      pw    = srcw
     756      wflag = .false.
     757      DO k  = srclev, nz
     758         IF ( p (k) <= plcl ) THEN
     759   
     760            !~ The first level after we pass the LCL, we're still going to
     761            !~ lift the parcel dry adiabatically, as we haven't added the
     762            !~ the required code to switch between the dry adiabatic and moist
     763            !~ adiabatic cooling.  Since the dry version results in a greater
     764            !~ temperature loss, doing that for the first step so we don't over
     765            !~ guesstimate the instability.
     766            !  ----------------------------------------------------------------
     767   
     768            IF ( wflag ) THEN
     769               flag  = .false.
     770   
     771               !~ Above the LCL, our parcel is now undergoing moist adiabatic
     772               !~ cooling.  Because of the latent heating being undergone as
     773               !~ the parcel rises above the LFC, must iterative solve for the
     774               !~ parcel temperature using equivalant potential temperature,
     775               !~ which is conserved during both dry adiabatic and
     776               !~ pseudoadiabatic displacements.
     777               !  --------------------------------------------------------------
     778               ptK   = The2T ( srcthetaeK, p(k), flag )
     779   
     780               !~ Calculate the parcel mixing ratio, which is now changing
     781               !~ as we condense moisture out of the parcel, and is equivalent
     782               !~ to the saturation mixing ratio, since we are, in theory, at
     783               !~ saturation.
     784               !  ------------------------------------------------------------
     785               pw = SaturationMixingRatio ( ptK, p(k) )
     786   
     787               !~ Now we can calculate the virtual temperature of the parcel
     788               !~ and the surrounding environment to assess the buoyancy.
     789               !  ----------------------------------------------------------
     790               ptvK  = VirtualTemperature ( ptK, pw )
     791               tvK   = VirtualTemperature ( tK (k), w (k) )
     792   
     793               !~ Modification to account for water loading
     794               !  -----------------------------------------
     795               freeze = 0.033 * ( 263.15 - pTvK )
     796               IF ( freeze > 1.0 ) freeze = 1.0
     797               IF ( freeze < 0.0 ) freeze = 0.0
     798   
     799               !~ Approximate how much of the water vapor has condensed out
     800               !~ of the parcel at this level
     801               !  ---------------------------------------------------------
     802               freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7
     803   
     804               pTvK = pTvK - pTvK * ( srcw - pw ) + freeze
     805               dTvK ( k ) = ptvK - tvK
     806               buoy ( k ) = g * ( dTvK ( k ) / tvK )
     807   
     808            ELSE
     809   
     810               !~ Since the theta remains constant whilst undergoing dry
     811               !~ adiabatic processes, can back out the parcel temperature
     812               !~ from potential temperature below the LCL
     813               !  --------------------------------------------------------
     814               ptK   = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
     815   
     816               !~ Grab the parcel virtual temperture, can use the source
     817               !~ mixing ratio since we are undergoing dry adiabatic cooling
     818               !  ----------------------------------------------------------
     819               ptvK  = VirtualTemperature ( ptK, srcw )
     820   
     821               !~ Virtual temperature of the environment
     822               !  --------------------------------------
     823               tvK   = VirtualTemperature ( tK (k), w (k) )
     824   
     825               !~ Buoyancy at this level
     826               !  ----------------------
     827               dTvK ( k ) = ptvK - tvK
     828               buoy ( k ) = g * ( dtvK ( k ) / tvK )
     829   
     830               wflag = .true.
     831   
     832            END IF
     833   
     834         ELSE
     835   
     836            !~ Since the theta remains constant whilst undergoing dry
     837            !~ adiabatic processes, can back out the parcel temperature
     838            !~ from potential temperature below the LCL
     839            !  --------------------------------------------------------
     840            ptK   = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
     841   
     842            !~ Grab the parcel virtual temperture, can use the source
     843            !~ mixing ratio since we are undergoing dry adiabatic cooling
     844            !  ----------------------------------------------------------
     845            ptvK  = VirtualTemperature ( ptK, srcw )
     846   
     847            !~ Virtual temperature of the environment
     848            !  --------------------------------------
     849            tvK   = VirtualTemperature ( tK (k), w (k) )
     850   
     851            !~ Buoyancy at this level
     852            !  ---------------------
     853            dTvK ( k ) = ptvK - tvK
     854            buoy ( k ) = g * ( dtvK ( k ) / tvK )
     855   
     856         END IF
     857
     858         !~ Chirp
     859         !  -----
     860  !          WRITE ( *,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw*1000.,ptvK,tvK,buoy(k)
     861   
     862      END DO
     863   
     864      !~ Add up the buoyancies, find the LFC
     865      !  -----------------------------------
     866      flag   = .false.
     867      lfclev = -1
     868      nbuoy  = REAL ( 0 )
     869      pbuoy = REAL ( 0 )
     870      DO k = sfc + 1, nz
     871         IF ( tK (k) < 253.15 ) EXIT
     872         CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
     873         CIN  = CIN  + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
     874   
     875         !~ If we've already passed the LFC
     876         !  -------------------------------
     877         IF ( flag .and. buoy (k) > REAL (0) ) THEN
     878            pbuoy = pbuoy + buoy (k)
     879         END IF
     880   
     881         !~ We are buoyant now - passed the LFC
     882         !  -----------------------------------
     883         IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN
     884            flag = .true.
     885            pbuoy = pbuoy + buoy (k)
     886            lfclev = k
     887         END IF
     888   
     889         !~ If we think we've passed the LFC, but encounter a negative layer
     890         !~ start adding it up.
     891         !  ----------------------------------------------------------------
     892         IF ( flag .and. buoy (k) < REAL (0) ) THEN
     893            nbuoy = nbuoy + buoy (k)
     894
     895            !~ If the accumulated negative buoyancy is greater than the
     896            !~ positive buoyancy, then we are capped off.  Got to go higher
     897            !~ to find the LFC. Reset positive and negative buoyancy summations
     898            !  ----------------------------------------------------------------
     899            IF ( ABS (nbuoy) > pbuoy ) THEN
     900               flag   = .false.
     901               nbuoy  = REAL ( 0 )
     902               pbuoy  = REAL ( 0 )
     903               lfclev = -1
     904            END IF
     905         END IF
     906
     907      END DO
     908
     909      !~ Calculate lifted index by interpolating difference between
     910      !~ parcel and ambient Tv to 500mb.
     911      !  ----------------------------------------------------------
     912      DO k = sfc + 1, nz
     913
     914         pm = 50000.
     915         pu = p ( k )
     916         pd = p ( k - 1 )
     917
     918         !~ If we're already above 500mb just set lifted index to 0.
     919         !~ --------------------------------------------------------
     920         IF ( pd .le. pm ) THEN
     921            lidx = zeroRK
     922            EXIT
     923   
     924         ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
     925
     926            !~ Found trapping pressure: up, middle, down.
     927            !~ We are doing first order interpolation. 
     928            !  ------------------------------------------
     929            lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
     930            lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
     931            lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
     932            EXIT
     933
     934         ENDIF
     935
     936      END DO
     937   
     938      !~ Assuming the the LFC is at a pressure level for now
     939      !  ---------------------------------------------------
     940      IF ( lfclev > zeroRK ) THEN
     941         PLFC = p   ( lfclev )
     942         ZLFC = hgt ( lfclev )
     943      END IF
     944   
     945      IF ( PLFC /= PLFC .OR. PLFC < zeroRK ) THEN
     946         PLFC = -oneRK
     947         ZLFC = -oneRK
     948      END IF
     949   
     950      IF ( CAPE /= CAPE ) cape = zeroRK
     951   
     952      IF ( CIN  /= CIN  ) cin  = zeroRK
     953
     954      !~ Chirp
     955      !  -----
     956  !       WRITE ( *,* ) ' CAPE: ', cape, ' CIN:  ', cin
     957  !       WRITE ( *,* ) ' LFC:  ', ZLFC, ' PLFC: ', PLFC
     958  !       WRITE ( *,* ) ''
     959  !       WRITE ( *,* ) ' Exiting buoyancy.'
     960  !       WRITE ( *,* ) ' ==================================== '
     961  !       WRITE ( *,* ) ''
     962   
     963    RETURN
     964
     965  END FUNCTION var_cape_afwa1D
     966
     967! ---- END modified from module_diag_afwa.F ---- !
     968
     969
    117970END MODULE module_ForDiagnosticsVars
Note: See TracChangeset for help on using the changeset viewer.