diff --git a/QuartoBook/AquaFortR_Codes/cape.f90 b/QuartoBook/AquaFortR_Codes/cape.f90 new file mode 100644 index 0000000..1cd9081 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/cape.f90 @@ -0,0 +1,331 @@ +module tools + implicit none + +contains + + function specific_heat_dry_air(temperature) result(C_pd) + ! Written by Klemens Barfus. Modified by Ahmed Homoudi + implicit none + double precision, intent(in) :: temperature + double precision :: t_temp, C_pd + ! temperature should be: -40°C < T < 40^C + ! output is in [J kg^-1 C^-1] + t_temp = temperature - 273.15 + C_pd = 1005.60 + 0.017211*t_temp + 0.000392*t_temp**2.0 + + end function specific_heat_dry_air + + function specific_heat_water_vapour(temperature) result(C_pv) + ! Written by Klemens Barfus. Modified by Ahmed Homoudi + implicit none + double precision, intent(in) :: temperature + double precision :: t_temp, C_pv + + ! Reid, R.C., J.M. Prausnitz, and B.E. Poling (1987) + ! The Properties of Gases and Liquids. 4th ed. McGraw-Hill, 741 pp. + ! output is in J kg^-1 K^-1 + + 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 + + end function + + function specific_heat_liquid_water(temperature) result(C_pw) + ! Written by Klemens Barfus. Modified by Ahmed Homoudi + implicit none + double precision, intent(in) :: temperature + double precision :: t_temp, C_pw + ! temperature [K] + ! output is in J kg^-1 K^-1 + + 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) + + end function + + function saturation_vapour_pressure(temperature, ice_in) result(es) + ! Written by Klemens Barfus. Modified by Ahmed Homoudi + implicit none + ! calculates the saturation vapour pressure in hPa using the Clausius-Claperon equation + double precision, intent(in):: temperature ! [K] + logical, intent(in), optional:: ice_in + + ! keyword ice, indicates if even in case of temperatures lower than 273.15 K es is calculated with + ! respect to liquid water (then ice must not been set) + ! output is in hPa + ! written by K.Barfus 12/2009 + + logical:: ice + + double precision, parameter:: e0 = 0.611 ! [kPa] + double precision, parameter:: T0 = 273.15 ! [K] + double precision, parameter:: Rv = 461.0 ! [J K**-1 kg**-1] gas constant for water vapour + double precision:: L, es + + if (present(ice_in)) then + ice = ice_in + else + ice = .false. + end if + + if (ice .eqv. (.true.)) then + if (temperature .gt. 273.15) then ! water + L = 2.5*10.0**6.0 ! J kg**-1 + else + L = 2.83*10.0**6.0 ! J kg**-1 + end if + else + L = 2.5*10.0**6.0 ! J kg**-1 + end if + + if (temperature .gt. 0) then + es = e0*exp((L/Rv)*(1.0/T0 - 1.0/temperature)) + es = es*10.0 + else + es = 0.0 + end if + + end function saturation_vapour_pressure + + function calc_rF1(pF1, p0) result(res) ! Frueh and Wirth, Eq. 4 + ! Written by Klemens Barfus. Modified by Ahmed Homoudi + implicit none + + double precision, intent(in):: pF1 ! saturation vapour pressure [hPa] + double precision, intent(in):: p0 ! partial pressure of dry air [hPa] + + double precision, parameter:: Rd = 287.058 ! gas constant for dry air [J * kg**-1 * K**-1] + double precision, parameter:: Rv = 461.5 ! gas constant for water vapour [J * kg**-1 * K**-1] + + double precision :: res + + res = (Rd*pF1)/(Rv*p0) + + end function + + function latent_heat_gas_to_liquid(temperature) result(res) + ! latent heat of condensation due to Rogers and Yau in J/kg + ! valid for 248.15 K < T < 313.15 K + ! Written by Klemens Barfus. Modified by Ahmed Homoudi + implicit none + double precision, intent(in):: temperature ! temperature [K] + + double precision:: t_temp, latent_heat, res + + 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 + res = latent_heat*1000.0 + + ! alternative approach + ! calculates the latent heat of condensation (gas -> liquid) due to + ! Fleagle, R.G. and J.A. Businger, (1980) + ! An Introduction to Atmospheric Physics. 2d ed. Academic Press, 432 pp. + ! input + ! T in K + ! output in J kg^-1 K^-1 + ! t_temp = T - 273.15 + ! Lv = (25.00 - 0.02274 * t_temp) * 10.0^5.0 + + end function + + function calc_dtdp_dry(temperature, pressure) result(dtdp) + ! Written by Klemens Barfus. Modified by Ahmed Homoudi + implicit none + double precision, intent(in):: temperature ! temperature [K] + double precision, intent(in):: pressure ! pressure [hPa] + ! output is: + ! [K/hPa] + + double precision, parameter:: Rd = 287.058 ! gas constant for dry air [J * kg**-1 * K**-1] + double precision:: cp0, dtdp + + cp0 = specific_heat_dry_air(temperature) + dtdp = (temperature*Rd)/(pressure*cp0) + + end function + + function calc_dtdp_wet(temperature, pressure, rF) result(dtdp) + ! Written by Klemens Barfus. Modified by Ahmed Homoudi + implicit none + ! not applying mixed-phase model ! + double precision, intent(in):: temperature ! temperature [K] + double precision, intent(in):: pressure ! pressure [hPa] + double precision, intent(in), optional:: rF ! liquid mixing ratio [kg/kg] + ! rF is liquid mixing ratio <- here 0.0 because of an irreversible process + ! output is: + ! [K/hPa] + double precision, parameter:: Rd = 287.058 ! gas constant for dry air [J * kg**-1 * K**-1] + double precision, parameter:: Rv = 461.5 ! gas constant for water vapour [J * kg**-1 * K**-1] + double precision:: pF1, rF1, lF1, LLF1 + double precision:: cp0, cp1, cp2, Cp, v, dtdp, p0 + + pF1 = saturation_vapour_pressure(temperature) ! hPa + p0 = pressure - pF1 ! hPa + rF1 = calc_rF1(pF1, p0) ! saturation mixing ratio in g/g + lF1 = latent_heat_gas_to_liquid(temperature) ! J/kg + LLF1 = pF1*lF1 ! hPa * (J/kg) + cp0 = specific_heat_dry_air(temperature) ! J/(kg*K) + cp1 = specific_heat_water_vapour(temperature) ! J/(kg*K) + cp2 = specific_heat_liquid_water(temperature) ! J/(kg*K) + Cp = cp0 + cp1*rF1 + cp2*rF ! J/(kg*K) + 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) + + end function + +end module tools +!-----------------------------------! +subroutine cape_f(t_parcel, dwpt_parcel, mr_parcel, & + nlevel, p_profile, t_profile, mr_profile, & + vtc, nresult, result) + use :: tools + use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_is_nan + implicit none + + double precision, intent(in) :: t_parcel, dwpt_parcel, mr_parcel + integer :: nlevel, nresult + double precision, intent(in), dimension(nlevel) :: p_profile, t_profile, mr_profile + integer :: vtc +! result CAPE, CIN, p_LCL, p_LFC + double precision, intent(out), dimension(nresult) :: result + +! Constants + double precision, parameter:: Rd = 287.058 ! gas constant for dry air [J * kg**-1 * K**-1] + double precision, parameter:: Rv = 461.5 ! gas constant for water vapour [J * kg**-1 * K**-1] + +! Dummy variables + integer :: ilevel + double precision :: dp, dt, pF, rF + double precision :: t_parcel_buoyancy, t_env_buoyancy, t_parcel_tmp + double precision :: CIN_add, CAPE_add, rF1, p_EL + logical, parameter :: ice = .false. +!double precision, parameter :: NaN = TRANSFER((/Z'00000000', Z'7FF80000'/), 1.0_8) + double precision :: NaN + +! Starting values + NaN = ieee_value(1., ieee_quiet_nan) + t_parcel_tmp = t_parcel + ilevel = 1 + result(1) = 0.0 ! CAPE + result(2) = 0.0 ! CIN + result(3) = NaN ! LCL + result(4) = NaN ! LFC + CIN_add = 0.0 + CAPE_add = 0.0 + +! Find LCL & Estimate CIN + do + if ((t_parcel_tmp .le. dwpt_parcel) .or. (ilevel .ge. nlevel)) exit +!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 .eq. 1) then +! Früh and Wirth Eq. 12: virtual temperature coorection 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) + end if +! Accumlate CIN + CIN_add = -1.0*Rd*((t_parcel_buoyancy - t_env_buoyancy)* & + (log(p_profile(ilevel)) - log(p_profile(ilevel + 1)))) + result(2) = result(2) + CIN_add + ilevel = ilevel + 1 +!print*, ilevel, t_parcel_tmp + end do + result(3) = p_profile(ilevel) ! LCL +! find LFC + if (t_parcel_tmp .le. t_profile(ilevel)) then + do + if ((t_parcel_tmp .ge. t_profile(ilevel)) .or. (ilevel .ge. nlevel)) exit + dp = p_profile(ilevel + 1) - p_profile(ilevel) +! + pF = saturation_vapour_pressure(t_parcel_tmp, ice) +! Saturation mixing ratio + rF = calc_rF1(pF, p_profile(ilevel) - pF) +! + dt = calc_dtdp_wet(t_parcel_tmp, p_profile(ilevel), rF)*dp +! new parcel temperature + t_parcel_tmp = t_parcel_tmp + dt + if (vtc .eq. 1) then +! Früh and Wirth Eq. 12: virtual temperature coorection 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) + end if +! Accumlate CIN + CIN_add = -1.0*Rd*((t_parcel_buoyancy - t_env_buoyancy)* & + (log(p_profile(ilevel)) - log(p_profile(ilevel + 1)))) + result(2) = result(2) + CIN_add + ilevel = ilevel + 1 +!print*, ilevel, t_parcel_tmp + end do + if (ilevel .lt. (nlevel - 1)) then + result(4) = p_profile(ilevel) + else + result(4) = NaN + end if + + else + result(4) = p_profile(ilevel) + end if + +! find EL and estimate CAPE + do + if ((t_parcel_tmp .le. t_profile(ilevel)) .or. (ilevel .ge. nlevel)) exit +! + pF = saturation_vapour_pressure(t_parcel_tmp, ice) +! Saturation mixing ratio + rF = calc_rF1(pF, p_profile(ilevel) - pF) + if (vtc .eq. 1) then +! Früh and Wirth Eq. 12: virtual temperature coorection 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) + end if + CAPE_add = -1.0*Rd*((t_parcel_buoyancy - t_env_buoyancy) & + *(log(p_profile(ilevel)) - log(p_profile(ilevel + 1)))) + + result(1) = result(1) + CAPE_add + dp = p_profile(ilevel + 1) - p_profile(ilevel) + rF1 = 0.0 + dt = calc_dtdp_wet(t_parcel_tmp, p_profile(ilevel), rF1)*dp + t_parcel_tmp = t_parcel_tmp + dt + ilevel = ilevel + 1 +!print*, ilevel, t_parcel_tmp + end do + + result(1) = abs(result(1)) +! + if (ilevel .eq. nlevel) then + p_EL = NaN + else + p_EL = p_profile(ilevel) + end if +! Calculate upper part above EL + do + if (ilevel .ge. nlevel) exit + rF1 = 0.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 +!print*, ilevel, t_parcel_tmp + end do +!print*, p_EL +end subroutine cape_f diff --git a/QuartoBook/AquaFortR_Codes/cape_f.R b/QuartoBook/AquaFortR_Codes/cape_f.R new file mode 100644 index 0000000..4567011 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/cape_f.R @@ -0,0 +1,25 @@ +cape_f0 <- function(t_parcel, dwpt_parcel, mr_parcel, + p_profile, t_profile, mr_profile, + vtc = TRUE) { + dyn.load("fortran/cape_f.so") + + nlevel <- length(p_profile) + nresult <- 4 + + result <- .Fortran("cape_f", + t_parcel = as.double(t_parcel), + dwpt_parcel = as.double(dwpt_parcel), + mr_parcel = as.double(mr_parcel), + nlevel = as.integer(nlevel), + p_profile = as.double(p_profile), + t_profile = as.double(t_profile), + mr_profile = as.double(mr_profile), + vtc = as.integer(vtc), + nresult = as.integer(nresult), + result = as.double(rep(0, nresult)) + )$result + + # END + names(result) <- c("CAPE", "CIN", "p_LCL", "p_LFC") + return(result) +} diff --git a/QuartoBook/AquaFortR_Codes/cape_r.R b/QuartoBook/AquaFortR_Codes/cape_r.R new file mode 100644 index 0000000..8b7fb81 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/cape_r.R @@ -0,0 +1,224 @@ +cape_r0 <- 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) +} + +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) +} diff --git a/QuartoBook/AquaFortR_Codes/conv2D.f90 b/QuartoBook/AquaFortR_Codes/conv2D.f90 new file mode 100644 index 0000000..2855fc9 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/conv2D.f90 @@ -0,0 +1,50 @@ +subroutine conv2d_f(m, n, p, q, k, l, a, b, conv) + implicit none + integer :: m, n, p, q, k, l, i, j + double precision, dimension(m, n) :: a + double precision, dimension(p, q) :: b, b_flipped + double precision, dimension(k, l) :: conv + ! dummy vars + integer :: min_row_shift, min_col_shift + integer :: max_row_shift, max_col_shift + integer :: rows_padded, cols_padded + integer :: icol, irow, iconv, icol2, irow2 + real, allocatable, dimension(:, :) :: padded_a, padded_b + + ! obtain possible shfits + min_row_shift = -1*(p - 1) + max_row_shift = m - 1 + min_col_shift = -1*(q - 1) + max_col_shift = n - 1 + + ! Flip the kernel i.e. B + b_flipped = 0.0 + do i = 1, p + do j = 1, q + b_flipped(p - i + 1, q - j + 1) = b(i, j) + end do + end do + + ! Padded arrray + rows_padded = abs(min_row_shift) + m + abs(max_row_shift) + cols_padded = abs(min_col_shift) + n + abs(max_col_shift) + ! A + allocate (padded_a(rows_padded, cols_padded)) + padded_a = 0.0 + padded_a((abs(min_row_shift) + 1):(abs(min_row_shift) + m), & + (abs(min_col_shift) + 1):(abs(min_col_shift) + n)) = a + + ! B + allocate (padded_b(rows_padded, cols_padded)) + padded_b = 0.0 + do icol = 1, l + do irow = 1, k + iconv = irow + ((icol - 1)*k) + icol2 = icol + q - 1 + irow2 = irow + p - 1 + padded_b(irow:irow2, icol:icol2) = b_flipped + conv(irow, icol) = sum(padded_a*padded_b) + padded_b = 0.0 + end do + end do +end subroutine conv2d_f diff --git a/QuartoBook/AquaFortR_Codes/conv2D_f.R b/QuartoBook/AquaFortR_Codes/conv2D_f.R new file mode 100644 index 0000000..b4c76b9 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/conv2D_f.R @@ -0,0 +1,24 @@ +conv2D_f0 <- function(a, b) { + dyn.load("fortran/conv2D.so") + + # the full convolution matrix + conv_row <- nrow(a) + nrow(b) - 1 + conv_col <- ncol(a) + ncol(b) - 1 + conv <- matrix(1:c(conv_row * conv_col), + byrow = FALSE, ncol = conv_col + ) + + result <- .Fortran("conv2d_f", + m = as.integer(dim(a)[1]), + n = as.integer(dim(a)[2]), + p = as.integer(dim(b)[1]), + q = as.integer(dim(b)[2]), + k = as.integer(conv_row), + l = as.integer(conv_col), + a = as.double(a), + b = as.double(b), + conv = as.double(conv) + )$conv + + return(matrix(result, nrow = conv_row, ncol = conv_col)) +} diff --git a/QuartoBook/AquaFortR_Codes/conv2D_r.R b/QuartoBook/AquaFortR_Codes/conv2D_r.R new file mode 100644 index 0000000..6fc90f3 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/conv2D_r.R @@ -0,0 +1,38 @@ +conv2D_r0 <- function(a, b) { + # the full convolution matrix + conv_row <- nrow(a) + nrow(b) - 1 + conv_col <- ncol(a) + ncol(b) - 1 + conv <- matrix(1:c(conv_row * conv_col), byrow = FALSE, ncol = conv_col) + + # obtain possible shifts + min_row_shift <- -(nrow(b) - 1) + max_row_shift <- (nrow(a) - 1) + min_col_shift <- -(ncol(b) - 1) + max_col_shift <- (ncol(a) - 1) + + # Padded matrix + rows_padded <- abs(min_row_shift) + nrow(a) + abs(max_row_shift) + cols_padded <- abs(min_col_shift) + ncol(a) + abs(max_col_shift) + # a + padded_a <- matrix(0, nrow = rows_padded, ncol = cols_padded) + padded_a[ + (abs(min_row_shift) + 1):(abs(min_row_shift) + nrow(a)), + (abs(min_col_shift) + 1):(abs(min_col_shift) + ncol(a)) + ] <- a + + for (icol in 1:conv_col) { + for (irow in 1:conv_row) { + iconv <- irow + ((icol - 1) * conv_row) + cols <- (icol):(icol + ncol(b) - 1) + rows <- (irow):(irow + nrow(b) - 1) + # b + padded_b <- array(0, dim = c(rows_padded, cols_padded)) + # flip the kernel i.e. b + padded_b[rows, cols] <- b[nrow(b):1, ncol(b):1] + + conv[irow, icol] <- sum(padded_a * padded_b) + } + } + + return(conv) +} diff --git a/QuartoBook/AquaFortR_Codes/xcorr2D.f90 b/QuartoBook/AquaFortR_Codes/xcorr2D.f90 new file mode 100644 index 0000000..8f66f40 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/xcorr2D.f90 @@ -0,0 +1,42 @@ +subroutine xcorr2d_f(m, n, p, q, k, l, a, b, cc) + implicit none + integer :: m, n, p, q, k, l + double precision, dimension(m, n) :: a + double precision, dimension(p, q) :: b + double precision, dimension(k, l) :: cc + ! dummy vars + integer :: min_row_shift, min_col_shift + integer :: max_row_shift, max_col_shift + integer :: rows_padded, cols_padded + integer :: icol, irow, icc, icol2, irow2 + real, allocatable, dimension(:, :) :: padded_a, padded_b + + ! obtain possible shfits + min_row_shift = -1*(p - 1) + max_row_shift = m - 1 + min_col_shift = -1*(q - 1) + max_col_shift = n - 1 + + ! Padded arrray + rows_padded = abs(min_row_shift) + m + abs(max_row_shift) + cols_padded = abs(min_col_shift) + n + abs(max_col_shift) + ! A + allocate (padded_a(rows_padded, cols_padded)) + padded_a = 0.0 + padded_a((abs(min_row_shift) + 1):(abs(min_row_shift) + m), & + (abs(min_col_shift) + 1):(abs(min_col_shift) + n)) = a + + ! B + allocate (padded_b(rows_padded, cols_padded)) + padded_b = 0.0 + do icol = 1, l + do irow = 1, k + icc = irow + ((icol - 1)*k) + icol2 = icol + q - 1 + irow2 = irow + p - 1 + padded_b(irow:irow2, icol:icol2) = b + cc(irow, icol) = sum(padded_a*padded_b) + padded_b = 0.0 + end do + end do +end subroutine xcorr2d_f diff --git a/QuartoBook/AquaFortR_Codes/xcorr2D_f.R b/QuartoBook/AquaFortR_Codes/xcorr2D_f.R new file mode 100644 index 0000000..b2c1ce6 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/xcorr2D_f.R @@ -0,0 +1,24 @@ +xcorr2D_f0 <- function(a, b) { + # Please adjust the path to your setup. In my machine, + # I have a folder called "fortran" containing all f90 files + dyn.load("fortran/xcorr2D.so") + + # the full CC matrix + cc_row <- nrow(a) + nrow(b) - 1 + cc_col <- ncol(a) + ncol(b) - 1 + cc <- matrix(1:c(cc_row * cc_col), byrow = FALSE, ncol = cc_col) + + result <- .Fortran("xcorr2d_f", + m = as.integer(dim(a)[1]), + n = as.integer(dim(a)[2]), + p = as.integer(dim(b)[1]), + q = as.integer(dim(b)[2]), + k = as.integer(cc_row), + l = as.integer(cc_row), + a = as.double(a), + b = as.double(b), + cc = as.double(cc) + )$cc + + return(matrix(result, nrow = nrow(cc), ncol = ncol(cc))) +} diff --git a/QuartoBook/AquaFortR_Codes/xcorr2D_r.R b/QuartoBook/AquaFortR_Codes/xcorr2D_r.R new file mode 100644 index 0000000..dec5969 --- /dev/null +++ b/QuartoBook/AquaFortR_Codes/xcorr2D_r.R @@ -0,0 +1,39 @@ +xcorr2D_r0 <- function(a, b) { + # the full CC matrix + cc_row <- nrow(a) + nrow(b) - 1 + cc_col <- ncol(a) + ncol(b) - 1 + cc <- matrix(1:c(cc_row * cc_col), + byrow = FALSE, ncol = cc_col + ) + + # obtain possible shifts + min_row_shift <- -(nrow(b) - 1) + max_row_shift <- (nrow(a) - 1) + min_col_shift <- -(ncol(b) - 1) + max_col_shift <- (ncol(a) - 1) + + # Padded matrix + rows_padded <- abs(min_row_shift) + nrow(a) + abs(max_row_shift) + cols_padded <- abs(min_col_shift) + ncol(a) + abs(max_col_shift) + # a + padded_a <- matrix(0, nrow = rows_padded, ncol = cols_padded) + padded_a[ + (abs(min_row_shift) + 1):(abs(min_row_shift) + nrow(a)), + (abs(min_col_shift) + 1):(abs(min_col_shift) + ncol(a)) + ] <- a + + for (icol in 1:cc_col) { + for (irow in 1:cc_row) { + icc <- irow + ((icol - 1) * cc_row) + cols <- (icol):(icol + ncol(b) - 1) + rows <- (irow):(irow + nrow(b) - 1) + # b + padded_b <- array(0, dim = c(rows_padded, cols_padded)) + padded_b[rows, cols] <- b + + cc[irow, icol] <- sum(padded_a * padded_b) + } + } + + return(cc) +} diff --git a/QuartoBook/QuartoBook.Rproj b/QuartoBook/QuartoBook.Rproj index 8e3c2eb..ef99afe 100644 --- a/QuartoBook/QuartoBook.Rproj +++ b/QuartoBook/QuartoBook.Rproj @@ -11,3 +11,9 @@ Encoding: UTF-8 RnwWeave: Sweave LaTeX: pdfLaTeX + +PythonType: system +PythonVersion: 3.11.8 +PythonPath: /usr/bin/python3.11 + +SpellingDictionary: en_GB diff --git a/QuartoBook/_quarto.yml b/QuartoBook/_quarto.yml index 04ccb13..9f9d297 100644 --- a/QuartoBook/_quarto.yml +++ b/QuartoBook/_quarto.yml @@ -1,27 +1,78 @@ project: type: book +execute: + cache: false + warning: false + book: - title: "QuartoBook" - author: "Ahmed Homoudi" + title: "AquaFortR: Streamlining Atmospheric Science, Oceanography, Climate, + and Water Research with Fortran-accelerated R" + author: + - name: "Ahmed Homoudi" + orcid: 0000-0001-6906-5986 + email: ahmed.homoudi@tu-dresden.de + affiliations: + - name: | + | TUD Dresden University of Technology + | Institute of Hydrology and Meteorology + date-format: "DD MMMM YYYY" date: "16/08/2023" + #cover-image: + downloads: pdf + reader-mode: true + site-url: "https://ahomoudi.github.io/AquaFortR/book/" + repo-url: "https://github.com/AHomoudi/AquaFortR/" + sharing: [twitter, facebook, linkedin] + open-graph: true + + chapters: - index.qmd + - preface.qmd - intro.qmd - scripts.qmd - package.qmd - summary.qmd - references.qmd + page-navigation: true + + bibliography: references.bib +csl: apa.csl format: html: - theme: cosmo - code-fold: true + theme: + light: + - flatly + dark: + - darkly + code-link: true pdf: - documentclass: scrreprt + documentclass: book + template-partials: + - title.tex toc-depth: 2 + geometry: + - paperwidth=21.0cm + - paperheight=29.7cm + - left=2.5cm + - right=2.5cm + - top=2cm + - bottom=2.5cm + - textheight=26.0cm + - marginparsep=0cm + - marginparwidth=0cm + - footskip=1cm + include-in-header: preamble.tex + linkcolor: blue + citecolor: blue + fontsize: 12pt + + + diff --git a/QuartoBook/apa.csl b/QuartoBook/apa.csl new file mode 100644 index 0000000..1b497e0 --- /dev/null +++ b/QuartoBook/apa.csl @@ -0,0 +1,1859 @@ + + diff --git a/QuartoBook/cover.png b/QuartoBook/cover.png deleted file mode 100644 index e1f5bc6..0000000 Binary files a/QuartoBook/cover.png and /dev/null differ diff --git a/QuartoBook/index.qmd b/QuartoBook/index.qmd index 377f434..78d1e37 100644 --- a/QuartoBook/index.qmd +++ b/QuartoBook/index.qmd @@ -1,8 +1,21 @@ -# Preface {.unnumbered} +#   {.unnumbered .unlisted} + Welcome to the enlightening journey through the pages of **"AquaFortR: Streamlining + Atmospheric Science, Oceanography, Climate, and Water Research with Fortran-accelerated R"**. + In this book, you will gain invaluable insights into seamlessly integrating Fortran + into R scripts by harnessing the power of Fortran. You will acquire + invaluable insights into how to speed up your package using simple C and Fortran codes. + Furthermore, you will accumulate tweaks to further accelerate your scripts or packages, + and further reading will prove to be both advantageous and highly beneficial for further + optimization, and efficiency. -## Acknowledgement {.unnumbered} -This work has been funded by the German Research Foundation (DFG) through the project NFDI4Earth (DFG project no. 460036893, ) within the German National Research Data Infrastructure (NFDI, ). +  -## License {.unnumbered} -This book is licensed under the Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0) License (). +This work has been funded by the German Research Foundation (DFG) through the +project NFDI4Earth (DFG project no. 460036893, ) +within the German National Research Data Infrastructure (NFDI, ). + +  + +This book is licensed under the Creative Commons Attribution-NonCommercial 4.0 +International (CC BY-NC 4.0) License (). diff --git a/QuartoBook/intro.qmd b/QuartoBook/intro.qmd index b4823ec..7abf09b 100644 --- a/QuartoBook/intro.qmd +++ b/QuartoBook/intro.qmd @@ -1,12 +1,158 @@ -# Introduction +# Introduction {#chap-intro} -This chapter briefly introduces the R programming language and how to install R on different operating systems. Furthermore, we will have a look at the installation of RStudio. Finally, we will learn briefly about Fortran. +This chapter briefly introduces the R programming language and how to install R +on different operating systems. Furthermore, we will have a look at the installation +of RStudio. Finally, we will learn briefly about Fortran. -The book is designed to cater to individuals with a foundational understanding of programming, irrespective of their familiarity with a specific programming language. +The book is designed to cater to individuals with a foundational understanding +of programming, irrespective of their familiarity with a specific programming language. + +## R + +R is a powerful and versatile open-source language and environment for statistical +computing and graphics [@RManual]. R was developed to facilitate data manipulation, +exploration, and visualisation, providing various statistical and graphical tools. + + +The R environment is designed for effective data handling, array and matrix +operations, and is a well-developed and simple programming language [@RManual]. +Its syntax is concise and expressive, making it an accessible language for +newbies and seasoned programmers. R can perform tasks either by executing scripts +or interactively. The latter is advantageous for beginners and during the first +stages of script development. The interactive environment is accessible through +command lines or various integrated development environments (IDEs), such [RStudio](https://posit.co/products/open-source/rstudio/). + + +Many factors play an important role in the popularity of R. For example, the +ability to produce graphics with a publication's quality. Simultaneously, the +users maintain full control over customising the plots according to their +preferences and intricate details. Another significant factor is its ability +to be extended to fit the user's demand. The Comprehensive R Archive Network +[(CRAN)](https://cran.r-project.org/) contains an extensive array of libraries +and packages to extend R's functionality for various tasks. With over +[20447](https://cran.r-project.org/web/packages/) available packages, the +CRAN package repository is a testament to R's versatility and adaptability. + + +[Tidyverse](https://www.tidyverse.org/) is the most popular bundle of R packages +for data science. It was developed by Hadley Wickham and his team. The common +shared design among tidyverse packages increases the consistency across functions +and makes each new function or package a little easier to comprehend and utilise [@Wickham2023]. + + +Many packages have been developed for spatial data analysis to address various +challenging tasks. [R-spatial](https://r-spatial.org) provides a rich set of +packages for handling spatial or spatio-temporal data. For example, [sf](https://cran.r-project.org/web/packages/sf/index.html) provides simple features +access in R, which is a standardized method for encoding spatial vector data. The [stars](https://cran.r-project.org/web/packages/stars/index.html) package aims +to handle spatiotemporal arrays. Additionally, the +[terra](https://cran.r-project.org/web/packages/terra/index.html) +package works with spatial data and has the ability to process large datasets on +the disk when loading into memory (RAM) is not feasible. + + +Furthermore, research produces a large data sets; therefore, it is essential to +store them according to the FAIR principles (**F**indability, **A**ccessibility, +**I**nteroperability, and **R**euse of digital assets; @Wilkinson2016). NetCDF +and HDF5 are among the most prominent scientific data formats owing to their +numerous capabilities. The R package +[ncdf4](https://cran.r-project.org/web/packages/ncdf4/index.html) delivers a +high-level R interface to data files written using Unidata's netCDF library. +Additionally, [rhdf5](https://bioconductor.org/packages/release/bioc/html/rhdf5.html) +provides an interface between HDF5 and R. + + + +Regarding the integration of other programming languages, R has diverse interfaces, +which are either a fundamental implementation of R or attainable via another R +package [@Chambers2016]. The fundamental interfaces are `.Call()`, `.C()`, and +`.Fortran()` to C and Fortran. The development of the +[Rcpp](https://cran.r-project.org/web/packages/Rcpp/index.html) package has +revolutionised seamless access to C++. Python is also accessible using the [reticulate](https://cran.r-project.org/web/packages/reticulate/index.html) package. +Finally, the Java interface was granted using the +[rJava](https://cran.r-project.org/web/packages/rJava/index.html) package. + + +In conclusion, R is a valuable tool in the Earth System, as it can effectively +tackle multifarious scientific tasks and address numerous outstanding research +inquiries. -## R Programming lanaguage -## Installing R -## Installing Rstudio ## Fortran -@intro2rbook +Fortran, short for **For**mula **Tran**slation, is one of the oldest high-level +programming languages. It was first developed in the 1950s, and is nonetheless widely +used in scientific and engineering applications. Key features of Fortran include +its ability to efficiently handle arrays and matrices, making it well-suited for +numerical computations. Additionally, Fortran has a simple and straightforward +syntax that makes it easy to learn and use, because it is possible to write +mathematical formulas almost as they written in mathematical texts [@metcalf2018modern]. + +Fortran has been revised multiple times, with the most recent iteration being Fortran [2023](https://wg5-fortran.org/f2023.html). Another important feature of Fortran +is its support for parallel programming, enabling developers to take advantage of +multicore processors and high-performance computing architectures. + +There are frequently numerous good reasons to integrate different programming languages +to achieve tasks. Interoperability with C programming language is a feature that +was introduced with Fortran 95 [@metcalf2018modern]. Given that C is widely used +for system-level programming, many of the other languages include +support for C. Therefore, the C Application Programming Interface (API) can also +be used to connect two non-C languages [@chirila2014introduction]. For instance, +in atmospheric modelling, Fortran is used for its high performance and capcity to +handle large data sets, while C is utilised for its efficiency and control over +memory usage. + +Noteworthy, developing software using Fortran necessitates utilisation of its primitive +procedures and developing from scratch. This is because Fortran is not similar to +scripting languages (i.e. R) that requires a special environment [@Masuda_2020]. +However, Fortran is privileged with its persistent backward compatibility, resulting +in the usability of countless (legacy) codes written decades ago. + + + +Since R was designed to streamline data analysis and statistics +[@wickham2015advanced] and Fortran is renowned for its high performance, it makes +seance to integrate the two languages. It should be noted that both programming +languages present arrays in column-major order, which makes it easier to bridge +without causing confusion. Additionally, R is developed using Fortran, C, and R +programming languages, and it features `.Call` and `.External()` functions that +allows users to utilise compiled code from other R packages. + +Despite the development of newer programming languages, Fortran remains a popular +choice for many scientists and engineers due to its reliability, efficiency, +and ability to handle large amounts of data. + + +## Installation + +**R** + +Installation of R differs according to the operating system: + + - The webpage [here](https://cran.r-project.org/bin/linux/) provides + information on installing R according to the Linux distribution. + - For Windows, The executable can be downloaded from + [here](https://cran.r-project.org/bin/windows/base/). + Additionally, [Previous releases](https://cran.r-project.org/bin/windows/base/old/) + of R for Windows exit. + - For macOS, various releases and versions are accessible + [here](https://cran.r-project.org/bin/macosx/). + +**RStudio** + +As mentioned earlier, RStudio is an IDE for R. Although it is the most prominent, +other IDEs exist, such as Jupyter Notebook and Visual Studio Code. In this book, +RStudio will be the main IDE. To install Rstudio, visit the posit +[page](https://posit.co/download/rstudio-desktop/) to download the suitable installers. + +**Fortran** + +Fortran doesn't require an explicit installation, unlike interpreted languages +such as R. The source code would be translated to the machine language using the +Fortran compiler, and then it can be executed. Therefore, it is important to make +sure that a Fortran compiler, i.e., the [GNU Fortran](https://gcc.gnu.org/fortran/) +(gfortran) compiler, exists in the working machine. + + - In the majority of Linux distributions, the GCC compilers, including gfortran, + come pre-installed. + - For Windows, installing RTools should ensure the existence of gfortran + - For macOS, binaires for gfortran are avaiable [here](https://gcc.gnu.org/wiki/GFortranBinaries). + Furthermore, more information is available on R for macOS [here](https://cran.r-project.org/bin/macosx/tools/) diff --git a/QuartoBook/package.qmd b/QuartoBook/package.qmd index 8cc823c..c841fe8 100644 --- a/QuartoBook/package.qmd +++ b/QuartoBook/package.qmd @@ -4,7 +4,29 @@ This chapter will focus on wrapping the developed routines from the previous cha ## Brief Introduction to R Packages +### devtools package + ## Developing AquaFortR Package -```{r} -1 + 1 -``` + + +### R function + +### Fortran subroutine + +- iso_c_binidng + +### C funtions + - C to Fortran + - R to C + +### Package's files + + + + + + + + + + diff --git a/QuartoBook/preamble.tex b/QuartoBook/preamble.tex new file mode 100644 index 0000000..b7362bb --- /dev/null +++ b/QuartoBook/preamble.tex @@ -0,0 +1,18 @@ +\usepackage[defaultsans]{opensans} + +\usepackage[noblocks]{authblk} + \renewcommand*{\Authsep}{, } + \renewcommand*{\Authand}{, } + \renewcommand*{\Authands}{, } + \renewcommand\Affilfont{\small} + +\usepackage{hyperref} % For hyperlinks in the PDF +\usepackage{xcolor} +\hypersetup{ + colorlinks, + linkcolor={blue!100!black}, + citecolor={blue!100!black}, + urlcolor={blue!100!black} +} +\usepackage{setspace} +\onehalfspacing \ No newline at end of file diff --git a/QuartoBook/preface.qmd b/QuartoBook/preface.qmd new file mode 100644 index 0000000..df301b1 --- /dev/null +++ b/QuartoBook/preface.qmd @@ -0,0 +1,9 @@ +# Preface {.unnumbered} + + + + + + + + diff --git a/QuartoBook/references.bib b/QuartoBook/references.bib index f7eea72..e232d39 100644 --- a/QuartoBook/references.bib +++ b/QuartoBook/references.bib @@ -24,3 +24,173 @@ @online{intro2rbook url = {https://intro2r.com/}, urldate = {2024-01-08} } +@misc{Wang2019, + author = {Chen Wang}, + title = {Kernel learning for visual perception}, + year = {2019}, + doi = {10.32657/10220/47835}, +} + +@book{metcalf2018modern, + title={Modern Fortran Explained: Incorporating Fortran 2018}, + author={Metcalf, M. and Reid, J. and Cohen, M.}, + isbn={9780192539878}, + series={Numerical Mathematics and Scientific Computation}, + url={https://books.google.de/books?id=sB1rDwAAQBAJ}, + year={2018}, + publisher={OUP Oxford} +} + +@article{Warren2016, +author = {Warren, M. A. and Quartly, G. D. and Shutler, J. D. and Miller, P. I. and Yoshikawa, Y.}, +title = {Estimation of ocean surface currents from maximum cross correlation applied to GOCI geostationary satellite remote sensing data over the Tsushima (Korea) Straits}, +journal = {Journal of Geophysical Research: Oceans}, +volume = {121}, +number = {9}, +pages = {6993-7009}, +keywords = {surface currents, remote sensing}, +doi = {https://doi.org/10.1002/2016JC011814}, +url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016JC011814}, +eprint = {https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2016JC011814}, +year = {2016} +} +@article{Seelig2021, +author = {Seelig, Torsten and Deneke, Hartwig and Quaas, Johannes and Tesche, Matthias}, +doi = {10.1029/2021JD035577}, +file = {:home/ahmed/Mendeley/Documents/2021/Seelig et al. - 2021 - Life Cycle of Shallow Marine Cumulus Clouds From Geostationary Satellite Observations(2).pdf:pdf}, +issn = {21698996}, +journal = {Journal of Geophysical Research: Atmospheres}, +keywords = {cloud physical properties,cloud tracking,geostationary observations,life cycle,shallow marine cumulus}, +mendeley-groups = {PhD/Tracking}, +number = {22}, +title = {{Life Cycle of Shallow Marine Cumulus Clouds From Geostationary Satellite Observations}}, +volume = {126}, +year = {2021} +} +@article{Willert1991, +author = {Willert, C. E. and Gharib, M.}, +doi = {10.1007/BF00190388}, +file = {:home/ahmed/Mendeley/Documents/1991/Willert, Gharib - 1991 - Digital particle image velocimetry.pdf:pdf}, +issn = {0723-4864}, +journal = {Experiments in Fluids}, +mendeley-groups = {PhD/Tracking}, +month = {jan}, +number = {4}, +pages = {181--193}, +title = {{Digital particle image velocimetry}}, +url = {http://link.springer.com/10.1007/BF00190388}, +volume = {10}, +year = {1991} +} + +@book{Goodfellow-et-al-2016, + title={Deep Learning}, + author={Ian Goodfellow and Yoshua Bengio and Aaron Courville}, + publisher={MIT Press}, + note={\url{http://www.deeplearningbook.org}}, + year={2016} +} + +@misc{Draelos_2019, +title={Convolution vs. cross-correlation}, +url={https://glassboxmedicine.com/2019/07/26/convolution-vs-cross-correlation/}, journal={Glass Box}, author={Draelos, Rachel}, +year={2019}, +month={Jul}} + +@book{Davis_2023, +place={New York}, +title={{Introduction to Environmental Data Science}}, +publisher={CRC Press}, +author={Davis, Jerry}, +doi={https://doi.org/10.1201/9781003317821}, +pages = 402, +year={2023}} + +@Manual{RManual, + title = {R: A Language and Environment for Statistical Computing}, + author = {{R Core Team}}, + organization = {R Foundation for Statistical Computing}, + address = {Vienna, Austria}, + year = {2023}, + url = {https://www.R-project.org/}, +} + +@book {Chambers2016, +author = { Chambers, John M. CRC Press Boca Raton, Fla }, +title = { Extending R }, +publisher = {CRC Press, Taylor & Francis Group}, +publisher = {}, +isbn = {9781498775717}, +keywords = { R Programm }, +year = {2016}, +abstract = {Literaturverzeichnis: Seite 355-357}, +booktitle = {The R series}, +booktitle = {A Chapman & Hall book}, +address = { Boca Raton, Fla. , }, +url = { http://slubdd.de/katalog?TN_libero_mab2 } +} + +@book {Wickham2023, +author = { Wickham, Hadley and Çetinkaya-Rundel, Mine and Grolemund, Garrett }, +title = { R for data science import, tidy, transform, visualize, and model data }, +edition = { 2nd edition } , +publisher = {O'Reilly}, +isbn = {9781492097402}, +keywords = { Data mining Computer programs , Information visualization Computer programs , R (Computer program language) , Big data , Databases , R Programm , Data Science }, +year = {2023}, +abstract = {Includes bibliographical references and index}, +address = { Beijing }, +url = { http://slubdd.de/katalog?TN_libero_mab2 } +} + +@book{sds, + author = {Edzer Pebesma and Roger Bivand}, + year = 2023, + title = {Spatial Data Science: With Applications in {R}}, + publisher = {Chapman and Hall/CRC}, + address = {Boca Raton}, + URL = {https://doi.org/10.1201/9780429459016}, + doi = {10.1201/9780429459016} +} + +@book{chirila2014introduction, + title={Introduction to Modern Fortran for the Earth System Sciences}, + author={Chirila, D.B. and Lohmann, G.}, + isbn={9783642370090}, + series={SpringerBriefs in earth systems sciences}, + url={https://books.google.de/books?id=bZCeBQAAQBAJ}, + year={2014}, + publisher={Springer Berlin Heidelberg} +} + +@article{Wilkinson2016, +author = {Wilkinson, Mark D and Dumontier, Michel and Aalbersberg, IJsbrand Jan and Appleton, Gabrielle and Axton, Myles and Baak, Arie and Blomberg, Niklas and Boiten, Jan-Willem and {da Silva Santos}, Luiz Bonino and Bourne, Philip E and Bouwman, Jildau and Brookes, Anthony J and Clark, Tim and Crosas, Merc{\`{e}} and Dillo, Ingrid and Dumon, Olivier and Edmunds, Scott and Evelo, Chris T and Finkers, Richard and Gonzalez-Beltran, Alejandra and Gray, Alasdair J.G. and Groth, Paul and Goble, Carole and Grethe, Jeffrey S and Heringa, Jaap and {'t Hoen}, Peter A.C and Hooft, Rob and Kuhn, Tobias and Kok, Ruben and Kok, Joost and Lusher, Scott J and Martone, Maryann E and Mons, Albert and Packer, Abel L and Persson, Bengt and Rocca-Serra, Philippe and Roos, Marco and van Schaik, Rene and Sansone, Susanna-Assunta and Schultes, Erik and Sengstag, Thierry and Slater, Ted and Strawn, George and Swertz, Morris A and Thompson, Mark and van der Lei, Johan and van Mulligen, Erik and Velterop, Jan and Waagmeester, Andra and Wittenburg, Peter and Wolstencroft, Katherine and Zhao, Jun and Mons, Barend}, +doi = {10.1038/sdata.2016.18}, +issn = {2052-4463}, +journal = {Scientific Data}, +month = {mar}, +number = {1}, +pages = {160018}, +title = {{The FAIR Guiding Principles for scientific data management and stewardship}}, +url = {https://doi.org/10.1038/sdata.2016.18 https://www.nature.com/articles/sdata201618}, +volume = {3}, +year = {2016} +} + +@online{Masuda_2020, +title={Modern Fortran Tutorial}, +url={https://masuday.github.io/fortran_tutorial/introduction.html}, +journal={Introduction}, +author={Masuda, Yutaka}, +year={2020}, +urldate = {04.03.2024}} + +@book{wickham2015advanced, + title={Advanced R}, + author={Wickham, H.}, + isbn={9781498759809}, + series={Chapman \& Hall/CRC The R Series}, + url={https://books.google.de/books?id=FfsYCwAAQBAJ}, + year={2015}, + publisher={CRC Press} +} \ No newline at end of file diff --git a/QuartoBook/scripts.qmd b/QuartoBook/scripts.qmd new file mode 100644 index 0000000..82bb4bf --- /dev/null +++ b/QuartoBook/scripts.qmd @@ -0,0 +1,524 @@ +--- +editor: + markdown: + wrap: 72 +--- + +# Accelerate R Scripts with Fortran {#chap-scripts} + +In this chapter, we will compare the efficiency of running three +computationally demanding examples in pure R script versus using another +R script with the core computations performed in Fortran. + +## 2D Coss-Correlation {#sec-2Dxcorr} + + + + + +In signal processing, cross-correlation is a measure of the similarity +between two signals as a function of the displacement of one relative to +the other [@Wang2019]. It can deliver information about the time lag +between the two signals. Cross-correlation is often applied in computer +vision for visual tracking. For example, it is used in template +matching, feature detection, and motion tracking. 2D cross-correlation +also plays an important role in convolutional networks and machine +learning. + + + +In atmospheric science, oceanography, climate, and water research, 2D +cross-correlation can be applied in various ways. For example, it can be +used to estimate ocean surface currents [@Warren2016], cloud tracking +using satellite imagery [@Seelig2021], and Particle Image Velocimetry +(PIV) in fluid dynamics applications [@Willert1991]. + +The 2D cross-correlation of an array $F$, with dimension $(M, N)$, and +array $G$, with dimension $(P, Q)$, can be given as the array $CC$ that +has a dimension of $(M+P-1, N+Q-1)$. + +$$ +CC_{(s,t)} = \sum_{m = 0}^{M-1}\sum_{n = 0}^{N-1} F_{(m,n)} G_{(m-s, n-t)} +$$ {#eq-2DCC} + +where $s$ varies between $-(P-1)$ and $(M-1)$ while $t$ varies between +$-(Q-1)$ and $(N-1)$. $F$ and $G \in R$. + +Now, let us define the `xcorr2D_r` function as shown in @lst-xcorrr. The +function receives two matrices or arrays `a` & `b` and return the full +cross-correlation plane `cc`. + +```{r} +#| lst-label: lst-xcorrr +#| lst-cap: Cross-correlation in R +xcorr2D_r0 <- function(a, b) { + # the full CC matrix + cc_row <- nrow(a) + nrow(b) - 1 + cc_col <- ncol(a) + ncol(b) - 1 + cc <- matrix(1:c(cc_row * cc_col), + byrow = FALSE, ncol = cc_col + ) + + # obtain possible shifts + min_row_shift <- -(nrow(b) - 1) + max_row_shift <- (nrow(a) - 1) + min_col_shift <- -(ncol(b) - 1) + max_col_shift <- (ncol(a) - 1) + + # Padded matrix + rows_padded <- abs(min_row_shift) + + nrow(a) + abs(max_row_shift) + cols_padded <- abs(min_col_shift) + + ncol(a) + abs(max_col_shift) + # a + padded_a <- matrix(0, + nrow = rows_padded, + ncol = cols_padded + ) + padded_a[ + (abs(min_row_shift) + 1):(abs(min_row_shift) + nrow(a)), + (abs(min_col_shift) + 1):(abs(min_col_shift) + ncol(a)) + ] <- a + + for (icol in 1:cc_col) { + for (irow in 1:cc_row) { + icc <- irow + ((icol - 1) * cc_row) + cols <- (icol):(icol + ncol(b) - 1) + rows <- (irow):(irow + nrow(b) - 1) + # b + padded_b <- array(0, + dim = c(rows_padded, cols_padded) + ) + padded_b[rows, cols] <- b + + cc[irow, icol] <- sum(padded_a * padded_b) + } + } + + return(cc) +} +``` + + + +Moving forward, we can define the `xcorr2d_f` subroutine in Fortran as +shown in @lst-xcorrf. Subroutines are generally the approach for integrating +Fortran in R. Function in Fortran return a single value with no option of +altering the input arguments, while subroutines have the +ability to perform complex tasks while altering input arguments. +This proofs to be helpful e.g., in solving equations system. + +Another imperative point is to define the dimension of the arrays when +passing them to Fortran (i.e. explicit-shape arrays). To illustrate, +`m, n, p, q, k, l` are the dimension of input arrays `a`and `b`, and +the out array `cc`. + +```{fortran, eval=FALSE} +#| lst-label: lst-xcorrf +#| lst-cap: Cross-correlation in Fortran +subroutine xcorr2d_f(m, n, p, q, k, l, a, b, cc) + implicit none + integer :: m, n, p, q, k, l + double precision, dimension(m, n) :: a + double precision, dimension(p, q) :: b + double precision, dimension(k, l) :: cc + ! dummy vars + integer :: min_row_shift, min_col_shift + integer :: max_row_shift, max_col_shift + integer :: rows_padded, cols_padded + integer :: icol, irow, icc, icol2, irow2 + real, allocatable, dimension(:, :) :: padded_a, padded_b + + ! obtain possible shfits + min_row_shift = -1*(p - 1) + max_row_shift = m - 1 + min_col_shift = -1*(q - 1) + max_col_shift = n - 1 + + ! Padded arrray + rows_padded = abs(min_row_shift) + m + abs(max_row_shift) + cols_padded = abs(min_col_shift) + n + abs(max_col_shift) + ! A + allocate (padded_a(rows_padded, cols_padded)) + padded_a = 0.0 + padded_a((abs(min_row_shift) + 1):(abs(min_row_shift) + m), & + (abs(min_col_shift) + 1):(abs(min_col_shift) + n)) = a + + ! B + allocate (padded_b(rows_padded, cols_padded)) + padded_b = 0.0 + do icol = 1, l + do irow = 1, k + icc = irow + ((icol - 1)*k) + icol2 = icol + q - 1 + irow2 = irow + p - 1 + padded_b(irow:irow2, icol:icol2) = b + cc(irow, icol) = sum(padded_a*padded_b) + padded_b = 0.0 + end do + end do +end subroutine xcorr2d_f +``` + +Since Fortran is a compiled language, we need to save the subroutine in +`xcorr2D.f90` file and compile it using: `R CMD SHLIB xcorr2D.f90` + +::: callout-note +Please use the terminal tab in Rstudio or open a new terminal using +`Alt+Shift+R` +::: + +As mentioned earlier, we need to pass the dimension of the arrays to +Fortran. Therefore, it would logical to write a wrapping function for +Fortran subroutine that provides other input arguments. + +In the wrapper function (@lst-xcorrfr), we initially require loading +the shared object, which is the compiled Fortran subroutine, as +`dyn.load("path/to/xcorr2D.so")`. Furthermore, it is important to prepare +other input variables for Fortran such as the dimensions of the input and +output arrays. Imperatively, data types should be approached carefully. Before +calling `.Fortran()`, all storage mode of the variables in R was +converted to the appropriate type using either `as.double()` or `as.integer()`. +If the wrong type is passed, it can result in a hard-to-catch error or +unexpected results^[Writing R Extensions, 5.2 Interface functions .C and .Fortran]. + + + +```{r} +#| lst-label: lst-xcorrfr +#| lst-cap: Cross-correlation wrapping function +xcorr2D_f0 <- function(a, b) { + # Please adjust the path to your setup. In my machine, + dyn.load("AquaFortR_Codes/xcorr2D.so") + + # the full CC matrix + cc_row <- nrow(a) + nrow(b) - 1 + cc_col <- ncol(a) + ncol(b) - 1 + cc <- matrix(1:c(cc_row * cc_col), byrow = FALSE, ncol = cc_col) + + result <- .Fortran("xcorr2d_f", + m = as.integer(dim(a)[1]), + n = as.integer(dim(a)[2]), + p = as.integer(dim(b)[1]), + q = as.integer(dim(b)[2]), + k = as.integer(cc_row), + l = as.integer(cc_row), + a = as.double(a), + b = as.double(b), + cc = as.double(cc) + )$cc + + return(result) +} +``` + +**Example** + +Now, we can use an example to compare the performance of the two functions. +In order to do so, `microbenchmark` package needs to be installed, and `ggplot2` +is required for plotting. + +The obtained benchmarking data allows (`mbm`) for a quantitative comparison of +the computational efficiency between the two methods. By printing "mbm" in the console +(`print(mbm)`) it is evident that Fortran outperforms the R implementation of 2D +cross-correlation by a factor of 10. The significance of leveraging Fortran becomes +evident in @fig-xcorr2D. + + +```{r} +#| label: fig-xcorr2D +#| fig-cap: "Performance comparison of 2D Cross-correlation in R and Fortran. Median is shown as red vertical line." +library(microbenchmark) +library(ggplot2) + +set.seed(72) +# Assume a +a <- structure(runif(64), dim = c(8L, 8L)) +# Assume b +b <- structure(runif(64), dim = c(8L, 8L)) +mbm <- microbenchmark( + xcorr2D_r = xcorr2D_r0(a, b), + xcorr2D_f = xcorr2D_f0(a, b) +) + +autoplot(mbm) + + stat_summary( + fun = "median", + geom = "crossbar", + width = 0.6, + colour = "red" + ) +``` + +## 2D Convolution + +Convolution and cross-correlation are both operations applied to images. +Cross-correlation means sliding a kernel (filter) across an image while +convolution means sliding a flipped kernel across an image. +[@Draelos_2019]. Therefore, the 2D convolution of an array $F$, with +dimension $(M, N)$, and array $G$, with dimension $(P, Q)$, can be given +as the array $CC$ that has a dimension of $(M+P-1, N+Q-1)$. $hv$ means +that $G$ is flipped. + +$$ +CC_{(s,t)} = \sum_{m = 0}^{M-1}\sum_{n = 0}^{N-1} F_{(m,n)} G_{(m-s, n-t)}^{hv} +$$ + +where $s$ varies between $-(P-1)$ and $(M-1)$ while $t$ varies between +$-(Q-1)$ and $(N-1)$. $F$ and $G \in R$. + +It is true that we can just flip the second array and use the functions from @sec-2Dxcorr. +However, we are interested in the complete workflow. + +@lst-convrr + +```{r} +#| lst-label: lst-convrr +#| lst-cap: Convolution in R +conv2D_r0 <- function(a, b) { + # the full convolution matrix + conv_row <- nrow(a) + nrow(b) - 1 + conv_col <- ncol(a) + ncol(b) - 1 + conv <- matrix(1:c(conv_row * conv_col), byrow = FALSE, ncol = conv_col) + + # obtain possible shifts + min_row_shift <- -(nrow(b) - 1) + max_row_shift <- (nrow(a) - 1) + min_col_shift <- -(ncol(b) - 1) + max_col_shift <- (ncol(a) - 1) + + # Padded matrix + rows_padded <- abs(min_row_shift) + nrow(a) + abs(max_row_shift) + cols_padded <- abs(min_col_shift) + ncol(a) + abs(max_col_shift) + # a + padded_a <- matrix(0, nrow = rows_padded, ncol = cols_padded) + padded_a[ + (abs(min_row_shift) + 1):(abs(min_row_shift) + nrow(a)), + (abs(min_col_shift) + 1):(abs(min_col_shift) + ncol(a)) + ] <- a + + for (icol in 1:conv_col) { + for (irow in 1:conv_row) { + iconv <- irow + ((icol - 1) * conv_row) + cols <- (icol):(icol + ncol(b) - 1) + rows <- (irow):(irow + nrow(b) - 1) + # b + padded_b <- array(0, dim = c(rows_padded, cols_padded)) + # flip the kernel i.e. b + padded_b[rows, cols] <- b[nrow(b):1, ncol(b):1] + + conv[irow, icol] <- sum(padded_a * padded_b) + } + } + + return(conv) +} +``` + +@lst-convrf + +```{fortran, eval=FALSE} +#| lst-label: lst-convrf +#| lst-cap: Convolution in Fortran +subroutine conv2d_f(m, n, p, q, k, l, a, b, conv) + implicit none + integer :: m, n, p, q, k, l, i, j + double precision, dimension(m, n) :: a + double precision, dimension(p, q) :: b, b_flipped + double precision, dimension(k, l) :: conv + ! dummy vars + integer :: min_row_shift, min_col_shift + integer :: max_row_shift, max_col_shift + integer :: rows_padded, cols_padded + integer :: icol, irow, iconv, icol2, irow2 + real, allocatable, dimension(:, :) :: padded_a, padded_b + + ! obtain possible shfits + min_row_shift = -1*(p - 1) + max_row_shift = m - 1 + min_col_shift = -1*(q - 1) + max_col_shift = n - 1 + + ! Flip the kernel i.e. B + b_flipped = 0.0 + do i = 1, p + do j = 1, q + b_flipped(p - i + 1, q - j + 1) = b(i, j) + end do + end do + + ! Padded arrray + rows_padded = abs(min_row_shift) + m + abs(max_row_shift) + cols_padded = abs(min_col_shift) + n + abs(max_col_shift) + ! A + allocate (padded_a(rows_padded, cols_padded)) + padded_a = 0.0 + padded_a((abs(min_row_shift) + 1):(abs(min_row_shift) + m), & + (abs(min_col_shift) + 1):(abs(min_col_shift) + n)) = a + + ! B + allocate (padded_b(rows_padded, cols_padded)) + padded_b = 0.0 + do icol = 1, l + do irow = 1, k + iconv = irow + ((icol - 1)*k) + icol2 = icol + q - 1 + irow2 = irow + p - 1 + padded_b(irow:irow2, icol:icol2) = b_flipped + conv(irow, icol) = sum(padded_a*padded_b) + padded_b = 0.0 + end do + end do +end subroutine conv2d_f +``` + +You might need some flags when compile Fortran subroutine, For example, to enable +the generation of run-time check `-fcheck=all`. The code below shows two options to +compile `conv2D.f90` either by R or gfortran compiler. + +```{bash} +#| eval: FALSE +# R +R CMD SHLIB conv2D.f90 + +# gfortran +gfortran -fpic -shared conv2D.f90 -o conv2D.so +``` + +@lst-convrfr + + +```{r} +#| lst-label: lst-convrfr +#| lst-cap: Convolution wrapping function +conv2D_f0 <- function(a, b) { + require(dotCall64) + dyn.load("AquaFortR_Codes/conv2D.so") + + # the full convolution matrix + conv_row <- nrow(a) + nrow(b) - 1 + conv_col <- ncol(a) + ncol(b) - 1 + conv <- matrix(0, byrow = FALSE, ncol = conv_col) + + result <- .C64("conv2d_f", + SIGNATURE = c(rep("integer",6), + rep("double",3)), + m = dim(a)[1], + n = dim(a)[2], + p = dim(b)[1], + q = dim(b)[2], + k = integer(conv_row), + l = integer(conv_col), + a = as.double(a), + b = as.double(b), + conv = as.double(conv))$conv + + return(result) +} +``` + +Similar to cross-correlation calculation, the Fortran implementation of convolution +outperforms the R one by + +```{r} +#| label: fig-conv2D +#| fig-cap: "Performance comparison of 2D Convolution in R and Fortran. Median is shown as red vertical line." +library(microbenchmark) +library(ggplot2) + + +set.seed(72) +# Assume a +a <- structure(runif(64), dim = c(8L, 8L)) +# Assume b +b <- structure(runif(64), dim = c(8L, 8L)) +mbm <- microbenchmark( + conv2D_r = conv2D_r0(a, b), + conv2D_f = conv2D_f0(a, b) +) + +autoplot(mbm) + + stat_summary( + fun = "mean", + geom = "crossbar", + width = 0.6, + colour = "red" + ) +``` + +::: {.callout-tip} +## Question + +After learning about `.Fortran()` and `.C64()`, you can use one of the two examples +above and compare the performance of the two interfaces using `microbenchmark()`. +Which function is faster? +::: + + +## Convective Available Potential Energy (CAPE) + + + + + + + +The maximum buoyancy of an undiluted air parcel, related to the +potential up-draft strength of thunderstorms. + +```{r} +#| warning: false +#| +if (!require(AquaFortR)) { + remotes::install_github("AHomoudi/AquaFortR", subdir = "RPackage") +} + +library(AquaFortR) +data("radiosonde") + +Temperature <- radiosonde$temp + 273.15 # K +Dewpoint <- radiosonde$dpt + 273.15 # K +Pressure <- radiosonde$pressure # hPa +# Mixing ratio +MixingRatio <- mixing_ratio_from_dewpoint(Dewpoint, Pressure) +t_parcel <- Temperature[1] +dwpt_parcel <- Dewpoint[1] +mr_parcel <- MixingRatio[1] + +source("AquaFortR_Codes/cape_r.R") +source("AquaFortR_Codes/cape_f.R") +``` + +```{r} +#| label: fig-cape +#| fig-cap: "Performance comparison of CAPE in R and Fortran" +library(microbenchmark) +library(ggplot2) + +mbm <- microbenchmark( + cape_r = cape_r0(t_parcel, dwpt_parcel, mr_parcel, + Pressure, Temperature, MixingRatio, + vtc = TRUE + ), + cape_f = cape_f0(t_parcel, dwpt_parcel, mr_parcel, + Pressure, Temperature, MixingRatio, + vtc = TRUE + ) +) + +autoplot(mbm) + + stat_summary( + fun = "mean", + geom = "crossbar", + width = 0.3, + colour = "green" + ) +``` + +https://www.f5wx.com/pages/CAPECalc.htm + +https://geo.libretexts.org/Bookshelves/Meteorology_and_Climate_Science/Practical_Meteorology\_(Stull)/14%3A_Thunderstorm_Fundamentals/14.03%3A_Section_4- + +https://journals.ametsoc.org/view/journals/clim/33/6/jcli-d-19-0461.1.xml + +https://www.weather.gov/source/zhu/ZHU_Training_Page/convective_parameters/Sounding_Stuff/MesoscaleParameters.html diff --git a/QuartoBook/summary.qmd b/QuartoBook/summary.qmd index 6767624..2c606b3 100644 --- a/QuartoBook/summary.qmd +++ b/QuartoBook/summary.qmd @@ -1,19 +1,39 @@ # Summary +- positive features of R +- positive features of Fortran + +- Why .Fortran() should not be used. dotCall64 is recommended for scripts. +- Mention the memory allocation issues related to .Fortran +- OpenMp implementation in Fortran, [Romp](https://github.com/wrathematics/Romp) +- Explain how beneficial R, Fortran with OpenMP is. For examples calculating cape for each vertical profile +in a high resolution (1-4km) climate simulation (years) output (nlat, lon, nlevels, ntime) can be accelerated +by passing the whole input arrays from R to Fortran, then loop over (nlat, lon, ntime) with OpenMP. + +- The cross-correlation can be speed up by using FFT. Also, in Fortran, the FFTW C library can be used. In summary, this book has no content whatsoever. - -```{r} -1 + 1 -``` - -Some issues on memory allocation https://masuday.github.io/fortran_tutorial/r.html - -OpenMP support https://cran.r-project.org/doc/manuals/R-exts.html#OpenMP-support - -using FFTW library - -The Need for Speed Part 1: Building an R Package with Fortran (or C) https://www.r-bloggers.com/2018/12/the-need-for-speed-part-1-building-an-r-package-with-fortran-or-c/ - -Advanced R by Hadley Wickham R’s C interface - -Conv and Xcorr shapes https://cran.r-project.org/web/packages/gsignal/gsignal.pdf +- Possibility of calling BLAS, LAPACK and LINPACK linear algebra functions Fortran using C (R_ext/BLAS.h, R_ext/Lapack.h and R_ext/Linpack.h) + +- Final words why packaging is better. + +## Further Reading + * [Fortran and R – Speed Things Up](https://www.r-bloggers.com/2014/04/fortran-and-r-speed-things-up/) + * [The Need for Speed Part 1: Building an R Package with Fortran (or C)](https://www.r-bloggers.com/2018/12/the-need-for-speed-part-1-building-an-r-package-with-fortran-or-c/) + * [The Need for Speed Part 2: C++ vs. Fortran vs. C](https://www.avrahamadler.com/2018/12/23/the-need-for-speed-part-2-c-vs-fortran-vs-c/) + * [Writing R Extensions](https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Writing-R-Extensions) + * [Extend R with Fortran](https://fortran-lang.discourse.group/t/extend-r-with-fortran/2386) + * [Advanced R by Hadley Wickham](http://adv-r.had.co.nz/) + * [Fortran 90 Tutorial](http://web.stanford.edu/class/me200c/tutorial_90/) + * [Modern Fortran Tutorial](https://masuday.github.io/fortran_tutorial/index.html) + * [Fortran Wiki - Libraries](https://fortranwiki.org/fortran/show/Libraries) + * [Fortran Best Practices](https://fortran-lang.org/en/learn/best_practices/#fortran-best-practices) + * [Hands-On Programming with R](https://rstudio-education.github.io/hopr/) + * [Spatial Data Science](https://rspatial.org/index.html) + * [R for Data Science (2e)](https://r4ds.hadley.nz/) + * [Introduction to Environmental Data Science](https://bookdown.org/igisc/EnvDataSci/) + + + + + + diff --git a/QuartoBook/title.tex b/QuartoBook/title.tex new file mode 100644 index 0000000..ae7b66c --- /dev/null +++ b/QuartoBook/title.tex @@ -0,0 +1,22 @@ +$if(title)$ +\title{$title$$if(thanks)$\thanks{$thanks$}$endif$} +$endif$ + +$if(subtitle)$ +\subtitle{$subtitle$} +$endif$ + +$for(by-author)$ + \author{$by-author.name.literal$} + $if(by-author.affiliations)$ + $for(by-author.affiliations)$ + \affil{% + $if(by-author.affiliations.name)$ + $by-author.affiliations.name$ + $endif$ + } + $endfor$ + $endif$ +$endfor$ + +\date{$date$}