Changeset 1769 in lmdz_wrf for trunk/tools/module_ForDiagnostics.f90


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

Re-ordering subroutines/functions to their apropriate module

File:
1 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,    &
Note: See TracChangeset for help on using the changeset viewer.