!! 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 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_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_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 END MODULE module_ForDiagnosticsVars