RPackage

RPackage/NAMESPACE

@@ -1,7 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
+export(cape_f)
+export(cape_r)
 export(conv2D_f)
 export(conv2D_r)
+export(mixing_ratio_from_dewpoint)
+export(saturation_vapour_pressure)
 export(xcorr2D_f)
 export(xcorr2D_r)
 useDynLib(AquaFortR)

RPackage/R/cape_f.R

@@ -0,0 +1,52 @@
+#' @title CAPE using Fortran
+#'
+#' @description  The function calculates the Convective Available Potential
+#'  Energy (CAPE) using R. It returns a vector of CAPE, convective inhibition
+#'  (CIN), pressure at lifted condensation level (LCL),and level of free
+#'  convection (LFC).
+#'
+#' @param t_parcel Parcel temperature [K].
+#' @param dwpt_parcel Parcel dew point temperature [K].
+#' @param mr_parcel Parcel mixing ratio [kg/kg].
+#' @param p_profile Pressure profile from sounding or modelling [hPa].
+#' @param t_profile Temperature profile from sounding or modelling [K].
+#' @param mr_profile Mixing ratio profile from sounding or modelling [kg/kg].
+#' @param vtc logical refers is virtual temperature correction due to Doswell and
+#' Rasmussen (1994).
+#' @return A vector containing CAPE, CIN, p_LCL, p_LFC
+#' @examples
+#' data("radiosonde")
+#' t_profile <- radiosonde$temp + 273.15 # K
+#' dwpt_profile <- radiosonde$dpt + 273.15 # K
+#' p_profile <- radiosonde$pressure # hPa
+#' # Mixing ratio
+#' mr_profile <- mixing_ratio_from_dewpoint(dwpt_profile, p_profile)
+#' cape_f(t_profile[1], dwpt_profile[1], mr_profile[1],
+#'   p_profile, t_profile, mr_profile,
+#'   vtc = TRUE
+#' )
+#' @author Klemens Barfus (Original in Fortran), Ahmed Homoudi (Integration in R)
+#' @useDynLib AquaFortR
+#' @export
+cape_f <- function(t_parcel, dwpt_parcel, mr_parcel,
+                   p_profile, t_profile, mr_profile,
+                   vtc = TRUE) {
+  DIM <- c(
+    length(p_profile),
+    4,
+    as.integer(vtc)
+  )
+  result <- .Call(
+    c_cape_f,
+    as.double(t_parcel),
+    as.double(dwpt_parcel),
+    as.double(mr_parcel),
+    as.double(p_profile),
+    as.double(t_profile),
+    as.double(mr_profile),
+    as.integer(DIM)
+  )
+
+  names(result) <- c("CAPE", "CIN", "p_LCL", "p_LFC")
+  return(result)
+}

RPackage/R/cape_r.R

@@ -0,0 +1,152 @@
+#' @title CAPE using R
+#'
+#' @description  The function calculates the Convective Available Potential
+#'  Energy (CAPE) using R. It returns a vector of CAPE, convective inhibition
+#'  (CIN), pressure at lifted condensation level (LCL),and level of free
+#'  convection (LFC).
+#'
+#' @param t_parcel Parcel temperature [K].
+#' @param dwpt_parcel Parcel dew point temperature [K].
+#' @param mr_parcel Parcel mixing ratio [kg/kg].
+#' @param p_profile Pressure profile from sounding or modelling [hPa].
+#' @param t_profile Temperature profile from sounding or modelling [K].
+#' @param mr_profile Mixing ratio profile from sounding or modelling [kg/kg].
+#' @param vtc logical refers is virtual temperature correction due to Doswell and
+#' Rasmussen (1994).
+#' @return A vector containing CAPE, CIN, p_LCL, p_LFC
+#' @examples
+#' data("radiosonde")
+#' t_profile <- radiosonde$temp + 273.15 # K
+#' dwpt_profile <- radiosonde$dpt + 273.15 # K
+#' p_profile <- radiosonde$pressure # hPa
+#' # Mixing ratio
+#' mr_profile <- mixing_ratio_from_dewpoint(dwpt_profile, p_profile)
+#' cape_r(t_profile[1], dwpt_profile[1], mr_profile[1],
+#'   p_profile, t_profile, mr_profile,
+#'   vtc = TRUE
+#' )
+#' @author Klemens Barfus (Original in Fortran), Ahmed Homoudi (Translation to R)
+#' @useDynLib AquaFortR
+#' @export
+cape_r <- function(t_parcel, dwpt_parcel, mr_parcel,
+                   p_profile, t_profile, mr_profile,
+                   vtc = TRUE) {
+  # Constants
+  Rd <- 287.058 # gas constant for dry air [J * kg**-1 * K**-1]
+  Rv <- 461.5 # gas constant for water vapour [J * kg**-1 * K**-1]
+
+
+  t_parcel_tmp <- t_parcel
+  ilevel <- 1
+  nlevel <- length(p_profile)
+  CAPE <- 0.0
+  CIN <- 0.0
+  p_LCL <- NaN
+  p_LFC <- NaN
+
+  # LCL and CIN
+  while ((t_parcel_tmp > dwpt_parcel) & ilevel < nlevel) {
+    # change in pressure
+    dp <- p_profile[ilevel + 1] - p_profile[ilevel]
+    #  change in temperature
+    dt <- calc_dtdp_dry(t_parcel_tmp, p_profile[ilevel]) * dp
+    # new parcel temperature
+    t_parcel_tmp <- t_parcel_tmp + dt
+    # liquid mixing ratio
+    rF <- mr_parcel
+    if (vtc) {
+      # Früh and Wirth Eq. 12: virtual temperature correction to calc CAPE
+      t_parcel_buoyancy <- t_parcel_tmp * (1.0 + (Rv / Rd) * rF) / (1.0 + rF)
+      t_env_buoyancy <- t_profile[ilevel] + (1.0 + (Rv / Rd) * mr_profile[ilevel]) /
+        (1.0 + mr_profile[ilevel])
+    } else {
+      t_parcel_buoyancy <- t_parcel_tmp
+      t_env_buoyancy <- t_profile[ilevel]
+    }
+    CIN_add <- -1.0 * Rd * ((t_parcel_buoyancy - t_env_buoyancy) *
+      (log(p_profile[ilevel]) - log(p_profile[ilevel + 1])))
+
+    CIN <- CIN + CIN_add
+    ilevel <- ilevel + 1
+  }
+  p_LCL <- p_profile[ilevel]
+
+  # LFC
+  if (t_parcel_tmp < t_profile[ilevel]) {
+    while (t_parcel_tmp < t_profile[ilevel] & ilevel < nlevel) {
+      # change in pressure
+      dp <- p_profile[ilevel + 1] - p_profile[ilevel]
+      #
+      pF <- saturation_vapour_pressure(t_parcel_tmp)
+      rF <- calc_rF1(pF, p_profile[ilevel] - pF)
+      #  change in temperature
+      dt <- calc_dtdp_wet(t_parcel_tmp, p_profile[ilevel], rF) * dp
+      # new parcel temperature
+      t_parcel_tmp <- t_parcel_tmp + dt
+      if (vtc) {
+        # Früh and Wirth Eq. 12: virtual temperature correction to calc CAPE
+        t_parcel_buoyancy <- t_parcel_tmp * (1.0 + (Rv / Rd) * rF) / (1.0 + rF)
+        t_env_buoyancy <- t_profile[ilevel] + (1.0 + (Rv / Rd) * mr_profile[ilevel]) /
+          (1.0 + mr_profile[ilevel])
+      } else {
+        t_parcel_buoyancy <- t_parcel_tmp
+        t_env_buoyancy <- t_profile[ilevel]
+      }
+      CIN_add <- -1.0 * Rd * ((t_parcel_buoyancy - t_env_buoyancy) *
+        (log(p_profile[ilevel]) - log(p_profile[ilevel + 1])))
+
+      CIN <- CIN + CIN_add
+      ilevel <- ilevel + 1
+    }
+    if (ilevel < (nlevel - 1)) {
+      p_LFC <- p_profile[ilevel]
+    } else {
+      p_LFC <- NaN
+    }
+  } else {
+    p_LFC <- p_profile[ilevel]
+  }
+
+  # EL and CAPE
+  while (t_parcel_tmp > t_profile[ilevel] & ilevel < nlevel) {
+    pF <- saturation_vapour_pressure(t_parcel_tmp)
+    rF <- calc_rF1(pF, p_profile[ilevel] - pF)
+    if (vtc) {
+      # Früh and Wirth Eq. 12: virtual temperature correction to calc CAPE
+      t_parcel_buoyancy <- t_parcel_tmp * (1.0 + (Rv / Rd) * rF) / (1.0 + rF)
+      t_env_buoyancy <- t_profile[ilevel] + (1.0 + (Rv / Rd) * mr_profile[ilevel]) /
+        (1.0 + mr_profile[ilevel])
+    } else {
+      t_parcel_buoyancy <- t_parcel_tmp
+      t_env_buoyancy <- t_profile[ilevel]
+    }
+    CAPE_add <- -1.0 * Rd * ((t_parcel_buoyancy - t_env_buoyancy) *
+      (log(p_profile[ilevel]) - log(p_profile[ilevel + 1])))
+    CAPE <- CAPE + CAPE_add
+    dp <- p_profile[ilevel + 1] - p_profile[ilevel]
+    rF1 <- 0
+    dt <- calc_dtdp_wet(t_parcel_tmp, p_profile[ilevel], rF1) * dp
+    t_parcel_tmp <- t_parcel_tmp + dt
+    ilevel <- ilevel + 1
+  }
+
+  CAPE <- abs(CAPE)
+  if (ilevel == nlevel) {
+    p_EL <- NaN
+  } else {
+    p_EL <- p_profile[ilevel]
+  }
+  while (ilevel < nlevel) {
+    rF1 <- 0
+    dp <- p_profile[ilevel + 1] - p_profile[ilevel]
+    dt <- calc_dtdp_wet(t_parcel_tmp, p_profile[ilevel], rF1) * dp
+    t_parcel_tmp <- t_parcel_tmp + dt
+    ilevel <- ilevel + 1
+  }
+
+  # END
+  result <- c(CAPE, CIN, p_LCL, p_LFC)
+  names(result) <- c("CAPE", "CIN", "p_LCL", "p_LFC")
+  return(result)
+}
+

RPackage/R/conv2D_f.R

@@ -15,7 +15,7 @@
 #' @useDynLib AquaFortR
 #' @export
 conv2D_f <- function(a, b) {
-  stopifnot(length(dim(a))==2 | length(dim(b))==2)
+  stopifnot(length(dim(a)) == 2 | length(dim(b)) == 2)
   result <- .Call(
     c_conv2d_f,
     as.integer(dim(a)),

RPackage/R/conv2D_r.R

@@ -15,7 +15,7 @@
 #' @useDynLib AquaFortR
 #' @export
 conv2D_r <- function(a, b) {
-  stopifnot(length(dim(a))==2 | length(dim(b))==2)
+  stopifnot(length(dim(a)) == 2 | length(dim(b)) == 2)
   # the full convolution matrix
   conv_row <- nrow(a) + nrow(b) - 1
   conv_col <- ncol(a) + ncol(b) - 1

RPackage/R/radiosonde.R

@@ -1,11 +1,9 @@
 #' @title Radiosonde data
 #' @description radiosonde example data to calculate CAPE.
-#' The data measured at Stuttgart radiosonde station on 2021-06-01 12:00 mid day.
-#' For more information about about the
-#' variables meaning, please refer to this sheet
-#' \url{https://opendata.dwd.de/climate_environment/CDC/help/Abkuerzung_Spaltenname_CDC_20180914.xlsx}
+#' The data measured at Pointe-Noire,  Republic of the Congo, radiosonde station
+#' on 2024.02.02 12Z. The radiosonde data are available online by the University
+#' of Wyoming \url{http://weather.uwyo.edu/upperair/sounding.html/}.
 #'
-#' @format A data frame with 14 columns and 2962 rows
-#' @source German Weather Service (DWD). (2024). Radiosondes data set.
-#' Climate Data Center (CDC). \url{https://opendata.dwd.de/}
+#'
+#' @format A data frame with 7 columns and 66 rows
 "radiosonde"

RPackage/R/utilities.R

@@ -0,0 +1,121 @@
+specific_heat_dry_air <- function(temperature) {
+  # Written by Klemens Barfus. Translated to R by Ahmed Homoudi
+  # temperature [K]
+  t_temp <- temperature - 273.15
+  C_pd <- 1005.60 + (0.017211 * t_temp) + (0.000392 * t_temp**2.0)
+  return(C_pd)
+}
+
+specific_heat_water_vapour <- function(temperature) {
+  # Written by Klemens Barfus. Translated to R by Ahmed Homoudi
+  # temperature [K]
+  t_temp <- temperature - 273.15
+  c_pv <- 1858.0 + (3.820 * 10.0**(-1.0) * t_temp) +
+    (4.220 * 10.0**(-4.0) * t_temp**2.0) -
+    (1.996 * 10.0**(-7.0) * temperature**3.0)
+  return(c_pv)
+}
+
+specific_heat_liquid_water <- function(temperature) {
+  # Written by Klemens Barfus. Translated to R by Ahmed Homoudi
+  # temperature [K]
+  t_temp <- temperature - 273.15
+  c_pw <- 4217.4 - (3.720283 * t_temp) +
+    (0.1412855 * t_temp**2.0) -
+    (2.654387 * 10.0**(-3.0) * t_temp**3.0) +
+    (2.093236 * 10.0**(-5.0) * t_temp**(4.0))
+  return(c_pw)
+}
+#' @title Saturation vapour pressure
+#' @description
+#' Estimation of Saturation vapour pressure [hPa] from temperature [k].
+#' @param temperature in [k]
+#' @param ice TRUE or FALSE, if to consider ice state or not.
+#' @author Klemens Barfus
+#' @export
+saturation_vapour_pressure <- function(temperature, ice = FALSE) {
+  e0 <- 0.611 # kPa
+  Rv <- 461.0 # J K**-1 kg**-1
+  T0 <- 273.15 # K
+  Lv <- 2.50e6 # J kg **-1
+  Ld <- 2.83e6 # J kg **-1
+
+  if (ice) {
+    if (temperature > T0) {
+      L <- Lv
+    } else {
+      L <- Ld
+    }
+  } else {
+    L <- Lv
+  }
+
+  es <- ifelse(temperature > 0,
+    e0 * exp((L / Rv) * ((1 / T0) - (1 / temperature))),
+    0
+  )
+  return(es*10)
+}
+
+calc_rF1 <- function(pF1, p0) {
+  # Written by Klemens Barfus. Translated to R by Ahmed Homoudi
+  # Constants
+  Rd <- 287.058 # gas constant for dry air [J * kg**-1 * K**-1]
+  Rv <- 461.5 # gas constant for water vapour [J * kg**-1 * K**-1]
+  return((Rd * pF1) / (Rv * p0))
+}
+latent_heat_gas_to_liquid <- function(temperature) {
+  # temperature [K]
+  t_temp <- temperature - 273.15
+  latent_heat <- 2500.8 - 2.36 * t_temp + 0.0016 * t_temp**2.0 - 0.00006 * t_temp**3.0
+  return(latent_heat * 1000.0)
+}
+calc_dtdp_dry <- function(temperature, pressure) {
+  # Written by Klemens Barfus. Translated to R by Ahmed Homoudi
+  # temperature [K]
+  # pressure [hPa]
+
+  Rd <- 287.058 # gas constant for dry air [J * kg**-1 * K**-1]
+
+  cp0 <- specific_heat_dry_air(temperature)
+  dtdp <- (temperature * Rd) / (pressure * cp0)
+  return(dtdp)
+}
+calc_dtdp_wet <- function(temperature, pressure, rF) {
+  # Written by Klemens Barfus. Translated to R by Ahmed Homoudi
+  # Constants
+  Rd <- 287.058 # gas constant for dry air [J * kg**-1 * K**-1]
+  Rv <- 461.5 # gas constant for water vapour [J * kg**-1 * K**-1]
+
+  pF1 <- saturation_vapour_pressure(temperature)
+  p0 <- pressure - pF1
+  rF1 <- calc_rF1(pF1, p0)
+  lF1 <- latent_heat_gas_to_liquid(temperature)
+  LLF1 <- pF1 * lF1
+  cp0 <- specific_heat_dry_air(temperature)
+  cp1 <- specific_heat_water_vapour(temperature)
+  cp2 <- specific_heat_liquid_water(temperature)
+  Cp <- cp0 + cp1 * rF1 + cp2 * rF
+  v <- (rF1 * lF1) / pF1 * (1.0 + (Rv / Rd) * rF1) * (LLF1 / (Rv * temperature**2.0))
+  dtdp <- ((rF1 * Rv * temperature / pF1) * (1.0 + (rF1 * lF1 / (Rd * temperature)))) / (Cp + v)
+  return(dtdp)
+}
+#' @title Mixing ratio
+#' @description
+#' Estimation of mixing ratio from pressure [hPa] and dewpoint temperature [K].
+#' @param dewpoint temperature [K]
+#' @param pressure  in [hPa]
+#' @author Ahmed Homoudi
+#' @export
+mixing_ratio_from_dewpoint <- function(dewpoint, pressure) {
+  eipslon <- 0.622 # g g**-1 or kg kg **-1
+  e_in <- saturation_vapour_pressure(dewpoint)
+  r <- (eipslon * e_in) / (pressure - e_in)
+  return(r)
+}
+# dimensionless
+mixing_ratio <- function(P_in, e_in) {
+  eipslon <- 0.622 # g g**-1 or kg kg **-1
+  r <- (eipslon * e_in) / (P_in - e_in)
+  return(r)
+}

RPackage/R/xcorr2D_f.R

@@ -15,7 +15,7 @@
 #' @useDynLib AquaFortR
 #' @export
 xcorr2D_f <- function(a, b) {
-  stopifnot(length(dim(a))==2 | length(dim(b))==2)
+  stopifnot(length(dim(a)) == 2 | length(dim(b)) == 2)
   result <- .Call(
     c_xcorr2d_f,
     as.integer(dim(a)),

RPackage/R/xcorr2D_r.R

@@ -15,7 +15,7 @@
 #' @useDynLib AquaFortR
 #' @export
 xcorr2D_r <- function(a, b) {
-  stopifnot(length(dim(a))==2 | length(dim(b))==2)
+  stopifnot(length(dim(a)) == 2 | length(dim(b)) == 2)
   # the full CC matrix
   cc_row <- nrow(a) + nrow(b) - 1
   cc_col <- ncol(a) + ncol(b) - 1