############################Radiosounding########################################
# C.Burckhardt
# 2021-03
##############################################################################

####packages####
library(readr)
library(lubridate)
library(DescTools)
library(RadioSonde)

####working directory####
wd<- "Z:/burc/Radiosound"
setwd(wd)

####functions####
matlab2POS = function(x, timez = "UTC") {
  days = x - 719529 	# 719529 = days from 1-1-0000 to 1-1-1970
  secs = days * 86400 # 86400 seconds in a day
  # This next string of functions is a complete disaster, but it works.
  # It tries to outsmart R by converting the secs value to a POSIXct value
  # in the UTC time zone, then converts that to a time/date string that 
  # should lose the time zone, and then it performs a second as.POSIXct()
  # conversion on the time/date string to get a POSIXct value in the user's 
  # specified timezone. Time zones are a goddamned nightmare.
  return(as.POSIXct(strftime(as.POSIXct(secs, origin = '1970-1-1', 
                                        tz = 'UTC'), format = '%Y-%m-%d %H:%M', 
                             tz = 'UTC', usetz = FALSE), tz = timez))
}

#	Calculates boundary layer height using Richardson bulk approach. 
#	Translated from FLEXPART
#
#' @param height	height above ground [m]
#' @param ust	friction velocity [m/s]
#' @param ttlev air temperature [K]; 2D vector with dimensions height, time
#' @param qvlev water vapour mixing ratio [kg/kg]; 2D vector with dimensions height, time

#' @param uulev west-east wind component [m/s]; 2D vector with dimensions height, time

#' @param vvlev south-north wind component [m/s]; 2D vector with dimensions height, time

#' @param plev  air pressure [Pa]; 2D vector with dimensions height, time

#' @param tt2 	2 m air temperature [K] (currently not used); vector over time
#' @param hf 	surface sensible heat flux [W/m2]; vector over time 
richardson = function(height, ust, ttlev, qvlev, uulev, vvlev, plev, tt2, hf){
  
  #	constants
  R.air = 287.058
  c.p = 1004.6
  ga = 9.81
  const = R.air / ga
  ric = 0.25 
  b = 100.
  bs = 8.5
  itmax = 3
  
  nn.lev = dim(ttlev)[1]
  nn.tim = dim(ttlev)[2]
  
  excess = rep(0.,nn.tim)
  wst = rep(0.,nn.tim)
  hmix = rep(0.,nn.tim)
  
  for (it in 1:nn.tim){
    iter = 0 
    #repeat{
    tvold = ttlev[1,it] * (1. + 0.608 * qvlev[1,it])
    zref = 2.0
    
    thetaref = tvold * (100000./(plev[1,it]))^(R.air/c.p) + excess[it]
    thetaold = thetaref
    
    kk = 2
    repeat{
      tv = ttlev[kk,it] * (1. + 0.608 * qvlev[kk,it])
      z = height[kk]
      
      theta = tv * (100000./plev[kk,it])^(R.air/c.p)
      
      ri = ga / thetaref * (theta - thetaref) * (z - zref) /
        max(c((uulev[kk,it]-uulev[1,it])^2 + (vvlev[kk,it]-vvlev[1,it])^2 +
                b*ust[it]^2, 0.1))
      
      if (ri > ric && thetaold <= theta) break
      
      tvold = tv
      thetaold = theta
      if (kk >= nn.lev) break
      kk = kk + 1
    }
    
    
    #	linear interpolation between levels to obtain more detailed mixing height
    zz = seq(height[kk-1], height[kk], length.out=20)
    uul = seq(uulev[kk-1,it], uulev[kk,it], length.out=20)
    vvl = seq(vvlev[kk-1,it], vvlev[kk,it], length.out=20)
    thetal = seq(thetaold, theta, length.out=20)
    tmp = ((uul - uulev[1,it])^2 + (vvl - vvlev[1,it])^2 + b * ust[it]^2)
    tmp[tmp<0.1] = 0.1
    ril = ga / thetaref * (thetal - thetaref) * (zz-zref) / tmp
    hmix[it] = zz[which(ril>ric)[1]]
    
    #	calculate excess temperature
    #if (hf[it] < 0){
    #	wst[it] = (-hmix[it] * ga / thetaref * hf[it]/c.p)^(1./3.)
    #	excess[it] = -bs * hf[it] / c.p / wst[it]
    #	if (iter >= itmax) break
    #} else {
    #	wst[it] = 0.
    #	break
    #}
    
    iter = iter + 1
    #}	#	end of iterative loop for excess temperature
    
  }
  
  return(list(hmix=hmix, T.excess=excess, wst=wst))
}

####read data####

RS_pay_2020 <- read_csv("Z:/burc/Data/Radiosounding/sounding_pay_2020.csv")

#convert date format
RS_pay_2020$date <- matlab2POS(RS_pay_2020$t)

####constants####
R.air = 287.058 #J/(kg*K)
c.p = 1004.6    #J/(kg*K)
####extract needed period#####
Date1="2020-07-13"

####Routine####
if(TRUE==TRUE){
Date2=as.character(as.Date(Date1)+1)
data<-RS_pay_2020[as.Date(RS_pay_2020$date)>=Date1 & as.Date(RS_pay_2020$date)<=Date2 &
                              hour(RS_pay_2020$date)==0,]
#only up to z=1000m
data <- data[data$z<=2000,]
#remove NaN to NA
data$T[data$T=="NaN"]<-NA
data$rh[data$rh=="NaN"]<-NA
data$u[data$u=="NaN"]<-NA
data$v[data$v=="NaN"]<-NA

####define data####
#Radiosounding:
data<-na.omit(data)

#meteorology:
dates <- unique(format(as.POSIXct(data$date),"%Y-%m-%d %H:%M:%S"))
data_met<-PAY_MET[as.character(PAY_MET$date) %in% dates,]

#soundingstimes
RS_no <- unique(data$date)

####prepare specific humidity####
#from Wallace Hobbs p.81:

#mixing ratio w [g kg^-1]
#specific humidity q [kg kg^-1]
#factor Mw/Md eps
eps=0.622

# saturation water vapor pressure in moist air [Pa] (Wang 2013, p.104):
data$esat_w <- exp(54.842763-6763.22/data$T-4.210*log(data$T)+0.000367*data$T+tanh(0.0415*(data$T-218.8))*
             (53.878-1331.22/data$T-9.44523*log(data$T)+0.014025*data$T))
#saturation mixing ratio:
data$ws <- eps*(data$esat_w/(data$p*100))
#mixing ratio:
data$w <- data$rh/100*data$ws
#specific humidity
data$q <- data$w/(1+data$w)

####approximate u*####
# roughnesslength z_0 from http://www.verwaltungsvorschriften-im-internet.de/pdf/BMU-IG-20020724-KF-A003.pdf 
# Payerne between grassland and city
z_0 = 0.5
z_10=10
kappa = 0.4
u_10m <- data_met$Windspeed
u_st <- (u_10m*kappa)/log(z_10/z_0)

####prepare temperature, u, v, pressure vectors####
#create new matrix
j=1
ttlev<-list()
plev<-list()
uulev<-list()
vvlev<-list()
qvlev<-list()
Ri_b<-data.frame(date=POSIXct(),PBL=numeric())
data$z<-data$z-data$z[1]

while(j <= length(RS_no)) {
  RS_Matrix <- matrix(data=NA, nrow=length(unique(data$z[data$date==RS_no[j]])),ncol = 1)
  rownames(RS_Matrix)<-sort(unique(data$z[data$date==RS_no[j]]))
  colnames(RS_Matrix)<-as.character(RS_no)[j]
  
  ##fill matrix with variables
  k=1
  while (k>0) {
    if(k==1){
      #write temperature 2D vector
      for (i in 1:dim(data)[1]) {
        RS_Matrix[rownames(RS_Matrix) == data$z[i],
                  colnames(RS_Matrix) == as.character(RS_no)[j]] <- data$T[i]
      }
      ttlev[[j]]<-RS_Matrix
      RS_Matrix[!is.na(RS_Matrix)]<-NA
      k=k+1
    } else if (k==2) {
      #write pressure 2D vector
      for (i in 1:dim(data)[1]) {
        RS_Matrix[rownames(RS_Matrix) == data$z[i],
                  colnames(RS_Matrix) == as.character(RS_no)[j]] <- data$p[i]
      }
      plev[[j]]<-RS_Matrix
      RS_Matrix[!is.na(RS_Matrix)]<-NA
      k=k+1
    } else if (k==3) {
      #write u 2D vector
      for (i in 1:dim(data)[1]) {
        RS_Matrix[rownames(RS_Matrix) == data$z[i],
                  colnames(RS_Matrix) == as.character(RS_no)[j]] <- data$u[i]
      }
      uulev[[j]]<-RS_Matrix
      RS_Matrix[!is.na(RS_Matrix)]<-NA
      k=k+1
    } else if (k==4) {
      #write v 2D vector
      for (i in 1:dim(data)[1]) {
        RS_Matrix[rownames(RS_Matrix) == data$z[i],
                  colnames(RS_Matrix) == as.character(RS_no)[j]] <- data$v[i]
      }
      vvlev[[j]]<-RS_Matrix
      RS_Matrix[!is.na(RS_Matrix)]<-NA
      k=k+1
    } else if (k==5){
      #write q 2D vector
      for (i in 1:dim(data)[1]) {
        RS_Matrix[rownames(RS_Matrix) == data$z[i],
                  colnames(RS_Matrix) == as.character(RS_no)[j]] <- data$q[i]
      }
      qvlev[[j]]<-RS_Matrix
      RS_Matrix[!is.na(RS_Matrix)]<-NA
      k=0
    }
  }
  ####calculate bulk Richardson####

  Ri_b[j,2]<- richardson(height = data$z[data$date==RS_no[j]], ust = u_st[j], ttlev = ttlev[[j]], qvlev = qvlev[[j]], uulev = uulev[[j]],
                    vvlev = vvlev[[j]], plev = plev[[j]])[1]
  Ri_b[j,1]<- RS_no[j]
  j=j+1
}

#write Richardson bulk
write.table(Ri_b,paste0("Z:/burc/Radiosound/Ri_b",substr(Date1,1,4),substr(Date1,6,7),substr(Date1,9,10),"_",substr(Date2,9,10),".csv"),sep=";",row.names = F)

####SBI and SBL pT####
RS_SBL=data.frame(date=as.POSIXct(unique(data$date)),SBI=NA,SBLpT=NA)
data$theta<-data$T*(1000/data$p)^(R.air/c.p)
#Date1
RS_del_T_date1<-diff(data$T[as.Date(data$date)==Date1])/diff(data$z[as.Date(data$date)==Date1])
RS_min<-which(RS_del_T_date1<0)
for (i in 3:length(RS_min)) {
  if(RS_min[i]==RS_min[i-2]+2){
    n=RS_min[i-2]
    break
  }
}
RS_SBL[1,2]<- data$z[n]
RS_del_theta_date1<-diff(data$theta[as.Date(data$date)==Date1])/diff(data$z[as.Date(data$date)==Date1])
RS_min<- which(!is.nan(RS_del_theta_date1) &abs(RS_del_theta_date1)<=0.002)
for (i in 6:length(RS_min)) {
  if(RS_min[i]==RS_min[i-5]+5){
    n=RS_min[i-5]
    break
  }
}
RS_SBL[1,3]<- data$z[n]
#Date2
RS_del_T_date2<-diff(data$T[as.Date(data$date)==Date2])/diff(data$z[as.Date(data$date)==Date2])
RS_min<-which(RS_del_T_date2<0)
for (i in 3:length(RS_min)) {
  if(RS_min[i]==RS_min[i-2]+2){
    n=RS_min[i-2]
    break
  }
}
RS_SBL[2,2]<- data$z[n]
RS_del_theta_date2<-diff(data$theta[as.Date(data$date)==Date2])/diff(data$z[as.Date(data$date)==Date2])
RS_min<- which(!is.nan(RS_del_theta_date2) &abs(RS_del_theta_date2)<=0.005)
for (i in 6:length(RS_min)) {
  if(RS_min[i]==RS_min[i-5]+5){
    n=RS_min[i-5]
    break
  }
}
RS_SBL[2,3]<- data$z[n+length(data$z[as.Date(data$date)==Date1])]
remove(RS_min)
write.table(RS_SBL,paste0("Z:/burc/Radiosound/RS_SBL",substr(Date1,1,4),substr(Date1,6,7),substr(Date1,9,10),"_",substr(Date2,9,10),".csv"),sep=";",row.names = F)
}

####plot sounding####

RS_plot<-data[as.Date(data$date)==Date2,]
RS_plot$dewpt<- RS_plot$T-((100-RS_plot$rh)/5)
RS_plot$dewpt<-RS_plot$dewpt*9/5-459.67
#RS_plot$dewpt <- RS_plot$dewpt-273.15
#RS_plot$temp<-RS_plot$T*9/5-459.67
#RS_plot$temp <- RS_plot$T-273.15
RS_plot<-RS_plot[,2:10]
colnames(RS_plot)<-c("height","T","rh","press","uwind","vwind","date","dewpt","temp")

plotsonde(RS_plot)

####SBI and SBL pT####
RS_SBL=data.frame(date=as.POSIXct(unique(data$date)),SBI=NA,SBLpT=NA)
data$theta<-data$T*(1000/data$p)^(R.air/c.p)
#Date1
RS_del_T_date1<-diff(data$T[as.Date(data$date)==Date1])/diff(data$z[as.Date(data$date)==Date1])
RS_SBL[1,2]<- data$z[min(which(RS_del_T_date1<0))]
RS_del_theta_date1<-diff(data$theta[as.Date(data$date)==Date1])/diff(data$z[as.Date(data$date)==Date1])
RS_SBL[1,3]<- data$z[min(which(!is.nan(RS_del_theta_date1) &RS_del_theta_date1<=0.005))]
#Date2
RS_del_T_date2<-diff(data$T[as.Date(data$date)==Date2])/diff(data$z[as.Date(data$date)==Date2])
RS_SBL[2,2]<- data$z[min(which(RS_del_T_date2<0))]
RS_del_theta_date2<-diff(data$theta[as.Date(data$date)==Date2])/diff(data$z[as.Date(data$date)==Date2])
RS_SBL[2,3]<- data$z[min(which(!is.nan(RS_del_theta_date2) &RS_del_theta_date2<=0.005))+length(data$z[as.Date(data$date)==Date1])]

write.table(RS_SBL,paste0("Z:/burc/Radiosound/RS_SBL",substr(Date1,1,4),substr(Date1,6,7),substr(Date1,9,10),"_",substr(Date2,9,10),".csv"),sep=";",row.names = F)