R/W_flux.R

In santiagoh719/ClimFunctions: Functions For Climate Data

Takaya_Nakamura_Flux <- function(GPH, U, V, lat, lon, lev,
                                 lat_lon = TRUE, spherical_x = FALSE,
                                 a = 6371000, rad = FALSE){
  # GPH anomaly
  # U and V climatology
  # level in which is working (in hPa)

  checks <- makeAssertCollection()
  assert_array(GPH,mode = 'numeric',d = 2, add = checks)
  assert_array(U,mode = 'numeric',d = 2, add = checks)
  assert_array(V,mode = 'numeric',d = 2, add = checks)

  assert_numeric(lat, lower = -90, upper = 90,add = checks)
  assert_numeric(lon, lower = -360, upper = 360,add = checks)
  assert_numeric(lev, lower = 10, upper = 1000,add = checks)
  assert_numeric(a, len = 1, finite = TRUE,add = checks)

  assert_logical(spherical_x, add = checks)

  assert_logical(lat_lon, add = checks)
  assert_logical(rad, add = checks)
  reportAssertions(checks)

  if( ! identical(dim(GPH),dim(U)) | ! identical(dim(V),dim(U)) ){
    stop('the dims of GPH, U and V must be equal')
  }

  if(lat_lon){
    if(dim(GPH)[1] != length(lat) | dim(GPH)[2] != length(lon)){
      stop('The dimetions of GPH do not match with the length of lat and lon, check if the data is order lat_lon or not')
    }
  }else{
    if(dim(GPH)[1] != length(lon) | dim(GPH)[2] != length(lat)){
      stop('The dimetions of GPH do not match with the length of lat and lon, check if the data is order lat_lon or not')
    }
    GPH <- aperm(a = GPH,perm = c(2,1))
    U <- aperm(a = U,perm = c(2,1))
    V <- aperm(a = V,perm = c(2,1))

  }

  if(!rad){
    lat <- lat / 180 * pi
    lon <- lon / 180 * pi
  }

  GPH <- GPH * 9.8 / array(data = coriolis_f(lat = lat, rad = TRUE),dim=dim(GPH))

  GPH_dx_dy <-   derivate(data_array = GPH, dx = TRUE, dy = TRUE,
                          lat = lat, lon = lon, lat_lon = TRUE,
                          spherical_x = spherical_x,
                          a = a, rad = TRUE)
  GPH_dxdx_dxdy <-   derivate(data_array = GPH_dx_dy$dx, dx = TRUE, dy = TRUE,
                              lat = lat, lon = lon, lat_lon = TRUE,
                              spherical_x = spherical_x,
                              a = a, rad = TRUE)
  GPH_dydy <-   derivate(data_array = GPH_dx_dy$dy, dx = FALSE, dy = TRUE,
                         lat = lat, lon = lon, lat_lon = TRUE,
                         spherical_x = spherical_x,
                         a = a, rad = TRUE)

  termxv <- GPH_dx_dy$dx * GPH_dx_dy$dy - GPH * GPH_dxdx_dxdy$dy
  termxu <- GPH_dx_dy$dx**2 - GPH * GPH_dxdx_dxdy$dx
  termyv <- GPH_dx_dy$dy**2 - GPH * GPH_dydy$dy

  coef <-((lev/1000)*cos(array(data = lat,dim=dim(GPH))))/(2*sqrt(U**2 + V**2))

  tx <- coef*( U*termxu + V*termxv)
  ty <- coef*( U*termxv + V*termyv)

  if(!lat_lon){
    tx <- aperm(a = tx, perm = 2:1)
    ty <- aperm(a = ty, perm = 2:1)
  }
  return(list(Flux_x = tx,
              Flux_y = ty))
}