!! Fortran version of different diagnostics ! L. Fita. LMD May 2016 ! gfortran module_generic.o -c module_ForDiagnostics.F90 ! ! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnostics.F90 MODULE module_ForDiagnosticsVars USE module_definitions USE module_generic USE module_scientific IMPLICIT NONE CONTAINS !!!!!!! Variables ! Cdrag_0: Fuction to compute a first order generic approximation of the drag coefficient ! compute_psl_ptarget4d2: Compute sea level pressure using a target pressure. Similar to the Benjamin ! and Miller (1990). Method found in p_interp.F90 ! compute_tv4d: 4D calculation of virtual temperaure ! SaturationMixingRatio: WRF's AFWA method to compute the saturation mixing ratio ! The2T: WRF's AFWA method to compute the temperature at any pressure level along a saturation adiabat ! by iteratively solving for it from the parcel thetae. ! Theta: WRF's AFWA method to compute potential temperature ! Thetae: WRF's AFWA method to compute equivalent potential temperature ! TLCL: WRF's AFWA method to compute the temperature of a parcel of air would have if lifed dry ! adiabatically to it's lifting condensation level (lcl) ! var_cape_afwa1D: WRF's AFWA method to compute cape, cin, fclp, fclz and li ! var_cllmh: low, medium, high-cloud [0,1] ! var_clt: total cloudiness [0,1] ! var_hur: relative humidity using August-Roche-Magnus approximation [1] ! var_fog_K84: fog and visibility following Kunkel, (1984) ! var_fog_RUC: fog and visibility following RUC method Smirnova, (2000) ! var_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010) ! var_front_R04: Subroutine to compute presence of a front following Rodrigues et al.(2004) ! var_potevap_orPM: potential evapotranspiration following Penman-Monteith formulation implemented in ORCHIDEE ! var_psl_ecmwf: sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa] ! var_range_faces: compute faces [uphill, valleys, downhill] of a monuntain range along a face ! var_rh: Subroutine to compute relative humidity following 'Tetens' equation (T,P) ...' ! var_tws_S11: Function to compute wet bulb temperature after Stull, 2011 ! var_zmla_generic: Subroutine to compute pbl-height following a generic method ! var_zwind: Subroutine to extrapolate the wind at a given height following the 'power law' methodology ! var_zwind_log: Subroutine to extrapolate the wind at a given height following the 'logarithmic law' methodology ! var_zwind_MOtheor: Subroutine of wind extrapolation following Moin-Obukhov theory ! VirtualTemp1D: Function for 1D calculation of virtual temperaure ! VirtualTemperature: WRF's AFWA method to compute virtual temperature ! stabfunc_businger: Fucntion of the stability function after Businger et al. (1971) !!!!!!! Calculations ! compute_clt: Computation of total cloudiness !!! ! Variables !!! FUNCTION var_cllmh(clfra, p, dz) ! Function to compute cllmh on a 1D column 1: low-cloud; 2: medium-cloud; 3: high-cloud [1] IMPLICIT NONE INTEGER, INTENT(in) :: dz REAL(r_k), DIMENSION(dz), INTENT(in) :: clfra, p REAL(r_k), DIMENSION(3) :: var_cllmh ! Local INTEGER :: iz REAL(r_k) :: zclearl, zcloudl, zclearm, zcloudm, & zclearh, zcloudh !!!!!!! Variables ! clfra: cloudfraction as 1D verical-column [1] ! p: pressure values of the column fname = 'var_cllmh' zclearl = oneRK zcloudl = zeroRK zclearm = oneRK zcloudm = zeroRK zclearh = oneRK zcloudh = zeroRK var_cllmh = oneRK DO iz=1, dz IF (p(iz) < prmhc) THEN var_cllmh(3) = var_cllmh(3)*(oneRK-MAX(clfra(iz),zcloudh))/(oneRK-MIN(zcloudh,oneRK-ZEPSEC)) zcloudh = clfra(iz) ELSE IF ( (p(iz) >= prmhc) .AND. (p(iz) < prmlc) ) THEN var_cllmh(2) = var_cllmh(2)*(oneRK-MAX(clfra(iz),zcloudm))/(oneRK-MIN(zcloudm,oneRK-ZEPSEC)) zcloudm = clfra(iz) ELSE IF (p(iz) >= prmlc) THEN var_cllmh(1) = var_cllmh(1)*(oneRK-MAX(clfra(iz),zcloudl))/(oneRK-MIN(zcloudl,oneRK-ZEPSEC)) zcloudl = clfra(iz) ELSE PRINT *,' ' // TRIM(fname) // ': This is weird, pressure:', p(iz), ' Pa fails out!!' PRINT *,' from high, low cloud pressures:', prmhc, ' ,', prmlc,' Pa at z-level:', iz PRINT *,' p_high > p:', prmhc,'> ',p(iz),' Pa' PRINT *,' p_low > p >= p_high:', prmlc,'> ',p(iz),' >=', prmhc,' Pa' PRINT *,' p_low >= p:', prmlc,'>= ',p(iz),' Pa' STOP END IF END DO var_cllmh = oneRK - var_cllmh RETURN END FUNCTION var_cllmh REAL(r_k) FUNCTION var_clt(clfra, dz) ! Function to compute the total cloud following 'newmicro.F90' from LMDZ using 1D vertical ! column values IMPLICIT NONE INTEGER, INTENT(in) :: dz REAL(r_k), DIMENSION(dz), INTENT(in) :: clfra ! Local INTEGER :: iz REAL(r_k) :: zclear, zcloud !!!!!!! Variables ! cfra: 1-column cloud fraction values fname = 'var_clt' zclear = oneRK zcloud = zeroRK DO iz=1,dz zclear = zclear*(oneRK-MAX(clfra(iz),zcloud))/(oneRK-MIN(zcloud,1.-ZEPSEC)) var_clt = oneRK - zclear zcloud = clfra(iz) END DO RETURN END FUNCTION var_clt SUBROUTINE var_psl_ecmwf(PRPRESS, hgt, PTB, PRESBH, PRESBF, psl) ! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier ! method found in LMDZ in phylmd/pppmer.F90 in combination with phylmd/ctsar.F90 ! IMPLICIT ARGUMENTS : CONSTANTS FROM YOMCST,YOMGEM,YOMSTA. ! -------------------- IMPLICIT NONE REAL, INTENT(in) :: PRPRESS, hgt, PTB, PRESBH, PRESBF REAL, INTENT(out) :: psl ! Local REAL :: ghgt, PTSTAR, PT0, ZTSTAR REAL :: ZALPHA, POROG REAL :: ZDTDZSG, ZOROG, ZT0, ZTX, ZTY, ZX, ZY, ZY2 REAL, PARAMETER :: RDTDZ1 = -gammav !!!!!!! Variables ! PRPRESS: Surface pressure [Pa] ! hgt: Terrain height [m] ! PTB: Temperature first half-level [K] ! PRESBH: Pressure first half-level [Pa] ! PRESBF: Pressure second full-level [Pa] ! psl: sea-level pressure fname = 'var_psl_ecmwf' ! Height by gravity POROG = hgt*g !* 1. COMPUTES SURFACE TEMPERATURE !* THEN STANDARD SURFACE TEMPERATURE. ZDTDZSG=-RDTDZ1/g ZALPHA=ZDTDZSG*r_d PTSTAR=PTB*(1.0+ZALPHA*(PRESBH/PRESBF-1.0)) PT0=PTSTAR+ZDTDZSG*POROG !* 2. POST-PROCESS MSL PRESSURE. ! -------------------------- !* 2.1 COMPUTATION OF MODIFIED ALPHA AND TSTAR. ZTX=290.5 ZTY=255.0 IF (PTSTAR < ZTY) THEN ZTSTAR=0.5*(ZTY+PTSTAR) ELSEIF (PTSTAR < ZTX) THEN ZTSTAR=PTSTAR ELSE ZTSTAR=0.5*(ZTX+PTSTAR) ENDIF ZT0=ZTSTAR+ZDTDZSG*POROG IF (ZTX > ZTSTAR .AND. ZT0 > ZTX) THEN ZT0=ZTX ELSEIF (ZTX <= ZTSTAR .AND. ZT0 > ZTSTAR) THEN ZT0=ZTSTAR ELSE ZT0=PT0 ENDIF ZOROG=SIGN(MAX(1.0,ABS(POROG)),POROG) ZALPHA=r_d*(ZT0-ZTSTAR)/ZOROG !* 2.2 COMPUTATION OF MSL PRESSURE. IF (ABS(POROG) >= 0.001) THEN ZX=POROG/(r_d*ZTSTAR) ZY=ZALPHA*ZX ZY2=ZY*ZY psl=PRPRESS*EXP(ZX*(1.0-0.5*ZY+1.0/3.*ZY2)) ELSE psl=PRPRESS ENDIF RETURN END SUBROUTINE var_psl_ecmwf SUBROUTINE compute_psl_ptarget4d2(press, ps, hgt, ta, qv, ptarget, psl, d1, d2, d3, d4) ! Subroutine to compute sea level pressure using a target pressure. Similar to the Benjamin ! and Miller (1990). Method found in p_interp.F90 IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: press, ta, qv REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt REAL(r_k), INTENT(in) :: ptarget REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl ! Local INTEGER :: i, j, it INTEGER :: kin INTEGER :: kupper REAL(r_k) :: dpmin, dp, tbotextrap, & tvbotextrap, virtual ! Exponential related to standard atmosphere lapse rate r_d*gammav/g REAL(r_k), PARAMETER :: expon=r_d*gammav/grav !!!!!!! Variables ! press: Atmospheric pressure [Pa] ! ps: surface pressure [Pa] ! hgt: surface height ! ta: temperature [K] ! qv: water vapor mixing ratio ! dz: number of vertical levels ! psl: sea-level pressure fname = 'compute_psl_ptarget4d2' ! Minimal distance between pressures [Pa] dpmin=1.e4 psl=0. DO i=1,d1 DO j=1,d2 IF (hgt(i,j) /= 0.) THEN DO it=1,d4 ! target pressure to be used for the extrapolation [Pa] (defined in namelist.input) ! ptarget = 70000. default value ! We are below both the ground and the lowest data level. ! First, find the model level that is closest to a "target" pressure ! level, where the "target" pressure is delta-p less that the local ! value of a horizontally smoothed surface pressure field. We use ! delta-p = 150 hPa here. A standard lapse rate temperature profile ! passing through the temperature at this model level will be used ! to define the temperature profile below ground. This is similar ! to the Benjamin and Miller (1990) method, using ! 700 hPa everywhere for the "target" pressure. kupper = 0 loop_kIN: DO kin=d3,1,-1 kupper = kin dp=abs( press(i,j,kin,it) - ptarget ) IF (dp .GT. dpmin) EXIT loop_kIN dpmin=min(dpmin,dp) ENDDO loop_kIN tbotextrap=ta(i,j,kupper,it)*(ps(i,j,it)/ptarget)**expon tvbotextrap=virtualTemp1D(tbotextrap,qv(i,j,kupper,it)) psl(i,j,it) = ps(i,j,it)*((tvbotextrap+gammav*hgt(i,j))/tvbotextrap)**(1/expon) END DO ELSE psl(i,j,:) = ps(i,j,:) END IF END DO END DO RETURN END SUBROUTINE compute_psl_ptarget4d2 SUBROUTINE compute_tv4d(ta,qv,tv,d1,d2,d3,d4) ! 4D calculation of virtual temperaure IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, qv REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out) :: tv ! Variables ! ta: temperature [K] ! qv: mixing ratio [kgkg-1] ! tv: virtual temperature tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv) END SUBROUTINE compute_tv4d FUNCTION VirtualTemp1D (ta,qv) result (tv) ! 1D calculation of virtual temperaure IMPLICIT NONE REAL(r_k), INTENT(in) :: ta, qv REAL(r_k) :: tv ! Variables ! ta: temperature [K] ! qv: mixing ratio [kgkg-1] tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv) END FUNCTION VirtualTemp1D ! ---- BEGIN modified from module_diag_afwa.F ---- ! !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! !~ !~ Name: !~ Theta !~ !~ Description: !~ This function calculates potential temperature as defined by !~ Poisson's equation, given temperature and pressure ( hPa ). !~ !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! FUNCTION Theta ( t, p ) IMPLICIT NONE !~ Variable declaration ! -------------------- REAL(r_k), INTENT ( IN ) :: t REAL(r_k), INTENT ( IN ) :: p REAL(r_k) :: theta ! Using WRF values !REAL :: Rd ! Dry gas constant !REAL :: Cp ! Specific heat of dry air at constant pressure !REAL :: p00 ! Standard pressure ( 1000 hPa ) REAL(r_k) :: Rd, p00 !Rd = 287.04 !Cp = 1004.67 !p00 = 1000.00 Rd = r_d p00 = p1000mb/100. !~ Poisson's equation ! ------------------ theta = t * ( (p00/p)**(Rd/Cp) ) END FUNCTION Theta !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! !~ !~ Name: !~ Thetae !~ !~ Description: !~ This function returns equivalent potential temperature using the !~ method described in Bolton 1980, Monthly Weather Review, equation 43. !~ !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! FUNCTION Thetae ( tK, p, rh, mixr ) IMPLICIT NONE !~ Variable Declarations ! --------------------- REAL(r_k) :: tK ! Temperature ( K ) REAL(r_k) :: p ! Pressure ( hPa ) REAL(r_k) :: rh ! Relative humidity REAL(r_k) :: mixr ! Mixing Ratio ( kg kg^-1) REAL(r_k) :: te ! Equivalent temperature ( K ) REAL(r_k) :: thetae ! Equivalent potential temperature ! Using WRF values !REAL, PARAMETER :: R = 287.04 ! Universal gas constant (J/deg kg) !REAL, PARAMETER :: P0 = 1000.0 ! Standard pressure at surface (hPa) REAL(r_k) :: R, p00, Lv !REAL, PARAMETER :: lv = 2.54*(10**6) ! Latent heat of vaporization ! (J kg^-1) !REAL, PARAMETER :: cp = 1004.67 ! Specific heat of dry air constant ! at pressure (J/deg kg) REAL(r_k) :: tlc ! LCL temperature R = r_d p00 = p1000mb/100. lv = XLV !~ Calculate the temperature of the LCL ! ------------------------------------ tlc = TLCL ( tK, rh ) !~ Calculate theta-e ! ----------------- thetae = (tK * (p00/p)**( (R/Cp)*(1.- ( (.28E-3)*mixr*1000.) ) ) )* & exp( (((3.376/tlc)-.00254))*& (mixr*1000.*(1.+(.81E-3)*mixr*1000.)) ) END FUNCTION Thetae !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! !~ !~ Name: !~ The2T.f90 !~ !~ Description: !~ This function returns the temperature at any pressure level along a !~ saturation adiabat by iteratively solving for it from the parcel !~ thetae. !~ !~ Dependencies: !~ function thetae.f90 !~ !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel ) IMPLICIT NONE !~ Variable Declaration ! -------------------- REAL(r_k), INTENT ( IN ) :: thetaeK REAL(r_k), INTENT ( IN ) :: pres LOGICAL, INTENT ( INOUT ) :: flag REAL(r_k) :: tparcel REAL(r_k) :: thetaK REAL(r_k) :: tovtheta REAL(r_k) :: tcheck REAL(r_k) :: svpr, svpr2 REAL(r_k) :: smixr, smixr2 REAL(r_k) :: thetae_check, thetae_check2 REAL(r_k) :: tguess_2, correction LOGICAL :: found INTEGER :: iter ! Using WRF values !REAL :: R ! Dry gas constant !REAL :: Cp ! Specific heat for dry air !REAL :: kappa ! Rd / Cp !REAL :: Lv ! Latent heat of vaporization at 0 deg. C REAL(r_k) :: R, kappa, Lv R = r_d Lv = XLV !R = 287.04 !Cp = 1004.67 Kappa = R/Cp !Lv = 2.500E+6 !~ Make initial guess for temperature of the parcel ! ------------------------------------------------ tovtheta = (pres/100000.0)**(r/cp) tparcel = thetaeK/exp(lv*.012/(cp*295.))*tovtheta iter = 1 found = .false. flag = .false. DO IF ( iter > 105 ) EXIT tguess_2 = tparcel + REAL ( 1 ) svpr = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) ) smixr = ( 0.622*svpr ) / ( (pres/100.0)-svpr ) svpr2 = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) ) smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 ) ! ------------------------------------------------------------------ ~! !~ When this function was orinially written, the final parcel ~! !~ temperature check was based off of the parcel temperature and ~! !~ not the theta-e it produced. As there are multiple temperature- ~! !~ mixing ratio combinations that can produce a single theta-e value, ~! !~ we change the check to be based off of the resultant theta-e ~! !~ value. This seems to be the most accurate way of backing out ~! !~ temperature from theta-e. ~! !~ ~! !~ Rentschler, April 2010 ~! ! ------------------------------------------------------------------ ! !~ Old way... !thetaK = thetaeK / EXP (lv * smixr /(cp*tparcel) ) !tcheck = thetaK * tovtheta !~ New way thetae_check = Thetae ( tparcel, pres/100., 100., smixr ) thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 ) !~ Whew doggies - that there is some accuracy... !IF ( ABS (tparcel-tcheck) < .05) THEN IF ( ABS (thetaeK-thetae_check) < .001) THEN found = .true. flag = .true. EXIT END IF !~ Old !tparcel = tparcel + (tcheck - tparcel)*.3 !~ New correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check ) tparcel = tparcel + correction iter = iter + 1 END DO !IF ( .not. found ) THEN ! print*, "Warning! Thetae to temperature calculation did not converge!" ! print*, "Thetae ", thetaeK, "Pressure ", pres !END IF END FUNCTION The2T !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! !~ !~ Name: !~ VirtualTemperature !~ !~ Description: !~ This function returns virtual temperature given temperature ( K ) !~ and mixing ratio. !~ !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! FUNCTION VirtualTemperature ( tK, w ) result ( Tv ) IMPLICIT NONE !~ Variable declaration real(r_k), intent ( in ) :: tK !~ Temperature real(r_k), intent ( in ) :: w !~ Mixing ratio ( kg kg^-1 ) real(r_k) :: Tv !~ Virtual temperature Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w ) END FUNCTION VirtualTemperature !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! !~ !~ Name: !~ SaturationMixingRatio !~ !~ Description: !~ This function calculates saturation mixing ratio given the !~ temperature ( K ) and the ambient pressure ( Pa ). Uses !~ approximation of saturation vapor pressure. !~ !~ References: !~ Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10 !~ !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! FUNCTION SaturationMixingRatio ( tK, p ) result ( ws ) IMPLICIT NONE REAL(r_k), INTENT ( IN ) :: tK REAL(r_k), INTENT ( IN ) :: p REAL(r_k) :: ws REAL(r_k) :: es es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) ) ws = ( 0.622*es ) / ( (p/100.0)-es ) END FUNCTION SaturationMixingRatio !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! !~ !~ Name: !~ tlcl !~ !~ Description: !~ This function calculates the temperature of a parcel of air would have !~ if lifed dry adiabatically to it's lifting condensation level (lcl). !~ !~ References: !~ Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22 !~ !!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!! FUNCTION TLCL ( tk, rh ) IMPLICIT NONE REAL(r_k), INTENT ( IN ) :: tK !~ Temperature ( K ) REAL(r_k), INTENT ( IN ) :: rh !~ Relative Humidity ( % ) REAL(r_k) :: tlcl REAL(r_k) :: denom, term1, term2 term1 = 1.0 / ( tK - 55.0 ) !! Lluis ! IF ( rh > REAL (0) ) THEN IF ( rh > zeroRK ) THEN term2 = ( LOG (rh/100.0) / 2840.0 ) ELSE term2 = ( LOG (0.001/oneRK) / 2840.0 ) END IF denom = term1 - term2 !! Lluis ! tlcl = ( 1.0 / denom ) + REAL ( 55 ) tlcl = ( oneRK / denom ) + 55*oneRK END FUNCTION TLCL FUNCTION var_cape_afwa1D(nz, tk, rhv, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, parcel) RESULT (ostat) ! Function to compute cape on a 1D column following implementation in phys/module_diag_afwa.F IMPLICIT NONE INTEGER, INTENT(in) :: nz, sfc REAL(r_k), DIMENSION(nz), INTENT(in) :: tk, rhv, p, hgt REAL(r_k), INTENT(out) :: cape, cin, zlfc, plfc, lidx INTEGER :: ostat INTEGER, INTENT(in) :: parcel ! Local !~ Derived profile variables ! ------------------------- REAL(r_k), DIMENSION(nz) :: rh, ws, w, dTvK, buoy REAL(r_k) :: tlclK, plcl, nbuoy, pbuoy !~ Source parcel information ! ------------------------- REAL(r_k) :: srctK, srcrh, srcws, srcw, srcp, & srctheta, srcthetaeK INTEGER :: srclev REAL(r_k) :: spdiff !~ Parcel variables ! ---------------- REAL(r_k) :: ptK, ptvK, tvK, pw !~ Other utility variables ! ----------------------- INTEGER :: i, j, k INTEGER :: lfclev INTEGER :: prcl INTEGER :: mlev INTEGER :: lyrcnt LOGICAL :: flag LOGICAL :: wflag REAL(r_k) :: freeze REAL(r_k) :: pdiff REAL(r_k) :: pm, pu, pd REAL(r_k) :: lidxu REAL(r_k) :: lidxd REAL(r_k), PARAMETER :: Rd = r_d REAL(r_k), PARAMETER :: RUNDEF = -9.999E30 !!!!!!! Variables ! nz: Number of vertical levels ! sfc: Surface level in the profile ! tk: Temperature profile [K] ! rhv: Relative Humidity profile [1] ! rh: Relative Humidity profile [%] ! p: Pressure profile [Pa] ! hgt: Geopotential height profile [gpm] ! cape: CAPE [Jkg-1] ! cin: CIN [Jkg-1] ! zlfc: LFC Height [gpm] ! plfc: LFC Pressure [Pa] ! lidx: Lifted index ! FROM: https://en.wikipedia.org/wiki/Lifted_index ! lidx >= 6: Very Stable Conditions ! 6 > lidx > 1: Stable Conditions, Thunderstorms Not Likely ! 0 > lidx > -2: Slightly Unstable, Thunderstorms Possible, With Lifting Mechanism (i.e., cold front, daytime heating, ...) ! -2 > lidx > -6: Unstable, Thunderstorms Likely, Some Severe With Lifting Mechanism ! -6 > lidx: Very Unstable, Severe Thunderstorms Likely With Lifting Mechanism ! ostat: Function return status (Nonzero is bad) ! parcel: ! Most Unstable = 1 (default) ! Mean layer = 2 ! Surface based = 3 !~ Derived profile variables ! ------------------------- ! ws: Saturation mixing ratio ! w: Mixing ratio ! dTvK: Parcel / ambient Tv difference ! buoy: Buoyancy ! tlclK: LCL temperature [K] ! plcl: LCL pressure [Pa] ! nbuoy: Negative buoyancy ! pbuoy: Positive buoyancy !~ Source parcel information ! ------------------------- ! srctK: Source parcel temperature [K] ! srcrh: Source parcel rh [%] ! srcws: Source parcel sat. mixing ratio ! srcw: Source parcel mixing ratio ! srcp: Source parcel pressure [Pa] ! srctheta: Source parcel theta [K] ! srcthetaeK: Source parcel theta-e [K] ! srclev: Level of the source parcel ! spdiff: Pressure difference !~ Parcel variables ! ---------------- ! ptK: Parcel temperature [K] ! ptvK: Parcel virtual temperature [K] ! tvK: Ambient virtual temperature [K] ! pw: Parcel mixing ratio !~ Other utility variables ! ----------------------- ! lfclev: Level of LFC ! prcl: Internal parcel type indicator ! mlev: Level for ML calculation ! lyrcnt: Number of layers in mean layer ! flag: Dummy flag ! wflag: Saturation flag ! freeze: Water loading multiplier ! pdiff: Pressure difference between levs ! pm, pu, pd: Middle, upper, lower pressures ! lidxu: Lifted index at upper level ! lidxd: Lifted index at lower level fname = 'var_cape_afwa' !~ Initialize variables ! -------------------- rh = rhv*100. ostat = 0 CAPE = zeroRK CIN = zeroRK ZLFC = RUNDEF PLFC = RUNDEF !~ Look for submitted parcel definition !~ 1 = Most unstable !~ 2 = Mean layer !~ 3 = Surface based ! ------------------------------------- IF ( parcel > 3 .or. parcel < 1 ) THEN prcl = 1 ELSE prcl = parcel END IF !~ Initalize our parcel to be (sort of) surface based. Because of !~ issues we've been observing in the WRF model, specifically with !~ excessive surface moisture values at the surface, using a true !~ surface based parcel is resulting a more unstable environment !~ than is actually occuring. To address this, our surface parcel !~ is now going to be defined as the parcel between 25-50 hPa !~ above the surface. UPDATE - now that this routine is in WRF, !~ going to trust surface info. GAC 20140415 ! ---------------------------------------------------------------- !~ Compute mixing ratio values for the layer ! ----------------------------------------- DO k = sfc, nz ws ( k ) = SaturationMixingRatio ( tK(k), p(k) ) w ( k ) = ( rh(k)/100.0 ) * ws ( k ) END DO srclev = sfc srctK = tK ( sfc ) srcrh = rh ( sfc ) srcp = p ( sfc ) srcws = ws ( sfc ) srcw = w ( sfc ) srctheta = Theta ( tK(sfc), p(sfc)/100.0 ) !~ Compute the profile mixing ratio. If the parcel is the MU parcel, !~ define our parcel to be the most unstable parcel within the lowest !~ 180 mb. ! ------------------------------------------------------------------- mlev = sfc + 1 DO k = sfc + 1, nz !~ Identify the last layer within 100 hPa of the surface ! ----------------------------------------------------- pdiff = ( p (sfc) - p (k) ) / REAL ( 100 ) IF ( pdiff <= REAL (100) ) mlev = k !~ If we've made it past the lowest 180 hPa, exit the loop ! ------------------------------------------------------- IF ( pdiff >= REAL (180) ) EXIT IF ( prcl == 1 ) THEN !IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN IF ( (w(k) > srcw) ) THEN srctheta = Theta ( tK(k), p(k)/100.0 ) srcw = w ( k ) srclev = k srctK = tK ( k ) srcrh = rh ( k ) srcp = p ( k ) END IF END IF END DO !~ If we want the mean layer parcel, compute the mean values in the !~ lowest 100 hPa. ! ---------------------------------------------------------------- lyrcnt = mlev - sfc + 1 IF ( prcl == 2 ) THEN srclev = sfc srctK = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt ) srcw = SUM ( w (sfc:mlev) ) / REAL ( lyrcnt ) srcrh = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt ) srcp = SUM ( p (sfc:mlev) ) / REAL ( lyrcnt ) srctheta = Theta ( srctK, srcp/100. ) END IF srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw ) !~ Calculate temperature and pressure of the LCL ! --------------------------------------------- tlclK = TLCL ( tK(srclev), rh(srclev) ) plcl = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) ) !~ Now lift the parcel ! ------------------- buoy = REAL ( 0 ) pw = srcw wflag = .false. DO k = srclev, nz IF ( p (k) <= plcl ) THEN !~ The first level after we pass the LCL, we're still going to !~ lift the parcel dry adiabatically, as we haven't added the !~ the required code to switch between the dry adiabatic and moist !~ adiabatic cooling. Since the dry version results in a greater !~ temperature loss, doing that for the first step so we don't over !~ guesstimate the instability. ! ---------------------------------------------------------------- IF ( wflag ) THEN flag = .false. !~ Above the LCL, our parcel is now undergoing moist adiabatic !~ cooling. Because of the latent heating being undergone as !~ the parcel rises above the LFC, must iterative solve for the !~ parcel temperature using equivalant potential temperature, !~ which is conserved during both dry adiabatic and !~ pseudoadiabatic displacements. ! -------------------------------------------------------------- ptK = The2T ( srcthetaeK, p(k), flag ) !~ Calculate the parcel mixing ratio, which is now changing !~ as we condense moisture out of the parcel, and is equivalent !~ to the saturation mixing ratio, since we are, in theory, at !~ saturation. ! ------------------------------------------------------------ pw = SaturationMixingRatio ( ptK, p(k) ) !~ Now we can calculate the virtual temperature of the parcel !~ and the surrounding environment to assess the buoyancy. ! ---------------------------------------------------------- ptvK = VirtualTemperature ( ptK, pw ) tvK = VirtualTemperature ( tK (k), w (k) ) !~ Modification to account for water loading ! ----------------------------------------- freeze = 0.033 * ( 263.15 - pTvK ) IF ( freeze > 1.0 ) freeze = 1.0 IF ( freeze < 0.0 ) freeze = 0.0 !~ Approximate how much of the water vapor has condensed out !~ of the parcel at this level ! --------------------------------------------------------- freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7 pTvK = pTvK - pTvK * ( srcw - pw ) + freeze dTvK ( k ) = ptvK - tvK buoy ( k ) = g * ( dTvK ( k ) / tvK ) ELSE !~ Since the theta remains constant whilst undergoing dry !~ adiabatic processes, can back out the parcel temperature !~ from potential temperature below the LCL ! -------------------------------------------------------- ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp) !~ Grab the parcel virtual temperture, can use the source !~ mixing ratio since we are undergoing dry adiabatic cooling ! ---------------------------------------------------------- ptvK = VirtualTemperature ( ptK, srcw ) !~ Virtual temperature of the environment ! -------------------------------------- tvK = VirtualTemperature ( tK (k), w (k) ) !~ Buoyancy at this level ! ---------------------- dTvK ( k ) = ptvK - tvK buoy ( k ) = g * ( dtvK ( k ) / tvK ) wflag = .true. END IF ELSE !~ Since the theta remains constant whilst undergoing dry !~ adiabatic processes, can back out the parcel temperature !~ from potential temperature below the LCL ! -------------------------------------------------------- ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp) !~ Grab the parcel virtual temperture, can use the source !~ mixing ratio since we are undergoing dry adiabatic cooling ! ---------------------------------------------------------- ptvK = VirtualTemperature ( ptK, srcw ) !~ Virtual temperature of the environment ! -------------------------------------- tvK = VirtualTemperature ( tK (k), w (k) ) !~ Buoyancy at this level ! --------------------- dTvK ( k ) = ptvK - tvK buoy ( k ) = g * ( dtvK ( k ) / tvK ) END IF !~ Chirp ! ----- ! WRITE ( *,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw*1000.,ptvK,tvK,buoy(k) END DO !~ Add up the buoyancies, find the LFC ! ----------------------------------- flag = .false. lfclev = -1 nbuoy = REAL ( 0 ) pbuoy = REAL ( 0 ) DO k = sfc + 1, nz IF ( tK (k) < 253.15 ) EXIT CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) ) CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) ) !~ If we've already passed the LFC ! ------------------------------- IF ( flag .and. buoy (k) > REAL (0) ) THEN pbuoy = pbuoy + buoy (k) END IF !~ We are buoyant now - passed the LFC ! ----------------------------------- IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN flag = .true. pbuoy = pbuoy + buoy (k) lfclev = k END IF !~ If we think we've passed the LFC, but encounter a negative layer !~ start adding it up. ! ---------------------------------------------------------------- IF ( flag .and. buoy (k) < REAL (0) ) THEN nbuoy = nbuoy + buoy (k) !~ If the accumulated negative buoyancy is greater than the !~ positive buoyancy, then we are capped off. Got to go higher !~ to find the LFC. Reset positive and negative buoyancy summations ! ---------------------------------------------------------------- IF ( ABS (nbuoy) > pbuoy ) THEN flag = .false. nbuoy = REAL ( 0 ) pbuoy = REAL ( 0 ) lfclev = -1 END IF END IF END DO !~ Calculate lifted index by interpolating difference between !~ parcel and ambient Tv to 500mb. ! ---------------------------------------------------------- DO k = sfc + 1, nz pm = 50000. pu = p ( k ) pd = p ( k - 1 ) !~ If we're already above 500mb just set lifted index to 0. !~ -------------------------------------------------------- IF ( pd .le. pm ) THEN lidx = zeroRK EXIT ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN !~ Found trapping pressure: up, middle, down. !~ We are doing first order interpolation. ! ------------------------------------------ lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp) lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp) lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd) EXIT ENDIF END DO !~ Assuming the the LFC is at a pressure level for now ! --------------------------------------------------- IF ( lfclev > zeroRK ) THEN PLFC = p ( lfclev ) ZLFC = hgt ( lfclev ) END IF IF ( PLFC /= PLFC .OR. PLFC < zeroRK ) THEN PLFC = -oneRK ZLFC = -oneRK END IF IF ( CAPE /= CAPE ) cape = zeroRK IF ( CIN /= CIN ) cin = zeroRK !~ Chirp ! ----- ! WRITE ( *,* ) ' CAPE: ', cape, ' CIN: ', cin ! WRITE ( *,* ) ' LFC: ', ZLFC, ' PLFC: ', PLFC ! WRITE ( *,* ) '' ! WRITE ( *,* ) ' Exiting buoyancy.' ! WRITE ( *,* ) ' ==================================== ' ! WRITE ( *,* ) '' RETURN END FUNCTION var_cape_afwa1D ! ---- END modified from module_diag_afwa.F ---- ! SUBROUTINE var_zmla_generic(dz, qv, tpot, z, topo, zmla) ! Subroutine to compute pbl-height following a generic method ! from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim. ! applied also in Garcia-Diez et al., 2013, QJRMS ! where ! "The technique identifies the ML height as a threshold increase of potential temperature from ! its minimum value within the boundary layer." ! here applied similarly to Garcia-Diez et al. where ! zmla = "...first level where potential temperature exceeds the minimum potential temperature ! reached in the mixed layer by more than 1.5 K" IMPLICIT NONE INTEGER, INTENT(in) :: dz REAL(r_k), DIMENSION(dz), INTENT(in) :: qv, tpot, z REAL(r_k), INTENT(in) :: topo REAL(r_k), INTENT(out) :: zmla ! Local INTEGER :: i INTEGER :: mldlev, bllev REAL(r_k) :: dqvar, tpotmin, refdt !!!!!!! Variables ! qv: water vapour mixing ratio ! tpot: potential temperature [K] ! z: height above sea level [m] ! topo: topographic height [m] ! zmla: boundary layer height [m] fname = 'var_zmla_generic' ! Pecentage of difference of mixing ratio used to determine Mixed layer depth dqvar = 0.1 ! MLD = Mixed layer with no substantial variation of mixing ratio /\qv < 10% ? !PRINT *,' Mixed layer mixing ratios qv[1] lev qv[lev] dqvar% _______' DO mldlev = 2, dz IF (ABS(qv(mldlev)-qv(1))/qv(1) > dqvar ) EXIT ! PRINT *,qv(1), mldlev, qv(mldlev), ABS(qv(mldlev)-qv(1))/qv(1) END DO ! Looking for the minimum potential temperature within the MLD [tpotmin = min(tpot)_0^MLD] tpotmin = MINVAL(tpot(1:mldlev)) ! Change in temperature to determine boundary layer height refdt = 1.5 ! Determine the first level where tpot > tpotmin + 1.5 K !PRINT *,' Mixed layer tpotmin lev tpotmin[lev] dtpot _______' DO bllev = 1, dz IF (ABS(tpot(bllev)-tpotmin) > refdt ) EXIT ! PRINT *,tpotmin, bllev, tpot(bllev), ABS(tpot(bllev)-tpotmin) END DO !PRINT *,' height end MLD:', z(mldlev) !PRINT *,' pbl height:', z(bllev) zmla = z(bllev) - topo RETURN END SUBROUTINE var_zmla_generic SUBROUTINE var_zwind(d1, u, v, z, u10, v10, sa, ca, newz, unewz, vnewz) ! Subroutine to extrapolate the wind at a given height following the 'power law' methodology ! wss[newz] = wss[z1]*(newz/z1)**alpha ! alpha = (ln(wss[z2])-ln(wss[z1]))/(ln(z2)-ln(z1)) ! AFTER: Phd Thesis: ! Benedicte Jourdier. Ressource eolienne en France metropolitaine : methodes d’evaluation du ! potentiel, variabilite et tendances. Climatologie. Ecole Doctorale Polytechnique, 2015. French ! IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL(r_k), DIMENSION(d1), INTENT(in) :: u,v,z REAL(r_k), INTENT(in) :: u10, v10, sa, ca, newz REAL(r_k), INTENT(out) :: unewz, vnewz ! Local INTEGER :: inear REAL(r_k) :: zaground REAL(r_k), DIMENSION(2) :: v1, v2, zz, alpha, uvnewz !!!!!!! Variables ! u,v: vertical wind components [ms-1] ! z: height above surface on half-mass levels [m] ! u10,v10: 10-m wind components [ms-1] ! sa, ca: local sine and cosine of map rotation [1.] ! newz: desired height above grpund of extrapolation ! unewz,vnewz: Wind compoonents at the given height [ms-1] fname = 'var_zwind' !PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______' IF (z(1) < newz ) THEN DO inear = 1,d1-2 ! L. Fita, CIMA. Feb. 2018 !! Choose between extra/inter-polate. Maybe better interpolate? ! Here we extrapolate from two closest lower levels !zaground = z(inear+2) zaground = z(inear+1) !PRINT *, inear, z(inear), newz, z(inear+1), z(inear+2) IF ( zaground >= newz) EXIT END DO ELSE !PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz' inear = d1 - 2 END IF IF (inear == d1-2) THEN ! No vertical pair of levels is below newz, using 10m wind as first value and the first level as the second v1(1) = u10 v1(2) = v10 v2(1) = u(1) v2(2) = v(1) zz(1) = 10. zz(2) = z(1) ELSE v1(1) = u(inear) v1(2) = v(inear) v2(1) = u(inear+1) v2(2) = v(inear+1) zz(1) = z(inear) zz(2) = z(inear+1) END IF ! Computing for each component alpha = (LOG(ABS(v2))-LOG(ABS(v1)))/(LOG(zz(2))-LOG(zz(1))) !PRINT *,' Computing with v1:', v1, ' ms-1 v2:', v2, ' ms-1' !PRINT *,' z1:', zz(1), 'm z2:', zz(2), ' m' !PRINT *,' alhpa u:', alpha(1), ' alpha 2:', alpha(2) uvnewz = v1*(newz/zz(1))**alpha ! Earth-rotation unewz = uvnewz(1)*ca - uvnewz(2)*sa vnewz = uvnewz(1)*sa + uvnewz(2)*ca !PRINT *,' result vz:', uvnewz !STOP RETURN END SUBROUTINE var_zwind SUBROUTINE var_zwind_log(d1, u, v, z, u10, v10, sa, ca, newz, unewz, vnewz) ! Subroutine to extrapolate the wind at a given height following the 'logarithmic law' methodology ! wsz = wss[z2]*(ln(newz)-ln(z0))(ln(z2)-ln(z0)) ! ln(z0) = (ws(z2)*ln(z1)-ws(z1)*ln(z2))/(ws(z2)-ws(z1)) ! AFTER: Phd Thesis: ! Benedicte Jourdier. Ressource eolienne en France metropolitaine : methodes d’evaluation du ! potentiel, variabilite et tendances. Climatologie. Ecole Doctorale Polytechnique, 2015. French ! IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL(r_k), DIMENSION(d1), INTENT(in) :: u,v,z REAL(r_k), INTENT(in) :: u10, v10, sa, ca, newz REAL(r_k), INTENT(out) :: unewz, vnewz ! Local INTEGER :: inear REAL(r_k) :: zaground REAL(r_k), DIMENSION(2) :: v1, v2, zz, logz0, uvnewz !!!!!!! Variables ! u,v: vertical wind components [ms-1] ! z: height above surface on half-mass levels [m] ! u10,v10: 10-m wind components [ms-1] ! sa, ca: local sine and cosine of map rotation [1.] ! newz: desired height above grpund of extrapolation ! unewz,vnewz: Wind compoonents at the given height [ms-1] fname = 'var_zwind_log' !PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______' IF (z(1) < newz ) THEN DO inear = 1,d1-2 ! L. Fita, CIMA. Feb. 2018 !! Choose between extra/inter-polate. Maybe better interpolate? ! Here we extrapolate from two closest lower levels !zaground = z(inear+2) zaground = z(inear+1) !PRINT *, inear, z(inear), newz, z(inear+1), z(inear+2) IF ( zaground >= newz) EXIT END DO ELSE !PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz' inear = d1 - 2 END IF IF (inear == d1-2) THEN ! No vertical pair of levels is below newz, using 10m wind as first value and the first level as the second v1(1) = u10 v1(2) = v10 v2(1) = u(1) v2(2) = v(1) zz(1) = 10. zz(2) = z(1) ELSE v1(1) = u(inear) v1(2) = v(inear) v2(1) = u(inear+1) v2(2) = v(inear+1) zz(1) = z(inear) zz(2) = z(inear+1) END IF ! Computing for each component logz0 = (v2*LOG(zz(1))-v1*LOG(zz(2)))/(v2-v1) uvnewz = v2*(LOG(newz)-logz0)/(LOG(zz(2))-logz0) ! Earth-rotation unewz = uvnewz(1)*ca - uvnewz(2)*sa vnewz = uvnewz(1)*sa + uvnewz(2)*ca !PRINT *,' result vz:', uvnewz !STOP RETURN END SUBROUTINE var_zwind_log SUBROUTINE var_zwind_MOtheor(ust, znt, rmol, u10, v10, sa, ca, newz, uznew, vznew) ! Subroutine of wind extrapolation following Moin-Obukhov theory R. B. Stull, 1988, ! Springer (p376-383) ! Only usefull for newz < 80. m ! Ackonwledgement: M. A. Jimenez, UIB IMPLICIT NONE REAL, INTENT(in) :: ust, znt, rmol, u10, v10, sa, ca REAL, INTENT(in) :: newz REAL, INTENT(out) :: uznew, vznew ! Local REAL :: OL REAL :: stability REAL :: wsz, alpha REAL, DIMENSION(2) :: uvnewz !!!!!!! Variables ! ust: u* in similarity theory [ms-1] ! z0: roughness length [m] !!! L. Fita, CIMA. Feb. 2018 !! NOT SURE if it should be z0 instead? ! znt: thermal time-varying roughness length [m] ! rmol: inverse of Obukhov length [m-1] ! u10: x-component 10-m wind speed [ms-1] ! v10: y-component 10-m wind speed [ms-1] ! sa, ca: local sine and cosine of map rotation [1.] ! fname = 'var_zwind_MOtheor' ! Obukhov Length (using the Boussinesq approximation giving Tv from t2) OL = 1/rmol ! Wind speed at desired height PRINT *,'ust:', ust, 'znt:', znt, 'OL:', OL CALL stabfunc_businger(newz,OL,stability) PRINT *,' z/L:', newz/OL,' stabfunc:', stability, 'log:', LOG(newz/znt), 'log+stability:', LOG(newz/znt) + stability PRINT *,' ust/karman:', ust/karman wsz = ust/karman*( LOG(newz/znt) + stability) PRINT *,' wsz:', wsz ! Without taking into account Ekcman pumping, etc... redistributed by components unsing 10-m wind ! as reference... alpha = ATAN2(v10,u10) uvnewz(1) = wsz*COS(alpha) uvnewz(2) = wsz*SIN(alpha) PRINT *,' uvnewz:', uvnewz ! Earth-rotation uznew = uvnewz(1)*ca - uvnewz(2)*sa vznew = uvnewz(1)*sa + uvnewz(2)*ca PRINT *,' uznew:', uznew, ' vznew:', vznew RETURN END SUBROUTINE var_zwind_MOtheor ! L. Fita, CIMA. Feb. 2018 ! WRF seems to have problems with my functions, let'suse subroutine instead !REAL FUNCTION stabfunc_businger(z,L) SUBROUTINE stabfunc_businger(z,L,stabfunc_busingerv) ! Fucntion of the stability function after Businger et al. (1971), JAS, 28(2), 181–189 IMPLICIT NONE REAL, INTENT(in) :: z,L REAL, INTENT(out) :: stabfunc_busingerv ! Local REAL :: zL, X !!!!!!! Variables ! z: height [m] ! L: Obukhov length [-] fname = 'stabfunc_businger' IF (L /= 0.) THEN zL = z/L ELSE ! Neutral zL = 0. END IF IF (zL > 0.) THEN ! Stable case stabfunc_busingerv = 4.7*z/L ELSE IF (zL < 0.) THEN ! unstable X = (1. - 15.*z/L)**(0.25) !stabfunc_busingerv = -2.*LOG((1.+X)/2.)-LOG((1.+X**2)/2.)+2.*ATAN(X)-piconst/2. stabfunc_busingerv = LOG( ((1.+X**2)/2.)*((1.+X)/2.)**2)-2.*ATAN(X)+piconst/2. ELSE stabfunc_busingerv = 0. END IF RETURN ! END FUNCTION stabfunc_businger END SUBROUTINE stabfunc_businger REAL(r_k) FUNCTION Cdrag_0(ust,uas,vas) ! Fuction to compute a first order generic approximation of the drag coefficient as ! CD = (ust/wss)**2 ! after, Garratt, J.R., 1992.: The Atmospheric Boundary Layer. Cambridge Univ. Press, ! Cambridge, U.K., 316 pp ! Ackonwledgement: M. A. Jimenez, UIB ! IMPLICIT NONE REAL(r_k), INTENT(in) :: ust, uas, vas !!!!!!! Variables ! ust: u* in similarity theory [ms-1] ! uas, vas: x/y-components of wind at 10 m fname = 'Cdrag_0' Cdrag_0 = ust**2/(uas**2+vas**2) END FUNCTION Cdrag_0 SUBROUTINE var_potevap_orPM(rho1, ust, uas, vas, tas, ps, qv1, potevap) ! Subroutine to compute the potential evapotranspiration following Penman-Monteith formulation ! implemented in ORCHIDEE ! potevap = dt*rho1*qc*(q2sat-qv1) IMPLICIT NONE REAL(r_k), INTENT(in) :: rho1, ust, uas, vas, tas, ps, qv1 REAL(r_k), INTENT(out) :: potevap ! Local REAL(r_k) :: q2sat, Cd, qc !!!!!!! Variables ! rho1: atsmophere density at the first layer [kgm-3] ! ust: u* in similarity theory [ms-1] ! uas, vas: x/y-components of 10-m wind [ms-1] ! tas: 2-m atmosphere temperature [K] ! ps: surface pressure [Pa] ! qv1: 1st layer atmospheric mixing ratio [kgkg-1] ! potevap: potential evapo transpiration [kgm-2s-1] fname = 'var_potevap_orPM' ! q2sat: Saturated air at 2m (can be assumed to be q2 == qsfc?) q2sat = SaturationMixingRatio(tas, ps) ! Cd: drag coeffiecient Cd = Cdrag_0(ust, uas, vas) ! qc: surface drag coefficient qc = SQRT(uas**2 + vas**2)*Cd potevap = MAX(zeroRK, rho1*qc*(q2sat - qv1)) END SUBROUTINE var_potevap_orPM SUBROUTINE var_fog_K84(qc, qi, fog, vis) ! Computation of fog (vis < 1km) only computed where qcloud, qice /= 0. ! And visibility following Kunkel, B. A., (1984): Parameterization of droplet terminal velocity and ! extinction coefficient in fog models. J. Climate Appl. Meteor., 23, 34–41. IMPLICIT NONE REAL(r_k), INTENT(in) :: qc, qi INTEGER, INTENT(out) :: fog REAL(r_k), INTENT(out) :: vis ! Local REAL(r_k) :: visc, visi !!!!!!! Variables ! qc: cloud mixing ratio [kgkg-1] ! qi, ice mixing ratio [kgkg-1] ! fog: presence of fog (1: yes, 0: no) ! vis: visibility within fog [km] fname = 'var_fog_K84' IF (qi > nullv .OR. qc > nullv) THEN visc = 100000.*oneRK visi = 100000.*oneRK ! From: Gultepe, 2006, JAM, 45, 1469-1480 IF (qc > nullv) visc = 0.027*(qc*1000.)**(-0.88) IF (qi > nullv) visi = 0.024*(qi*1000.)**(-1.0) vis = MINVAL((/visc, visi/)) IF (vis <= oneRK) THEN fog = 1 ELSE fog = 0 vis = -oneRK END IF ELSE fog = 0 vis = -oneRK END IF END SUBROUTINE var_fog_K84 SUBROUTINE var_fog_RUC(qv, ta, pres, fog, vis) ! Computation of fog (vis < 1km) only computed where qcloud, qice /= 0. ! And visibility following RUC method Smirnova, T. G., S. G. Benjamin, and J. M. Brown, 2000: Case ! study verification of RUC/MAPS fog and visibility forecasts. Preprints, 9 th Conference on ! Aviation, Range, and Aerospace Meteorlogy, AMS, Orlando, FL, Sep. 2000. Paper#2.3, 6 pp. IMPLICIT NONE REAL(r_k), INTENT(in) :: qv, ta, pres INTEGER, INTENT(out) :: fog REAL(r_k), INTENT(out) :: vis ! Local REAL(r_k) :: rh !!!!!!! Variables ! qc: cloud mixing ratio [kgkg-1] ! qi, ice mixing ratio [kgkg-1] ! fog: presence of fog (1: yes, 0: no) ! vis: visibility within fog [km] fname = 'var_fog_RUC' CALL var_hur(ta, pres, qv, rh) ! Avoiding supersaturation rh = MINVAL((/1.,rh/)) IF (rh > 0.3) THEN ! From: Gultepe, I., and G. Isaac, 2006: Visbility versus precipitation rate and relative ! humidity. Preprints, 12th Cloud Physics Conf, Madison, WI, Amer. Meteor. Soc., P2.55. ! [Available online at http://ams.confex.com/ams/Madison2006/techprogram/paper_l13177.htm] vis = 60.*EXP(-2.5*(rh*100.-15.)/80.) IF (vis <= oneRK) THEN fog = 1 ELSE fog = 0 vis = -oneRK END IF ELSE fog = 0 vis = -oneRK END IF END SUBROUTINE var_fog_RUC SUBROUTINE var_fog_FRAML50(qv, ta, pres, fog, vis) ! Computation of fog (vis < 1km) ! And visibility following Gultepe, I. and J.A. Milbrandt, 2010: Probabilistic Parameterizations ! of Visibility Using Observations of Rain Precipitation Rate, Relative Humidity, and Visibility. ! J. Appl. Meteor. Climatol., 49, 36-46, https://doi.org/10.1175/2009JAMC1927.1 ! Interest is focused on a 'general' fog/visibilty approach, thus the fit at 50 % of probability ! is chosen ! Effects from precipitation are not considered IMPLICIT NONE REAL(r_k), INTENT(in) :: qv, ta, pres INTEGER, INTENT(out) :: fog REAL(r_k), INTENT(out) :: vis ! Local REAL(r_k) :: rh !!!!!!! Variables ! qv: mixing ratio in [kgkg-1] ! ta: temperature [K] ! pres: pressure field [Pa] ! rh: relative humidity [1] ! fog: presence of fog (1: yes, 0: no) ! vis: visibility within fog [km] fname = 'var_fog_FRAML50' CALL var_hur(ta, pres, qv, rh) ! Avoiding supersaturation rh = MINVAL((/1.,rh/)) IF (rh > 0.3) THEN vis = -5.19*10.**(-10)*(rh*100.)**5.44+40.10 ! Fog definition (vis <= 1. km) IF (vis <= oneRK) THEN fog = 1 ELSE vis = -oneRK fog = 0 END IF ELSE vis = -oneRK fog = 0 END IF END SUBROUTINE var_fog_FRAML50 SUBROUTINE var_range_faces(d, lon, lat, hgt, dist, filt, newrange, hvalleyrange, hgtmax, ihgtmax, & dhgt, Npeaks, ipeaks, Nvalleys, ivalleys, faces0, faces, Nranges, ranges, rangeshgtmax, & irangeshgtmax) ! Subroutine to compute faces [uphill, valleys, downhill] of a monuntain range along a face IMPLICIT NONE INTEGER, INTENT(in) :: d REAL(r_k), INTENT(in) :: filt, newrange, hvalleyrange REAL(r_k), DIMENSION(d), INTENT(in) :: lon, lat, hgt, dist INTEGER, INTENT(out) :: ihgtmax, Npeaks, Nvalleys, Nranges INTEGER, DIMENSION(d), INTENT(out) :: ipeaks, ivalleys, faces0, faces, & irangeshgtmax INTEGER, DIMENSION(2,d), INTENT(out) :: ranges REAL(r_k), INTENT(out) :: hgtmax REAL(r_k), DIMENSION(d), INTENT(out) :: dhgt, rangeshgtmax ! Local INTEGER :: i, j, j1, k, l, m INTEGER :: iface INTEGER :: Nfaces, Nfaces1, Npeaks1, Nvalleys1 INTEGER :: fbeg, fend REAL(r_k) :: dd INTEGER, DIMENSION(1) :: ihmax INTEGER, DIMENSION(d) :: ddhgt, Ndhgt INTEGER, DIMENSION(d) :: faces1, Ndhgt1, ipeaks1, ivalleys1 REAL(r_k), DIMENSION(d) :: dLl, peaks, valleys, peaks1, valleys1 REAL(r_k), DIMENSION(d) :: sortedpeaks INTEGER, DIMENSION(d) :: isortedpeaks LOGICAL :: firstvalley, rangewithin, peakwithin, & valleywithin !!!!!!! Variables ! lon: longitude [degrees east] ! lat: latitude [degrees north] ! hgt: topograpical height [m] ! dist: distance between grid points [m] ! filt: distance to use to filter the topographic values [m]. used to define: ! the minimum length of a given section ! newrange: distance to start a new mountain range [m] ! hvalleyrange: maximum height of a valley to mark change of range [m] fname = 'var_range_faces' !PRINT *, 'heights:', hgt ! Looking for the maximum height hgtmax = MAXVAL(hgt) ihmax = MAXLOC(hgt) ihgtmax = ihmax(1) !PRINT *, 'height max:', hgtmax, 'location:', ihmax ! Filtering height [NOT needed] ! CALL runmean_F1D(d, hgt, filt, 'progressfill', hgtrunmean) ! range slope dhgt(1:d) = hgt(2:d) - hgt(1:d-1) dhgt = dhgt/dist(1:d) !PRINT *, 'slope:', dhgt ! Initial classification Npeaks = 0 Nvalleys = 0 DO i=1, d-1 IF (dhgt(i) > 0.) THEN IF (hgt(i) < hvalleyrange) THEN faces0(i) = fillValI ELSE faces0(i) = 2 END IF ELSE IF (dhgt(i) < 0.) THEN IF (hgt(i) < hvalleyrange) THEN faces0(i) = fillValI ELSE faces0(i) = -2 END IF ELSE IF (hgt(i) == zeroRK .AND. hgt(i+1) == zeroRK) THEN faces0(i) = fillValI ELSE faces0(i) = 0 END IF END IF ! peaks IF (dhgt(i) > 0. .AND. dhgt(i+1) < 0.) THEN Npeaks = Npeaks + 1 peaks(Npeaks) = hgt(i+1) ipeaks(Npeaks) = i+1 END IF ! valleys IF (dhgt(i) < 0. .AND. dhgt(i+1) > 0.) THEN Nvalleys = Nvalleys + 1 valleys(Nvalleys) = hgt(i+1) ivalleys(Nvalleys) = i+1 END IF END DO !PRINT *, 'faces:', faces0 !PRINT *, Npeaks, ' peaks:', peaks(1:Npeaks), ' ipeak:', ipeaks(1:Npeaks) !PRINT *, Nvalleys, ' valleys:', valleys(1:Nvalleys), ' ivalley:', ivalleys(1:Nvalleys) ! tendency of the slope ddhgt(1) = 1 Nfaces = 1 Ndhgt(Nfaces) = 1 DO i=2, d-1 IF (faces0(i) /= faces0(i-1)) THEN Nfaces = Nfaces + 1 Ndhgt(Nfaces) = 1 ddhgt(i) = ddhgt(i-1) + 1 ELSE Ndhgt(Nfaces) = Ndhgt(Nfaces) + 1 ddhgt(i) = ddhgt(i-1) END IF END DO !PRINT *, Nfaces, ' length faces:', Ndhgt(1:Nfaces) ! Defining quantitiy of ranges within the face ! sorting peaks within the face and defining ranges as maximum peaks distanced > newrage !sortedpeaks = peaks !CALL SortR_K(sortedpeaks, Npeaks) !isortedpeaks = 0 !DO i=1, Npeaks ! DO j=1, Npeaks ! IF (peaks(j) == sortedpeaks(i)) isortedpeaks(i) = ipeaks(j) ! END DO !END DO ranges = -1 Nranges = 0 l = 0 DO i=1,d IF (hgt(i) >= hvalleyrange) THEN IF (Nranges == 0) THEN Nranges = 1 ranges(1,Nranges) = i ELSE IF (ranges(2,Nranges) /= -1) THEN Nranges = Nranges + 1 ranges(1,Nranges) = i END IF END IF ELSE IF ((ranges(1,Nranges) /= -1) .AND. (ranges(2,Nranges) == -1)) ranges(2,Nranges) = i-1 END IF END DO IF (ranges(2,Nranges) == -1) ranges(2,Nranges) = d ! Getting characteristics of ranges irangeshgtmax = -1 rangeshgtmax = -1000. k = 1 DO i=1, Nranges DO j=ranges(1,i), ranges(2,i) IF (hgt(j) > rangeshgtmax(i)) THEN irangeshgtmax(i) = j rangeshgtmax(i) = hgt(j) END IF END DO END DO !PRINT *, Nranges, ' _______' !DO i=1, Nranges ! PRINT *,i,':', ranges(:,i),' max:', rangeshgtmax(i), ' ,', irangeshgtmax(i) !END DO ! Defining valleys as that consecutive grid points below surrounding peaks from the same range and ! below range max k = 1 IF (Npeaks > 1) THEN j = 1 j1 = 2 DO i=2, d IF (i >= ipeaks(j) .AND. i < ipeaks(j1)) THEN IF (hgt(i) < peaks(j) .AND. hgt(i) < peaks(j1)) faces0(i) = 0 IF (i == ipeaks(j) .AND. hgt(i) < rangeshgtmax(k)) faces0(i) = 0 ELSE IF (i == ipeaks(j1)) THEN j = j1 j1 = j1 + 1 IF (j1 > Npeaks) EXIT END IF END DO IF (ipeaks(j1) > ranges(2,k)) k = k+1 END IF ! Using defined ranges for the last refinement. ! Uphills after range maximum peak as valley: if dhgt[i] > 0. and i > rangemax --> face[i] = 0 ! Downhills before range maximum peak as valley: if dhgt[i] < 0. and i < rangemax --> face[i] = 0 ! face values before first range, removed: if dhgt[i] != 0. and i < range[1,1] --> face[i] = fillvalueI ! face values after last range, removed: if dhgt[i] != 0. and i < range[2,Nranges] --> face[i] = fillvalueI ! Uphill valleys as 0.5: if face[i] == 0. and i < rangemax --> face[i] = 0.5 ! downhill valleys as 0.5: if face[i] == 0. and i > rangemax --> face[i] = -0.5 ! range peaks as 0. face[rangemax] = 0. DO i=1,Nranges DO j=ranges(1,i), ranges(2,i) IF (dhgt(j) > 0. .AND. j > irangeshgtmax(i)) faces0(j) = 0 IF (dhgt(j) < 0. .AND. j < irangeshgtmax(i)) faces0(j) = 0 IF (faces0(j) == 0. .AND. j < irangeshgtmax(i)) faces0(j) = 1 IF (faces0(j) == 0. .AND. j > irangeshgtmax(i)) faces0(j) = -1 END DO faces0(irangeshgtmax(i)) = 0 END DO IF (dhgt(j) /= 0. .AND. j < ranges(1,1)) faces0(j) = fillValI IF (dhgt(j) /= 0. .AND. j > ranges(2,Nranges)) faces0(j) = fillValI ! Removing all that grid points outside any range DO j=1, d IF (j < ranges(1,1)) THEN faces0(j) = fillValI ELSE IF (j > ranges(2,Nranges)) THEN faces0(j) = fillValI ELSE rangewithin = .FALSE. DO i=1,Nranges IF (j >= ranges(1,i) .AND. j <= ranges(2,i)) THEN rangewithin = .TRUE. EXIT END IF END DO IF (.NOT.rangewithin) faces0(j) = fillValI END IF END DO ! classification using filt. ! On that section of the hill smaller than the selected filter. Replace values by the ones from ! the previous section (except if it is already valley!) ! e.g.: filt= 3 ! faces0 = 0, 0, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1 -1, -1, 0, 0 ! faces = 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1 -1, -1, 0, 0 !!faces1 = faces0 !!Nfaces1 = 1 !!Npeaks1 = 0 !!Nvalleys1 = 0 !!Ndhgt1(1) = Ndhgt(1) !!DO i=2, Nfaces !! fbeg = SUM(Ndhgt(1:i-1)) !! fend = fbeg + Ndhgt(i) !! peakwithin = .FALSE. !! valleywithin = .FALSE. !! DO j=1, Npeaks !! IF (ipeaks(j) == fbeg) THEN !! k = j !! peakwithin = .TRUE. !! EXIT !! END IF !! END DO !! DO j=1, Nvalleys !! IF (ivalleys(j) == fbeg) THEN !! m = j !! valleywithin = .TRUE. !! EXIT !! END IF !! END DO !! dd = SUM(dist(fbeg:fend)) !! IF (Ndhgt(i) < filt) THEN !! PRINT *, 'Lluis replacing !!!', i, ':', fbeg, fend, '<>', faces0(fbeg-1) !! IF (faces0(fbeg) /= 0) faces1(fbeg:fend) = faces0(fbeg-1) !! Ndhgt1(Nfaces1) = Ndhgt1(Nfaces1) + fend - fbeg + 1 !! peakwithin = .FALSE. !! valleywithin = .FALSE. !! ELSE !! Nfaces1 = Nfaces1 + 1 !! Ndhgt1(Nfaces1) = Ndhgt(i) !! END IF !! IF (peakwithin) THEN !! Npeaks1 = Npeaks1 + 1 !! peaks1(Npeaks1) = peaks(k) !! END IF !! IF (valleywithin) THEN !! Nvalleys1 = Nvalleys1 + 1 !! valleys1(Nvalleys1) = valleys(m) !! END IF !!END DO ! Introducing valleys between peaks lower than the highest peak in the section faces = faces1 IF (Npeaks1 > 1) THEN DO i=1,Npeaks1,2 fbeg = ipeaks1(i) fend = ipeaks1(i+1) faces(fbeg:fend) = 0 END DO END IF ! Removing above sea points DO i=1, d-1 IF (hgt(i) == zeroRK .AND. hgt(i+1) == zeroRK) dhgt(i) = fillVal64 END DO !PRINT *, '# lon[0] lat[1] heights[2] slope[3] faces0[4] faces1[5] faces1[6] ddhgt[7]' !DO i=1, d ! PRINT *, lon(i), lat(i), hgt(i), dhgt(i), faces0(i), faces1(i), faces(i), ddhgt(i) !END DO RETURN END SUBROUTINE var_range_faces SUBROUTINE var_hur(t, press, qv, hur) ! Subroutine to compute relative humidity using August-Roche-Magnus approximation [1] IMPLICIT NONE REAL, INTENT(in) :: t, press, qv REAL, INTENT(out) :: hur ! Local REAL :: tC, es, ws !!!!!!! Variables ! t: temperature [K] ! press: pressure [Pa] ! q: mixing ratio [kgkg-1] ! dz: vertical extension ! hur: relative humidity [1] fname = 'var_hur' ! August - Roche - Magnus formula (Avoiding overflow on last level) tC = t - SVPT0 es = ARM1 * exp(ARM2*tC/(tC+ARM3)) ! Saturated mixing ratio ws = mol_watdry*es/(0.01*press-es) ! Relative humidity hur = qv / ws RETURN END SUBROUTINE var_hur REAL(r_k) FUNCTION var_tws_S11(ta0, hur0) ! Function to compute wet bulb temperature using equation after: ! Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1 IMPLICIT NONE REAL(r_k), INTENT(in) :: ta0, hur0 ! Local REAL :: ta, hur !!!!!!! Variables ! ta0: temperature [K] (does it only make sense if it is at 2 m?) ! hur0: relative humidity [1] (does it only make sense if it is at 2 m?) ! tws: wet bulb temperature [C] fname = 'var_tws_S11' ta = ta0 - SVPT0 hur = hur0*100.*oneRK var_tws_S11 = ta * ATAN(0.151977*SQRT(hur+8.313659)) + ATAN(ta+hur) - ATAN(hur-1.676331) + & 0.00391838*(hur)**(1.5)*ATAN(0.023101*hur) - 4.686035 RETURN END FUNCTION var_tws_S11 SUBROUTINE Svar_tws_S11(ta0, hur0, tws) ! Subroutine to compute wet bulb temperature using equation after: ! Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1 IMPLICIT NONE REAL(r_k), INTENT(in) :: ta0, hur0 REAL(r_k), INTENT(out) :: tws ! Local REAL :: ta, hur !!!!!!! Variables ! ta0: temperature [K] (does it only make sense if it is at 2 m?) ! hur0: relative humidity [1] (does it only make sense if it is at 2 m?) ! tws: wet bulb temperature [C] fname = 'var_tws_S11' ta = ta0 - SVPT0 hur = hur0*100.*oneRK tws = ta * ATAN(0.151977*SQRT(hur+8.313659)) + ATAN(ta+hur) - ATAN(hur-1.676331) + & 0.00391838*(hur)**(1.5)*ATAN(0.023101*hur) - 4.686035 RETURN END SUBROUTINE Svar_tws_S11 SUBROUTINE var_front_R04(dx, dy, dt, tas, uas, vas, ddx, ddy, front) ! Subroutine to compute presence of a front following ! Rodrigues et al.(2004), Rev. Bras. Geofis. 22, 135-151, doi: 10.1590/S0102-261X2004000200004 IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dt REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ddx, ddy REAL(r_k), DIMENSION(dx,dy,dt), INTENT(in) :: tas, uas, vas INTEGER, DIMENSION(dx,dy,dt), INTENT(out) :: front ! Local INTEGER :: i, j, l REAL, DIMENSION(dx,dy,dt) :: dt1uas, dt1vas, dt1tas, dc1wss, d2tas !!!!!!! Variables ! tas: 2-m temperature [K] ! dd[x/y]: real distance between grid points along x and y axes [m] ! uas: 10m eastward wind speed [ms-1] ! vas: 10m northward wind speed [ms-1] ! front: presence of a front in the grid point [0: no, 1: yes] fname = 'var_front_R04' dt1uas = zeroRK dt1vas = zeroRK dt1tas = zeroRK dc1wss = zeroRK d2tas = zeroRK ! 1-time-step derivatives DO l=1,dt-1 dt1tas(:,:,l) = ABS(tas(:,:,l+1) - tas(:,:,l)) dt1vas(:,:,l) = ABS(vas(:,:,l+1) - vas(:,:,l)) END DO ! First order curl DO l=1,dt-1 CALL curl2D_1o(dx,dy,ddx,ddy,uas(:,:,l),vas(:,:,l),dc1wss(:,:,l)) END DO ! 2-time-step centered derivatives DO l=2,dt-1 d2tas(:,:,l) = ABS(tas(:,:,l+1) - tas(:,:,l-1)) END DO front = 0 DO l=1,dt DO i=1,dx DO j=1,dy IF ( (dc1wss(i,j,l) /= zeroRK) .AND. (d2tas(i,j,l) > 0.5) ) front(i,j,l) = 1 END DO END DO END DO RETURN END SUBROUTINE var_front_R04 END MODULE module_ForDiagnosticsVars