diff --git a/DESCRIPTION b/DESCRIPTION index b23e150d6eab012afa734c3f3fd5d095706010a3..2f5c132250a2ae2d4498ef534675b84b1da24c27 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -79,4 +79,4 @@ VignetteBuilder: knitr License: Apache License 2.0 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.0.2 +RoxygenNote: 7.1.1 diff --git a/NAMESPACE b/NAMESPACE index 6c95ac0826159823d987706fb8d064f9d53bda39..1153623fd38fe20490d13a0858e13240f1643154 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ export(Analogs) export(BEI_PDFBest) export(BEI_Weights) export(CST_Analogs) +export(CST_AnalogsPredictors) export(CST_Anomaly) export(CST_BEI_Weighting) export(CST_BiasCorrection) @@ -46,6 +47,7 @@ export(SplitDim) export(WeatherRegime) export(as.s2dv_cube) export(s2dv_cube) +export(training_analogs) import(abind) import(ggplot2) import(multiApply) @@ -108,5 +110,7 @@ importFrom(s2dverification,Mean1Dim) importFrom(s2dverification,Subset) importFrom(utils,glob2rx) importFrom(utils,head) +importFrom(utils,read.table) importFrom(utils,tail) importFrom(verification,verify) +useDynLib(CSTools, .registration=TRUE) diff --git a/R/AnalogsPred_train.R b/R/AnalogsPred_train.R new file mode 100644 index 0000000000000000000000000000000000000000..c68c48b05cd5991de93ff5f1282fa761148693a1 --- /dev/null +++ b/R/AnalogsPred_train.R @@ -0,0 +1,534 @@ +#' AEMET Training +#' Training method (pre-downscaling) based on analogs: +#' synoptic situations and significant predictors. +#' +#'@author Marta Dominguez Alonso - AEMET, \email{mdomingueza@aemet.es} +#'@author Nuria Perez-Zanon - BSC, \email{nuria.perez@bsc.es} +#' +#'@description This function caracterizes the synoptic situations in a past period based on +#' low resolution reanalysis data (e.g, ERAInterim 1.5º x 1.5º) and an observational high +#' resolution (HR) dataset (AEMET 5 km gridded daily precipitation and maximum and +#' minimum temperature) (Peral et al., 2017)). +#' The method uses three domains: +#' - peninsular Spain and Balearic Islands domain (5 km resolution): HR domain +#' - synoptic domain (low resolution): it should be centered over Iberian Peninsula and +#' cover enough extension to detect as much synoptic situations as possible. +#' - extended domain (low resolution): it is an extension of the synoptic +#' domain. It is used for 'slp_ext' parameter (see 'slp_lon' and 'slp_lat' below). +#'@param pred List of matrix reanalysis data in a synoptic domain. The list +#' has to contain reanalysis atmospheric variables (instantaneous 12h data) +#' that must be indentify by parenthesis name. +#' For precipitation: +#' - u component of wind at 500 hPa (u500) in m/s +#' - v component of wind at 500 hPa (v500) in m/s +#' - temperature at 500 hPa (t500) in K +#' - temperature at 850 hPa (t850) in K +#' - temperature at 1000 hPa (t1000) in K +#' - geopotential height at 500 hPa (z500) in m +#' - geopotential height at 1000 hPa (z1000) in m +#' - sea level pressure (slp) in hPa +#' - specific humidity at 700 hPa (q700) in g/kg +#' For maximum and minimum temperature: +#' - temperature at 1000 hPa (t1000) in K +#' - sea level pressure (slp) in hPa +#' All matrix must have [time,gridpoint] dimensions. +#' (time = number of training days, gridpoint = number of synoptic gridpoints). +#'@param slp_ext Matrix with atmospheric reanalysis sea level pressure +#' (instantaneous 12h data)(hPa). It has the same resolution as 'pred' parameter +#' but with an extended domain. This domain contains extra degrees (most in the +#' north and west part) compare to synoptic domain. The matrix must have +#' [time,gridpoint] dimensions. +#' (time = number of training days, gridpoint = number of extended gridpoints). +#'@param lon Vector of the synoptic longitude (from (-180º) to 180º), +#' The vector must go from west to east. +#'@param lat Vector of the synoptic latitude. The vector must go from north to south. +#'@param slp_lon Vector of the extended longitude (from (-180º) to 180º) +#' The vector must go from west to east. +#'@param slp_lat Vector of the extended latitude. The vector must go from north to south. +#'@param var Variable name to downscale. There are two options: 'prec' for +#' precipitation and 'temp' for maximum and minimum temperature. +#'@param HR_path Local path of HR observational files (maestro and pcp/tmx-tmn). +#' For precipitation can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/v2/Serie_AEMET_v2_pcp_1951a202006_txt.tar.gz +#' For maximum and minimum temperature can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmax/Serie_AEMET_v1_tmax_1951a202006_txt.tar.gz and http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmin/Serie_AEMET_v1_tmin_1951a202006_txt.tar.gz respetively. +#' Maestro file (maestro_red_hr_SPAIN.txt) has gridpoint (nptos), longitude (lon), latitude (lat) and +#' altitude (alt) in columns (vector structure). +#' Data file (pcp/tmx/tmn_red_SPAIN_1951-201903.txt) includes 5km resolution spanish daily data +#' (precipitation or maximum and minimum temperature from january 1951 to june 2020. See README +#' file for more information. +#' IMPORTANT!: HR observational period must be the same as for reanalysis variables. +#' It is assumed that the training period is smaller than the HR original one (1951-2020), so it is +#' needed to make a new ascii file with the new period and the same structure as original, +#' specifying the training dates ('tdates' parameter) in the name +#' (e.g. 'pcp_red_SPAIN_19810101-19961231.txt' for '19810101-19961231' period). +#'@param tdates Training period dates in format YYYYMMDD(start)-YYYYMMDD(end) (e.g. 19810101-19961231). +#'@return matrix list (e.g. restrain) as a result of characterize the past synoptic +#' situations and the significant predictors needed to downscale seasonal forecast variables. +#' For precipitation the output includes: +#' um: u component of geostrophic wind in all period (numeric matrix with [time,gridpoint] dimensions) +#' vm: v component of geostrophic wind in all period (numeric matrix with [time,gridpoint] dimensions) +#' nger: number of synoptic situations (integer) +#' gu92: u component of geostrophic wind for each synoptic situation (numeric matrix with +#' [nger,gridpoint] dimensions) +#' gv92: v component of geostrophic wind for each synoptic situation (numeric matrix with +#' [nger,gridpoint] dimensions) +#' gu52: u component of wind at 500 hPa for each synotic situation (numeric matrix with +#' [nger,gridpoint] dimensions) +#' gv52: v component of wind at 500 hPa for each synotic situation (numeric matrix with +#' [nger,gridpoint] dimensions) +#' neni: number of reference centers where predictors are calculated (integer) +#' vdmin: minimum distances between each HR gridpoint and the four nearest synoptic +#' gridpoints (numeric matrix with [nptos,4] dimensions) (nptos = number of HR gridpoints) +#' vref: four nearest synoptic gridpoints to each HR gridpoint (integer matrix with +#' [nptos,4] dimensions) +#' ccm: multiple correlation coeficients (numeric matrix with [nger,nptos] dimensions) +#' indices: +#' - lab_pred: numeric labels of selected predictors (integer matrix +#' with [nger,nptos,11,1] dimensions) +#' - cor_pred: partial correlation of selected predictors (numeric matrix with +#' [nger,nptos,11,2] dimensions) +#' For maximum and minimum temperature the output includes: +#' um: u component of geostrophic wind in all training period (numeric matrix with [time,gridpoint] dimensions) +#' vm: v component of geostrophic wind in all training period (numeric matrix with [time,gridpoint] dimensions) +#' insol: insolation in all training period (numeric vector with [time] dimension) +#' neni: number of reference centers where predictors are calculated (integer) +#' vdmin: minimum distances between each HR gridpoint and the four nearest synoptic +#' gridpoints (numeric matrix with [nptos,4] dimensions) (nptos = number of HR gridpoints) +#' vref: four nearest synoptic gridpoints to each HR gridpoint (integer matrix with +#' [nptos,4] dimensions) +#' +#' The output can directly use as argument to 'CST_AnalogsPredictors' function +#' (e.g. resdowns <- CST_AnalogsPredictors(...,restrain)) +#' +#'@importFrom utils read.table +#' +#'@useDynLib CSTools +#' +#'@export + +training_analogs <- function(pred, + slp_ext, + lon, + lat, + slp_lon, + slp_lat, + var, + HR_path, + tdates) { + +if (!is.list(pred)) { + stop("Parameter 'pred' must be a list of 'matrix' objects") + } + +if (!(all(sapply(pred, inherits, 'matrix')))) { + stop("Elements of the list in parameter 'pred' must be of the class ", + "'matrix'.") + } + +if (var == "prec") { + if (length(pred) != 9) { + stop("Parameter 'pred' must be a length of 9.") + } else { + if (is.null(names(dim(pred[[1]]))) || + is.null(names(dim(pred[[2]]))) || + is.null(names(dim(pred[[3]]))) || + is.null(names(dim(pred[[4]]))) || + is.null(names(dim(pred[[5]]))) || + is.null(names(dim(pred[[6]]))) || + is.null(names(dim(pred[[7]]))) || + is.null(names(dim(pred[[8]]))) || + is.null(names(dim(pred[[9]])))) { + stop("Parameter 'pred' should have dimmension names.") + } + if (!(any(names(pred) %in% "u500"))) { + stop("Variable 'u500' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "v500"))) { + stop("Variable 'v500' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "t500"))) { + stop("Variable 't500' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "t850"))) { + stop("Variable 't850' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "t1000"))) { + stop("Variable 't1000' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "z500"))) { + stop("Variable 'z500' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "z1000"))) { + stop("Variable 'z1000' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "slp"))) { + stop("Variable 'slp' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "q700"))) { + stop("Variable 'q700' in pred parameter is missed.") + } + } +} else { + if (length(pred) != 2) { + stop("Parameter 'pred' must be a length of 2.") + } else { + if (is.null(names(dim(pred[[1]]))) || + is.null(names(dim(pred[[2]])))) { + stop("Parameter 'pred' should have dimmension names.") + } + if (!(any(names(pred) %in% "t1000"))) { + stop("Variable 't1000' in pred parameter is missed.") + } else if (!(any(names(pred) %in% "slp"))) { + stop("Variable 'slp' in pred parameter is missed.") + } + } +} + + if (all((sapply(pred,dim))==dim(pred[[1]])) & + all((sapply(pred,function(pred){names(dim(pred))}))==names(dim(pred[[1]])))) { + dim_pred <- dim(pred[[1]]) + if (!(any(names(dim_pred) %in% "time"))) { + stop("Dimension 'time' in pred parameter is missed.") + } + if (!(any(names(dim_pred) %in% "gridpoint"))) { + stop("Dimension 'gridpoint' in pred parameter is missed.") + } + if (names(dim_pred)[1] == "gridpoint") { + pred <- lapply(pred,aperm) + } else { + pred <- pred + } + } else { + stop("All 'pred' variables must have the same dimensions and name dimensions.") + } + +if (!is.vector(lon) || !is.numeric(lon)) { + stop("Parameter 'lon' must be a numeric vector") +} else { + if (is.unsorted(lon)) { + lon <- sort(lon) + warning("'lon' vector has been sorted in increasing order") + } +} + +if (!is.vector(lat) || !is.numeric(lat)) { + stop("Parameter 'lat' must be a numeric vector") +} else { + if (!is.unsorted(lat)) { + lat <- sort(lat, decreasing = TRUE) + warning("'lat' vector has been sorted in decreasing order") + } +} + +if (!is.character(HR_path)) { + stop("Parameter 'HR_path' must be a character.") +} else { + if (!dir.exists(HR_path)) { + stop("'HR_path' directory does not exist") + } +} + +if (!is.character(tdates)) { + stop("Parameter 'tdates' must be a character.") +} else { + if (nchar(tdates) != "17") { + stop("Parameter 'tdates' must be a string with 17 charecters.") + } else { + dateini <- as.Date(substr(tdates,start=1,stop=8),format="%Y%m%d") + dateend <- as.Date(substr(tdates,start=10,stop=18),format="%Y%m%d") + if (dateend <= dateini) { + stop("Parameter 'tdates' must be at least of one day") + } + } +} + +#! REANALYSIS GRID PARAMETERS + + rlon <- c(lon, NA) - c(NA, lon) + rlon <- rlon[!is.na(rlon)] + if (!all(rlon == rlon[1])) { + stop("Parameter 'lon' must be in regular grid.") + } else { + rlon <- rlon[1] + } + + rlat <- c(lat, NA) - c(NA, lat) + rlat <- rlat[!is.na(rlat)] + if (!all(rlat == rlat[1])) { + stop("Parameter 'lat' must be in regular grid.") + } else { + rlat <- rlat[1] + } + + if (rlon != (-rlat)) { + stop("Parameters 'lon' and 'lat' must have the same resolution.") + } else { + res <- rlon + } + + nlat <- ((lat[length(lat)] - lat[1]) / rlat) + 1 + nlon <- ((lon[length(lon)] - lon[1]) / rlon) + 1 + + ic <- nlat * nlon +# + slp_rlon <- c(slp_lon, NA) - c(NA, slp_lon) + slp_rlon <- slp_rlon[!is.na(slp_rlon)] + if (!all(slp_rlon == slp_rlon[1])) { + stop("Parameter 'slp_lon' must be in regular grid.") + } else { + slp_rlon <- slp_rlon[1] + } + + slp_rlat <- c(slp_lat, NA) - c(NA, slp_lat) + slp_rlat <- slp_rlat[!is.na(slp_rlat)] + if (!all(slp_rlat == slp_rlat[1])) { + stop("Parameter 'slp_lat' must be in regular grid.") + } else { + slp_rlat <- slp_rlat[1] + } + + if (slp_rlon != (-slp_rlat)) { + stop("Parameters 'slp_lon' and 'slp_lat' must have the same resolution.") + } else { + slp_res <- slp_rlon + } + + nlatt <- ((slp_lat[length(slp_lat)] - slp_lat[1]) / slp_rlat) + 1 + nlont <- ((slp_lon[length(slp_lon)] - slp_lon[1]) / slp_rlon) + 1 + + id <- nlatt * nlont + + slat <- max(lat) + slon <- min(c(lon[which(lon > 180)] - 360, + lon[which(lon <= 180)])) + + slatt <- max(slp_lat) + slont <- min(c(slp_lon[which(slp_lon > 180)] - 360, + slp_lon[which(slp_lon <= 180)])) + + ngridd <- ((2*nlatt)-1)*((2*nlont)-1) + + if (all((sapply(pred,nrow))==nrow(pred[[1]]))){ + nd <- nrow(pred[[1]]) + } else { + stop("All 'pred' variables must be in the same period.") + } + +#!!!!! COMPROBAR QUE SLP TAMBIEN TIENE EL MISMO NROW + + seqdates <- seq(as.Date(substr(tdates,start=1,stop=8),format="%Y%m%d"), + as.Date(substr(tdates,start=10,stop=18),format="%Y%m%d"),by="days") + month <- format(seqdates,format="%m") + day <- format(seqdates,format="%d") + +#! TRAINING REANALYSIS VARIABLES +t1000 <- pred[['t1000']] +msl_si <- pred[['slp']] +msl_lr <- slp_ext + +if (var == "prec") { +u500 <- pred[['u500']] +v500 <- pred[['v500']] +t500 <- pred[['t500']] +t850 <- pred[['t850']] +z500 <- pred[['z500']] +z1000 <- pred[['z1000']] +q700 <- pred[['q700']] +} + +#! HIGH-RESOLUTION (HR) OBSERVATIONAL DATASET +maestro_hr_file <- paste(HR_path, "maestro_red_hr_SPAIN.txt",sep="") +if (!file.exists(maestro_hr_file)) { + stop("'maestro_red_hr_SPAIN.txt' does not exist.") +} else { + maestro <- read.table(maestro_hr_file) + lon_hr <- unlist(maestro[2]) + lat_hr <- unlist(maestro[3]) + nptos <- length(readLines(maestro_hr_file)) +} + +if (var == "prec") { + prec_hr_file <- paste(HR_path, "pcp_red_SPAIN_",tdates,".txt",sep="") + if (!file.exists(prec_hr_file)) { + stop(sprintf("precipitation HR file for %s does not exist.",tdates)) + } else { + nd_hr <- length(readLines(prec_hr_file)) + preprec_hr <- matrix(scan(prec_hr_file), nrow=nd_hr ,ncol= nptos+1, byrow=TRUE) + prec_hr <- preprec_hr[1:nd_hr,-c(1)] + } +} else { + tmx_hr_file <- paste(HR_path, "tmx_red_SPAIN_",tdates,".txt",sep="") + tmn_hr_file <- paste(HR_path, "tmn_red_SPAIN_",tdates,".txt",sep="") + if (!file.exists(tmx_hr_file)) { + stop(sprintf("maximum temperature HR file for %s does not exist.",tdates)) + } else if (!file.exists(tmn_hr_file)) { + stop(sprintf("minimum temperature HR file for %s does not exist.",tdates)) + } else if (length(readLines(tmx_hr_file)) != length(readLines(tmn_hr_file))) { + stop("maximum and minimum temperature HR observation files must have the same period.") + } else { + nd_hr <- length(readLines(tmx_hr_file)) + pretmx_hr <- matrix(scan(tmx_hr_file), nrow=nd_hr ,ncol= nptos+1, byrow=TRUE) + tmx_hr <- pretmx_hr[1:nd_hr,-c(1)] + pretmn_hr <- matrix(scan(tmn_hr_file), nrow=nd_hr ,ncol= nptos+1, byrow=TRUE) + tmn_hr <- pretmn_hr[1:nd_hr,-c(1)] + } +} + if (nd_hr != nd) { + stop("Reanalysis variables and HR observations must have the same period.") + } + +#! OTHER PARAMETERS that should not be changed +#! Number of analog situations to consider +nanx <- 155 +#! Number of predictors +npx <- 11 + +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +if (var == "prec") { + + prePro <- .Fortran("training_part1_prec", + u500 = as.numeric(u500), + v500 = as.numeric(v500), + t1000 = as.numeric(t1000), + z500 = as.numeric(z500), + z1000 = as.numeric(z1000), + msl_si = as.numeric(msl_si), + msl_lr = as.numeric(msl_lr), + ngridd = as.integer(ngridd), + nlat = as.integer(nlat), + nlon = as.integer(nlon), + ic = as.integer(ic), + nlatt = as.integer(nlatt), + nlont = as.integer(nlont), + id = as.integer(id), + slat = as.numeric(slat), + slon = as.numeric(slon), + rlat = as.numeric(rlat), + rlon = as.numeric(rlon), + slatt = as.numeric(slatt), + slont = as.numeric(slont), + nd = as.integer(nd), + um = matrix(as.double(seq(1,nd*ic)),c(nd,ic)), + vm = matrix(as.double(seq(1,nd*ic)),c(nd,ic)), + gu92 = matrix(as.double(seq(1,nd*ic)),c(nd,ic)), + gv92 = matrix(as.double(seq(1,nd*ic)),c(nd,ic)), + gu52 = matrix(as.double(seq(1,nd*ic)),c(nd,ic)), + gv52 = matrix(as.double(seq(1,nd*ic)),c(nd,ic)), + nger = as.integer(1), + PACKAGE = 'CSTools') + + a <- prePro$um + b <- prePro$vm + c <- prePro$gu92[1:prePro$nger,] + d <- prePro$gv92[1:prePro$nger,] + e <- prePro$gu52[1:prePro$nger,] + f <- prePro$gv52[1:prePro$nger,] + + g <- prePro$nger + + predSig <- .Fortran("training_part2_prec", + u500 = as.numeric(u500), + v500 = as.numeric(v500), + t500 = as.numeric(t500), + t850 = as.numeric(t850), + msl_si = as.numeric(msl_si), + q700 = as.numeric(q700), + lon_hr = as.numeric(lon_hr), + lat_hr = as.numeric(lat_hr), + prec_hr = as.numeric(prec_hr), + nanx = as.integer(nanx), + nlat = as.integer(nlat), + nlon = as.integer(nlon), + ic = as.integer(ic), + nlatt = as.integer(nlatt), + nlont = as.integer(nlont), + id = as.integer(id), + slat = as.numeric(slat), + slon = as.numeric(slon), + rlat = as.numeric(rlat), + rlon = as.numeric(rlon), + slatt = as.numeric(slatt), + slont = as.numeric(slont), + nd = as.integer(nd), + um = as.double(a), + vm = as.double(b), + gu92 = as.double(c), + gv92 = as.double(d), + gu52 = as.double(e), + gv52 = as.double(f), + nger = as.integer(g), + vdmin = matrix(as.double(seq(1,nptos*4)),c(nptos,4)), + vref = matrix(as.integer(seq(1,nptos*4)),c(nptos,4)), + neni = as.integer(1), + mi = matrix(as.integer(seq(1,prePro$nger*nptos)),c(prePro$nger,nptos)), + ccm = matrix(as.double(seq(1,prePro$nger*nptos)),c(prePro$nger,nptos)), + lab_pred = matrix(as.integer(seq(1,prePro$nger*nptos*npx)),c(prePro$nger,nptos,npx)), + cor_pred = matrix(as.double(seq(1,prePro$nger*nptos*npx)),c(prePro$nger,nptos,npx)), + PACKAGE = 'CSTools') + + h <- predSig$mi + i <- predSig$ccm + j <- predSig$lab_pred + k <- predSig$cor_pred + l <- predSig$vdmin + m <- predSig$vref + + indices <- array(c(j,k),c(g,nptos,npx,2)) + dimnames(indices)[[4]] <- c("lab_pred","cor_pred") + + output <- list("um" = a, + "vm" = b, + "nger" = g, + "gu92" = c, + "gv92" = d, + "gu52" = e, + "gv52" = f, + "neni" = predSig$neni, + "vdmin" = l, + "vref" = m, + "ccm" = i, + "indices" = indices) +} else { + + prePro <- .Fortran("training_temp", + t1000 = as.numeric(t1000), + msl_si = as.numeric(msl_si), + msl_lr = as.numeric(msl_lr), + lon_hr = as.numeric(lon_hr), + lat_hr = as.numeric(lat_hr), + ngridd = as.integer(ngridd), + nlat = as.integer(nlat), + nlon = as.integer(nlon), + ic = as.integer(ic), + nlatt = as.integer(nlatt), + nlont = as.integer(nlont), + id = as.integer(id), + slat = as.numeric(slat), + slon = as.numeric(slon), + rlat = as.numeric(rlat), + rlon = as.numeric(rlon), + slatt = as.numeric(slatt), + slont = as.numeric(slont), + nd = as.integer(nd), + day = as.integer(day), + month = as.integer(month), + um = matrix(as.double(seq(1,nd*ic)),c(nd,ic)), + vm = matrix(as.double(seq(1,nd*ic)),c(nd,ic)), + insol = vector(mode="double",length=nd), + vdmin = matrix(as.double(seq(1,nptos*4)),c(nptos,4)), + vref = matrix(as.integer(seq(1,nptos*4)),c(nptos,4)), + neni = as.integer(1), + PACKAGE = 'CSTools') + + a <- prePro$um + b <- prePro$vm + c <- prePro$insol + d <- prePro$vdmin + e <- prePro$vref + f <- prePro$neni + + output <- list("um" = a, + "vm" = b, + "insol" = c, + "vdmin" = d, + "vref" = e, + "neni" = f) + +} + + return(output) + +} + diff --git a/R/CST_AnalogsPredictors.R b/R/CST_AnalogsPredictors.R new file mode 100644 index 0000000000000000000000000000000000000000..a15a4c0cab79b383293dcbd4c34de41b95a8a19c --- /dev/null +++ b/R/CST_AnalogsPredictors.R @@ -0,0 +1,847 @@ +#' AEMET Downscaling +#' Precipitation and maximum and minimum temperature downscaling method +#' based on analogs: synoptic situations and significant predictors. +#' +#'@author Marta Dominguez Alonso - AEMET, \email{mdomingueza@aemet.es} +#'@author Nuria Perez-Zanon - BSC, \email{nuria.perez@bsc.es} +#' +#'@description This function downscales low resolution precipitation data (e.g. from +#' Seasonal Forecast Models) through the association with an observational high +#' resolution (HR) dataset (AEMET 5 km gridded data of daily precipitation (Peral et al., 2017)) +#' and a collection of predictors and past synoptic situations similar to estimated day. +#' The method uses three domains: +#' - peninsular Spain and Balearic Islands domain (5 km resolution): HR precipitation +#' and the downscaling result domain. +#' - synoptic domain (low resolution, e.g. 1.5º x 1.5º): it should be centered over Iberian Peninsula +#' and cover enough extension to detect as much synoptic situations as possible. +#' - extended domain (low resolution, e.g. 1.5º x 1.5º): it should have the same resolution +#' as synoptic domain. It is used for SLP Seasonal Forecast Models. +#'@param exp List of arrays with downscaled period seasonal forecast data. The list +#' has to contain model atmospheric variables (instantaneous 12h data) that must +#' be indentify by parenthesis name. +#' For precipitation: +#' - u component of wind at 500 hPa (u500_mod) in m/s +#' - v component of wind at 500 hPa (v500_mod) in m/s +#' - temperature at 500 hPa (t500_mod) in K +#' - temperature at 850 hPa (t850_mod) in K +#' - specific humidity at 700 hPa (q700_mod) in g/kg +#' For temperature: +#' - u component of wind at 500 hPa (u500_mod) in m/s +#' - v component of wind at 500 hPa (v500_mod) in m/s +#' - temperature at 500 hPa (t500_mod) in K +#' - temperature at 700 hPa (t700_mod) in K +#' - temperature at 850 hPa (t850_mod) in K +#' - specific humidity at 700 hPa (q700_mod) in g/kg +#' - 2 meters temperature (tm2m_mod) in K +#' The arrays must have at least three dimensions with names 'lon', 'lat' and 'time'. +#' (lon = gridpoints of longitude, lat = gridpoints of latitude, time = number of downscaling days) +#' Seasonal forecast variables must have the same resolution and +#' domain as reanalysis variables ('obs' parameter, below). +#'@param slp Array with atmospheric seasonal forecast model sea level pressure +#' (instantaneous 12h data) that must be indentify as 'slp' (hPa). It has the same +#' resolution as 'exp' and 'obs' paremeters but with an extended domain. +#' This domain contains extra degrees (most in the north and west part) compare to +#' synoptic domain. The array must have at least three dimensions +#' with names 'lon', 'lat' and 'time'. +#'@param obs List of arrays with training period reanalysis data. +#' The list has to contain reanalysis atmospheric variables (instantaneous +#' 12h data) that must be indentify by parenthesis name. +#' For precipitation: +#' - u component of wind at 500 hPa (u500) in m/s +#' - v component of wind at 500 hPa (v500) in m/s +#' - temperature at 500 hPa (t500) in K +#' - temperature at 850 hPa (t850) in K +#' - sea level pressure (slp) in hPa +#' - specific humidity at 700 hPa (q700) in g/kg +#' For maximum and minimum temperature: +#' - u component of wind at 500 hPa (u500) in m/s +#' - v component of wind at 500 hPa (v500) in m/s +#' - temperature at 500 hPa (t500) in K +#' - temperature at 700 hPa (t700) in K +#' - temperature at 850 hPa (t850) in K +#' - sea level pressure (slp) in hPa +#' - specific humidity at 700 hPa (q700) in g/kg +#' - 2 meters temperature (tm2m) in K +#' The arrays must have at least three dimensions with names 'lon', 'lat' and 'time'. +#'@param lon Vector of the synoptic longitude (from (-180º) to 180º), +#' The vector must go from west to east. The same as for the training function. +#'@param lat Vector of the synoptic latitude. The vector must go from north to south. +#' The same as for the training function. +#'@param slp_lon Vector of the extended longitude (from (-180º) to 180º), +#' The vector must go from west to east. The same as for the training function. +#'@param slp_lat Vector of the extended latitude. The vector must go from north to south. +#' The same as for the training function. +#'@param var_name Variable name to downscale. There are two options: 'prec' for +#' precipitation and 'temp' for maximum and minimum temperature. +#'@param hr_obs Local path of HR observational files (maestro and pcp/tmx-tmn). +#' For precipitation can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/v2/Serie_AEMET_v2_pcp_1951a202006_txt.tar.gz +#' For maximum and minimum temperature can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmax/Serie_AEMET_v1_tmax_1951a202006_txt.tar.gz and http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmin/Serie_AEMET_v1_tmin_1951a202006_txt.tar.gz respetively. +#' Maestro file (maestro_red_hr_SPAIN.txt) has gridpoint (nptos), longitude (lon), latitude (lat) and +#' altitude (alt) in columns (vector structure). +#' Data file (pcp/tmx/tmn_red_SPAIN_1951-201903.txt) includes 5km resolution spanish daily data +#' (precipitation or maximum and minimum temperature from january 1951 to june 2020. See README +#' file for more information. +#' IMPORTANT!: HR observational period must be the same as for reanalysis variables. +#' It is assumed that the training period is smaller than the HR original one (1951-2019), so it is +#' needed to make a new ascii file with the new period and the same structure as original, +#' specifying the training dates in the name (e.g. 'pcp_red_SPAIN_19810101-19961231.txt' for +#' '19810101-19961231' period). +#'@param tdates Training period dates in format YYYYMMDD(start)-YYYYMMDD(end) +#' (e.g. 19810101-20181231). +#'@param ddates Downscaling period dates in format YYYYMMDD(start)-YYYYMMDD(end) (e.g. 20191001-20200331). +#'@param restrain Output (list of matrix) obtained from 'training_analogs' function. +#' For precipitation, 'restrain' object must contains um, vm, nger, gu92, gv92, +#' gu52, gv52, neni, vdmin, vref, ccm, lab_pred and cor_pred variables. +#' For maximum and minimum temperature, 'restrain' object must contains um, vm, +#' insol, neni, vdmin y vref. See 'AnalogsPred_train.R' for more information. +#'@param dim_name_longitude A character string indicating the name of the longitude +#'dimension, by default 'longitude'. +#'@param dim_name_latitude A character string indicating the name of the latitude +#'dimension, by default 'latitude'. +#'@param dim_name_time A character string indicating the name of the time +#'dimension, by default 'time'. +#'@return Matrix with seasonal forecast precipitation (mm) or +#' maximum and minimum temperature (dozens of ºC) in a 5km x 5km regular grid +#' over peninsular Spain and Balearic Islands. The resulted matrices have two +#' dimensions ('ddates' x 'nptos').(ddates = number of downscaling days +#' and nptos = number of 'hr_obs' gridpoints). +#' +#'@useDynLib CSTools +#' +#'@export +#' +CST_AnalogsPredictors <- function(exp, + slp, + obs, + lon, + lat, + slp_lon, + slp_lat, + var_name, + hr_obs, + tdates, + ddates, + restrain, + dim_name_longitude = "lon", + dim_name_latitude = "lat", + dim_name_time = "time") { + +if (!is.list(exp)) { + stop("Parameter 'exp' must be a list of 'array' objects") + } + +if (!(all(sapply(exp, inherits, 'array')))) { + stop("Elements of the list in parameter 'exp' must be of the class ", + "'array'.") + } + +if (!is.array(slp)) { + stop("Parameter 'slp' must be of the class 'array'.") + } + +if (!is.list(obs)) { + stop("Parameter 'obs' must be a list of 'array' objects") + } + +if (!(all(sapply(obs, inherits, 'array')))) { + stop("Elements of the list in parameter 'obs' must be of the class ", + "'array'.") + } + +if (var_name == "prec") { + if (length(exp) != 5) { + stop("Parameter 'exp' must be a length of 5.") + } else { + if (!(any(names(exp) %in% "u500_mod"))) { + stop("Variable 'u500_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "v500_mod"))) { + stop("Variable 'v500_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "t500_mod"))) { + stop("Variable 't500_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "t850_mod"))) { + stop("Variable 't850_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "q700_mod"))) { + stop("Variable 'q700_mod' in 'exp' parameter is missed.") + } + } + if (length(obs) != 6) { + stop("Parameter 'obs' must be a length of 6.") + } else { + if (!(any(names(obs) %in% "u500"))) { + stop("Variable 'u500' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "v500"))) { + stop("Variable 'v500' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "t500"))) { + stop("Variable 't500' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "t850"))) { + stop("Variable 't850' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "slp"))) { + stop("Variable 'slp' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "q700"))) { + stop("Variable 'q700' in 'obs' parameter is missed.") + } + } +} else { + if (length(exp) != 7) { + stop("Parameter 'exp' must be a length of 7.") + } else { + if (!(any(names(exp) %in% "u500_mod"))) { + stop("Variable 'u500_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "v500_mod"))) { + stop("Variable 'v500_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "t500_mod"))) { + stop("Variable 't500_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "t700_mod"))) { + stop("Variable 't700_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "t850_mod"))) { + stop("Variable 't850_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "q700_mod"))) { + stop("Variable 'q700_mod' in 'exp' parameter is missed.") + } else if (!(any(names(exp) %in% "tm2m_mod"))) { + stop("Variable 'tm2m_mod' in 'exp' parameter is missed.") + } + } + if (length(obs) != 8) { + stop("Parameter 'obs' must be a length of 8.") + } else { + if (!(any(names(obs) %in% "u500"))) { + stop("Variable 'u500' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "v500"))) { + stop("Variable 'v500' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "t500"))) { + stop("Variable 't500' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "t700"))) { + stop("Variable 't700' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "t850"))) { + stop("Variable 't850' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "slp"))) { + stop("Variable 'slp' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "q700"))) { + stop("Variable 'q700' in 'obs' parameter is missed.") + } else if (!(any(names(obs) %in% "tm2m"))) { + stop("Variable 'tm2m' in 'obs' parameter is missed.") + } + } +} + +if (all((sapply(exp,dim))==dim(exp[[1]]))) { + dim_exp <- dim(exp[[1]]) + if (!(any(names(dim_exp) %in% dim_name_longitude))) { + stop("Dimension 'lon' in exp parameter is missed.") + } + if (!(any(names(dim_exp) %in% dim_name_latitude))) { + stop("Dimension 'lat' in exp parameter is missed.") + } + if (!(any(names(dim_exp) %in% dim_name_time))) { + stop("Dimension 'time' in exp parameter is missed.") + } +} else { + stop("All 'exp' variables must have the same dimensions.") +} + +dim_slp <- dim(slp) +if (!(any(names(dim_slp) %in% dim_name_longitude))) { + stop("Dimension 'lon' in slp parameter is missed.") +} +if (!(any(names(dim_slp) %in% dim_name_latitude))) { + stop("Dimension 'lat' in slp parameter is missed.") +} +if (!(any(names(dim_slp) %in% dim_name_time))) { + stop("Dimension 'time' in slp parameter is missed.") +} + +if (all((sapply(obs,dim))==dim(obs[[1]]))) { + dim_obs <- dim(obs[[1]]) + if (!(any(names(dim_obs) %in% dim_name_longitude))) { + stop("Dimension 'lon' in obs parameter is missed.") + } + if (!(any(names(dim_obs) %in% dim_name_latitude))) { + stop("Dimension 'lat' in obs parameter is missed.") + } + if (!(any(names(dim_obs) %in% dim_name_time))) { + stop("Dimension 'time' in obs parameter is missed.") + } +} else { + stop("All 'obs' variables must have the same dimensions.") +} + +if (!is.vector(lon) || !is.numeric(lon)) { + stop("Parameter 'lon' must be a numeric vector") +} else { + if (is.unsorted(lon)) { + lon <- sort(lon) + warning("'lon' vector has been sorted in increasing order") + } +} + +if (!is.vector(lat) || !is.numeric(lat)) { + stop("Parameter 'lat' must be a numeric vector") +} else { + if (!is.unsorted(lat)) { + lat <- sort(lat, decreasing = TRUE) + warning("'lat' vector has been sorted in decreasing order") + } +} + +if (!is.vector(slp_lon) || !is.numeric(slp_lon)) { + stop("Parameter 'slp_lon' must be a numeric vector") +} else { + if (is.unsorted(slp_lon)) { + lon <- sort(slp_lon) + warning("'slp_lon' vector has been sorted in increasing order") + } +} + +if (!is.vector(slp_lat) || !is.numeric(slp_lat)) { + stop("Parameter 'slp_lat' must be a numeric vector") +} else { + if (!is.unsorted(slp_lat)) { + lat <- sort(slp_lat, decreasing = TRUE) + warning("'slp_lat' vector has been sorted in decreasing order") + } +} + +if (!is.character(hr_obs)){ + stop("Parameter 'hr_obs' must be a character.") +} else { + if (!dir.exists(hr_obs)) { + stop("'hr_obs' directory does not exist") + } +} + +if (!is.character(tdates)) { + stop("Parameter 'tdates' must be a character.") +} else { + if (nchar(tdates) != "17") { + stop("Parameter 'tdates' must be a string with 17 charecters.") + } else { + dateini <- as.Date(substr(tdates,start=1,stop=8),format="%Y%m%d") + dateend <- as.Date(substr(tdates,start=10,stop=18),format="%Y%m%d") + if (dateend <= dateini) { + stop("Parameter 'tdates' must be at least of one day") + } + } +} + +if (!is.character(ddates)) { + stop("Parameter 'ddates' must be a character.") +} else { + if (nchar(ddates) != "17") { + stop("Parameter 'ddates' must be a string with 17 charecters.") + } else { + dateini <- as.Date(substr(ddates,start=1,stop=8),format="%Y%m%d") + dateend <- as.Date(substr(ddates,start=10,stop=18),format="%Y%m%d") + if (dateend <= dateini) { + stop("Parameter 'ddates' must be at least of one day") + } + } +} + +# + +if (names(dim(exp[[1]]))[1] == "lon" & names(dim(exp[[1]]))[2] == "lat" + || names(dim(exp[[1]]))[2] == "lon" & names(dim(exp[[1]]))[3] == "lat") { + texp2D <- lapply(exp, MergeDims, merge_dims = c('lon', 'lat'), + rename_dim = 'gridpoint') +} else if (names(dim(exp[[1]]))[1] == "lat" & names(dim(exp[[1]]))[2] == "lon" + || names(dim(exp[[1]]))[2] == "lat" & names(dim(exp[[1]]))[3] == "lon") { + texp2D <- lapply(exp, MergeDims, merge_dims = c('lat', 'lon'), + rename_dim = 'gridpoint') +} + +if (names(dim(slp))[1] == "lon" & names(dim(slp))[2] == "lat" + || names(dim(slp))[2] == "lon" & names(dim(slp))[3] == "lat") { + tslp2D <- MergeDims(slp,merge_dims = c('lon', 'lat'), + rename_dim = 'gridpoint') +} else if (names(dim(slp))[1] == "lat" & names(dim(slp))[2] == "lon" + || names(dim(slp))[2] == "lat" & names(dim(slp))[3] == "lon") { + tslp2D <- MergeDims(slp,merge_dims = c('lat', 'lon'), + rename_dim = 'gridpoint') +} + +if (names(dim(obs[[1]]))[1] == "lon" & names(dim(obs[[1]]))[2] == "lat" + || names(dim(obs[[1]]))[2] == "lon" & names(dim(obs[[1]]))[3] == "lat") { + tobs2D <- lapply(obs, MergeDims, merge_dims = c('lon', 'lat'), + rename_dim = 'gridpoint') +} else if (names(dim(obs[[1]]))[1] == "lat" & names(dim(obs[[1]]))[2] == "lon" + || names(dim(obs[[1]]))[2] == "lat" & names(dim(obs[[1]]))[3] == "lon") { + tobs2D <- lapply(obs, MergeDims, merge_dims = c('lat', 'lon'), + rename_dim = 'gridpoint') +} + +if (names(dim(texp2D[[1]]))[1] == "gridpoint") { + exp2D <- lapply(texp2D,aperm) +} else { + exp2D <- texp2D +} + +if (names(dim(tslp2D))[1] == "gridpoint") { + slp2D <- aperm(tslp2D) +} else { + slp2D <- tslp2D +} + +if (names(dim(tobs2D[[1]]))[1] == "gridpoint") { + obs2D <- lapply(tobs2D,aperm) +} else { + obs2D <- tobs2D +} + + downres <- .analogspred(exp2D, + slp2D, + obs2D, + lon, + lat, + slp_lon, + slp_lat, + var_name, + hr_obs, + tdates, + ddates, + restrain) + +} + +#' Atomic .analogspred function +#' +#'@author Marta Dom\'inguez Alonso - AEMET, \email{mdomingueza@aemet.es} +#' +#' This function works with lists of matrix from reanalysis and seasonal +#' forecast data and uses a Fortran interface (.Fortran) to run an +#' analogs method developed in AEMET. +#'@param pred_mod List of matrix with downscaled period seasonal forecast data. The list +#' has to contain model atmospheric variables (instantaneous 12h data) that must +#' be indentify by parenthesis name. +#' For precipitation: +#' - u component of wind at 500 hPa (u500_mod) in m/s +#' - v component of wind at 500 hPa (v500_mod) in m/s +#' - temperature at 500 hPa (t500_mod) in K +#' - temperature at 850 hPa (t850_mod) in K +#' - specific humidity at 700 hPa (q700_mod) in g/kg +#' For temperature: +#' - u component of wind at 500 hPa (u500_mod) in m/s +#' - v component of wind at 500 hPa (v500_mod) in m/s +#' - temperature at 500 hPa (t500_mod) in K +#' - temperature at 700 hPa (t500_mod) in K +#' - temperature at 850 hPa (t850_mod) in K +#' - specific humidity at 700 hPa (q700_mod) in g/kg +#' - 2 meters temperature (tm2m_mod) in K +#' Seasonal forecast variables must have the same resolution and +#' domain as 'pred_rea' parameter. +#' All matrices must have two dimensions with names 'time' and 'gridpoint'. +#'@param pred_slp Matrix with atmospheric seasonal forecast model sea level pressure +#' (instantaneous 12h data) that must be indentify as 'slp'. It has the same +#' resolution as 'pred_mod' paremeter but with an extended domain. This domain contains +#' extra degrees (most in the north and west part) compare to synoptic domain. +#' The matrix must have two dimensions with names 'time' and 'gridpoint'. +#'@param pred_rea List of matrix with training period reanalysis data. +#' The list has to contain reanalysis atmospheric variables (instantaneous +#' 12h data) that must be indentify by parenthesis name. +#' For precipitation: +#' - u component of wind at 500 hPa (u500) in m/s +#' - v component of wind at 500 hPa (v500) in m/s +#' - temperature at 500 hPa (t500) in K +#' - temperature at 850 hPa (t850) in K +#' - sea level pressure (slp) in hPa +#' - specific humidity at 700 hPa (q700) in g/kg +#' For maximum and minimum temperature: +#' - u component of wind at 500 hPa (u500) in m/s +#' - v component of wind at 500 hPa (v500) in m/s +#' - temperature at 500 hPa (t500) in K +#' - temperature at 700 hPa (t500) in K +#' - temperature at 850 hPa (t850) in K +#' - sea level pressure (slp) in hPa +#' - specific humidity at 700 hPa (q700) in g/kg +#' - 2 meters temperature (tm2m) in K +#' All matrices must have two dimensions with names 'ddates' and 'gridpoint'. +#'@param lon Vector of the synoptic longitude (from (-180º) to 180º), +#' The vector must go from west to east. +#'@param lat Vector of the synoptic latitude. The vector must go from north to south. +#'@param slp_lon Vector of the extended longitude (from (-180º) to 180º), +#' The vector must go from west to east. +#'@param slp_lat Vector of the extended latitude. The vector must go from north to south. +#'@param var Variable name to downscale. There are two options: 'prec' for +#' precipitation and 'temp' for maximum and minimum temperature. +#'@param HR_path Local path of HR observational files (maestro and pcp/tmx-tmn). +#' For precipitation can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/v2/Serie_AEMET_v2_pcp_1951a201903_txt.tar.gz +#' For maximum and minimum temperature can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmax/Serie_AEMET_v1_tmax_1951a202006_txt.tar.gz and http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmin/Serie_AEMET_v1_tmin_1951a202006_txt.tar.gz respetively. +#' Maestro file (maestro_red_hr_SPAIN.txt) has gridpoint (nptos), longitude (lon), latitude (lat) and +#' altitude (alt) in columns (vector structure). +#' Data file (pcp/tmx/tmn_red_SPAIN_1951-201903.txt) includes 5km resolution spanish daily data +#' (precipitation or maximum and minimum temperature from january 1951 to march 2019. See README +#' file for more information. +#' IMPORTANT!: HR observational period must be the same as for reanalysis variables +#' ('pred_rea' parameter). +#' It is assumed that the training period is smaller than the HR original one (1951-2019), so it is +#' needed to make a new ascii file with the new period and the same structure as original, +#' specifying the training dates in the name (e.g. 'pcp_red_SPAIN_19810101-19961231.txt' for +#' '19810101-19961231' period). +#'@param tdates Training period dates in format YYYYMMDD(start)-YYYYMMDD(end) +#' (e.g. 19810101-20181231). The same as for the training function. +#'@param ddates Downscaling period dates in format YYYYMMDD(start)-YYYYMMDD(end) (e.g. 20191001-20200331). +#'@param restrain Output (list of matrix) obtained from 'training_analogs' function. +#' For precipitation, 'restrain' object must contains um, vm, nger, gu92, gv92, +#' gu52, gv52, neni, vdmin, vref, ccm, lab_pred and cor_pred variables. +#' For maximum and minimum temperature, 'restrain' object must contains um, vm, +#' insol, neni, vdmin y vref. See 'AnalogsPred_train.R' for more information. +#'@return .analogspred returns seasonal forecast precipitation (mm) or +#' maximum and minimum temperature (dozens of ºC) in a 5km x 5km regular grid over +#' peninsular Spain and Balearic Islands. Each matrix of the list has two dimensions +#' ('ddates' x 'nptos'). +#' +#'@importFrom utils read.table +#' +#'@useDynLib CSTools +#'@noRd + + +.analogspred <- function(pred_mod, + pred_slp, + pred_rea, + lon, + lat, + slp_lon, + slp_lat, + var, + HR_path, + tdates, + ddates, + restrain) { + + +if (!is.list(pred_mod)) { + stop("Parameter 'pred_mod' must be a list of 'matrix' objects") + } + +if (!(all(sapply(pred_mod, inherits, 'matrix')))) { + stop("Elements of the list in parameter 'pred_mod' must be of the class ", + "'matrix'.") + } + +if (!is.matrix(pred_slp)) { + stop("Parameter 'pred_slp' must be of the class 'matrix'.") + } + +if (!is.list(pred_rea)) { + stop("Parameter 'pred_rea' must be a list of 'matrix' objects") + } + +if (!(all(sapply(pred_rea, inherits, 'matrix')))) { + stop("Elements of the list in parameter 'pred_rea' must be of the class ", + "'matrix'.") + } + +if (var == "prec") { + if (length(pred_rea) != 6) { + stop("Parameter 'pred_rea' must be a length of 6.") + } + if (length(pred_mod) != 5) { + stop("Parameter 'pred_mod' must be a length of 5.") + } +} else { + if (length(pred_rea) != 8) { + stop("Parameter 'pred_rea' must be a length of 8.") + } + if (length(pred_mod) != 7) { + stop("Parameter 'pred_mod' must be a length of 7.") + } +} + +if (!is.vector(lon) || !is.numeric(lon)) { + stop("Parameter 'lon' must be a numeric vector") + } + +if (!is.vector(lat) || !is.numeric(lat)) { + stop("Parameter 'lat' must be a numeric vector") + } + +if (!is.vector(slp_lon) || !is.numeric(slp_lon)) { + stop("Parameter 'slp_lon' must be a numeric vector") + } + +if (!is.vector(slp_lat) || !is.numeric(slp_lat)) { + stop("Parameter 'slp_lat' must be a numeric vector") + } + +if (!is.character(HR_path)){ + stop("Parameter 'HR_path' must be a character.") + } + +if (!is.character(tdates)) { + stop("Parameter 'tdates' must be a character.") + } + +if (!is.character(ddates)) { + stop("Parameter 'ddates' must be a character.") + } + +if (!is.list(restrain)) { + stop("Parameter 'restrain' must be a list of 'matrix' and 'parameter' objects") + } + +#! REANALYSIS GRID PARAMETERS + + rlon <- c(lon, NA) - c(NA, lon) + rlon <- rlon[!is.na(rlon)] + if (!all(rlon == rlon[1])) { + stop("Parameter 'lon' must be in regular grid.") + } else { + rlon <- rlon[1] + } + + rlat <- c(lat, NA) - c(NA, lat) + rlat <- rlat[!is.na(rlat)] + if (!all(rlat == rlat[1])) { + stop("Parameter 'lat' must be in regular grid.") + } else { + rlat <- rlat[1] + } + + if (rlon != (-rlat)) { + stop("Parameters 'lon' and 'lat' must have the same resolution.") + } else { + res <- rlon + } + + nlat <- ((lat[length(lat)] - lat[1]) / rlat) + 1 + nlon <- ((lon[length(lon)] - lon[1]) / rlon) + 1 + + ic <- nlat * nlon +# + slp_rlon <- c(slp_lon, NA) - c(NA, slp_lon) + slp_rlon <- slp_rlon[!is.na(slp_rlon)] + if (!all(slp_rlon == slp_rlon[1])) { + stop("Parameter 'slp_lon' must be in regular grid.") + } else { + slp_rlon <- slp_rlon[1] + } + + slp_rlat <- c(slp_lat, NA) - c(NA, slp_lat) + slp_rlat <- slp_rlat[!is.na(slp_rlat)] + if (!all(slp_rlat == slp_rlat[1])) { + stop("Parameter 'slp_lat' must be in regular grid.") + } else { + slp_rlat <- slp_rlat[1] + } + + if (slp_rlon != (-slp_rlat)) { + stop("Parameters 'slp_lon' and 'slp_lat' must have the same resolution.") + } else { + slp_res <- slp_rlon + } + + nlatt <- ((slp_lat[length(slp_lat)] - slp_lat[1]) / slp_rlat) + 1 + nlont <- ((slp_lon[length(slp_lon)] - slp_lon[1]) / slp_rlon) + 1 + + id <- nlatt * nlont + + slat <- max(lat) + slon <- min(c(lon[which(lon > 180)] - 360, + lon[which(lon <= 180)])) + + slatt <- max(slp_lat) + slont <- min(c(slp_lon[which(slp_lon > 180)] - 360, + slp_lon[which(slp_lon <= 180)])) + + ngridd <- ((2*nlatt)-1)*((2*nlont)-1) + + if (all((sapply(pred_rea,nrow))==nrow(pred_rea[[1]]))){ + nd <- nrow(pred_rea[[1]]) + } else { + stop("All 'pred_rea' variables must have the same period.") + } + + if (all((sapply(pred_mod,nrow))==nrow(pred_mod[[1]]))){ + nm <- nrow(pred_mod[[1]]) + } else { + stop("All 'pred_mod' variables must have the same period.") + } + + seqdates <- seq(as.Date(substr(ddates,start=1,stop=8),format="%Y%m%d"),as.Date(substr(ddates,start=10,stop=18),format="%Y%m%d"),by="days") + month <- format(seqdates,format="%m") + day <- format(seqdates,format="%d") + +#! TRAINING REANALYSIS VARIABLES +u500 <- pred_rea[['u500']] +v500 <- pred_rea[['v500']] +t500 <- pred_rea[['t500']] +t850 <- pred_rea[['t850']] +msl_si <- pred_rea[['slp']] +q700 <- pred_rea[['q700']] + +if (var == "temp") { +t700 <- pred_rea[['t700']] +tm2m <- pred_rea[['tm2m']] +} + +#! SEASONAL FORECAST MODEL VARIABLES +u500_mod <- pred_mod[['u500_mod']] +v500_mod <- pred_mod[['v500_mod']] +t500_mod <- pred_mod[['t500_mod']] +t850_mod <- pred_mod[['t850_mod']] +msl_lr_mod <- pred_slp +q700_mod <- pred_mod[['q700_mod']] + +if (var == "temp") { +t700_mod <- pred_mod[['t700_mod']] +tm2m_mod <- pred_mod[['tm2m_mod']] +} + +#! HIGH-RESOLUTION (HR) OBSERVATIONAL DATASET +maestro_hr_file <- paste(HR_path, "maestro_red_hr_SPAIN.txt",sep="") +if (!file.exists(maestro_hr_file)) { + stop("'maestro_red_hr_SPAIN.txt' does not exist.") +} else { + maestro <- read.table(maestro_hr_file) + lon_hr <- unlist(maestro[2]) + lat_hr <- unlist(maestro[3]) + nptos <- length(readLines(maestro_hr_file)) +} + +if (var == "prec") { + prec_hr_file <- paste(HR_path, "pcp_red_SPAIN_",tdates,".txt",sep="") + if (!file.exists(prec_hr_file)) { + stop(sprintf("precipitation HR file for %s does not exist.",tdates)) + } else { + nd_hr <- length(readLines(prec_hr_file)) + preprec_hr <- matrix(scan(prec_hr_file), nrow=nd_hr ,ncol= nptos+1, byrow=TRUE) + prec_hr <- preprec_hr[1:nd_hr,-c(1)] + } +} else { + tmx_hr_file <- paste(HR_path, "tmx_red_SPAIN_",tdates,".txt",sep="") + tmn_hr_file <- paste(HR_path, "tmn_red_SPAIN_",tdates,".txt",sep="") + if (!file.exists(tmx_hr_file)) { + stop(sprintf("maximum temperature HR file for %s does not exist.",tdates)) + } else if (!file.exists(tmn_hr_file)) { + stop(sprintf("minimum temperature HR file for %s does not exist.",tdates)) + } else if (length(readLines(tmx_hr_file)) != length(readLines(tmn_hr_file))) { + stop("maximum and minimum temperature HR observation files must have the same period.") + } else { + nd_hr <- length(readLines(tmx_hr_file)) + pretmx_hr <- matrix(scan(tmx_hr_file), nrow=nd_hr ,ncol= nptos+1, byrow=TRUE) + tmx_hr <- pretmx_hr[1:nd_hr,-c(1)] + pretmn_hr <- matrix(scan(tmn_hr_file), nrow=nd_hr ,ncol= nptos+1, byrow=TRUE) + tmn_hr <- pretmn_hr[1:nd_hr,-c(1)] + } +} + + if (nd_hr != nd) { + stop("Reanalysis variables and HR observations must have the same period.") + } + +#! OTHER PARAMETERS that should not be changed +#! Number of analog situations to consider +nanx <- 155 +#! Number of temperature predictors +nvar <- 7 + +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +if (var == "prec") { + + downs <- .Fortran("down_prec", + ic = as.integer(ic), + id = as.integer(id), + nd = as.integer(nd), + nm = as.integer(nm), + nlat = as.integer(nlat), + nlon = as.integer(nlon), + nlatt = as.integer(nlatt), + nlont = as.integer(nlont), + slat = as.numeric(slat), + slon = as.numeric(slon), + rlat = as.numeric(rlat), + rlon = as.numeric(rlon), + slatt = as.numeric(slatt), + slont = as.numeric(slont), + ngridd = as.integer(ngridd), + u500 = as.numeric(u500), + v500 = as.numeric(v500), + t500 = as.numeric(t500), + t850 = as.numeric(t850), + msl_si = as.numeric(msl_si), + q700 = as.numeric(q700), + prec_hr = as.numeric(prec_hr), + nanx = as.integer(nanx), + restrain$um, + restrain$vm, + restrain$nger, + restrain$gu92, + restrain$gv92, + restrain$gu52, + restrain$gv52, + restrain$neni, + restrain$vdmin, + restrain$vref, + restrain$ccm, + restrain$indices[,,,1],#lab_pred + restrain$indices[,,,2],#cor_pred + u500_mod = as.numeric(u500_mod), + v500_mod = as.numeric(v500_mod), + t500_mod = as.numeric(t500_mod), + t850_mod = as.numeric(t850_mod), + msl_lr_mod = as.numeric(msl_lr_mod), + q700_mod = as.numeric(q700_mod), + pp=matrix(as.double(seq(1,nm*nptos)),c(nm,nptos)), + PACKAGE = 'CSTools') + + output <- downs$pp + +} else { + + downs <- .Fortran("down_temp", + ic = as.integer(ic), + id = as.integer(id), + nd = as.integer(nd), + nm = as.integer(nm), + nlat = as.integer(nlat), + nlon = as.integer(nlon), + nlatt = as.integer(nlatt), + nlont = as.integer(nlont), + slat = as.numeric(slat), + slon = as.numeric(slon), + rlat = as.numeric(rlat), + rlon = as.numeric(rlon), + slatt = as.numeric(slatt), + slont = as.numeric(slont), + ngridd = as.integer(ngridd), + u500 = as.numeric(u500), + v500 = as.numeric(v500), + t500 = as.numeric(t500), + t850 = as.numeric(t850), + msl_si = as.numeric(msl_si), + q700 = as.numeric(q700), + t700 = as.numeric(t700), + tm2m = as.numeric(tm2m), + tmx_hr = as.numeric(tmx_hr), + tmn_hr = as.numeric(tmn_hr), + nanx = as.integer(nanx), + nvar = as.integer(nvar), + day = as.integer(day), + month = as.integer(month), + restrain$um, + restrain$vm, + restrain$insol, + restrain$neni, + restrain$vdmin, + restrain$vref, + u500_mod = as.numeric(u500_mod), + v500_mod = as.numeric(v500_mod), + t500_mod = as.numeric(t500_mod), + t850_mod = as.numeric(t850_mod), + msl_lr_mod = as.numeric(msl_lr_mod), + q700_mod = as.numeric(q700_mod), + t700_mod = as.numeric(t700_mod), + tm2m_mod = as.numeric(tm2m_mod), + tmx=matrix(as.double(seq(1,nm*nptos)),c(nm,nptos)), + tmn=matrix(as.double(seq(1,nm*nptos)),c(nm,nptos)), + PACKAGE = 'CSTools') + + output <- list("tmax" = downs$tmx, + "tmin" = downs$tmn) + +} + return(output) +} + + diff --git a/man/CST_AnalogsPredictors.Rd b/man/CST_AnalogsPredictors.Rd new file mode 100644 index 0000000000000000000000000000000000000000..152b0c8a7060026f06ff1c72b04b0a9b2ae12015 --- /dev/null +++ b/man/CST_AnalogsPredictors.Rd @@ -0,0 +1,151 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_AnalogsPredictors.R +\name{CST_AnalogsPredictors} +\alias{CST_AnalogsPredictors} +\title{AEMET Downscaling +Precipitation and maximum and minimum temperature downscaling method +based on analogs: synoptic situations and significant predictors.} +\usage{ +CST_AnalogsPredictors( + exp, + slp, + obs, + lon, + lat, + slp_lon, + slp_lat, + var_name, + hr_obs, + tdates, + ddates, + restrain, + dim_name_longitude = "lon", + dim_name_latitude = "lat", + dim_name_time = "time" +) +} +\arguments{ +\item{exp}{List of arrays with downscaled period seasonal forecast data. The list +has to contain model atmospheric variables (instantaneous 12h data) that must +be indentify by parenthesis name. +For precipitation: +- u component of wind at 500 hPa (u500_mod) in m/s +- v component of wind at 500 hPa (v500_mod) in m/s +- temperature at 500 hPa (t500_mod) in K +- temperature at 850 hPa (t850_mod) in K +- specific humidity at 700 hPa (q700_mod) in g/kg +For temperature: +- u component of wind at 500 hPa (u500_mod) in m/s +- v component of wind at 500 hPa (v500_mod) in m/s +- temperature at 500 hPa (t500_mod) in K +- temperature at 700 hPa (t700_mod) in K +- temperature at 850 hPa (t850_mod) in K +- specific humidity at 700 hPa (q700_mod) in g/kg +- 2 meters temperature (tm2m_mod) in K +The arrays must have at least three dimensions with names 'lon', 'lat' and 'time'. +(lon = gridpoints of longitude, lat = gridpoints of latitude, time = number of downscaling days) +Seasonal forecast variables must have the same resolution and +domain as reanalysis variables ('obs' parameter, below).} + +\item{slp}{Array with atmospheric seasonal forecast model sea level pressure +(instantaneous 12h data) that must be indentify as 'slp' (hPa). It has the same +resolution as 'exp' and 'obs' paremeters but with an extended domain. +This domain contains extra degrees (most in the north and west part) compare to +synoptic domain. The array must have at least three dimensions +with names 'lon', 'lat' and 'time'.} + +\item{obs}{List of arrays with training period reanalysis data. +The list has to contain reanalysis atmospheric variables (instantaneous +12h data) that must be indentify by parenthesis name. +For precipitation: +- u component of wind at 500 hPa (u500) in m/s +- v component of wind at 500 hPa (v500) in m/s +- temperature at 500 hPa (t500) in K +- temperature at 850 hPa (t850) in K +- sea level pressure (slp) in hPa +- specific humidity at 700 hPa (q700) in g/kg +For maximum and minimum temperature: +- u component of wind at 500 hPa (u500) in m/s +- v component of wind at 500 hPa (v500) in m/s +- temperature at 500 hPa (t500) in K +- temperature at 700 hPa (t700) in K +- temperature at 850 hPa (t850) in K +- sea level pressure (slp) in hPa +- specific humidity at 700 hPa (q700) in g/kg +- 2 meters temperature (tm2m) in K +The arrays must have at least three dimensions with names 'lon', 'lat' and 'time'.} + +\item{lon}{Vector of the synoptic longitude (from (-180º) to 180º), +The vector must go from west to east. The same as for the training function.} + +\item{lat}{Vector of the synoptic latitude. The vector must go from north to south. +The same as for the training function.} + +\item{slp_lon}{Vector of the extended longitude (from (-180º) to 180º), +The vector must go from west to east. The same as for the training function.} + +\item{slp_lat}{Vector of the extended latitude. The vector must go from north to south. +The same as for the training function.} + +\item{var_name}{Variable name to downscale. There are two options: 'prec' for +precipitation and 'temp' for maximum and minimum temperature.} + +\item{hr_obs}{Local path of HR observational files (maestro and pcp/tmx-tmn). +For precipitation can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/v2/Serie_AEMET_v2_pcp_1951a202006_txt.tar.gz +For maximum and minimum temperature can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmax/Serie_AEMET_v1_tmax_1951a202006_txt.tar.gz and http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmin/Serie_AEMET_v1_tmin_1951a202006_txt.tar.gz respetively. +Maestro file (maestro_red_hr_SPAIN.txt) has gridpoint (nptos), longitude (lon), latitude (lat) and +altitude (alt) in columns (vector structure). +Data file (pcp/tmx/tmn_red_SPAIN_1951-201903.txt) includes 5km resolution spanish daily data +(precipitation or maximum and minimum temperature from january 1951 to june 2020. See README +file for more information. +IMPORTANT!: HR observational period must be the same as for reanalysis variables. +It is assumed that the training period is smaller than the HR original one (1951-2019), so it is +needed to make a new ascii file with the new period and the same structure as original, +specifying the training dates in the name (e.g. 'pcp_red_SPAIN_19810101-19961231.txt' for +'19810101-19961231' period).} + +\item{tdates}{Training period dates in format YYYYMMDD(start)-YYYYMMDD(end) +(e.g. 19810101-20181231).} + +\item{ddates}{Downscaling period dates in format YYYYMMDD(start)-YYYYMMDD(end) (e.g. 20191001-20200331).} + +\item{restrain}{Output (list of matrix) obtained from 'training_analogs' function. +For precipitation, 'restrain' object must contains um, vm, nger, gu92, gv92, +gu52, gv52, neni, vdmin, vref, ccm, lab_pred and cor_pred variables. +For maximum and minimum temperature, 'restrain' object must contains um, vm, +insol, neni, vdmin y vref. See 'AnalogsPred_train.R' for more information.} + +\item{dim_name_longitude}{A character string indicating the name of the longitude +dimension, by default 'longitude'.} + +\item{dim_name_latitude}{A character string indicating the name of the latitude +dimension, by default 'latitude'.} + +\item{dim_name_time}{A character string indicating the name of the time +dimension, by default 'time'.} +} +\value{ +Matrix with seasonal forecast precipitation (mm) or +maximum and minimum temperature (dozens of ºC) in a 5km x 5km regular grid +over peninsular Spain and Balearic Islands. The resulted matrices have two +dimensions ('ddates' x 'nptos').(ddates = number of downscaling days +and nptos = number of 'hr_obs' gridpoints). +} +\description{ +This function downscales low resolution precipitation data (e.g. from +Seasonal Forecast Models) through the association with an observational high +resolution (HR) dataset (AEMET 5 km gridded data of daily precipitation (Peral et al., 2017)) +and a collection of predictors and past synoptic situations similar to estimated day. +The method uses three domains: +- peninsular Spain and Balearic Islands domain (5 km resolution): HR precipitation + and the downscaling result domain. +- synoptic domain (low resolution, e.g. 1.5º x 1.5º): it should be centered over Iberian Peninsula + and cover enough extension to detect as much synoptic situations as possible. +- extended domain (low resolution, e.g. 1.5º x 1.5º): it should have the same resolution +as synoptic domain. It is used for SLP Seasonal Forecast Models. +} +\author{ +Marta Dominguez Alonso - AEMET, \email{mdomingueza@aemet.es} + +Nuria Perez-Zanon - BSC, \email{nuria.perez@bsc.es} +} diff --git a/man/training_analogs.Rd b/man/training_analogs.Rd new file mode 100644 index 0000000000000000000000000000000000000000..447f8b0c27ea0f3a2b1c5feef79c75165a3dd1a2 --- /dev/null +++ b/man/training_analogs.Rd @@ -0,0 +1,132 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AnalogsPred_train.R +\name{training_analogs} +\alias{training_analogs} +\title{AEMET Training +Training method (pre-downscaling) based on analogs: +synoptic situations and significant predictors.} +\usage{ +training_analogs( + pred, + slp_ext, + lon, + lat, + slp_lon, + slp_lat, + var, + HR_path, + tdates +) +} +\arguments{ +\item{pred}{List of matrix reanalysis data in a synoptic domain. The list +has to contain reanalysis atmospheric variables (instantaneous 12h data) +that must be indentify by parenthesis name. +For precipitation: +- u component of wind at 500 hPa (u500) in m/s +- v component of wind at 500 hPa (v500) in m/s +- temperature at 500 hPa (t500) in K +- temperature at 850 hPa (t850) in K +- temperature at 1000 hPa (t1000) in K +- geopotential height at 500 hPa (z500) in m +- geopotential height at 1000 hPa (z1000) in m +- sea level pressure (slp) in hPa +- specific humidity at 700 hPa (q700) in g/kg +For maximum and minimum temperature: +- temperature at 1000 hPa (t1000) in K +- sea level pressure (slp) in hPa +All matrix must have [time,gridpoint] dimensions. +(time = number of training days, gridpoint = number of synoptic gridpoints).} + +\item{slp_ext}{Matrix with atmospheric reanalysis sea level pressure +(instantaneous 12h data)(hPa). It has the same resolution as 'pred' parameter +but with an extended domain. This domain contains extra degrees (most in the +north and west part) compare to synoptic domain. The matrix must have +[time,gridpoint] dimensions. +(time = number of training days, gridpoint = number of extended gridpoints).} + +\item{lon}{Vector of the synoptic longitude (from (-180º) to 180º), +The vector must go from west to east.} + +\item{lat}{Vector of the synoptic latitude. The vector must go from north to south.} + +\item{slp_lon}{Vector of the extended longitude (from (-180º) to 180º) +The vector must go from west to east.} + +\item{slp_lat}{Vector of the extended latitude. The vector must go from north to south.} + +\item{var}{Variable name to downscale. There are two options: 'prec' for +precipitation and 'temp' for maximum and minimum temperature.} + +\item{HR_path}{Local path of HR observational files (maestro and pcp/tmx-tmn). +For precipitation can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/v2/Serie_AEMET_v2_pcp_1951a202006_txt.tar.gz +For maximum and minimum temperature can be downloaded from http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmax/Serie_AEMET_v1_tmax_1951a202006_txt.tar.gz and http://www.aemet.es/documentos/es/serviciosclimaticos/cambio_climat/datos_diarios/dato_observacional/rejilla_5km/temperatura/v1/tmin/Serie_AEMET_v1_tmin_1951a202006_txt.tar.gz respetively. +Maestro file (maestro_red_hr_SPAIN.txt) has gridpoint (nptos), longitude (lon), latitude (lat) and +altitude (alt) in columns (vector structure). +Data file (pcp/tmx/tmn_red_SPAIN_1951-201903.txt) includes 5km resolution spanish daily data +(precipitation or maximum and minimum temperature from january 1951 to june 2020. See README +file for more information. +IMPORTANT!: HR observational period must be the same as for reanalysis variables. +It is assumed that the training period is smaller than the HR original one (1951-2020), so it is +needed to make a new ascii file with the new period and the same structure as original, +specifying the training dates ('tdates' parameter) in the name +(e.g. 'pcp_red_SPAIN_19810101-19961231.txt' for '19810101-19961231' period).} + +\item{tdates}{Training period dates in format YYYYMMDD(start)-YYYYMMDD(end) (e.g. 19810101-19961231).} +} +\value{ +matrix list (e.g. restrain) as a result of characterize the past synoptic +situations and the significant predictors needed to downscale seasonal forecast variables. +For precipitation the output includes: +um: u component of geostrophic wind in all period (numeric matrix with [time,gridpoint] dimensions) +vm: v component of geostrophic wind in all period (numeric matrix with [time,gridpoint] dimensions) +nger: number of synoptic situations (integer) +gu92: u component of geostrophic wind for each synoptic situation (numeric matrix with + [nger,gridpoint] dimensions) +gv92: v component of geostrophic wind for each synoptic situation (numeric matrix with + [nger,gridpoint] dimensions) +gu52: u component of wind at 500 hPa for each synotic situation (numeric matrix with + [nger,gridpoint] dimensions) +gv52: v component of wind at 500 hPa for each synotic situation (numeric matrix with + [nger,gridpoint] dimensions) +neni: number of reference centers where predictors are calculated (integer) +vdmin: minimum distances between each HR gridpoint and the four nearest synoptic + gridpoints (numeric matrix with [nptos,4] dimensions) (nptos = number of HR gridpoints) +vref: four nearest synoptic gridpoints to each HR gridpoint (integer matrix with + [nptos,4] dimensions) +ccm: multiple correlation coeficients (numeric matrix with [nger,nptos] dimensions) +indices: + - lab_pred: numeric labels of selected predictors (integer matrix + with [nger,nptos,11,1] dimensions) + - cor_pred: partial correlation of selected predictors (numeric matrix with + [nger,nptos,11,2] dimensions) +For maximum and minimum temperature the output includes: +um: u component of geostrophic wind in all training period (numeric matrix with [time,gridpoint] dimensions) +vm: v component of geostrophic wind in all training period (numeric matrix with [time,gridpoint] dimensions) +insol: insolation in all training period (numeric vector with [time] dimension) +neni: number of reference centers where predictors are calculated (integer) +vdmin: minimum distances between each HR gridpoint and the four nearest synoptic + gridpoints (numeric matrix with [nptos,4] dimensions) (nptos = number of HR gridpoints) +vref: four nearest synoptic gridpoints to each HR gridpoint (integer matrix with + [nptos,4] dimensions) + +The output can directly use as argument to 'CST_AnalogsPredictors' function +(e.g. resdowns <- CST_AnalogsPredictors(...,restrain)) +} +\description{ +This function caracterizes the synoptic situations in a past period based on +low resolution reanalysis data (e.g, ERAInterim 1.5º x 1.5º) and an observational high +resolution (HR) dataset (AEMET 5 km gridded daily precipitation and maximum and +minimum temperature) (Peral et al., 2017)). +The method uses three domains: +- peninsular Spain and Balearic Islands domain (5 km resolution): HR domain +- synoptic domain (low resolution): it should be centered over Iberian Peninsula and + cover enough extension to detect as much synoptic situations as possible. +- extended domain (low resolution): it is an extension of the synoptic + domain. It is used for 'slp_ext' parameter (see 'slp_lon' and 'slp_lat' below). +} +\author{ +Marta Dominguez Alonso - AEMET, \email{mdomingueza@aemet.es} + +Nuria Perez-Zanon - BSC, \email{nuria.perez@bsc.es} +} diff --git a/src/CSTools.so b/src/CSTools.so new file mode 100755 index 0000000000000000000000000000000000000000..013c91ca1300cce79ee5a744faa387225f9fd38f Binary files /dev/null and b/src/CSTools.so differ diff --git a/src/Makevars b/src/Makevars new file mode 100755 index 0000000000000000000000000000000000000000..42ca55ca811d36a314de163b040aa50b3200711c --- /dev/null +++ b/src/Makevars @@ -0,0 +1,27 @@ + +FC = gfortran +MOD_OBJS = mod_csts.o mod_funcs.o +FF_OBJS = calc_geoswind.o calculo_tempes_densi_sealevel.o calc_utmcandelasgrid.o calc_utm_rej_era5_penin.o clasif_era_pen_kmeans.o insol.o predictores_significativos.o pts_ref_est_pen_4int.o pts_ref_est_pen.o +F_OBJS = training_part1_prec.o training_part2_prec.o training_temp.o downscaling_prec.o downscaling_temp.o + +all: + @$(MAKE) $(SHLIB) + @rm -f *.mod *.o + +$(SHLIB): $(MOD_OBJS) $(FF_OBJS) $(F_OBJS) + +calc_geoswind.o: mod_csts.o mod_funcs.o +calculo_tempes_densi_sealevel.o: calculo_tempes_densi_sealevel.f90 mod_csts.o +calc_utmcandelasgrid.o: calc_utmcandelasgrid.f90 mod_csts.o mod_funcs.o +calc_utm_rej_era5_penin.o: calc_utm_rej_era5_penin.f90 mod_csts.o mod_funcs.o +clasif_era_pen_kmeans.o: clasif_era_pen_kmeans.f90 mod_csts.o mod_funcs.o +insol.o: insol.f90 mod_csts.o mod_funcs.o +predictores_significativos.o: predictores_significativos.f90 mod_csts.o mod_funcs.o +pts_ref_est_pen_4int.o: pts_ref_est_pen_4int.f90 mod_csts.o +pts_ref_est_pen.o: pts_ref_est_pen.f90 mod_csts.o +training_part1_prec.o: training_part1_prec.f90 mod_csts.o mod_funcs.o calculo_tempes_densi_sealevel.o calc_geoswind.o clasif_era_pen_kmeans.o +training_part2_prec.o: training_part2_prec.f90 calc_utmcandelasgrid.o calc_utm_rej_era5_penin.o pts_ref_est_pen_4int.o pts_ref_est_pen.o predictores_significativos.o +training_temp.o: training_temp.f90 mod_csts.o mod_funcs.o calculo_tempes_densi_sealevel.o calc_geoswind.o insol.o calc_utm_rej_era5_penin.o calc_utmcandelasgrid.o pts_ref_est_pen_4int.o +downscaling_prec.o: downscaling_prec.f90 mod_csts.o mod_funcs.o +downscaling_temp.o: downscaling_temp.f90 mod_csts.o mod_funcs.o + diff --git a/src/calc_geoswind.f90 b/src/calc_geoswind.f90 new file mode 100755 index 0000000000000000000000000000000000000000..d6b4645476b0f01ada5ea29fb8c171dc6ca24dd7 --- /dev/null +++ b/src/calc_geoswind.f90 @@ -0,0 +1,65 @@ + +! Program to calculate the geostrophic wind based on mean +! sea level pressure (msl) and sea level density + +SUBROUTINE geos(ic,nd,id,slat,slon,slats,slons,rlat,& + rlon,rlats,rlons,nlat,nlon,nlats,nlons,den,msl_lr,& + ngridd,um,vm) + +!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc +!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +USE MOD_CSTS +USE MOD_FUNCS, ONLY : geostrofico,dobla,bessel + + Implicit none + + INTEGER, INTENT(IN) :: nd + INTEGER, INTENT(IN) :: ic + INTEGER, INTENT(IN) :: id + + INTEGER, INTENT(IN) :: nlat + INTEGER, INTENT(IN) :: nlon + INTEGER, INTENT(IN) :: nlats + INTEGER, INTENT(IN) :: nlons + DOUBLE PRECISION, INTENT(IN) :: slat + DOUBLE PRECISION, INTENT(IN) :: slon + DOUBLE PRECISION, INTENT(IN) :: slats + DOUBLE PRECISION, INTENT(IN) :: slons + DOUBLE PRECISION, INTENT(IN) :: rlat + DOUBLE PRECISION, INTENT(IN) :: rlon + DOUBLE PRECISION, INTENT(IN) :: rlats + DOUBLE PRECISION, INTENT(IN) :: rlons + REAL, INTENT(IN) :: den(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: msl_lr(nd,id) + INTEGER, INTENT(IN) :: ngridd + DOUBLE PRECISION, INTENT(OUT) :: um(nd,ic) + DOUBLE PRECISION, INTENT(OUT) :: vm(nd,ic) + + integer i + real psl(id),umint(ic),vmint(ic) + character for10*50 + +! print*,"program 2: geostrophic wind" + + do 1000 i=1,nd + +! Pressure data + psl(:)=msl_lr(i,:)*100./g + +! Calculate the geostrophic wind components (umint,vmint) + call geostrofico(psl,umint,vmint,id,ic,slat,slon,slats,slons,& + rlat,rlon,rlats,rlons,nlat,nlon,nlats,nlons,ngridd) + +! do j=1,ic +! um(i,j)=umint(j)/den(i,j) +! vm(i,j)=vmint(j)/den(i,j) +! enddo + + um(i,:)=umint(:)/den(i,:) + vm(i,:)=vmint(:)/den(i,:) + + 1000 continue + +END SUBROUTINE geos + diff --git a/src/calc_utm_rej_era5_penin.f90 b/src/calc_utm_rej_era5_penin.f90 new file mode 100755 index 0000000000000000000000000000000000000000..75fee505cbef04745464de4a2607f25f1d7b4ef0 --- /dev/null +++ b/src/calc_utm_rej_era5_penin.f90 @@ -0,0 +1,52 @@ + +! The utm_ERA program calculates the Reanalysis UTM coordinates +! with a time zone of 30. + +SUBROUTINE utm_ERA(ic,nlat,nlon,slat,slon,rlat,rlon,x,y) + +USE MOD_CSTS +USE MOD_FUNCS + + Implicit none + + INTEGER, INTENT(IN) :: ic + + INTEGER, INTENT(IN) :: nlat + INTEGER, INTENT(IN) :: nlon + DOUBLE PRECISION, INTENT(IN) :: slat + DOUBLE PRECISION, INTENT(IN) :: slon + DOUBLE PRECISION, INTENT(IN) :: rlat + DOUBLE PRECISION, INTENT(IN) :: rlon + REAL, INTENT(OUT) :: x(ic) + REAL, INTENT(OUT) :: y(ic) + + integer j + integer igrad,imin,rseg,igrad1,imin1 + real rseg1 + real rlt(ic),rln(ic) +! real*8 r1,r2,r3,r4,r5,r6,rad,rad1,xint,yint + double precision r1,r2,r3,r4,r5,r6,rad,rad1,xint,yint + +! print*,"program 4: reanalysis UTM coordinates" + +! Calculation of geostrophic coordinates in each synoptic grid points + + do j=1,ic + rlt(j)=slat+(((j-1)/nlon)*rlat) + rln(j)=slon+((mod(j-1,nlon)+1-1)*rlon) + enddo + +! calculation of UTM coordinates + + do j=1,ic + rad1=rln(j) + rad=rlt(j) + + call geoutm(rad1,rad,huso,xint,yint) + + x(j)=xint + y(j)=yint + + enddo + +END SUBROUTINE utm_ERA diff --git a/src/calc_utmcandelasgrid.f90 b/src/calc_utmcandelasgrid.f90 new file mode 100755 index 0000000000000000000000000000000000000000..1acf8fb05212a64ea3605d42707c0c8510290373 --- /dev/null +++ b/src/calc_utmcandelasgrid.f90 @@ -0,0 +1,41 @@ + +! The utm_obs program calculates the UTM coordinates of high resolution +! (5km x 5km) observational database created by AEMET (Peral et al., 2017). +! +SUBROUTINE utm_obs(lon_hr,lat_hr,xcand,ycand) + +USE MOD_CSTS +USE MOD_FUNCS, ONLY : geoutm + + Implicit none + + DOUBLE PRECISION, INTENT(IN) :: lon_hr(nptos) + DOUBLE PRECISION, INTENT(IN) :: lat_hr(nptos) + REAL, INTENT(OUT) :: xcand(nptos) + REAL, INTENT(OUT) :: ycand(nptos) + + integer n + integer i +! real*8 rad,rad1,xint,yint + double precision rad,rad1,xint,yint + +! print*,"program 5: UTM coordinates high resolution observational database" + + n=nptos + + do i=1,n + rad1=lon_hr(i) + rad=lat_hr(i) +! Calculation of UTM coordinates + call geoutm(rad1,rad,huso,xint,yint) + + xcand(i)=xint + ycand(i)=yint + + enddo + +END SUBROUTINE utm_obs + +!Peral, C., Navascués, B. and Ramos, P.: Serie de precipitación diaria en +!rejilla con fines climáticos. Nota Técnica nº 24, AEMET, +!http://hdl.handle.net/20.500.11765/7573, 2017. diff --git a/src/calculo_tempes_densi_sealevel.f90 b/src/calculo_tempes_densi_sealevel.f90 new file mode 100755 index 0000000000000000000000000000000000000000..362c071442188dc6deef2a4306b9cfdc5ca9cd67 --- /dev/null +++ b/src/calculo_tempes_densi_sealevel.f90 @@ -0,0 +1,37 @@ + +! Program to calculate sea level temperature and density from 1000Mb +! temperature and mean sea level pressure (msl) + +SUBROUTINE calc_tempes_densi_sealev(ic,nd,msl_si,t1000,den) + +USE MOD_CSTS + + IMPLICIT NONE + +!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + INTEGER, INTENT(IN) :: nd + INTEGER, INTENT(IN) :: ic + DOUBLE PRECISION, INTENT(IN) :: msl_si(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: t1000(nd,ic) + REAL, INTENT(OUT) :: den(nd,ic) + + real c,yy + real psl(ic),tmil(ic),tsl(ic) + integer i,j + + c=r*a/g + + do i=1,nd + psl(:)=msl_si(i,:) + tmil(:)=t1000(i,:) + do j=1,ic + yy=log(tmil(j))-c*log(1000./psl(j)) + tsl(j)=exp(yy) +! Air density equation + den(i,j)=(psl(j)*100.)/(r*tsl(j)) + enddo + enddo + +END SUBROUTINE calc_tempes_densi_sealev + diff --git a/src/clasif_era_pen_kmeans.f90 b/src/clasif_era_pen_kmeans.f90 new file mode 100755 index 0000000000000000000000000000000000000000..3db1e7e3d3576099831d5af9adec167e372ae762 --- /dev/null +++ b/src/clasif_era_pen_kmeans.f90 @@ -0,0 +1,273 @@ + +! Program to do synoptic clasification from reanalysis: geostrophic wind +! components, 500Mb wind components and 1000Mb and 500Mb geopotential + +SUBROUTINE clasif(ic,nd,nlon,nlat,slat,slon,rlat,rlon,um,vm,u500,v500,& + z1000,z500,nger,gu92,gv92,gu52,gv52) + +USE MOD_CSTS +USE MOD_FUNCS + + Implicit none + + INTEGER, INTENT(IN) :: ic + INTEGER, INTENT(IN) :: nd + + INTEGER, INTENT(IN) :: nlat + INTEGER, INTENT(IN) :: nlon + DOUBLE PRECISION, INTENT(IN) :: slat + DOUBLE PRECISION, INTENT(IN) :: slon + DOUBLE PRECISION, INTENT(IN) :: rlat + DOUBLE PRECISION, INTENT(IN) :: rlon + DOUBLE PRECISION, INTENT(IN) :: um(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: vm(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: u500(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: v500(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: z1000(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: z500(nd,ic) + INTEGER, INTENT(OUT) :: nger + DOUBLE PRECISION, INTENT(OUT) :: gu92(nd,ic) + DOUBLE PRECISION, INTENT(OUT) :: gv92(nd,ic) + DOUBLE PRECISION, INTENT(OUT) :: gu52(nd,ic) + DOUBLE PRECISION, INTENT(OUT) :: gv52(nd,ic) + + integer n,nm + integer i,j,k,i1,i2,iel,iger,igmin,ipos,iter + real dis,disu5,disv5,disu9,disv9,dmin + real u9(nd,ic),v9(nd,ic),u5(nd,ic),v5(nd,ic) + real gu91(nd,ic),gv91(nd,ic),gu51(nd,ic),gv51(nd,ic) + real ger9(nd,ic),ger5(nd,ic),geo9(nd,ic),geo5(nd,ic) + real ger9_newUnd(nd,ic) + real ser(nd),md,sg,mu9(ic),su9(ic),mv9(ic),sv9(ic) + real mu5(ic),su5(ic),mv5(ic),sv5(ic) + + integer cl(nd,umb_ger+1),ger(nd+1) + real p9(ic),p5(ic),rlt(ic),rln(ic) + + n=nd + nm=umb_ger + +! print*,'program 3: sinoptic clasification' + +! Calculation to assign the weights to each grid point. + + do j=1,ic + rlt(j)=slat+(((j-1)/nlon)*rlat) + rln(j)=slon+((mod(j-1,nlon)+1-1)*rlon) + enddo + + p9=0. + p5=1. + + do i1=1,ic + if((rlt(i1).le.fnor2).and.(rlt(i1).ge.fsur2)) then + if((rln(i1).ge.foes2).and.(rln(i1).le.fest2)) then + p9(i1)=1. + p5(i1)=4. + endif + endif + enddo + do i1=1,ic + if((rlt(i1).le.fnor1).and.(rlt(i1).ge.fsur1)) then + if((rln(i1).ge.foes1).and.(rln(i1).le.fest1)) then + p9(i1)=2. + p5(i1)=8. + endif + endif + enddo +! +! REANALYSIS VARIABLES + + u5(:,:)=u500(:,:) + v5(:,:)=v500(:,:) + geo9(:,:)=z1000(:,:) + geo5(:,:)=z500(:,:) + +! Mean and standard deviation of reference synoptic fields + + do j=1,ic + do i=1,n + ser(i)=um(i,j) + enddo + call estadis(ser,md,sg,n) + mu9(j)=md + su9(j)=sg + do i=1,n + ser(i)=vm(i,j) + enddo + call estadis(ser,md,sg,n) + mv9(j)=md + sv9(j)=sg + do i=1,n + ser(i)=u5(i,j) + enddo + call estadis(ser,md,sg,n) + mu5(j)=md + su5(j)=sg + do i=1,n + ser(i)=v5(i,j) + enddo + call estadis(ser,md,sg,n) + mv5(j)=md + sv5(j)=sg + enddo + + +! Geostrophic wind components standatization + + do i=1,n + do j=1,ic + u9(i,j)=(um(i,j)-mu9(j))/su9(j) + v9(i,j)=(vm(i,j)-mv9(j))/sv9(j) + u5(i,j)=(u5(i,j)-mu5(j))/su5(j) + v5(i,j)=(v5(i,j)-mv5(j))/sv5(j) + enddo + enddo + +! Finding the cluster centers + + ger(n+1)=1 + ger(1)=1 + do 200 i=2,n + do 210 j=1,ger(n+1) + iger=ger(j) + call distan9(u9,n,ic,i,iger,p9,disu9) + call distan9(v9,n,ic,i,iger,p9,disv9) + call distan5(u5,n,ic,i,iger,p5,disu5) + call distan5(v5,n,ic,i,iger,p5,disv5) + dis=(disu9+disv9+disu5+disv5)/4. + if(dis.lt.umb) go to 200 + 210 continue + ger(n+1)=ger(n+1)+1 + ipos=ger(n+1) + ger(ipos)=i + 200 continue + + do k=1,ger(n+1) + enddo + +! K-means method: weather types + + nger=ger(n+1) + +! print*,' number of synoptic types = ', nger + + do k=1,nger + iger=ger(k) + do j=1,ic + gu92(k,j)=u9(iger,j) + gv92(k,j)=v9(iger,j) + gu52(k,j)=u5(iger,j) + gv52(k,j)=v5(iger,j) + gu91(k,j)=u9(iger,j) + gv91(k,j)=v9(iger,j) + gu51(k,j)=u5(iger,j) + gv51(k,j)=v5(iger,j) + enddo + enddo + + iter=0 + 251 continue + cl=0 + iter=iter+1 + do 300 i1=1,n + dmin=1000. + igmin=0 + do 310 i2=1,nger + call distancia9(u9,n,gu92,n,i1,i2,p9,disu9,ic) + call distancia9(v9,n,gv92,n,i1,i2,p9,disv9,ic) + call distancia5(u5,n,gu52,n,i1,i2,p5,disu5,ic) + call distancia5(v5,n,gv52,n,i1,i2,p5,disv5,ic) + dis=(disu9+disv9+disu5+disv5)/4. + if(dis.lt.dmin) then + dmin=dis + igmin=i2 + endif + 310 continue + cl(igmin,nm+1)=cl(igmin,nm+1)+1 + ipos=cl(igmin,nm+1) + cl(igmin,ipos)=i1 + 300 continue + + ger9=0. + ger5=0. + gu92=0. + gv92=0. + gu52=0. + gv52=0. + + do i=1,nger + do j=1,cl(i,nm+1) + iel=cl(i,j) + do k=1,ic + ger9(i,k)=ger9(i,k)+geo9(iel,k) + ger5(i,k)=ger5(i,k)+geo5(iel,k) + gu92(i,k)=gu92(i,k)+u9(iel,k) + gv92(i,k)=gv92(i,k)+v9(iel,k) + gu52(i,k)=gu52(i,k)+u5(iel,k) + gv52(i,k)=gv52(i,k)+v5(iel,k) + enddo + enddo + do k=1,ic + gu92(i,k)=gu92(i,k)/real(cl(i,nm+1)) + gv92(i,k)=gv92(i,k)/real(cl(i,nm+1)) + gu52(i,k)=gu52(i,k)/real(cl(i,nm+1)) + gv52(i,k)=gv52(i,k)/real(cl(i,nm+1)) + enddo + + enddo + + do i=1,nger + call distancia9(gu91,n,gu92,n,i,i,p9,disu9,ic) + call distancia9(gv91,n,gv92,n,i,i,p9,disv9,ic) + call distancia5(gu51,n,gu52,n,i,i,p5,disu5,ic) + call distancia5(gv51,n,gv52,n,i,i,p5,disv5,ic) + dis=(disu9+disv9+disu5+disv5)/4. + if(dis.ge.0.10) go to 250 + enddo + + go to 252 + + 250 continue + + do i=1,nger + do j=1,ic + gu91(i,j)=gu92(i,j) + gv91(i,j)=gv92(i,j) + gu51(i,j)=gu52(i,j) + gv51(i,j)=gv52(i,j) + enddo + enddo + go to 251 + + 252 continue + +!cccccccccccccccccccccccccccccccccccccccccc + + do i=1,nger + do j=1,ic + gu92(i,j)=(gu92(i,j)*su9(j))+mu9(j) + gv92(i,j)=(gv92(i,j)*sv9(j))+mv9(j) + gu52(i,j)=(gu52(i,j)*su5(j))+mu5(j) + gv52(i,j)=(gv52(i,j)*sv5(j))+mv5(j) + enddo + enddo + +! These variables are not going to be used but they should not be delated + do 401 i=1,nger + do k=1,ic + ger9(i,k)=ger9(i,k)/real(cl(i,nm+1)) + ger9_newUnd(i,k)=1000.+(ger9(i,k)/8.) + ger5(i,k)=ger5(i,k)/real(cl(i,nm+1)) + enddo + + 401 continue + +END SUBROUTINE clasif + + + + + + + diff --git a/src/downscaling_prec.f90 b/src/downscaling_prec.f90 new file mode 100755 index 0000000000000000000000000000000000000000..c061488ee806e81b528589e7044d03b01c545de7 --- /dev/null +++ b/src/downscaling_prec.f90 @@ -0,0 +1,943 @@ +! Program to downscale precipitation based on analogs method +! for Iberian Peninsula and Balearic Islands (Autor: Petisco de Lara) +! ****************************************************** + +SUBROUTINE down_prec(ic,id,nd,nm,nlat,nlon,nlatt,nlont,slat,slon,rlat,rlon,& + slatt,slont,ngridd,u500,v500,t500,t850,msl_si,q700,& + prec_hr,nanx,um,vm,nger,gu92,gv92,gu52,gv52,& + neni,vdmin,vref4,new_ccm,new_kvars,new_corrpar,u500e,& + v500e,t500e,t850e,msle,q700e,pp) + + +USE MOD_CSTS +USE MOD_FUNCS + + IMPLICIT NONE + +!cccccccccccccccccccccccccccccccccc +!cccccccccccccccccccccccccccccccccc + + INTEGER, INTENT(IN) :: ic + INTEGER, INTENT(IN) :: id + INTEGER, INTENT(IN) :: nd + INTEGER, INTENT(IN) :: nm + + INTEGER, INTENT(IN) :: nlat + INTEGER, INTENT(IN) :: nlon + INTEGER, INTENT(IN) :: nlatt + INTEGER, INTENT(IN) :: nlont + DOUBLE PRECISION, INTENT(IN) :: slat + DOUBLE PRECISION, INTENT(IN) :: slon + DOUBLE PRECISION, INTENT(IN) :: rlat + DOUBLE PRECISION, INTENT(IN) :: rlon + DOUBLE PRECISION, INTENT(IN) :: slatt + DOUBLE PRECISION, INTENT(IN) :: slont + INTEGER, INTENT(IN) :: ngridd + + DOUBLE PRECISION, INTENT(IN) :: u500(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: v500(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: msl_si(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: q700(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: t500(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: t850(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: prec_hr(nd,nptos) + + INTEGER, INTENT(IN) :: nanx + + DOUBLE PRECISION, INTENT(IN) :: um(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: vm(nd,ic) + INTEGER, INTENT(IN) :: nger + DOUBLE PRECISION, INTENT(IN) :: gu92(nger,ic) + DOUBLE PRECISION, INTENT(IN) :: gv92(nger,ic) + DOUBLE PRECISION, INTENT(IN) :: gu52(nger,ic) + DOUBLE PRECISION, INTENT(IN) :: gv52(nger,ic) + + INTEGER, INTENT(IN) :: neni + DOUBLE PRECISION, INTENT(IN) :: vdmin(nptos,4) + INTEGER, INTENT(IN) :: vref4(nptos,4) + + DOUBLE PRECISION, INTENT(IN) :: new_ccm(nger,nptos) + INTEGER, INTENT(IN) :: new_kvars(nger,nptos,npx) + DOUBLE PRECISION, INTENT(IN) :: new_corrpar(nger,nptos,npx) + + DOUBLE PRECISION, INTENT(IN) :: u500e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: v500e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: t500e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: t850e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: msle(nm,id) + DOUBLE PRECISION, INTENT(IN) :: q700e(nm,ic) + + DOUBLE PRECISION, INTENT(OUT) :: pp(nm,nptos) + +!cccccccccccccccccccccccccccccccccc + integer m,n + integer i,j + integer nvar + integer ii + integer jt,jp + real rlx,rly + real sp,dd +!cccccccccccccccccccccccccccccccccc +!cccccccccccccccccccccccccccccccccc + integer i1,i2,i3,i4,i7,iana,ice,ien,ipos,ips,ipu,ir,ire,iti,jk,k + integer mesa,nan,nan2,nanf,nanv,nen,nmm,np + real disu5,disu9,disv5,disv9,dmin,dt,du5,du9,dv5,dv9,supo + real vorm,vorz + +!***************************************************************** + + integer anai(nanx),ana(nanx),puh(ic) + real u9(nd,ic),v9(nd,ic),u5(nd,ic),v5(nd,ic),he7(nd,ic),he7m(ic) + real psl(nd,ic),ut9(nger,ic),vt9(nger,ic),ut5(nger,ic),vt5(nger,ic) + real t5(nd,ic),t8(nd,ic),tm5(ic),tm8(ic) + real pslm(ic),um9(ic),vm9(ic),um5(ic),vm5(ic),pslma(ic) + real bdlon,bilat + real pres(id),bar(id) + real pred1(npx,nd,neni),pred1m(npx,neni) + character sc*8,pt*9,nomeb*90 + + integer nor(nanx),prs(nger,nptos,npx+7) + integer annoa + integer ior(nd),eqc(nptos) + integer ref(nptos),puce(neni),puen(neni,5001) + + integer est(nptos) + real prec(nd,nptos) + real p9(ic),p5(ic) + real dista(nd),prees(nm,nptos) + real rlt(ic),rln(ic),rltt(id),rlnt(id) + real dist(nanx),dist1(npx,nanx),serin(nanx) + real aaa(nanx) + real ser(nd),media(npx,neni),sigma(npx,neni) + real md,sg + real mu9(ic),su9(ic),mv9(ic),sv9(ic),copar(nger,nptos,npx) + real mu5(ic),su5(ic),mv5(ic),sv5(ic),ccm(nger,nptos) + +!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +!c Variables nuevas por la nueva interpolacion de los predictores + + integer Vref(nptos,4) + real Vdis(nptos,4) + integer iii + real distancion1, distancion2, distancion3, distancion4 + real peso1, peso2, peso3, peso4 + real calculin_modelo, calculin_calibracion + integer ien1, ien2, ien3, ien4 + integer ik + +!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + m=nger + n=nd + +!******************************** +! 1. Synoptic latitude and longitude calculation and assignment of +! weights windows +! + do j=1,ic + rlt(j)=slat+(((j-1)/nlon)*rlat) + rln(j)=slon+((mod(j-1,nlon)+1-1)*rlon) + enddo + p9=0. + p5=1. + do i1=1,ic + if((rlt(i1).le.fnor2).and.(rlt(i1).ge.fsur2)) then + if((rln(i1).ge.foes2).and.(rln(i1).le.fest2)) then + p9(i1)=1. + p5(i1)=4. + endif + endif + enddo + do i1=1,ic + if((rlt(i1).le.fnor1).and.(rlt(i1).ge.fsur1)) then + if((rln(i1).ge.foes1).and.(rln(i1).le.fest1)) then + p9(i1)=2. + p5(i1)=8. + endif + endif + enddo +! +! Latitude and longitude calculation in the extended domain (called low +! resolution) + + do j=1,id + rltt(j)=slatt+(((j-1)/nlont)*rlat) + rlnt(j)=slont+((mod(j-1,nlont)+1-1)*rlon) + enddo +! +!**************************************************************** +! TRAINING REANALYSIS VARIABLES + + u5(:,:)=u500(:,:) + v5(:,:)=v500(:,:) + psl(:,:)=msl_si(:,:) + he7(:,:)=q700(:,:) + t5(:,:)=t500(:,:) + t8(:,:)=t850(:,:) + +! HIGH RESOLUTION (5KM) OBSERVATIONS +! It is neccesary to convert to tenths of mm (multiplying by 10). + + prec(:,:)=prec_hr(:,:)*10. + +! Mean and standard deviation of reference synoptic fields (wind components). + + do j=1,ic + do i=1,n + ser(i)=um(i,j) + enddo + md=0. + call estadis(ser,md,sg,n) + mu9(j)=md + su9(j)=sg + do i=1,n + ser(i)=vm(i,j) + enddo + call estadis(ser,md,sg,n) + mv9(j)=md + sv9(j)=sg + do i=1,n + ser(i)=u5(i,j) + enddo + call estadis(ser,md,sg,n) + mu5(j)=md + su5(j)=sg + do i=1,n + ser(i)=v5(i,j) + enddo + call estadis(ser,md,sg,n) + mv5(j)=md + sv5(j)=sg + + enddo + +! A reference centers (matching points between synoptic and high +! resolution grids) are define to know where the predictor must be +! calculated. + + Vref(:,:)=vref4(:,:) + Vdis(:,:)=vdmin(:,:) + + nen=1 + puce(1)=Vref(1,1) + + do iii=1,4 + do j=1,nptos + do k=1,nen + if (Vref(j,iii).eq.puce(k)) go to 101 + enddo + nen=nen+1 + ipos=nen + puce(ipos)=Vref(j,iii) + 101 continue + enddo + enddo + +! Each reference point have associated to a group of high resolution grids. + puen=0 + + do k=1,nen + do j=1,nptos + do iii=1,4 + if(Vref(j,iii).eq.puce(k)) then + puen(k,5001)=puen(k,5001)+1 + ipos=puen(k,5001) + puen(k,ipos)=j + endif + enddo + enddo + enddo + +! The predictors are obtained and normalized +! OBTAINING THE SEA LEVEL PRESSURE (PREDICTOR 1) IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(1,i,j)=psl(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(1,i,j) + enddo + call estadis(ser,md,sg,n) + media(1,j)=md + sigma(1,j)=sg + do i=1,n + pred1(1,i,j)=(pred1(1,i,j)-media(1,j))/sigma(1,j) + enddo + enddo + +! OBTAINING THE TREND (PREDICTOR 11) IN THE REFERENCE CENTERS + + do j=1,nen + pred1(11,1,j)=0. + enddo + + do i=2,n + do j=1,nen + ice=puce(j) + pred1(11,i,j)=psl(i,ice)-psl((i-1),ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(11,i,j) + enddo + call estadis(ser,md,sg,n) + media(11,j)=md + sigma(11,j)=sg + do i=1,n + pred1(11,i,j)=(pred1(11,i,j)-media(11,j))/sigma(11,j) + enddo + enddo + +! OBTAINING THE VERTICAL THERMAL GRADIENT(PREDICTOR 3) +! IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(3,i,j)=t8(i,ice)-t5(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(3,i,j) + enddo + call estadis(ser,md,sg,n) + media(3,j)=md + sigma(3,j)=sg + do i=1,n + pred1(3,i,j)=(pred1(3,i,j)-media(3,j))/sigma(3,j) + enddo + enddo + +! OBTAINING THE 500 hPa TEMPERATURE (PREDICTOR 2) +! IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(2,i,j)=t5(i,ice) + enddo + enddo + + do j=1,nen + do i=1,n + ser(i)=pred1(2,i,j) + enddo + call estadis(ser,md,sg,n) + media(2,j)=md + sigma(2,j)=sg + do i=1,n + pred1(2,i,j)=(pred1(2,i,j)-media(2,j))/sigma(2,j) + enddo + enddo + +! OBTAINING THE VORTICITY (PREDICTOR 4) IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + rlx=rt*cos(rlt(ice)*pi/180.)*pi*rlon/180. + rly=rt*abs(rlat)*pi/180. + vorm=um(i,ice-nlon)-um(i,ice+nlon) + vorm=vorm/(2.*rly) + vorz=vm(i,ice+1)-vm(i,ice-1) + vorz=vorz/(2.*rlx) + pred1(4,i,j)=vorz-vorm + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(4,i,j) + enddo + call estadis(ser,md,sg,n) + media(4,j)=md + sigma(4,j)=sg + do i=1,n + pred1(4,i,j)=(pred1(4,i,j)-media(4,j))/sigma(4,j) + enddo + enddo + +! OBTAINING THE GEOSTROPHIC U/V COMPONENTS (PREDICTORS 5 AND 6) IN THE REFERENCE +! CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(5,i,j)=um(i,ice) + pred1(6,i,j)=vm(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(5,i,j) + enddo + call estadis(ser,md,sg,n) + media(5,j)=md + sigma(5,j)=sg + do i=1,n + pred1(5,i,j)=(pred1(5,i,j)-media(5,j))/sigma(5,j) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(6,i,j) + enddo + call estadis(ser,md,sg,n) + media(6,j)=md + sigma(6,j)=sg + do i=1,n + pred1(6,i,j)=(pred1(6,i,j)-media(6,j))/sigma(6,j) + enddo + enddo + +! OBTAINING THE VORTICITY IN 500 hPa (PREDICTOR 7) IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + rlx=rt*cos(rlt(ice)*pi/180.)*pi*rlon/180. + rly=rt*abs(rlat)*pi/180. + vorm=u5(i,ice-nlon)-u5(i,ice+nlon) + vorm=vorm/(2.*rly) + vorz=v5(i,ice+1)-v5(i,ice-1) + vorz=vorz/(2.*rlx) + pred1(7,i,j)=vorz-vorm + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(7,i,j) + enddo + call estadis(ser,md,sg,n) + media(7,j)=md + sigma(7,j)=sg + do i=1,n + pred1(7,i,j)=(pred1(7,i,j)-media(7,j))/sigma(7,j) + enddo + enddo + +! OBTAINING THE GEOSTROPHIC U/V COMPONENTS IN 500 hPa (PREDICTORS 8 AND 9) +! IN THE RERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(8,i,j)=u5(i,ice) + pred1(9,i,j)=v5(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(8,i,j) + enddo + call estadis(ser,md,sg,n) + media(8,j)=md + sigma(8,j)=sg + do i=1,n + pred1(8,i,j)=(pred1(8,i,j)-media(8,j))/sigma(8,j) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(9,i,j) + enddo + call estadis(ser,md,sg,n) + media(9,j)=md + sigma(9,j)=sg + do i=1,n + pred1(9,i,j)=(pred1(9,i,j)-media(9,j))/sigma(9,j) + enddo + enddo + +! OBTAINING THE ESPECIFIC HUMIDITY IN 700 hPa (PREDICTOR 10) IN THE REFERENCE +! CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(10,i,j)=he7(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(10,i,j) + enddo + call estadis(ser,md,sg,n) + media(10,j)=md + sigma(10,j)=sg + do i=1,n + pred1(10,i,j)=(pred1(10,i,j)-media(10,j))/sigma(10,j) + enddo + enddo + +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! ESTANDARIZATION OF REFERENCE WINDS (SYNOPTIC WINDS ALSO) + + do i=1,n + do j=1,ic + u9(i,j)=(um(i,j)-mu9(j))/su9(j) + v9(i,j)=(vm(i,j)-mv9(j))/sv9(j) + u5(i,j)=(u5(i,j)-mu5(j))/su5(j) + v5(i,j)=(v5(i,j)-mv5(j))/sv5(j) + enddo + enddo + + do i=1,m + do j=1,ic + ut9(i,j)=(gu92(i,j)-mu9(j))/su9(j) + vt9(i,j)=(gv92(i,j)-mv9(j))/sv9(j) + ut5(i,j)=(gu52(i,j)-mu5(j))/su5(j) + vt5(i,j)=(gv52(i,j)-mv5(j))/sv5(j) + enddo + enddo + +! SIGNIFICANT PREDICTORS FOR EACH SYNOPTIC TYPE IN EACH HIGH +! RESOLUTION GRID POINT. + + ccm(:,:)=new_ccm(:,:) + prs(:,:,:)=new_kvars(:,:,:) + copar(:,:,:)=new_corrpar(:,:,:) + +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!************************************************************** +!************************************************************** +! +! DOWNSCALING BEGINS (ESTIMATING THE PROBLEM DAYS PRECIPITATION +! IN EACH HIGH RESOLUTION GRID POINT) +! +!************************************************************** +!************************************************************** +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +! print*,'Downscaling begins...' +! print*,'estimated day... ' + + nmm=0 + do 1000 i=1,nm ! Estimated days loop + +! ESTIMATED REANALYSIS VARIABLES + + um5(:)=u500e(i,:) + vm5(:)=v500e(i,:) + pres(:)=msle(i,:) + he7m(:)=q700e(i,:) + tm5(:)=t500e(i,:) + tm8(:)=t850e(i,:) + +! print*,' ',i + nmm=nmm+1 + +! Pressure synoptic grid calculation from low resolution one + + bilat=slat+(nlat-1)*rlat + bdlon=slon+(nlon-1)*rlon + + ire=0 + do 111 j=1,id + if((rltt(j).gt.slat).or.(rltt(j).lt.bilat)) go to 111 + if((rlnt(j).lt.slon).or.(rlnt(j).gt.bdlon)) go to 111 + ire=ire+1 + pslm(ire)=pres(j) + 111 continue + +! Geostrophic wind components at sea level to the estimated +! day (pressure in Pa). +! "/g" to invalidate gravity acceleration. + + bar=pres*100./g + + call geostrofico(bar,um9,vm9,id,ic,slatt,slont,slat,slon,& + rlat,rlon,rlat,rlon,nlatt,nlont,nlat,nlon,ngridd) + +! It is divided by density at sea level (standard atmosphere) to obtain +! the geostrophic wind components. + + do j=1,ic + um9(j)=um9(j)/1.225 + vm9(j)=vm9(j)/1.225 + enddo + +! The estimated predictors are obtained and normalized + +! OBTAINING THE 500 hPa TEMPERATURE (PREDICTOR 2) +! IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(2,j)=tm5(ice) + pred1m(2,j)=(pred1m(2,j)-media(2,j))/sigma(2,j) + enddo + +! OBTAINING THE SEA LEVEL PRESSURE (PREDICTOR 1) IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(1,j)=pslm(ice) + pred1m(1,j)=(pred1m(1,j)-media(1,j))/sigma(1,j) + enddo + +! OBTAINING THE TREND (PREDICTOR 11) IN THE REFERENCE CENTERS + + if(i.eq.1) then + do j=1,nen + pred1m(11,j)=0. + enddo + else + do j=1,nen + ice=puce(j) + pred1m(11,j)=pslm(ice)-pslma(ice) + pred1m(11,j)=(pred1m(11,j)-media(11,j))/sigma(11,j) + enddo + endif + + pslma=pslm + +! OBTAINING THE VERTICAL THERMAL GRADIENT(PREDICTOR 3) +! IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(3,j)=tm8(ice)-tm5(ice) + pred1m(3,j)=(pred1m(3,j)-media(3,j))/sigma(3,j) + enddo + +! OBTAINING THE VORTICITY (PREDICTOR 4) IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + rlx=rt*cos(rlt(ice)*pi/180.)*pi*rlon/180. + rly=rt*abs(rlat)*pi/180. + vorm=um9(ice-nlon)-um9(ice+nlon) + vorm=vorm/(2.*rly) + vorz=vm9(ice+1)-vm9(ice-1) + vorz=vorz/(2.*rlx) + pred1m(4,j)=vorz-vorm + pred1m(4,j)=(pred1m(4,j)-media(4,j))/sigma(4,j) + enddo + +! OBTAINING THE GEOSTROPHIC U/V COMPONENTS (PREDICTORS 5 AND 6) IN THE REFERENCE +! CENTERS + + do j=1,nen + ice=puce(j) + pred1m(5,j)=um9(ice) + pred1m(5,j)=(pred1m(5,j)-media(5,j))/sigma(5,j) + pred1m(6,j)=vm9(ice) + pred1m(6,j)=(pred1m(6,j)-media(6,j))/sigma(6,j) + enddo + +! OBTAINING THE VORTICITY IN 500 hPa (PREDICTOR 7) IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + rlx=rt*cos(rlt(ice)*pi/180.)*pi*rlon/180. + rly=rt*abs(rlat)*pi/180. + vorm=um5(ice-nlon)-um5(ice+nlon) + vorm=vorm/(2.*rly) + vorz=vm5(ice+1)-vm5(ice-1) + vorz=vorz/(2.*rlx) + pred1m(7,j)=vorz-vorm + pred1m(7,j)=(pred1m(7,j)-media(7,j))/sigma(7,j) + enddo + +! OBTAINING THE GEOSTROPHIC U/V COMPONENTS IN 500 hPa (PREDICTORS 8 AND 9) +! IN THE RERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(8,j)=um5(ice) + pred1m(8,j)=(pred1m(8,j)-media(8,j))/sigma(8,j) + pred1m(9,j)=vm5(ice) + pred1m(9,j)=(pred1m(9,j)-media(9,j))/sigma(9,j) + enddo + +! OBTAINING THE ESPECIFIC HUMIDITY IN 700 hPa (PREDICTOR 10) IN THE REFERENCE +! CENTERS + + do j=1,nen + ice=puce(j) + pred1m(10,j)=he7m(ice) + pred1m(10,j)=(pred1m(10,j)-media(10,j))/sigma(10,j) + enddo + +! ESTANDARIZATION OF REFERENCE WINDS + + do j=1,ic + um9(j)=(um9(j)-mu9(j))/su9(j) + vm9(j)=(vm9(j)-mv9(j))/sv9(j) + um5(j)=(um5(j)-mu5(j))/su5(j) + vm5(j)=(vm5(j)-mv5(j))/sv5(j) + enddo + +! Synoptic type determination to which the estimated day belongs. + + dmin=99999. + iti=0 + do j=1,m + call distan9_2(um9,ut9,m,j,p9,du9,ic) + call distan9_2(vm9,vt9,m,j,p9,dv9,ic) + call distan5_2(um5,ut5,m,j,p5,du5,ic) + call distan5_2(vm5,vt5,m,j,p5,dv5,ic) + dd=(du9+dv9+du5+dv5)/4. + if(dd.lt.dmin) then + dmin=dd + iti=j + endif + enddo + +! Determine the "nanx" reference alements more similar to each synoptic type +! and the corresponding distances + + do i1=1,n + ior(i1)=i1 + dista(i1)=9999. + enddo + do 113 i1=1,n + call distan9_2(um9,u9,n,i1,p9,disu9,ic) + call distan9_2(vm9,v9,n,i1,p9,disv9,ic) + call distan5_2(um5,u5,n,i1,p5,disu5,ic) + call distan5_2(vm5,v5,n,i1,p5,disv5,ic) + dista(i1)=(disu9+disv9+disu5+disv5)/4. + 113 continue + call burbuja1(dista,ior,n,nanx) + do i1=1,nanx + anai(i1)=ior(i1) + enddo + +!******************************************************************* +!******************************************************************* + + do 1100 ien=1,nen + do 1200 i2=1,puen(ien,5001) + ipu=puen(ien,i2) + +!**************** +! An analog (nanf) have synoptic similarity regarding estimated day +! when it has value in a point and presents lower distance than a +! given threshold. + + nan=0 + nanf=0 + + do i3=1,nanx + iana=anai(i3) + if(prec(iana,ipu).ne.-999.) then + nan=nan+1 + ana(nan)=iana + dist(nan)=dista(i3) + if(dist(nan).eq.0.0) dist(nan)=0.1 + if(dist(nan).le.umb) nanf=nanf+1 + endif + enddo + + if(nan.lt.nmin) then +! print*,i,ipu,' there are not enough analogs ' + goto 1200 + endif + + if(nanf.le.nmin) nanf=nmin + +! Significant predictors for the synoptic type of the estimated day +! in each HR grid point. + + np=prs(iti,ipu,npx+7) + if(ccm(iti,ipu).lt.ccmi) np=0 + +! If no significant predictors, just a synoptic similarity is taken account. + + if(np.eq.0) then + if(nanf.gt.nmax) then + nan2=nmax + else + nan2=nanf + endif + +!!!!!!!!!!!!!!!!!!!!!!! + + prees(i,ipu)=0. + sp=0. + do i3=1,nan2 + iana=ana(i3) + dt=dist(i3) + sp=sp+1./dt + prees(i,ipu)=prees(i,ipu)+prec(iana,ipu)*(1./dt) + enddo + prees(i,ipu)=prees(i,ipu)/sp + go to 1200 + endif + +! If there are significant predictors: + + do ik=1,nen + ice=puce(ik) + + if (Vref(ipu,1).eq.ice) then + ien1=ik + distancion1=Vdis(ipu,1) + peso1=1/distancion1 + go to 251 + endif + enddo + 251 continue + + + do ik=1,nen + ice=puce(ik) + if (Vref(ipu,2).eq.ice) then + ien2=ik + distancion2=Vdis(ipu,2) + peso2=1/distancion2 + go to 252 + endif + enddo + 252 continue + + do ik=1,nen + ice=puce(ik) + if (Vref(ipu,3).eq.ice) then + ien3=ik + distancion3=Vdis(ipu,3) + peso3=1/distancion3 + go to 253 + endif + enddo + 253 continue + + + do ik=1,nen + ice=puce(ik) + if (Vref(ipu,4).eq.ice) then + ien4=ik + distancion4=Vdis(ipu,4) + peso4=1/distancion4 + go to 254 + endif + enddo + 254 continue + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + do 1250 i4=1,np + ips=prs(iti,ipu,i4) + do i1=1,nanf + iana=ana(i1) + calculin_modelo = pred1m(ips,ien1)*peso1+pred1m(ips,ien2)* & + peso2+pred1m(ips,ien3)*peso3+pred1m(ips,ien4)*peso4/(peso1+peso2+ & + peso3+peso4) + calculin_calibracion = pred1(ips,iana,ien1)*peso1+ & + pred1(ips,iana,ien2)*peso2+pred1(ips,iana,ien3)*peso3+ & + pred1(ips,iana,ien4)*peso4/(peso1+peso2+peso3+peso4) + + dist1(ips,i1)=(calculin_modelo - calculin_calibracion)**2 + + enddo + 1250 continue + +! The "nanf" analogs are sorted from higher to lower similarity (taken account +! both synoptic similarity and significant predictors) + + do ii=1,nanf + aaa(ii)=0. + nor(ii)=ii + enddo + nanv=0 + do ii=1,nanf + supo=0. + do i7=1,np + ips=prs(iti,ipu,i7) + aaa(ii)= aaa(ii)+dist1(ips,nor(ii))*copar(iti,ipu,i7) + supo=supo+copar(iti,ipu,i7) + enddo + aaa(ii)=aaa(ii)/supo + if(aaa(ii).eq.0.) aaa(ii)=0.1 + if(aaa(ii).le.umbl) nanv=nanv+1 + serin(ii)=(aaa(ii)+dist(nor(ii)))/2. + enddo + + call burbuja(serin,nor,nanf,nanf,nanf) + +!!!!!!!!!!!!!!!! + + if(nanv.lt.nmin) go to 1998 + prees(i,ipu)=0. + sp=0. + nan2=0 + do 8888 ii=1,nanf + if(aaa(nor(ii)).gt.umbl) go to 8888 + iana=ana(nor(ii)) + dt=serin(ii) + sp=sp+(1./dt) + prees(i,ipu)=prees(i,ipu)+prec(iana,ipu)*(1./dt) + nan2=nan2+1 + if(nan2.eq.nmax) go to 1995 + 8888 continue + 1995 continue + prees(i,ipu)=prees(i,ipu)/sp + + go to 1200 + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + 1998 continue + prees(i,ipu)=0. + sp=0. + nan2=0 +! 1) If there is some local similarity analog, it is taken account + if(nanv.gt.0) then + do 8889 ii=1,nanf + if(aaa(nor(ii)).gt.umbl) go to 8889 + nan2=nan2+1 + iana=ana(nor(ii)) + dt=serin(ii) + sp=sp+(1./dt) + prees(i,ipu)=prees(i,ipu)+prec(iana,ipu)*(1./dt) + 8889 continue +! and it is completed with the rest of analogs in order of +! total similarity until the minimum number of analogs are completed. + do 8890 ii=1,nanf + if(aaa(nor(ii)).le.umbl) go to 8890 + nan2=nan2+1 + iana=ana(nor(ii)) + dt=serin(ii) + sp=sp+(1./dt) + prees(i,ipu)=prees(i,ipu)+prec(iana,ipu)*(1./dt) + if(nan2.eq.nmin) go to 1997 + 8890 continue + 1997 continue + prees(i,ipu)=prees(i,ipu)/sp + go to 1200 +! 2)If no analogs with local similarity, analogs with total similarity +! are taken (synoptic+local predictors) until the minimum number of +! analogs are completed. + elseif(nanv.eq.0) then + do ii=1,nmin + nan2=nan2+1 + iana=ana(nor(ii)) + dt=serin(ii) + sp=sp+(1./dt) + prees(i,ipu)=prees(i,ipu)+prec(iana,ipu)*(1./dt) + enddo + prees(i,ipu)=prees(i,ipu)/sp + go to 1200 + endif + + + 1200 continue + 1100 continue + +!!!!!!!!!!!!!!!!!!!!!!!!!!!! + 1000 continue !End of estimated days loop +!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + pp(:,:)=prees(:,:) + +!++++++++++++++++++++++++++++++++++++++++++++++++++ + +END SUBROUTINE down_prec + diff --git a/src/downscaling_temp.f90 b/src/downscaling_temp.f90 new file mode 100755 index 0000000000000000000000000000000000000000..fb47b6bbacbcf7fb1a70a7dbeff1f4c8e64b65b8 --- /dev/null +++ b/src/downscaling_temp.f90 @@ -0,0 +1,664 @@ +! Program to downscale maximum and minimum temperature based on analogs method +! for Iberian Peninsula and Balearic Islands (Autor: Petisco de Lara) +! ****************************************************** + +SUBROUTINE down_temp(ic,id,nd,nm,nlat,nlon,nlatt,nlont,slat,slon,rlat,rlon,& + slatt,slont,ngridd,u500,v500,t500,t850,msl_si,q700,& + t700,tm2m,tmx_hr,tmn_hr,nanx,nvar,dia,mes,um,vm,& + insol,neni,vdmin,vref4,u500e,v500e,t500e,t850e,& + msle,q700e,t700e,tm2me,tmax,tmin) + + +USE MOD_CSTS +USE MOD_FUNCS + + IMPLICIT NONE + +!cccccccccccccccccccccccccccccccccc +!cccccccccccccccccccccccccccccccccc + + INTEGER, INTENT(IN) :: ic + INTEGER, INTENT(IN) :: id + INTEGER, INTENT(IN) :: nd + INTEGER, INTENT(IN) :: nm + + INTEGER, INTENT(IN) :: nlat + INTEGER, INTENT(IN) :: nlon + INTEGER, INTENT(IN) :: nlatt + INTEGER, INTENT(IN) :: nlont + DOUBLE PRECISION, INTENT(IN) :: slat + DOUBLE PRECISION, INTENT(IN) :: slon + DOUBLE PRECISION, INTENT(IN) :: rlat + DOUBLE PRECISION, INTENT(IN) :: rlon + DOUBLE PRECISION, INTENT(IN) :: slatt + DOUBLE PRECISION, INTENT(IN) :: slont + INTEGER, INTENT(IN) :: ngridd + + DOUBLE PRECISION, INTENT(IN) :: u500(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: v500(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: msl_si(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: q700(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: t500(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: t850(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: t700(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: tm2m(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: tmx_hr(nd,nptos) + DOUBLE PRECISION, INTENT(IN) :: tmn_hr(nd,nptos) + + INTEGER, INTENT(IN) :: nanx + INTEGER, INTENT(IN) :: nvar + INTEGER, INTENT(IN) :: dia(nm) + INTEGER, INTENT(IN) :: mes(nm) + + DOUBLE PRECISION, INTENT(IN) :: um(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: vm(nd,ic) + DOUBLE PRECISION, INTENT(IN) :: insol(nd) + + INTEGER, INTENT(IN) :: neni + DOUBLE PRECISION, INTENT(IN) :: vdmin(nptos,4) + INTEGER, INTENT(IN) :: vref4(nptos,4) + + DOUBLE PRECISION, INTENT(IN) :: u500e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: v500e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: t500e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: t850e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: msle(nm,id) + DOUBLE PRECISION, INTENT(IN) :: q700e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: t700e(nm,ic) + DOUBLE PRECISION, INTENT(IN) :: tm2me(nm,ic) + + DOUBLE PRECISION, INTENT(OUT) :: tmax(nm,nptos) + DOUBLE PRECISION, INTENT(OUT) :: tmin(nm,nptos) + +!cccccccccccccccccccccccccccccccccc + integer m,n + integer i,j + integer ii,jp + real sp,dd + real aaa,bdlon,bilat,dim,dift,ccm + integer jv,kk,kki,ik,it,kp,imes,ida,ida2,idia,mi,mm +! + integer i1,i2,i3,i4,i7,iana,ice,ien,ipos,ips,ipu,ir,ire,iti,jk,k + integer mesa,nan,nan2,nanf,nanv,nen,nmm,np + real disu5,disu9,disv5,disv9,dmin,dt,du5,du9,dv5,dv9,supo + real vorm,vorz +! + character sc*8,pt*9 + real u9(nd,ic),v9(nd,ic),u5(nd,ic),v5(nd,ic) + real mu9(ic),su9(ic),mv9(ic),sv9(ic) + real mu5(ic),su5(ic),mv5(ic),sv5(ic) + real p9(ic),p5(ic),rlt(ic),rln(ic),rltt(id),rlnt(id) + real inso(nd),t8(nd,ic),t7(nd,ic),t5(nd,ic),he7(nd,ic) + real efan(nd,ic),psl(nd,ic),pslm(ic),t2(nd,ic),t2m(ic) + real pres(id),bar(id),ser(nd),md,sg + integer ior(nd),ref(nptos),puce(neni) + + real um9(ic),vm9(ic),um5(ic),vm5(ic),he7m(ic) + real t8m(ic),t7m(ic),t5m(ic) + real insom,pred1(nps,nd,neni),pred1m(nps,neni),pred(nvar) + integer ana(nanx),anai(nanx) + real tmxr(nd,nptos),tmir(nd,nptos),dis(nd) + real dato1(nvar,nanx),tempx(nanx),tempi(nanx) + real ccmux(nptos),ccmui(nptos) + integer kvars(nvar) + real corrpar(nvar) + real tmxes(nm,nptos),tmies(nm,nptos) + real coe(nvar),con +!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + integer Vref(nptos,4) + real Vdis(nptos,4) + integer iii + real distancion1, distancion2, distancion3, distancion4 + real peso1, peso2, peso3, peso4, calculin + integer ien1, ien2, ien3, ien4 + +!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + m=nm + n=nd + +!******************************** +! 1. Synoptic latitude and longitude calculation and assignment of +! weights windows +! + do j=1,ic + rlt(j)=slat+(((j-1)/nlon)*rlat) + rln(j)=slon+((mod(j-1,nlon)+1-1)*rlon) + enddo + + p9=0. + p5=1. + do i1=1,ic + if((rlt(i1).le.fnor2).and.(rlt(i1).ge.fsur2)) then + if((rln(i1).ge.foes2).and.(rln(i1).le.fest2)) then + p9(i1)=1. + p5(i1)=4. + endif + endif + enddo + do i1=1,ic + if((rlt(i1).le.fnor1).and.(rlt(i1).ge.fsur1)) then + if((rln(i1).ge.foes1).and.(rln(i1).le.fest1)) then + p9(i1)=2. + p5(i1)=8. + endif + endif + enddo + +! +! Latitude and longitude calculation in the extended domain (called low +! resolution) + + do j=1,id + rltt(j)=slatt+(((j-1)/nlont)*rlat) + rlnt(j)=slont+((mod(j-1,nlont)+1-1)*rlon) + enddo + +! A reference centers (matching points between synoptic and high +! resolution grids) are define to know where the predictor must be +! calculated. + + Vref(:,:)=vref4(:,:) + Vdis(:,:)=vdmin(:,:) + + nen=1 + puce(1)=Vref(1,1) + + do iii=1,4 + do j=1,nptos + do k=1,nen + if (Vref(j,iii).eq.puce(k)) go to 101 + enddo + nen=nen+1 + ipos=nen + puce(ipos)=Vref(j,iii) + 101 continue + enddo + enddo + +! +!**************************************************************** +! TRAINING REANALYSIS VARIABLES + + u9(:,:)=um(:,:) + v9(:,:)=vm(:,:) + u5(:,:)=u500(:,:) + v5(:,:)=v500(:,:) + psl(:,:)=msl_si(:,:) + he7(:,:)=q700(:,:) + t5(:,:)=t500(:,:) + t8(:,:)=t850(:,:) + t7(:,:)=t700(:,:) + t2(:,:)=tm2m(:,:) + +! INSOLATION PARAMETER + + inso(:)=insol(:) + +! The predictors are obtained +! OBTAINING THE GEOSTROPHIC U/V COMPONENTS IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(2,i,j)=u9(i,ice) + pred1(3,i,j)=v9(i,ice) + enddo + enddo + +! OBTAINING THE SEA LEVEL PRESSURE IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(4,i,j)=psl(i,ice) + enddo + enddo + +! OBTAINING THE 850 hPa TEMPERATURE IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(1,i,j)=t8(i,ice) + enddo + enddo + +! OBTAINING THE ESPECIFIC HUMIDITY IN 700 hPa IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(5,i,j)=he7(i,ice) + enddo + enddo + +! OBTAINING THE 2 METERS TEMPERATURE IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(8,i,j)=t2(i,ice) + enddo + enddo + +! OBTAINING THE 700 hPa TEMPERATURE IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(6,i,j)=t7(i,ice) + enddo + enddo + +! OBTAINING THE 500 hPa TEMPERATURE IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(7,i,j)=t5(i,ice) + enddo + enddo + +! MEAN AND DEVIATION OF REFERENCE FIELDS (WINDS) + + do j=1,ic + do i=1,n + ser(i)=u9(i,j) + enddo + call estadis(ser,md,sg,n) + mu9(j)=md + su9(j)=sg + do i=1,n + ser(i)=v9(i,j) + enddo + call estadis(ser,md,sg,n) + mv9(j)=md + sv9(j)=sg + do i=1,n + ser(i)=u5(i,j) + enddo + call estadis(ser,md,sg,n) + mu5(j)=md + su5(j)=sg + do i=1,n + ser(i)=v5(i,j) + enddo + call estadis(ser,md,sg,n) + mv5(j)=md + sv5(j)=sg + enddo + +! REFERENCE WINDS ARE NORMALIZED + + do i=1,n + do j=1,ic + u9(i,j)=(u9(i,j)-mu9(j))/su9(j) + v9(i,j)=(v9(i,j)-mv9(j))/sv9(j) + u5(i,j)=(u5(i,j)-mu5(j))/su5(j) + v5(i,j)=(v5(i,j)-mv5(j))/sv5(j) + enddo + enddo + +! HIGH RESOLUTION (5KM) MAXIMUM AND MINIMUM OBSERVED TEMPERATURE + + tmxr(:,:)=tmx_hr(:,:) + tmir(:,:)=tmn_hr(:,:) + +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!************************************************************** +!************************************************************** +! +! DOWNSCALING BEGINS (ESTIMATING THE PROBLEM DAYS MAXIMUM AND +! MINIMUM TEMPERATURES IN EACH HIGH RESOLUTION GRID POINT) +! +!************************************************************** +!************************************************************** +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +! print*,'Downscaling begins...' +! print*,'estimated day... ' + + mm=0 + do 1000 i=1,m + +! ESTIMATED REANALYSIS VARIABLES + + um5(:)=u500e(i,:) + vm5(:)=v500e(i,:) + t8m(:)=t850e(i,:) + t7m(:)=t700e(i,:) + t5m(:)=t500e(i,:) + he7m(:)=q700e(i,:) + t2m(:)=tm2me(i,:) + pres(:)=msle(i,:) + +! print*,' ',i + mm=mm+1 + +! Pressure synoptic grid calculation from low resolution one + + bilat=slat+(nlat-1)*rlat + bdlon=slon+(nlon-1)*rlon + + ire=0 + do 111 j=1,id + if((rltt(j).gt.slat).or.(rltt(j).lt.bilat)) go to 111 + if((rlnt(j).lt.slon).or.(rlnt(j).gt.bdlon)) go to 111 + ire=ire+1 + pslm(ire)=pres(j) + 111 continue + +! Geostrophic wind components at sea level to the estimated +! day (pressure in Pa). +! "/g" to invalidate gravity acceleration. + + bar=pres*100./9.81 + + call geostrofico(bar,um9,vm9,id,ic,slatt,slont,slat,slon,& + rlat,rlon,rlat,rlon,nlatt,nlont,nlat,nlon,ngridd) + +! It is divided by density at sea level (standard atmosphere) to obtain +! the geostrophic wind components. + + do j=1,ic + um9(j)=um9(j)/1.225 + vm9(j)=vm9(j)/1.225 + enddo + +! The estimated predictors are obtained +! OBTAINING THE GEOSTROPHIC U/V COMPONENTS IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(2,j)=um9(ice) + pred1m(3,j)=vm9(ice) + enddo + +! OBTAINING THE SEA LEVEL PRESSURE IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(4,j)=pslm(ice) + enddo + +! OBTAINING THE 850 hPa TEMPERATURE IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(1,j)=t8m(ice) + enddo + +! OBTAINING THE ESPECIFIC HUMIDITY IN 700 hPa IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(5,j)=he7m(ice) + enddo + +! OBTAINING THE 2 METERS TEMPERATURE IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(8,j)=t2m(ice) + enddo + +! OBTAINING THE 700 hPa TEMPERATURE IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(6,j)=t7m(ice) + enddo + +! OBTAINING THE 500 hPa TEMPERATURE IN THE REFERENCE CENTERS + + do j=1,nen + ice=puce(j) + pred1m(7,j)=t5m(ice) + enddo + +! REFERENCE WINDS ARE NORMALIZED + + do j=1,ic + um9(j)=(um9(j)-mu9(j))/su9(j) + vm9(j)=(vm9(j)-mv9(j))/sv9(j) + um5(j)=(um5(j)-mu5(j))/su5(j) + vm5(j)=(vm5(j)-mv5(j))/sv5(j) + enddo + +! INSOLATION PARAMETER ARE CALCULATED + + idia=dia(i) + imes=mes(i) + call fechanno(idia,imes,ida) + ida2=ida-80 + if(ida2.le.0) ida2=ida2+365 + aaa=2.*pi*float(ida2)/365. + insom=sin(aaa) + + +! Synoptic type determination: the "nanx" reference alements +! more similar to each synoptic type and the corresponding +! distances. + + do k=1,n + ior(k)=k + dis(k)=9999. + enddo + + do 110 j=1,n + call distan9_2(um9,u9,n,j,p9,disu9,ic) + call distan9_2(vm9,v9,n,j,p9,disv9,ic) + call distan5_2(um5,u5,n,j,p5,disu5,ic) + call distan5_2(vm5,v5,n,j,p5,disv5,ic) + dim=(disu9+disv9+disu5+disv5)/4. + dis(j)=dim + 110 continue + call burbuja1(dis,ior,n,nanx) + do j=1,nanx + anai(j)=ior(j) + enddo + + do 1200 ipu=1,nptos + +! Reference environment + + do ik=1,nen + ice=puce(ik) + + if (Vref(ipu,1).eq.ice) then + ien1=ik + distancion1=Vdis(ipu,1) + peso1=1/distancion1 + go to 251 + endif + enddo + 251 continue + + do ik=1,nen + ice=puce(ik) + if (Vref(ipu,2).eq.ice) then + ien2=ik + distancion2=Vdis(ipu,2) + peso2=1/distancion2 + go to 252 + endif + enddo + 252 continue + + do ik=1,nen + ice=puce(ik) + if (Vref(ipu,3).eq.ice) then + ien3=ik + distancion3=Vdis(ipu,3) + peso3=1/distancion3 + go to 253 + endif + enddo + 253 continue + + do ik=1,nen + ice=puce(ik) + if (Vref(ipu,4).eq.ice) then + ien4=ik + distancion4=Vdis(ipu,4) + peso4=1/distancion4 + go to 254 + endif + enddo + 254 continue + + +! Predictors for the estimated day +! INSOLATION PREDICTOR + + pred(1)=insom + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + do jp=1,nps + kp=jp+1 + calculin = pred1m(jp,ien1)*peso1+pred1m(jp,ien2)*peso2+ & + pred1m(jp,ien3)*peso3+pred1m(jp,ien4)*peso4 + pred(kp)= calculin/(peso1+peso2+peso3+peso4) + enddo + +!!!!!!!!!!!!!!!!!!!!!!!!! +! FOR MAXIMUM TEMPERATURE +!!!!!!!!!!!!!!!!!!!!!!!!! + + nan=0 + do 1201 i2=1,nanx + iana=anai(i2) + if(tmxr(iana,ipu).ne.-999.) then + nan=nan+1 + ana(nan)=iana + + dato1(1,nan)=inso(iana) + do jp=1,nps + kp=jp+1 + dato1(kp,nan) = (pred1(jp,iana,ien1)*peso1+ & + pred1(jp,iana,ien2)*peso2 + pred1(jp,iana,ien3)*peso3+ & + pred1(jp,iana,ien4)*peso4 ) / (peso1+peso2+peso3+peso4) + enddo + + tempx(nan)=tmxr(iana,ipu) + else + go to 1201 + endif + 1201 continue + + if(nan.gt.150) nan=150 + +! Calculation of significant predictors, their coeficients and their +! multiple and partial correlation coeficients to estimate the +! maximum temperature +! +! mi: number of selected predictors +! ccm: multiple correlation coeficient +! kvars: selected predictors labels (vector) +! corrpar: partial correlation of selected predictors (vector) +! coe: regression coeficients associated to each predictor (vector). +! (value = 0 when there is no selected predictor). +! con: Y-intercept (independent equation term) +! tol: tolerance to select predictors + + call stepregrs & + (tempx,dato1,nanx,nvar,nan,mi,ccm,kvars,corrpar,coe,con,tol) + +! Maximum temperature estimation. When there are no significant predictors, +! estimated temperature is the temperature of the analog that has the 2 +! meters temperature more similar to the estimated day. + + if(mi.eq.0) then + dift=999999. + do kk=1,nan + if(abs(pred(9)-dato1(9,kk)).lt.dift) then + kki=kk + dift=abs(pred(9)-dato1(9,kk)) + endif + enddo + tmxes(i,ipu)=tmxr(ana(kki),ipu) + else + tmxes(i,ipu)=con + do jv=1,nvar + tmxes(i,ipu)=coe(jv)*pred(jv)+tmxes(i,ipu) + enddo + endif + + 1203 CONTINUE + +!!!!!!!!!!!!!!!!!!!!!!!!! +! FOR MINIMUM TEMPERATURE +!!!!!!!!!!!!!!!!!!!!!!!!! + + nan=0 + do 1202 i2=1,nanx + iana=anai(i2) + if(tmir(iana,ipu).ne.-999.) then + nan=nan+1 + ana(nan)=iana +!!!! +! With NAN observed tmin data, next lines should be included +!c dato1(1,nan)=inso(iana) +!c do jp=1,nps +!c kp=jp+1 +!c dato1(kp,nan)=pred1(jp,iana,ien) +!c enddo +!!!! + tempi(nan)=tmir(iana,ipu) + else + go to 1202 + endif + 1202 continue + + if(nan.gt.150) nan=150 + +! Calculation of significant predictors, their coeficients and their +! multiple and partial correlation coeficients to estimate the +! minimum temperature +! +! mi: number of selected predictors +! ccm: multiple correlation coeficient +! kvars: selected predictors labels (vector) +! corrpar: partial correlation of selected predictors (vector) +! coe: regression coeficients associated to each predictor (vector). +! (value = 0 when there is no selected predictor). +! con: Y-intercept (independent equation term) +! tol: tolerance to select predictors + + call stepregrs & + (tempi,dato1,nanx,nvar,nan,mi,ccm,kvars,corrpar,coe,con,tol) + +! Minimum temperature estimation. When there are no significant predictors, +! estimated temperature is the temperature of the analog that has the 2 +! meters temperature more similar to the estimated day. + + if(mi.eq.0) then + dift=999999. + do kk=1,nan + if(abs(pred(9)-dato1(9,kk)).lt.dift) kki=kk + enddo + tmies(i,ipu)=tmir(ana(kki),ipu) + else + tmies(i,ipu)=con + do jv=1,nvar + tmies(i,ipu)=coe(jv)*pred(jv)+tmies(i,ipu) + enddo + endif + + 1200 continue + +!!!!!!!!!!!!!!!!!!!!!!!!!!!! + 1000 continue !End of estimated days loop +!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + tmax(:,:)=nint(10.*tmxes(:,:)) + tmin(:,:)=nint(10.*tmies(:,:)) + +!++++++++++++++++++++++++++++++++++++++++++++++++++ + +END SUBROUTINE down_temp + + diff --git a/src/insol.f90 b/src/insol.f90 new file mode 100755 index 0000000000000000000000000000000000000000..df4717377d0e2b683bee38dad5aa9dc5b226ca38 --- /dev/null +++ b/src/insol.f90 @@ -0,0 +1,38 @@ +! The insol program calculates insolation of 'nd' period + +SUBROUTINE insolation(nd,day,month,insol) + +USE MOD_CSTS +USE MOD_FUNCS, ONLY : fechanno + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nd + INTEGER, INTENT(IN) :: day(nd) + INTEGER, INTENT(IN) :: month(nd) + + DOUBLE PRECISION, INTENT(OUT) :: insol(nd) +! REAL, INTENT(OUT) :: insol(nd) + + integer i,ida,ida2 + integer dd,mm + real aaa + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! print*,'program 3: insolation' + + do 1000 i=1,nd + dd=day(i) + mm=month(i) + call fechanno(dd,mm,ida) + + ida2=ida-80 + if(ida2.le.0) ida2=ida2+365 + aaa=2.*pi*float(ida2)/365. + insol(i)=sin(aaa) + 1000 continue + +END SUBROUTINE insolation + + diff --git a/src/mod_csts.f90 b/src/mod_csts.f90 new file mode 100755 index 0000000000000000000000000000000000000000..eb98684877e3a6eaf943144fd1d981060d091a5a --- /dev/null +++ b/src/mod_csts.f90 @@ -0,0 +1,67 @@ +! !!!!!!!!!!!!!!! + MODULE MOD_CSTS +! !!!!!!!!!!!!!!! +!* 0. DECLARATIONS +! ------------ +! +IMPLICIT NONE + +!CHARACTER(LEN=100),SAVE :: pathSI_estim="/home/sclim/cle/regionalizacion/asun/eraint/datosAscii_SI_estim/" +!CHARACTER(LEN=100),SAVE :: pathLR_estim="/home/sclim/cle/regionalizacion/asun/eraint/datosAscii_LR_estim/" +!*********************************************** +! CONSTANTS +!*********************************************** +! +CHARACTER(LEN=5),SAVE :: mt="aemet" +INTEGER,SAVE :: nvar_p=11 ! total number of predictors for precipitation +INTEGER,SAVE :: nvar_t=9 ! total number of predictors for temperature +INTEGER,SAVE :: nps=8 ! number of predictors from the reference grid point of the model +INTEGER, SAVE :: npx=11 +REAL,SAVE :: umb=0.75 +REAL,SAVE :: umbl=0.75 +INTEGER,SAVE :: nmin=10 +INTEGER,SAVE :: nmax=30 +REAL,SAVE :: ccmi=0.30 +!----------------------------------------------- +REAL,SAVE :: r=287.05 +REAL,SAVE :: a=0.0065 +REAL,SAVE :: g=9.81 +REAL,SAVE :: rt=6.37e6 +REAL,SAVE :: pi=3.14159265 +REAL,SAVE :: romega=7.292E-5 +!----------------------------------------------- +INTEGER,SAVE :: umb_ger=10000 +!----------------------------------------------- +!REAL,SAVE :: huso=30. +DOUBLE PRECISION,SAVE :: huso=30. +!----------------------------------------------- +REAL,SAVE :: tol=0.020 +! +!*********************************************** +! DOMAIN +!************************************************ +!------------------------------------------------ +! Limits (lat-lon) of strips to assign weights +! +REAL,SAVE :: fnor1=45.0 +REAL,SAVE :: fnor2=47.5 +REAL,SAVE :: fsur1=35.0 +REAL,SAVE :: fsur2=32.5 +REAL,SAVE :: foes1=-12.5 +REAL,SAVE :: foes2=-17.5 +REAL,SAVE :: fest1=5.0 +REAL,SAVE :: fest2=7.5 +!------------------------------------------------ +! High resolution (5km x 5km) observed AEMET grid +! +! For precipitation: the variable is multiplied to 10, to +! obtain it in decimes of mm, as the program requaire. +INTEGER, PARAMETER :: nptos=20945 ! number of point where make the estimations +!------------------------------------------------ +INTEGER,SAVE :: icl=24 +INTEGER,SAVE :: ipui=99 +INTEGER,SAVE :: ila=4 +INTEGER,SAVE :: ilo=6 + +END MODULE MOD_CSTS + diff --git a/src/mod_funcs.f90 b/src/mod_funcs.f90 new file mode 100755 index 0000000000000000000000000000000000000000..442750d49de1a35263314449fcd94cc47f70ea95 --- /dev/null +++ b/src/mod_funcs.f90 @@ -0,0 +1,1302 @@ + +MODULE MOD_FUNCS + +CONTAINS + + SUBROUTINE ESTADIS(SER,MEDIA,SIGMA,N) +! CALCULA LA MEDIA Y LA DESVIACION TIPO DE UNA SERIE DE DATOS + REAL SER(N),MEDIA,SIGMA + MEDIA=0. + DO I=1,N + MEDIA=MEDIA+SER(I) + ENDDO + MEDIA=MEDIA/REAL(N) + SIGMA=0. + DO I=1,N + SIGMA=SIGMA+(SER(I)-MEDIA)**2 + ENDDO + SIGMA=SIGMA/REAL(N) + SIGMA=SQRT(SIGMA) + RETURN + END SUBROUTINE + + SUBROUTINE BURBUJA(A,NOR,NAN,NAN1,IAN) + REAL A(NAN) + INTEGER NOR(NAN) + DO 100 I=1,IAN + DO 110 J=I+1,NAN1 + IF(A(I).GT.A(J)) THEN + TEM=A(J) + ITEM=NOR(J) + A(J)=A(I) + NOR(J)=NOR(I) + A(I)=TEM + NOR(I)=ITEM + ENDIF + 110 CONTINUE + 100 CONTINUE + RETURN + END SUBROUTINE + + subroutine burbuja1(a,nor,n,nan) + real a(n) + integer nor(n) + do 100 i=1,nan + do 110 j=i+1,n + if(a(i).gt.a(j)) then + tem=a(j) + item=nor(j) + a(j)=a(i) + nor(j)=nor(i) + a(i)=tem + nor(i)=item + endif + 110 continue + 100 continue + return + end subroutine + + subroutine distan9(ca,n,ic,i,nr,p,dis) + real ca(n,ic),p(ic) + dis=0. + sp=0. + do 10 k=1,ic + if(p(k).eq.0.) go to 10 + dis=dis+p(k)*(ca(i,k)-ca(nr,k))**2 + sp=sp+p(k) + 10 continue + dis=dis/sp + return + end subroutine + + subroutine distan5(ca,n,ic,i,nr,p,dis) + real ca(n,ic),p(ic) + dis=0. + sp=0. + do k=1,ic + dis=dis+p(k)*(ca(i,k)-ca(nr,k))**2 + sp=sp+p(k) + enddo + dis=dis/sp + return + end subroutine + +! + subroutine distancia9(ca,n,cg,m,i,nr,p,dis,ic) +! implicit none +! real ca(n,ic),cg(m,ic),p(ic) + real ca(n,ic),p(ic) +! real cg(m,ic) + double precision cg(m,ic) + dis=0. + sp=0. + do 100 k=1,ic + if(p(k).eq.0.) go to 100 + dis=dis+p(k)*(ca(i,k)-cg(nr,k))**2 + sp=sp+p(k) + 100 continue + dis=dis/sp + return + end subroutine + + subroutine distancia5(ca,n,cg,m,i,nr,p,dis,ic) +! implicit none +! real ca(n,ic),cg(m,ic),p(ic) + real ca(n,ic),p(ic) + double precision cg(m,ic) + dis=0. + sp=0. + do k=1,ic + dis=dis+p(k)*(ca(i,k)-cg(nr,k))**2 + sp=sp+p(k) + enddo + dis=dis/sp + return + end subroutine + + subroutine distan9_2(cb,ca,n,nr,p,dis,ic) + real cb(ic),ca(n,ic),p(ic) + dis=0. + sp=0. + do 100 k=1,ic + if(p(k).eq.0.) go to 100 + dis=dis+p(k)*(cb(k)-ca(nr,k))**2 + sp=sp+p(k) + 100 continue + dis=dis/sp + return + end subroutine + + subroutine distan5_2(cb,ca,n,nr,p,dis,ic) + real cb(ic),ca(n,ic),p(ic) + dis=0. + sp=0. + do k=1,ic + dis=dis+p(k)*(cb(k)-ca(nr,k))**2 + sp=sp+p(k) + enddo + dis=dis/sp + return + end subroutine + + subroutine distancia9_2(ca,n,cg,m,i,nr,p,dis,ic) + real ca(n,ic),cg(m,ic),p(ic) + dis=0. + sp=0. + do 100 k=1,ic + if(p(k).eq.0.) go to 100 + dis=dis+p(k)*(ca(i,k)-cg(nr,k))**2 + sp=sp+p(k) + 100 continue + dis=dis/sp + return + end subroutine + + subroutine distancia5_2(ca,n,cg,m,i,nr,p,dis,ic) + real ca(n,ic),cg(m,ic),p(ic) + dis=0. + sp=0. + do k=1,ic + dis=dis+p(k)*(ca(i,k)-cg(nr,k))**2 + sp=sp+p(k) + enddo + dis=dis/sp + return + end subroutine + + SUBROUTINE STEPREGRS & + (YI,XI,NX,NVARX,N,MI,CCM,IVAR,COPA,COE,CON,TOL) + +!EFECTUA UNA REGRESION LINEAL MULTIPLE POR ETAPAS +!MEDIANTE LA TECNICA 'PASO A PASO' INTRODUCIENDO +!EN CADA PASO COMO NUEVA VARIABLE LA DE MEJOR +!CORRELACION PARCIAL CON LA VARIABLE DEPENDIENTE +!Y ELIMINANDO AQUELLAS QUE DESPUES DE CADA +!REGRESION NO SEAN SIGNIFICATIVAS + +!NX=numero maximo de datos posibles +!N=numero de datos actuales a usar +!NVAR= numero total de variables posibles +! de regresion + +!LA VARIABLE YI(NX) CONTIENE LOS VALORES DEL PREDICTANDO +!LA VARIABLE XI(NVARX,NX) CONTIENE LOS VALORES DE LOS PREDICTORES + +!LA VARIABLE YY(N) CONTIENE LOS VALORES DEL PREDICTANDO +!LA VARIABLE XX(NVARX,N) CONTIENE LOS VALORES DE LOS PREDICTORES + +! MI es el numero de variables o predictores seleccionados + +! CCM es el coeficiente de correlacion multiple + +! IVAR(NVARX) contiene los numeros de etiqueta de los predictores +! seleccionados en la regresion + +!COPA(NVARX) contiene las correlaciones parciales de +! los predictores seleccionados. + +! COE(NVARX) contiene los coeficientes de regresion beta(i) de las +! variables seleccionadas + +! CON contiene la constante de la regresion (beta0) + +! DATO1(NVARX,N) contiene los datos de los predictores +! que se meten en cada paso de regresion + +! COEF(0:NVARX) contiene la constante y los coeficientes +! de regresion de las variables introducidas en cada paso + +! YYES(N) contiene valores del predictando estimados por la regresion + +! SST es la variabilidad total, +! SSE es la variabilidad residual no explicada por la regresion + +! CDET es el coeficiente de determinacion en el paso actual +! CDATA es el coeficiente de determinacion en el paso anterior + +! CDETP(NVARX) contiene los coeficientes de determinacion +! cuando se considera cada variable como introducida en el paso +! actual. Se utiliza como base para eliminar variables que +! se consideran no significativas + +! TOL representa el minimo incremento de variabilidad explicada +! por la introduccion de una variable para que esta se +! considere significativa + + + real yi(nx),xi(nvarx,nx) + real yy(n),xx(nvarx,n),res1(n),res(n),ser1(n),ser2(n),aa(n) + real yyr(n),cdet1,cp,ay(n),cormax,cor,ccm + real copa(nvarx),dato1(nvarx,n) + real coef(0:nvarx),coe(nvarx),con + real yyes(n),sst,sse,cdet,cdeta,cdetp(nvarx),myy,incr + + character var(nvarx)*5 + integer ivar(nvarx),ivar1(nvarx) + +!TRASPASAMOS LOS DATOS INICIALES ENVIADOS POR EL PROGRAMA PRINCIPAL +!DESDE LAS VARIABLES DE DIMENSION MAXIMA A LAS VARIABLES CON LA +!DIMENSION AJUSTADA AL NUMERO ACTUAL DE DATOS UTILIZADOS + + do i=1,n + yy(i)=yi(i) + do k=1,nvarx + xx(k,i)=xi(k,i) + enddo + enddo + + +!CALCULO DE LA MEDIA Y DE LA VARIABILIDAD TOTAL DEL PREDICTANDO + + myy=0. + do i=1,n + myy=myy+yy(i) + enddo + myy=myy/real(n) + + sst=0. + do i=1,n + sst=sst+(yy(i)-myy)**2 + enddo + + + + + +!******************************************************** +!INICIALIZACION DEL CONTROL +!DE LAS VARIABLES INTRODUCIDAS EN EL MODELO +!EN CADA PASO + + do j=1,nvarx + var(j)='nosel' + enddo + +!**************************************************** + +!BUSQUEDA DE LA PRIMERA VARIABLE DEL MODELO +! (UNICA VARIABLE EN LA PRIMERA ETAPA) + + cdeta=0. + cormax=-2.0 + nvx=0 + do j=1,nvarx + do i=1,n + ser2(i)=xx(j,i) + enddo + call corr1(yy,ser2,n,cor) + if(abs(cor).gt.cormax) then + cormax=abs(cor) + nvx=j + endif + enddo + var(nvx)='sisel' + nvult=nvx + + +!**************************************************** +!PREPARACION DE LA MATRIZ DE DATOS PARA EL +!CALCULO DE LA REGRESION DE LAS VARIABLES +!INDEPENDIENTES SELECCIONADAS CON LA VARIABLE +!DEPENDIENTE + + 222 nuvar=0 + dato1=0. + ivar=0 + do 100 j=1,nvarx + if(var(j).ne.'sisel') go to 100 + nuvar=nuvar+1 + ivar(nuvar)=j + do i=1,n + dato1(nuvar,i)=xx(j,i) + enddo + 100 continue + +!SE CALCULA LA REGRESION CON LAS VARIABLES SELECCIONADAS + + yyr=yy + call regr(yyr,dato1,nvarx,nuvar,n,coef) + + + +!SE CALCULA EL COEFICIENTE DE DETERMINACION (esta subrutina +! devuelve los residuos de la regresion y el coeficiente de +! determinacion) + + call coedet(yy,xx,n,nvarx,ivar,nuvar,coef,res1,sst,cdet) + + +!COMPROBAMOS SI EL COEFICIENTE DE DETERMINACION SE HA INCREMENTADO +!SUFICIENTEMENTE COMO PARA CONSIDERAR SIGNIFICATIVA LA ULTIMA VARIABLE +!INTRODUCIDA. SI NO LO ES SE ACABA EL PROCESO PASO A PASO Y OBTENEMOS +!LA REGRESION DEFINITIVA + + incr=cdet-cdeta + if(incr.lt.tol) then + if(nuvar.eq.1) then + mi=0 + ccm=-8.88 + go to 555 + else + var(nvult)='elimi' + go to 444 + endif + endif + + cdeta=cdet + +!SE COMPRUEBA SI ALGUNA DE LAS VARIABLES SELECCIONADAS RESULTA NO +!SIGNIFICATIVA. PARA ELLO SE COMPARAN LOS COEFICIENTES DE DETERMINACION +!OBTENIDOS QUITANDO CADA VARIABLE, CON EL OBTENIDO SIN QUITAR NINGUNA +!DE LAS QUE YA TENEMOS, SI PARA ALGUNA VARIABLE EL INCREMENTO NO +!SUPERA EL MINIMO LA VARIABLE SE ELIMINA DEFINITIVAMENTE + + if(nuvar.eq.1) go to 333 + +! Quitamos una variable cada vez + + do 200 k=1,nuvar + dato1=0. + nivar=0 + ivar1=0 + do 210 k1=1,nuvar + if(k1.eq.k) go to 210 + nivar=nivar+1 + ivar1(nivar)=ivar(k1) + do i=1,n + dato1(nivar,i)=xx(ivar1(nivar),i) + enddo + 210 continue + +! Se calcula la regresion con la variable quitada + + yyr=yy + call regr(yyr,dato1,nvarx,nivar,n,coef) + +! Se calcula el coeficiente de determinacion + + call coedet(yy,xx,n,nvarx,ivar1,nivar,coef,res,sst,cdet1) + +! Si la diferencia entre el coeficiente de determinacion +! con todas las variables y el mismo con la variable +! quitada es menor que el umbral, la variable se considera +! no significativa y se elimina definitivamente. + + if((cdet-cdet1).lt.tol) then + var(ivar(k))='elimi' + endif + 200 continue + + do k=1,nuvar + if(var(ivar(k)).eq.'elimi') go to 332 + enddo + go to 333 + + + +!ELIMINADAS LAS VARIABLES NO SIGNIFICATIVAS SE CALCULA DE NUEVO +!LA REGRESION CON LAS VARIABLES QUE HAN QUEDADO + + 332 continue + nuvar=0 + dato1=0. + ivar=0 + do 220 j=1,nvarx + if(var(j).ne.'sisel') go to 220 + nuvar=nuvar+1 + ivar(nuvar)=j + do i=1,n + dato1(nuvar,i)=xx(j,i) + enddo + 220 continue + + + yyr=yy + call regr(yyr,dato1,nvarx,nuvar,n,coef) + call coedet(yy,xx,n,nvarx,ivar,nuvar,coef,res1,sst,cdet) + + cdeta=cdet + + 333 continue + +!SE COMPRUEBA SI HAY AUN VARIABLES QUE PUEDAN SER SELECCIONADAS +!SI HAY SE TRATA DE BUSCAR UNA NUEVA, SI NO HAY SE TERMINA + + do j=1,nvarx + if(var(j).eq.'nosel') go to 334 + enddo + go to 444 + + 334 continue + +!SE BUSCA UNA NUEVA VARIABLE TOMANDO LA QUE TENGA MAYOR CORRELACION +!PARCIAL CON EL PREDICTANDO + +! Se construye matriz de datos con variables ya seleccionadas + + dato1=0. + nivar=0 + do j=1,nvarx + if(var(j).eq.'sisel') then + nivar=nivar+1 + do i=1,n + dato1(nivar,i)=xx(j,i) + enddo + endif + enddo + +! Se busca nueva variable + + cormax=-2.0 + nvx=0 + do 230 j=1,nvarx + if(var(j).ne.'nosel') go to 230 + do i=1,n + aa(i)=xx(j,i) + enddo + call corpar(res1,n,dato1,nvarx,nivar,aa,cp) + if(abs(cp).gt.cormax) then + cormax=abs(cp) + nvx=j + endif + 230 continue + if (nvx.gt.0) then + var(nvx)='sisel' + nvult=nvx + endif + go to 222 + + 444 continue + +! REGRESION DEFINITIVA + +! PREPARACION DE MATRIZ DE DATOS CON VARIABLES DEFINITIVAS + + nuvar=0 + dato1=0. + ivar=0 + do 250 j=1,nvarx + if(var(j).ne.'sisel') go to 250 + nuvar=nuvar+1 + ivar(nuvar)=j + do i=1,n + dato1(nuvar,i)=xx(j,i) + enddo + 250 continue + +! CALCULO DE LA REGRESION + + + yyr=yy + call regr(yyr,dato1,nvarx,nuvar,n,coef) + +! CALCULO DEL COEFICIENTE DE DETERMINACION Y DE LOS RESIDUOS + + call coedet(yy,xx,n,nvarx,ivar,nuvar,coef,res1,sst,cdet) + + + +! RESULTADOS FINALES + +! COEFICIENTES Y DEMAS DATOS DE LA REGRESION + + mi=nuvar + ccm=sqrt(cdet) + + con=coef(0) + coe=0. + do k=1,nuvar + coe(ivar(k))=coef(k) + enddo + + +! COEFICIENTES DE CORRELACION PARCIAL DE LAS VARIABLES +! SELECCIONADAS CON LA VARIABLE DEPENDIENTE + + copa=-1. + + do 300 j=1,nuvar + do i=1,n + aa(i)=xx(ivar(j),i) + ay(i)=yy(i) + enddo + nivar=0 + dato1=0. + do k=1,nuvar + if(k.ne.j) then + nivar=nivar+1 + do i=1,n + dato1(nivar,i)=xx(ivar(k),i) + enddo + endif + enddo + call corpar1(ay,n,dato1,nvarx,nivar,aa,cp) + copa(ivar(j))=abs(cp) + 300 continue + + 555 continue + + return + end subroutine + + + SUBROUTINE CORR1(CENT,COMP,IC,CORRE1) + REAL SUM1,SUM2,MED1,MED2 + REAL CENT(IC),COMP(IC),SUMC1,SUMC2,SUMCR + REAL COV,VAR1,VAR2,CORRE1 + SUM1=0.0 + SUM2=0.0 + DO 100 I=1,IC + SUM1=SUM1+CENT(I) + SUM2=SUM2+COMP(I) +100 CONTINUE + C=REAL(IC) + MED1=SUM1/C + MED2=SUM2/C + SUMC1=0.0 + SUMC2=0.0 + SUMCR=0.0 + DO 200 J=1,IC + SUMCR=SUMCR+((CENT(J)-MED1)*(COMP(J)-MED2)) + SUMC1=SUMC1+(CENT(J)-MED1)**2 + SUMC2=SUMC2+(COMP(J)-MED2)**2 +200 CONTINUE + COV=SUMCR/C + VAR1=SUMC1/C + VAR2=SUMC2/C + CORRE1=COV/SQRT(VAR1*VAR2) + RETURN + END SUBROUTINE + + + SUBROUTINE REGR(aa,bb,nvarx,nvar,ndat,creg) + + +!CALCULA LA ECUACION DE REGRESION A PARTIR DE UNA MUESTRA DE DATOS + +!TRABAJA CON LAS DESVIACIONES RESPECTO A LA MEDIA PARA MINIMIZAR +!LOS ERRORES DE REDONDEO POR LO QUE AL FINAL HAY QUE CALCULAR +!APARTE EL TERMINO INDEPENDIENTE DE LA ECUACION DE REGRESION +!(ver libro de D. Penna capitulo regresion multiple) + + +! ndat: numero de datos de la muestra +! nvar: numero de variables independientes +! yy(ndat): contiene las desviaciones del predictando +! xx(nvar,ndat): contiene las desviaciones de los predictores +! myy,mxx(nvar): contiene las medias de predictandos y predictores + +! Elementos de las ecuaciones normales + +!nn(nvar,nvar): Matriz de los coeficientes de las incognitas +! (coeficientes de regresion salvo termino +! independiente beta0) +!b(nvar) : En entrada contiene los terminos independientes +! de las ecuaciones normales que se pasan a +! a las subrutinas que resuelven el sistema +! En salida contiene los coeficientes de regresion +! (no el termino independiente beta0) +!creg(0:nvarx) : Contiene la salida al programa principal de los +! coeficientes de regresion y del termino +! independiente +!sxx(nvar): Suma de los valores de los predictores de todos los +! datos de la muestra +!syy: Suma de los valores de los predictandos +!syyxx(nvar): Suma de productos predictando-predictores +!sxxxx(nvar): Suma de productos predictores-predictores + + + real yy(ndat) + real xx(nvar,ndat) + real aa(ndat),bb(nvarx,ndat) + real myy,mxx(nvar) + real b(nvar),nn(nvar,nvar),creg(0:nvarx) + real sxx(nvar),syy,syyxx(nvar),sxxxx(nvar,nvar),d + integer indx(nvar) + + + +! SE CALCULAN LAS MEDIAS DE LOS VALORES DE PREDICTANDOS Y PREDICTORES + + myy=0. + do i=1,ndat + myy=myy+aa(i) + enddo + myy=myy/real(ndat) + + mxx=0. + do j=1,nvar + do i=1,ndat + mxx(j)=mxx(j)+bb(j,i) + enddo + mxx(j)=mxx(j)/real(ndat) + enddo + +! SE SUSTITUYEN LOS DATOS ORIGINALES POR SUS DESVIACIONES RESPECTO +! A LA MEDIA + + do i=1,ndat + yy(i)=aa(i)-myy + enddo + + do j=1,nvar + do i=1,ndat + xx(j,i)=bb(j,i)-mxx(j) + enddo + enddo + +! CALCULO DE LA SUMA DE VALORES DE PREDICTANDO Y PREDICTORES +! DE TODOS LOS DATOS DE LA MUESTRA ASI COMO LA DE PRODUCTOS +! CRUZADOS (utiliza ya como variables las desviaciones respecto +! a las medias) +! (En realidad utilizando el modelo de regresion en desviaciones +! no utilizamos las sumas de las variables aunque las calculamos) + + syy=0. + do i=1,ndat + syy=syy+yy(i) + enddo + + + sxx=0. + syyxx=0. + do j=1,nvar + do i=1,ndat + sxx(j)=sxx(j)+xx(j,i) + syyxx(j)=syyxx(j)+yy(i)*xx(j,i) + enddo + enddo + + sxxxx=0. + do j=1,nvar + do k=j,nvar + do i=1,ndat + sxxxx(j,k)=sxxxx(j,k)+xx(j,i)*xx(k,i) + enddo + if(j.ne.k) sxxxx(k,j)=sxxxx(j,k) + enddo + enddo + + +! CONSTRUYE LA MATRIZ DE LOS COEFICIENTES DE LAS ECUACIONES NORMALES + + do j=1,nvar + do k=1,nvar + nn(j,k)=sxxxx(j,k) + enddo + enddo + + +! CONSTRUYE EL VECTOR DE TERMINOS INDEPENDIENTES DE LAS ECUACIONES +! NORMALES. EN LA SALIDA CONTENDRA LOS VALORES DE LOS COEFICIENTES +! DE REGRESION + + + do j=1,nvar + b(j)=syyxx(j) + enddo + + +! SE RESUELVE EL SISTEMA DE ECUACIONES NORMALES Y SE OBTIENEN +! LOS COEFICIENTES DE REGRESION DE CADA VARIABLE EN LA ECUACION +! DE REGRESION + + call ludcmp(nn,nvar,nvar,indx,d) + call lubksb(nn,nvar,nvar,indx,b) + + do j=1,nvar + creg(j)=b(j) + enddo + +! SE CALCULA EL TERMINO INDEPENDIENTE DE LA ECUACION DE REGRESION + + creg(0)=myy + do j=1,nvar + creg(0)=creg(0)-b(j)*mxx(j) + enddo + + + return + end subroutine + + SUBROUTINE lubksb(a,n,np,indx,b) + INTEGER n,np,indx(nP) + REAL a(np,np),b(np) + INTEGER i,ii,j,ll + REAL sum + ii=0 + do 12 i=1,n + ll=indx(i) + sum=b(ll) + b(ll)=b(i) + if (ii.ne.0)then + do 11 j=ii,i-1 + sum=sum-a(i,j)*b(j) +11 continue + else if (sum.ne.0.) then + ii=i + endif + b(i)=sum +12 continue + do 14 i=n,1,-1 + sum=b(i) + do 13 j=i+1,n + sum=sum-a(i,j)*b(j) +13 continue + b(i)=sum/a(i,i) +14 continue + return + END SUBROUTINE +! (C) Copr. 1986-92 Numerical Recipes Software !)#. +!********************************************************** +! + SUBROUTINE ludcmp(a,n,np,indx,d) + INTEGER n,np,indx(nP),NMAX + REAL d,a(np,np),TINY + PARAMETER (NMAX=500,TINY=1.0e-20) + INTEGER i,imax,j,k + REAL aamax,dum,sum,vv(NMAX) + d=1. + do 12 i=1,n + aamax=0. + do 11 j=1,n + if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) +11 continue + if (aamax .eq. 0.) then +! write (*,*) 'singular matrix in ludcmp' + else + vv(i)=1./aamax + endif +12 continue + do 19 j=1,n + do 14 i=1,j-1 + sum=a(i,j) + do 13 k=1,i-1 + sum=sum-a(i,k)*a(k,j) +13 continue + a(i,j)=sum +14 continue + aamax=0. + do 16 i=j,n + sum=a(i,j) + do 15 k=1,j-1 + sum=sum-a(i,k)*a(k,j) +15 continue + a(i,j)=sum + dum=vv(i)*abs(sum) + if (dum.ge.aamax) then + imax=i + aamax=dum + endif +16 continue + if (j.ne.imax)then + do 17 k=1,n + dum=a(imax,k) + a(imax,k)=a(j,k) + a(j,k)=dum +17 continue + d=-d + vv(imax)=vv(j) + endif + indx(j)=imax + if(a(j,j).eq.0.)a(j,j)=TINY + if(j.ne.n)then + dum=1./a(j,j) + do 18 i=j+1,n + a(i,j)=a(i,j)*dum +18 continue + endif +19 continue + return + END SUBROUTINE + + SUBROUTINE COEDET(yy,xx,n,nvarx,ivar1,nivar,coef,res,sst,cdet1) + real yy(n),yyes(n),xx(nvarx,n),res(n),sst,cdet1 + real sse + real coef(0:nvarx) + integer ivar1(nvarx) + +! ESTA SUBRUTINA DEVUELVE LOS RESIDUOS DE LA REGRESION +! Y EL COEFICIENTE DE DETERMINACION + +!SE CALCULAN LOS VALORES ESTIMADOS DEL PREDICTANDO ESTIMADOS +!CON LA REGRESION + + + do i=1,n + yyes(i)=coef(0) + do k=1,nivar + yyes(i)=yyes(i)+coef(k)*xx(ivar1(k),i) + enddo + enddo + +!SE CALCULAN LOS RESIDUOS DE LA REGRESION Y LA VARIABILIDAD +!NO EXPLICADA + + sse=0. + do i=1,n + res(i)=yy(i)-yyes(i) + sse=sse+res(i)**2 + enddo + + +!SE CALCULA EL COEFICIENTE DE DETERMINACION + + cdet1=sse/sst + cdet1=1.-cdet1 + + return + end subroutine + + + SUBROUTINE CORPAR(res1,n,dato1,nvarx,nivar,aa,cp) + real res1(n),res2(n),dato1(nvarx,n),aa(n) + real coef(0:nvarx) + real aaes(n),aar(n),cp + +! SE OBTIENE LA REGRESION DE LA VARIABLE CUYA CORRELACION +! PARCIAL SE CALCULA, CON LAS OTRAS VARIABLES PRESENTES + + aar=aa + call regr(aar,dato1,nvarx,nivar,n,coef) + +! Y SE OBTIENEN LOS VALORES ESTIMADOS POR ESTA REGRESION + + do i=1,n + aaes(i)=coef(0) + do k=1,nivar + aaes(i)=aaes(i)+coef(k)*dato1(k,i) + enddo + enddo + +! SE OBTIENEN LOS RESIDUOS CORRESPONDIENTES + + do i=1,n + res2(i)=aa(i)-aaes(i) + enddo + +! SE CALCULA LA CORRELACION PARCIAL + + call corr1(res1,res2,n,cp) + + + return + end subroutine + + SUBROUTINE CORPAR1(ay,n,dato1,nvarx,nivar,aa,cp) + real ay(n),res1(n),res2(n),dato1(nvarx,n),aa(n) + real ayes(n),aaes(n),ayr(n),aar(n),cp + real coef(0:nvarx),coefy(0:nvarx) + + +! SE OBTIENE LA REGRESION DE LA VARIABLE DEPENDIENTE CON LAS +! OTRAS VARIABLES PRESENTES DISTINTAS DE AQUELLAS CUYA CORRELACION +! PARCIAL SE QUIERE CALCULAR + + ayr=ay + call regr(ayr,dato1,nvarx,nivar,n,coefy) + +! Y SE OBTIENEN LOS VALORES ESTIMADOS POR ESTA REGRESION + + do i=1,n + ayes(i)=coefy(0) + do k=1,nivar + ayes(i)=ayes(i)+coefy(k)*dato1(k,i) + enddo + enddo + +! SE OBTIENEN LOS RESIDUOS CORRESPONDIENTES + + do i=1,n + res1(i)=ay(i)-ayes(i) + enddo + + + +! SE OBTIENE LA REGRESION DE LA VARIABLE CUYA CORRELACION +! PARCIAL SE CALCULA, CON LAS OTRAS VARIABLES PRESENTES + + aar=aa + call regr(aar,dato1,nvarx,nivar,n,coef) + +! Y SE OBTIENEN LOS VALORES ESTIMADOS POR ESTA REGRESION + + do i=1,n + aaes(i)=coef(0) + do k=1,nivar + aaes(i)=aaes(i)+coef(k)*dato1(k,i) + enddo + enddo + +! SE OBTIENEN LOS RESIDUOS CORRESPONDIENTES + + do i=1,n + res2(i)=aa(i)-aaes(i) + enddo + +! SE CALCULA LA CORRELACION PARCIAL + + call corr1(res1,res2,n,cp) + return + end subroutine + + SUBROUTINE GEOSTROFICO(& + Z,U,V,NGRID,NGRIDS,SLAT,SLON,SLATS,SLONS,RLAT,RLON,RLATS,RLONS,& + NLAT,NLON,NLATS,NLONS,NGRIDD) + +USE MOD_CSTS + +! CALCULA VIENTO GEOSTROFICO A PARTIR DE CAMPOS DE DOBLE RESOLUCION, +! CON CALCULO CENTRADO, Y DEVUELVE VALORES EN LOS PUNTOS DE GRID DE +! RESOLUCION NORMAL + + IMPLICIT INTEGER(K) + +! --------- Parametros de la REJILLA DE BAJA RESOLUCION ----- +! NGRID es el npuntos de la rejilla de baja resolucion +! NGRIDD se deja igual +! SLAT es la latitud de la rejilla de baja resolucion -latitud superior +! izda, SLON es la longitud ...........de la rej baja resol-longitud +! superior oeste NLAT es el numero de latitudes + +! PARAMETER (NGRID=${id},NGRIDD=${ngridd},ROMEGA=${romega}) +! PARAMETER (SLAT=${slatt},SLON=${slont},RLAT=${rlat}, +! $ RLON=${rlon},NLAT=${nlatt},NLON=${nlont}) +! ---------------------------------------------------------------------- +! ------------- Parametros de la REJILLA SINOPTICA --------- +! PARAMETER (NGRIDS=${ic},NLATS=${nlat},NLONS=${nlon}) !TERMINADO EN S, +! GRID DE SALIDA +! PARAMETER (SLATS=${slat},SLONS=${slon},RLATS=${rlat}, +! $ RLONS=${rlon}) +! PARAMETER (GR=${g},RT=${rt},R=${r},PI=${pi}) +! ---------------------------------------------------------------------- + INTEGER, INTENT(IN) :: ngrid + INTEGER, INTENT(IN) :: ngrids + + INTEGER, INTENT(IN) :: nlat + INTEGER, INTENT(IN) :: nlon + INTEGER, INTENT(IN) :: nlats + INTEGER, INTENT(IN) :: nlons + DOUBLE PRECISION, INTENT(IN) :: slat + DOUBLE PRECISION, INTENT(IN) :: slon + DOUBLE PRECISION, INTENT(IN) :: slats + DOUBLE PRECISION, INTENT(IN) :: slons + DOUBLE PRECISION, INTENT(IN) :: rlat + DOUBLE PRECISION, INTENT(IN) :: rlon + DOUBLE PRECISION, INTENT(IN) :: rlats + DOUBLE PRECISION, INTENT(IN) :: rlons + INTEGER, INTENT(IN) :: ngridd + + REAL GR + +! RLX ES LA LONGITUD DEL PASO DE REJILLA SOBRE EL PARALELO CORRESPONDIENTE +! RLY ES LA LONGITUD DEL PASO DE REJILLA SOBRE EL MERIDIANO CORRESPONDIENTE + +! OJO, RLY Y RLX CORRESPONDEN A RESOL NORMAL, Y SON 2Ay Y 2Ax DE LA RES DOBLE + REAL, INTENT(IN) :: Z(NGRID) + REAL, INTENT(OUT) :: U(NGRIDS),V(NGRIDS) + + REAL RLX(NGRID),RLY + REAL F(NGRIDD) + REAL GG(NGRID),GD(NGRIDD) + REAL RLT(NGRID),RLN(NGRID) + REAL RLTS(NGRIDS),RLNS(NGRIDS) + + GR=g +! print*,'En GEOSTROFICO, g,ngrid,ngrids= ',g,ngrid,ngrids +! print*,'En GEOSTROFICO, ngrid= ',ngrid +! print*,'En GEOSTROFICO, ngrids= ',ngrids +! print*,'En GEOSTROFICO, ngridd= ',ngridd +! print*,'En GEOSTROFICO, nlat= ',nlat +! print*,'En GEOSTROFICO, nlon= ',nlon +! print*,'En GEOSTROFICO, nlats= ',nlats +! print*,'En GEOSTROFICO, nlons= ',nlons +! print*,'En GEOSTROFICO, slat= ',slat +! print*,'En GEOSTROFICO, slon= ',slon +! print*,'En GEOSTROFICO, slats= ',slats +! print*,'En GEOSTROFICO, slons= ',slons +! print*,'En GEOSTROFICO, rlat= ',rlat +! print*,'En GEOSTROFICO, rlon= ',rlon +! print*,'En GEOSTROFICO, rlats= ',rlats +! print*,'En GEOSTROFICO, rlons= ',rlons + +! +! CALCULA LATITUD Y LONGITUD DE CADA PUNTO J + DO J=1,NGRID + RLT(J)=SLAT+(((J-1)/NLON)*RLAT) + RLN(J)=SLON+((MOD(J-1,NLON)+1-1)*RLON) + +! IF(J.GE.1.AND.J.LE.50) THEN +! print*,"J, RLT = ",J,RLT(J) +! print*," RLN = ",RLN(J) +! ENDIF + ENDDO +! print*,"fuera del bucle de calculo" +! print*,"RLT=",RLT(1:50) +! print*,"RLN",RLN(1:50) + NLOND=(NLON*2)-1 +! +! CALCULA LOS VALORES DE RLX Y F,LEE P, Y ACTUALIZA KCOD + DO J=1,NGRID + RLX(J)=(2.*PI*RT*COS(RLT(J)*PI/180.))/(360./RLON) + F(J)=2.*ROMEGA*SIN(RLT(J)*PI/180.) + ENDDO + RLY=2.*PI*RT*ABS(RLAT)/360. + K0=0 +! +! print*,"antes de dobla" +! print*,"RLT=",RLT(1:50) +! print*,"RLN",RLN(1:50) +! +! TRANFORMA ALTURA GEOPOTENCIAL EN GEOPOTENCIAL PHI=Z*g Y CALCULA LOS VALORES +! EN DOBLE RESOLUCION + + DO IG=1,NGRID + GG(IG)=Z(IG)*GR + ENDDO + +! print*,"antes de DOBLA" + CALL DOBLA(SLAT,SLON,RLAT,RLON,NLAT,NLON,GG,GD) +! print*,"despues de DOBLA" +! +! print*,"despues de dobla" +! print*,"RLT=",RLT(1:50) +! print*,"RLN",RLN(1:50) +! +! CALCULO DEL VIENTO GEOSTROFICO + JS=0 + DO 17 J=1,NGRID + +! SI NO PERTENECE A LA VENTANA DE SALIDA, SALTA + IF(RLT(J).GT.SLATS.OR.RLT(J).LT.SLATS+((NLATS-1)*RLATS))cycle +! print*,"entro en el primer IF",j,RLT(J),SLATS,SLATS+((NLATS-1)& +! *RLATS) +! print*,"entro en el primer IF, j, RLT =",j,RLT(1:50) +! stop + +! cycle +! endif + IF(RLN(J).LT.SLONS.OR.RLN(J).GT.SLONS+((NLONS-1)*RLONS))cycle +! print*,"entro en el segundo IF",j +! cycle +! endif + + JS=JS+1 +! print*,"JS=",JS + + JD=((MOD(J-1,NLON)+1)*2)-1+(((((J-1)/NLON)*2)+1-1)*NLOND) +! POS EN DOBLE= PTO EN ESA LAT + NUM LATD PASADAS* NLOND + U(JS)=-(GD(JD-NLOND)-GD(JD+NLOND))/(RLY*F(J)) + V(JS)=(GD(JD+1)-GD(JD-1))/(RLX(J)*F(J)) + RLTS(JS)=RLT(J) + RLNS(JS)=RLN(J) + 17 ENDDO +! +! +! + RETURN + END SUBROUTINE + + SUBROUTINE DOBLA(SLAT,SLON,RLAT,RLON,NLAT,NLON,GG,GD) +! +! SLAT, LATITUD DEL LIMITE MERIDIONAL DEL CAMPO DE ENTRADA +! SLATS, LATITUD DEL LIMITE MERIDIONAL DEL CAMPO DE SALIDA +! SLON, LONGITUD DEL LIMITE OCCIDENTAL DEL CAMPO DE ENTRADA +! SLONS, LONGITUD DEL LIMITE OCCIDENTAL DEL CAMPO DE SALIDA + IMPLICIT INTEGER (K) + INTEGER,INTENT(IN) :: NLON,NLAT +! REAL,INTENT(IN) :: SLAT,SLON,RLAT,RLON + DOUBLE PRECISION,INTENT(IN) :: SLAT,SLON,RLAT,RLON + REAL,INTENT(IN) :: GG(NLON*NLAT) + REAL,INTENT(OUT) :: GD(((NLON*2)-1)*((NLAT*2)-1)) + INTEGER IGA(((NLON*2)-1)*((NLAT*2)-1)) + REAL A(NLON,NLAT),S(((NLON-2)*2)-1,((NLAT-2)*2)-1) +! NLON2 y NLAT2 son nlont y nlatt que son NLON y NLAT de +! GEOSTROFICO que llama a DOBLA + NLON2=NLON + NLAT2=NLAT +! +! IF(NLAT.NE.NLAT2 .OR. NLON.NE.NLON2)THEN +! PRINT*,' CAMBIAR NLAT2 Y NLON2 EN SUBRUTINA DOBLA',& +! NLAT,NLAT2,NLON,NLON2 +! STOP +! ENDIF +! + XLATMIH=SLAT+(RLAT*(NLAT-1)) + XLONMIH=SLON + XLATMIR=XLATMIH+ABS(RLAT) + XLONMIR=XLONMIH+RLON + DLATH=ABS(RLAT) + DLONH=ABS(RLON) + DLATR=ABS(RLAT/2.) + DLONR=ABS(RLON/2.) + NLONH=NLON + NLATH=NLAT + NLONR=((NLON-2)*2)-1 + NLATR=((NLAT-2)*2)-1 + ICA=0 + DO J=NLAT,1,-1 + DO I=1,NLON + ICA=ICA+1 + A(I,J)=GG(ICA) + ENDDO + ENDDO + CALL BESSEL(XLATMIH,XLONMIH,XLATMIR,XLONMIR,DLATH,DLONH,& + DLATR,DLONR,NLONH,NLATH,NLONR,NLATR,A,S) +! + NLONS=(NLON*2)-1 + NLATS=(NLAT*2)-1 + IGA=0 + ICA=0 + DO 10 IG=1,NLONS*NLATS + IF (MOD(((IG-1)/NLONS),2).EQ.1) cycle !LATITUDES PARES + IF (MOD(MOD((IG-1),NLONS)+1,2).EQ.0) cycle !LONGITUDES PARES + ICA=ICA+1 + IGA(IG)=ICA + 10 ENDDO +! ESCRIBE LOS PUNTOS DE LA REJILLA ORIGINAL + DO IG=1,NLONS*NLATS + IF(IGA(IG).NE.0) GD(IG)=GG(IGA(IG)) + ENDDO +! + DO IG=1,NLONS*NLATS +! SOBREESCRIBE INTERPOLACIONES HECHAS POR BESSEL + IF(IG.GT.NLONS*2 .AND. IG.LT.NLONS*(NLATS-2) .AND.& + MOD((IG-1),NLONS)+1.GT.2 .AND. MOD((IG-1),NLONS)+1.LT.NLONS-2)& + THEN !INT POR BESSEL, TODAS MENOS LAS DOS PRIMERAS Y ULTIMAS +! LATITUDES Y LONGITUDES + ILAT=NLATR-((((IG-1)/NLONS)-1)-1) + ILON=(MOD(IG-1,NLONS)+1)-2 + GD(IG)=S(ILON,ILAT) +! INTERPOLA PARA LA FRONTERA DE LA REJILLA + ELSE + IF(IGA(IG).EQ.0)THEN + IF(MOD(((IG-1)/NLONS),2).EQ.1 .AND.& + MOD(MOD((IG-1),NLONS)+1,2).EQ.0)THEN !FILA PAR, COLUMNA PAR, INT +!4 PTOS + GD(IG)=(GD(IG-NLONS-1)+GD(IG-NLONS+1)+GD(IG+NLONS-1)+& + GD(IG+NLONS+1))/4. + ELSEIF(MOD(((IG-1)/NLONS),2).EQ.1)THEN !FILA PAR, COLUMNA IMPAR, +!INT 2 PTOS + GD(IG)=(GD(IG-NLONS)+GD(IG+NLONS))/2. + ELSEIF(MOD(MOD((IG-1),NLONS)+1,2).EQ.0)THEN !FILA IMPAR, COLUMNA +!PAR, INT 2 PTOS + GD(IG)=(GD(IG-1)+GD(IG+1))/2. + ENDIF + ENDIF + ENDIF +! + ENDDO +! + RETURN + END SUBROUTINE +! + SUBROUTINE BESSEL(XLATMIH,XLONMIH,XLATMIR,XLONMIR,DLATH,DLONH,& + DLATR,DLONR,NLONH,NLATH,NLONR,NLATR,A,E) +! + INTEGER,INTENT(IN) :: NLONH,NLATH,NLONR,NLATR + REAL,INTENT(IN) :: DLONH,DLONR,DLATH,DLATR,XLATMIH,XLATMIR,XLONMIH,XLONMIR + REAL,INTENT(IN) :: A(NLONH,NLATH) + REAL,INTENT(OUT) :: E(NLONR,NLATR) +! +! COMPRUEBA QUE LOS LIMITES SON CORRECTOS + XLATMAH=XLATMIH+((NLATH-1)*DLATH) + XLONMAH=XLONMIH+((NLONH-1)*DLONH) + XLATMAR=XLATMIR+((NLATR-1)*DLATR) + XLONMAR=XLONMIR+((NLONR-1)*DLONR) +! IF(XLATMIR.LT.XLATMIH+DLATH .OR. XLONMIR.LT.XLONMIH+DLONH& +! .OR. XLATMAR.GT.XLATMAH-DLATH .OR. XLONMAR.GT.XLONMAH-DLONH)THEN +! PRINT*,' ERROR EN LIMITES DE REJILLA ESTIMADA:SLATE,ELATE,SLATS,E& +! LATS, Y LON RESPECTIVOS:',XLATMIH,XLATMAH,XLATMIR,XLATMAR,& +! XLONMIH,XLONMAH,XLONMIR,XLONMAR +! STOP +! ENDIF +! +! HAZ LA INTERPOLACION PARA CADA PUNTO DE LA REJILLA DE SALIDA +3 DO J=1,NLATR + DO I = 1,NLONR +! DETERMINA LA POSICION DEL PUNTO DE LA REJILLA DE SALIDA EN LAS COORDENADAS DE +! LA REJILLA DE ENTRADA + XX= (((XLONMIR+DLONR*(I-1)) - XLONMIH) / DLONH ) +1. + YY= (((XLATMIR+DLATR*(J-1)) - XLATMIH) / DLATH ) +1. + M = XX + N = YY + DX = XX - M + DY = YY - N +! APLICA EL ESQUEMA DE INTERPOLACI\324N DE 16 PT DE BESSEL + DXX = .25 *(DX - 1.) + DYY = .25 *(DY - 1.) + AA = A(M,N-1) + DX *(A(M+1,N-1) - A(M,N-1) + DXX *& + (A(MIN(M+2,NLONH),N-1) - A(M+1,N-1) + A(M-1,N-1) - A(M,N-1))) + AB = A(M,N) + DX*(A(M+1,N) - A(M,N) + DXX *(A(MIN(M+2,NLONH),N)& + - A(M+1,N) + A(M-1,N) - A(M,N))) + AC = A(M,N+1) + DX *(A(M+1,N+1) - A(M,N+1) + DXX *& + (A(MIN(M+2,NLONH),N+1) - A(M+1,N+1) + A(M-1,N+1) - A(M,N+1))) + AD = A(M,MIN(N+2,NLATH)) + DX *(A(M+1,MIN(N+2,NLATH)) -& + A(M,MIN(N+2,NLATH)) + DXX *(A(MIN(M+2,NLONH),MIN(N+2,NLATH))& + - A(M+1,MIN(N+2,NLATH)) + A(M-1,MIN(N+2,NLATH)) -& + A(M,MIN(N+2,NLATH)))) + E(I,J) = AB + DY *(AC - AB + DYY *(AD - AC + AA - AB)) + ENDDO + ENDDO + RETURN + END SUBROUTINE + + subroutine radian(t1,t2,t3,sol) + implicit real(a-h,o-z) +! implicit none +! double precision (a-h,o-z) + pi=3.14159265358979d0 + sol=t1+t2/60.d0+t3/3600.d0 + return + end subroutine + + SUBROUTINE GEOUTM (FLON, FLAT, HUSO, X, Y) +! IMPLICIT REAL (A-Z) + IMPLICIT DOUBLE PRECISION (A-Z) +! IMPLICIT NONE +! DOUBLE PRECISION (A-Z) + PI = 3.14159265 + RG = 180. / PI + E2 = 0.6722670E-02 + EP2 = 0.6768170E-02 + A0 = 0.998317208055891 + A2 = 5.050503255106305E-03 + A4 = 5.323041134969273E-06 + A6 = 6.981680670962105E-09 + A8 = 9.931708438892222E-12 + A10 = 1.44222427482031E-14 + RA = 6378388.0 + XM = (6. * HUSO - 183.) / RG + LOI = FLON / RG - XM + LAT = FLAT / RG + B = RA*(A0 * LAT - 0.5 * (A2 * SIN(2. * LAT) - A4 * SIN(4. * LAT)& + + A6*SIN(6. * LAT) - A8 * SIN(8. * LAT) + A10 * SIN(10. * LAT))) + PSI = LOI * COS(LAT) + W = SQRT(1. - E2 * (SIN(LAT) ** 2)) + HN = RA / W + V2 = 1. + EP2 * (COS(LAT) ** 2) + TF2 = TAN(LAT) ** 2 + C2 = (V2 - TF2) / 6 + C3 = V2 / 24. + V2 ** 2 / 6. - TF2 / 24. + C4 = (V2 * (14. - 58. * TF2) + 40. * TF2 + TF2 ** 2 - 9.) / 120. + C5 = (61.- 58.*TF2 + TF2**2 + (V2 - 1.)*(270. - 330.*TF2)) / 720. + X = 500000. + HN * PSI * (1. + C2 * PSI**2 + C4 * PSI**4)*0.9996 + Y = (B + HN * TAN(LAT) * (0.5 * PSI ** 2 + C3 * PSI ** 4 +& + C5 * PSI ** 6)) * 0.9996 + RETURN + END SUBROUTINE + + SUBROUTINE FECHANNO(DIA,MES,IDA) + INTEGER,INTENT(IN) :: DIA,MES + INTEGER,INTENT(OUT) :: IDA + INTEGER NORMAL(12) + DATA NORMAL/0,31,59,90,120,151,181,212,243,273,304,334/ + IDA=NORMAL(MES)+DIA + IF(MES.EQ.2 .AND. DIA.GT.28)IDA=60 + RETURN + END SUBROUTINE + + +END MODULE MOD_FUNCS + diff --git a/src/predictores_significativos.f90 b/src/predictores_significativos.f90 new file mode 100755 index 0000000000000000000000000000000000000000..9dbe189304a47c55c67c371ded4dbdfc012efede --- /dev/null +++ b/src/predictores_significativos.f90 @@ -0,0 +1,582 @@ + +! sig_predic program selects significance predictor from +! the finded collection +SUBROUTINE sig_predic(nlat,nlon,nlatt,nlont,slat,slon,rlat,rlon,slatt,& + slont,n,ic,id,prec_hr,nger,um,vm,gu92,gv92,gu52,& + gv52,iri,u500,v500,msl_si,q700,t500,t850,& + nanx,neni,new_mi,new_ccm,new_kvars,new_corrpar) + +USE MOD_CSTS +USE MOD_FUNCS + + IMPLICIT NONE + +! 0.1 Declarations of arguments +! ------------------------- + + INTEGER, INTENT(IN) :: nlat + INTEGER, INTENT(IN) :: nlon + INTEGER, INTENT(IN) :: nlatt + INTEGER, INTENT(IN) :: nlont + DOUBLE PRECISION, INTENT(IN) :: slat + DOUBLE PRECISION, INTENT(IN) :: slon + DOUBLE PRECISION, INTENT(IN) :: rlat + DOUBLE PRECISION, INTENT(IN) :: rlon + DOUBLE PRECISION, INTENT(IN) :: slatt + DOUBLE PRECISION, INTENT(IN) :: slont + INTEGER, INTENT(IN) :: n + INTEGER, INTENT(IN) :: ic + INTEGER, INTENT(IN) :: id + DOUBLE PRECISION, INTENT(IN) :: prec_hr(n,nptos) + INTEGER, INTENT(IN) :: nger + DOUBLE PRECISION, INTENT(IN) :: um(n,ic) + DOUBLE PRECISION, INTENT(IN) :: vm(n,ic) + DOUBLE PRECISION, INTENT(IN) :: gu92(n,ic) + DOUBLE PRECISION, INTENT(IN) :: gv92(n,ic) + DOUBLE PRECISION, INTENT(IN) :: gu52(n,ic) + DOUBLE PRECISION, INTENT(IN) :: gv52(n,ic) + INTEGER, INTENT(IN) :: iri(nptos) + DOUBLE PRECISION, INTENT(IN) :: u500(n,ic) + DOUBLE PRECISION, INTENT(IN) :: v500(n,ic) + DOUBLE PRECISION, INTENT(IN) :: msl_si(n,ic) + DOUBLE PRECISION, INTENT(IN) :: q700(n,ic) + DOUBLE PRECISION, INTENT(IN) :: t500(n,ic) + DOUBLE PRECISION, INTENT(IN) :: t850(n,ic) + INTEGER, INTENT(IN) :: nanx + INTEGER, INTENT(IN) :: neni + + INTEGER, INTENT(OUT) :: new_mi(nger,nptos) + DOUBLE PRECISION, INTENT(OUT) :: new_ccm(nger,nptos) + INTEGER, INTENT(OUT) :: new_kvars(nger,nptos,npx) + DOUBLE PRECISION, INTENT(OUT) :: new_corrpar(nger,nptos,npx) + + integer nvar,m + + integer nulon,nulat,nulev,nudays,ideb,ifin,ip + integer i,j,tt,vv + integer is + +!***************************************************************** + integer mi + real ccm + character mdl*20,sc*8,pt*9,nomeb*90,nomef*90,ta*3,nta*1 + real he7(n,ic),he7m(ic) + double precision u9(n,ic),v9(n,ic),u5(n,ic),v5(n,ic) + real psl(n,ic),ut9(nger,ic),vt9(nger,ic),ut5(nger,ic),vt5(nger,ic),pseal(id) + real presor(n,id) + real xlat(ic),xlon(ic) + real t5(n,ic),t8(n,ic),tm5(ic),tm8(ic) + real pslm(ic),um9(ic),vm9(ic),um5(ic),vm5(ic),pslma(ic) + real ue5(id),ve5(id),he7ms(ic),he7mr(id) +! + character (len=6) :: he7ca(id) + character (len=6) :: t8ca(id) + real te8(id),te5(id) + real he7me(24) + real pres(id),bar(id),den(ic) + real pred1(npx,n,neni),pred1m(npx,neni),predh(n,neni),predhm(neni) +! + integer anai(nanx),ana(nanx) + integer kvars(npx) + integer nor(nanx) + integer indi1(ic),indi2(ic) + integer annor(n),mesr(n),diar(n) + integer ior(n),anno,mes,dia,eqc(nptos) + integer ref(nptos),puce(neni),puen(neni,5001) + + integer i1,i2,i3,i4,iana,ice,ien,ipos,ipu,ir,iv,jk,k,nan,ndcp + integer nen,rlx,rly,vorm,vorz + real prec(n,nptos),dis(n) + real p9(ic),p5(ic) + real dato1(npx,nanx),pr(nanx) + real corrpar(npx) + real coe(npx),con + real rlt(ic),rln(ic),rltt(id),rlnt(id) + real dist(nanx),dist1(npx,nanx),serin(nanx) + real aaa(nanx) + real ser(n),media(npx,neni),sigma(npx,neni) + real md,sg,medh(neni),sigh(neni) + real mu9(ic),su9(ic),mv9(ic),sv9(ic) + real mu5(ic),su5(ic),mv5(ic),sv5(ic) + real disu5,disu9,disv5,disv9 +!******************************************************* + +! print*,"program 7: significant predictors" + + nvar=npx + m=nger + +!********************************* +! 1. Sinoptic latitude and longitude calculation and assignment of +! weights windows +! + do j=1,ic + rlt(j)=slat+(((j-1)/nlon)*rlat) + rln(j)=slon+((mod(j-1,nlon)+1-1)*rlon) + enddo + p9=0. + p5=1. + do i1=1,ic + if((rlt(i1).le.fnor2).and.(rlt(i1).ge.fsur2)) then + if((rln(i1).ge.foes2).and.(rln(i1).le.fest2)) then + p9(i1)=1. + p5(i1)=4. + endif + endif + enddo + do i1=1,ic + if((rlt(i1).le.fnor1).and.(rlt(i1).ge.fsur1)) then + if((rln(i1).ge.foes1).and.(rln(i1).le.fest1)) then + p9(i1)=2. + p5(i1)=8. + endif + endif + enddo +! +! Latitude and longitude calculation in the extended domain (called low +! resolution) + + do j=1,id + rltt(j)=slatt+(((j-1)/nlont)*rlat) + rlnt(j)=slont+((mod(j-1,nlont)+1-1)*rlon) + enddo + +!*********************************** +! REANALYSIS VARIABLES + + u5(:,:)=u500(:,:) + v5(:,:)=v500(:,:) + psl(:,:)=msl_si(:,:) + he7(:,:)=q700(:,:) + t5(:,:)=t500(:,:) + t8(:,:)=t850(:,:) + +! HIGH RESOLUTION (5KM) OBSERVATIONS +! It is neccesary to convert to tenths of mm (multiplying by 10). + + prec(:,:)=prec_hr(:,:)*10. + +! Mean and standard deviation of reference synoptic fields. + + do j=1,ic + do i=1,n + ser(i)=um(i,j) + enddo + call estadis(ser,md,sg,n) + mu9(j)=md + su9(j)=sg + do i=1,n + ser(i)=vm(i,j) + enddo + call estadis(ser,md,sg,n) + mv9(j)=md + sv9(j)=sg + do i=1,n + ser(i)=u5(i,j) + enddo + call estadis(ser,md,sg,n) + mu5(j)=md + su5(j)=sg + do i=1,n + ser(i)=v5(i,j) + enddo + call estadis(ser,md,sg,n) + mv5(j)=md + sv5(j)=sg + enddo + +! A reference centers (matching points between sinoptic and high +! resolution grids) are define to know where the predictor must be +! calculated. + + nen=1 + ref=iri + + puce(1)=ref(1) + do 101 j=1,nptos + do k=1,nen + if(ref(j).eq.puce(k)) go to 101 + enddo + nen=nen+1 + ipos=nen + puce(ipos)=ref(j) + 101 continue + +! Each reference point have associated a group of high resolution grids. + puen=0 + do k=1,nen + do j=1,nptos + if(ref(j).eq.puce(k)) then + puen(k,5001)=puen(k,5001)+1 + ipos=puen(k,5001) + puen(k,ipos)=j + endif + enddo + enddo + +! The predictors are obtained and normalized + +! OBTAINING THE SEA LEVEL PRESSURE (PREDICTOR 1) IN THE REFERENCE CENTERS + do i=1,n + do j=1,nen + ice=puce(j) + pred1(1,i,j)=psl(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(1,i,j) + enddo + call estadis(ser,md,sg,n) + media(1,j)=md + sigma(1,j)=sg + do i=1,n + pred1(1,i,j)=(pred1(1,i,j)-media(1,j))/sigma(1,j) + enddo + enddo + +! OBTAINING THE TREND (PREDICTOR 11) IN THE REFERENCE CENTERS + + do j=1,nen + pred1(11,1,j)=0. + enddo + + do i=2,n + do j=1,nen + ice=puce(j) + pred1(11,i,j)=psl(i,ice)-psl((i-1),ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(11,i,j) + enddo + call estadis(ser,md,sg,n) + media(11,j)=md + sigma(11,j)=sg + do i=1,n + pred1(11,i,j)=(pred1(11,i,j)-media(11,j))/sigma(11,j) + enddo + enddo + + +! OBTAINING THE VERTICAL THERMAL GRADIENT(PREDICTOR 3) +! IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(3,i,j)=t8(i,ice)-t5(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(3,i,j) + enddo + call estadis(ser,md,sg,n) + media(3,j)=md + sigma(3,j)=sg + do i=1,n + pred1(3,i,j)=(pred1(3,i,j)-media(3,j))/sigma(3,j) + enddo + enddo + +! OBTAINING THE 500 hPa TEMPERATURE (PREDICTOR 2) +! IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(2,i,j)=t5(i,ice) + enddo + enddo + + do j=1,nen + do i=1,n + ser(i)=pred1(2,i,j) + enddo + call estadis(ser,md,sg,n) + media(2,j)=md + sigma(2,j)=sg + do i=1,n + pred1(2,i,j)=(pred1(2,i,j)-media(2,j))/sigma(2,j) + enddo + enddo + +! OBTAINING THE VORTICITY (PREDICTOR 4) IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + rlx=rt*cos(rlt(ice)*pi/180.)*pi*rlon/180. + rly=rt*abs(rlat)*pi/180. + vorm=um(i,ice-nlon)-um(i,ice+nlon) + vorm=vorm/(2.*rly) + vorz=vm(i,ice+1)-vm(i,ice-1) + vorz=vorz/(2.*rlx) + pred1(4,i,j)=vorz-vorm + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(4,i,j) + enddo + call estadis(ser,md,sg,n) + media(4,j)=md + sigma(4,j)=sg + do i=1,n + pred1(4,i,j)=(pred1(4,i,j)-media(4,j))/sigma(4,j) + enddo + enddo + +! OBTAINING THE GEOSTROPHIC U/V COMPONENTS (PREDICTORS 5 AND 6) IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(5,i,j)=um(i,ice) + pred1(6,i,j)=vm(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(5,i,j) + enddo + call estadis(ser,md,sg,n) + media(5,j)=md + sigma(5,j)=sg + do i=1,n + pred1(5,i,j)=(pred1(5,i,j)-media(5,j))/sigma(5,j) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(6,i,j) + enddo + call estadis(ser,md,sg,n) + media(6,j)=md + sigma(6,j)=sg + do i=1,n + pred1(6,i,j)=(pred1(6,i,j)-media(6,j))/sigma(6,j) + enddo + enddo + +! OBTAINING THE VORTICITY IN 500 hPa (PREDICTOR 7) IN THE REFERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + rlx=rt*cos(rlt(ice)*pi/180.)*pi*rlon/180. + rly=rt*abs(rlat)*pi/180. + vorm=u5(i,ice-nlon)-u5(i,ice+nlon) + vorm=vorm/(2.*rly) + vorz=v5(i,ice+1)-v5(i,ice-1) + vorz=vorz/(2.*rlx) + pred1(7,i,j)=vorz-vorm + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(7,i,j) + enddo + call estadis(ser,md,sg,n) + media(7,j)=md + sigma(7,j)=sg + do i=1,n + pred1(7,i,j)=(pred1(7,i,j)-media(7,j))/sigma(7,j) + enddo + enddo + + +! OBTAINING THE GEOSTROPHIC U/V COMPONENTS IN 500 hPa (PREDICTORS 8 AND 9) +! IN THE RERENCE CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(8,i,j)=u5(i,ice) + pred1(9,i,j)=v5(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(8,i,j) + enddo + call estadis(ser,md,sg,n) + media(8,j)=md + sigma(8,j)=sg + do i=1,n + pred1(8,i,j)=(pred1(8,i,j)-media(8,j))/sigma(8,j) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(9,i,j) + enddo + call estadis(ser,md,sg,n) + media(9,j)=md + sigma(9,j)=sg + do i=1,n + pred1(9,i,j)=(pred1(9,i,j)-media(9,j))/sigma(9,j) + enddo + enddo + +! OBTAINING THE ESPECIFIC HUMIDITY IN 700 hPa (PREDICTOR 10) IN THE REFERENCE +! CENTERS + + do i=1,n + do j=1,nen + ice=puce(j) + pred1(10,i,j)=he7(i,ice) + enddo + enddo + do j=1,nen + do i=1,n + ser(i)=pred1(10,i,j) + enddo + call estadis(ser,md,sg,n) + media(10,j)=md + sigma(10,j)=sg + do i=1,n + pred1(10,i,j)=(pred1(10,i,j)-media(10,j))/sigma(10,j) + enddo + enddo +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! ESTANDARIZATION OF REFERENCE WINDS (SINOPTIC WINDS ALSO) + + do i=1,n + do j=1,ic + u9(i,j)=(um(i,j)-mu9(j))/su9(j) + v9(i,j)=(vm(i,j)-mv9(j))/sv9(j) + u5(i,j)=(u5(i,j)-mu5(j))/su5(j) + v5(i,j)=(v5(i,j)-mv5(j))/sv5(j) + enddo + enddo + + do i=1,m + do j=1,ic + ut9(i,j)=(gu92(i,j)-mu9(j))/su9(j) + vt9(i,j)=(gv92(i,j)-mv9(j))/sv9(j) + ut5(i,j)=(gu52(i,j)-mu5(j))/su5(j) + vt5(i,j)=(gv52(i,j)-mv5(j))/sv5(j) + enddo + enddo + +! OBTAINING SIGNIFICANT PREDICTORS FOR EACH SINOPTIC TYPE IN EACH HIGH +! RESOLUTION GRID POINT. + + do 1000 i=1,m + +! print*,i + + +! Determine the "nanx" reference alements more similar to each sinoptic type +! and the corresponding distances + + do i1=1,n + ior(i1)=i1 + dis(i1)=9999. + enddo + do 113 i1=1,n + call distancia9(ut9,m,u9,n,i,i1,p9,disu9,ic) + call distancia9(vt9,m,v9,n,i,i1,p9,disv9,ic) + call distancia5(ut5,m,u5,n,i,i1,p5,disu5,ic) + call distancia5(vt5,m,v5,n,i,i1,p5,disv5,ic) + dis(i1)=(disu9+disv9+disu5+disv5)/4. + + 113 continue + + call burbuja1(dis,ior,n,nanx) + + do i1=1,nanx + anai(i1)=ior(i1) + enddo + +! Consider all high resolution grid points associated with the low resolution +! ones + + do 1100 ien=1,nen + do 1200 i2=1,puen(ien,5001) + ipu=puen(ien,i2) + +! Consider predictand values (precipitation) and predictors from the analogs + + nan=0 + ndcp=0 + do i3=1,nanx + iana=anai(i3) + + if(prec(iana,ipu).ne.-999.) then + nan=nan+1 + ana(nan)=iana + do i4=1,nvar + dato1(i4,nan)=pred1(i4,iana,ien) + enddo + pr(nan)=prec(iana,ipu) + if(pr(nan).eq.-3.) pr(nan)=1. + if(pr(nan).ge.1.) ndcp=ndcp+1 + endif + enddo + if(nan.le.30) then + mi=0 + ccm=-7.77 + go to 1199 + endif + + if(ndcp.le.30) then + mi=0 + ccm=-9.99 + go to 1199 + endif + + if(nan.gt.150) nan=150 + +! Calculation of significant predictors, their coeficients and their +! multiple and partial correlation coeficients to estimate the +! precipitation +! +! mi: number of selected predictors +! ccm: multiple correlation coeficient +! kvars: selected predictors labels (vector) +! corrpar: partial correlation of selected predictors (vector) +! coe: regression coeficients associated to each predictor (vector). +! (value = 0 when there is no selected predictor). +! con: Y-intercept (independent equation term) +! tol: tolerance to select predictors + + call stepregrs& + (pr,dato1,nanx,nvar,nan,mi,ccm,kvars,corrpar,coe,con,tol) + + 1199 continue + + new_mi(i,ipu)=0 + new_ccm(i,ipu)=0. + new_kvars(i,ipu,:)=0 + new_corrpar(i,ipu,:)=0. + + new_mi(i,ipu)=mi + new_ccm(i,ipu)=ccm + + if (mi.ne.0) then + vv=mi + do tt=1,mi + new_kvars(i,ipu,tt)=kvars(tt) + new_corrpar(i,ipu,tt)=corrpar(kvars(tt)) + end do + else + vv=1 + tt=vv + new_kvars(i,ipu,tt)=0 + new_corrpar(i,ipu,tt)=0. + end if + + + 1200 continue + + 1100 continue + + 1000 continue + +END SUBROUTINE sig_predic + diff --git a/src/pts_ref_est_pen.f90 b/src/pts_ref_est_pen.f90 new file mode 100755 index 0000000000000000000000000000000000000000..39f504c5b59d417e3fb2d7b2d6e3e9099b2b6a81 --- /dev/null +++ b/src/pts_ref_est_pen.f90 @@ -0,0 +1,59 @@ + +! The ptos_ref program links Reanalysis grid and +! observed grid using the nearest neighbor + +SUBROUTINE ptos_ref(ic,x,y,xcand,ycand,iri,ipos) + +USE MOD_CSTS + + Implicit none + + INTEGER, INTENT(IN) :: ic + REAL, INTENT(IN) :: x(ic) + REAL, INTENT(IN) :: y(ic) + REAL, INTENT(IN) :: xcand(nptos) + REAL, INTENT(IN) :: ycand(nptos) + INTEGER, INTENT(OUT) :: iri(nptos) + INTEGER, INTENT(OUT) :: ipos + + integer np,i,j,k + real xe(nptos),ye(nptos),xr(ic),yr(ic),dis,dmin(nptos) + integer valores_unicos(nptos+1) + + np=nptos + +! print*,"program 6: reference points" + + xr=x/1000. + yr=y/1000. + + xe=xcand/1000. + ye=ycand/1000. + + valores_unicos=0 + + do 100 i=1,np + + dmin(i)=1600000000. + + do 110 j=1,ic + dis=(xe(i)-xr(j))**2+(ye(i)-yr(j))**2 + if(dis.lt.dmin(i)) then + dmin(i)=dis + iri(i)=j + endif + 110 continue + + dmin(i)=sqrt(dmin(i)) + + do k=1,valores_unicos(np+1) + if(valores_unicos(k).eq.iri(i)) go to 100 + enddo + valores_unicos(np+1)=valores_unicos(np+1)+1 + ipos=valores_unicos(np+1) + valores_unicos(ipos)=iri(i) + + 100 continue + +END SUBROUTINE ptos_ref + diff --git a/src/pts_ref_est_pen_4int.f90 b/src/pts_ref_est_pen_4int.f90 new file mode 100755 index 0000000000000000000000000000000000000000..0e44c895704255abd2fd6d31718bb89435d79efe --- /dev/null +++ b/src/pts_ref_est_pen_4int.f90 @@ -0,0 +1,84 @@ + +! The ptos_ref_4 program links Reanalysis grid and +! observed grid using 4 nearest points interpolation +! (bilineal interpolation approach) + + +SUBROUTINE ptos_ref_4(ic,x,y,xcand,ycand,Vdmin,Vref,ipos) + +USE MOD_CSTS + + Implicit none + + INTEGER, INTENT(IN) :: ic + REAL, INTENT(IN) :: x(ic) + REAL, INTENT(IN) :: y(ic) + REAL, INTENT(IN) :: xcand(nptos) + REAL, INTENT(IN) :: ycand(nptos) + DOUBLE PRECISION, INTENT(OUT) :: Vdmin(nptos,4) + INTEGER, INTENT(OUT) :: Vref(nptos,4) + INTEGER, INTENT(OUT) :: ipos + + integer iri + real dmin + integer np,i,j,k + real xe(nptos),ye(nptos),xr(ic),yr(ic),dis + real copiaXr(nptos),copiaYr(nptos) + integer valores_unicos(nptos+1) + integer cont + +! print*,"program 6: 4 nearest points of reference" + + np=nptos + + xr=x/1000. + yr=y/1000. + + xe=xcand/1000. + ye=ycand/1000. + + + valores_unicos=0 + + do i=1,np + + copiaXr=xr + copiaYr=yr + + do cont=1,4 !4 nearest pts loop + + dmin=1600000000. + + do 110 j=1,ic + dis=(xe(i)-copiaXr(j))**2+(ye(i)-copiaYr(j))**2 + if(dis.lt.dmin) then + dmin=dis + iri=j + endif + 110 continue + + Vdmin(i,cont)=sqrt(dmin) + if (Vdmin(i,cont) .lt. 0.1) then + Vdmin(i,cont)=0.1 + endif + + Vref(i,cont)=iri + copiaXr(iri)=99999999. + copiaYr(iri)=99999999. + + + do k=1,valores_unicos(np+1) + if(valores_unicos(k).eq.iri) go to 100 + enddo + + valores_unicos(np+1)=valores_unicos(np+1)+1 + ipos=valores_unicos(np+1) + valores_unicos(ipos)=iri + + 100 continue + enddo !4 rearest pts loop + + enddo + +END SUBROUTINE ptos_ref_4 + diff --git a/src/registerDynamicSymbol.c b/src/registerDynamicSymbol.c new file mode 100644 index 0000000000000000000000000000000000000000..ad203591419a8e391ea22bd7b2a496f6f94c0052 --- /dev/null +++ b/src/registerDynamicSymbol.c @@ -0,0 +1,10 @@ +// RegisteringDynamic Symbols + +#include +#include +#include + +void R_init_markovchain(DllInfo* info) { + R_registerRoutines(info, NULL, NULL, NULL, NULL); + R_useDynamicSymbols(info, TRUE); +} diff --git a/src/training_part1_prec.f90 b/src/training_part1_prec.f90 new file mode 100755 index 0000000000000000000000000000000000000000..37dc8f3b2653a7460577c8e1ebbfd44a5e9e89e7 --- /dev/null +++ b/src/training_part1_prec.f90 @@ -0,0 +1,89 @@ +! !!!!!!!!!!!!!!! + subroutine training_part1(u500,v500,t1000,z500,z1000,& + msl_si,msl_lr,ngridd,nlat,nlon,ic,nlatt,nlont,& + id,slat,slon,rlat,rlon,slatt,slont,nd,& + um,vm,insol,gu92,gv92,gu52,gv52,nger) + +! !!!!!!!!!!!!!!! +!* 0. DECLARATIONS +! ------------ +! MODULES with constants and functions +use mod_csts +use mod_funcs + +implicit none + +!!!!!!!!!!!!!!!!!!!!!!!! +! INPUT ARGUMENTS +!!!!!!!!!!!!!!!!!!!!!!!! +integer, intent(in) :: ngridd +!*********************************************** +! DOMAIN variables +!************************************************ +! sinoptic grid +integer, intent(in) :: nlat +integer, intent(in) :: nlon +integer, intent(in) :: ic +!----------------------------------------------- +! low resolution grid +integer, intent(in) :: nlatt +integer, intent(in) :: nlont +integer, intent(in) :: id +!------------------------------------------------ +double precision, intent(in) :: slat +double precision, intent(in) :: slon +double precision, intent(in) :: rlat +double precision, intent(in) :: rlon +!------------------------------------------------ +double precision, intent(in) :: slatt +double precision, intent(in) :: slont +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!! +! TIME variables +integer, intent(in) :: nd +!------------------------------------------------ +! Reanlysis fields +double precision, intent(in) :: u500(nd,ic) +double precision, intent(in) :: v500(nd,ic) +double precision, intent(in) :: t1000(nd,ic) +double precision, intent(in) :: z500(nd,ic) +double precision, intent(in) :: z1000(nd,ic) +double precision, intent(in) :: msl_si(nd,ic) +double precision, intent(in) :: msl_lr(nd,id) +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!! +! OUTPUT ARGUMENTS +!!!!!!!!!!!!!!!!!!!!!!!! +double precision, intent(out) :: um(nd,ic) +double precision, intent(out) :: vm(nd,ic) +double precision, intent(out) :: insol(nd) +!------------------------------------------------ +double precision, intent(out) :: gu92(nd,ic) +double precision, intent(out) :: gv92(nd,ic) +double precision, intent(out) :: gu52(nd,ic) +double precision, intent(out) :: gv52(nd,ic) +!------------------------------------------------ +integer, intent(out) :: nger +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! INNER FORTRAN VARIABLES +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!------------------------------------------------ +real :: den(nd,ic) +!------------------------------------------------ + +!print*,"" +!print*,"*** TRAINING PROCESS ***" +!print*,"" + +call calc_tempes_densi_sealev(ic,nd,msl_si,t1000,den) + +call geos(ic,nd,id,slatt,slont,slat,slon,rlat,& + rlon,rlat,rlon,nlatt,nlont,nlat,nlon,den,msl_lr,ngridd,& + um,vm) + +call clasif(ic,nd,nlon,nlat,slat,slon,rlat,rlon,um,vm,u500,v500,z1000,& + z500,nger,gu92,gv92,gu52,gv52) + +end subroutine + diff --git a/src/training_part2_prec.f90 b/src/training_part2_prec.f90 new file mode 100755 index 0000000000000000000000000000000000000000..733f1d881427b0180895638a58637dbebf552a66 --- /dev/null +++ b/src/training_part2_prec.f90 @@ -0,0 +1,112 @@ +! !!!!!!!!!!!!!!! + subroutine training_part2(u500,v500,t500,t850,msl_si,q700,& + lon_hr,lat_hr,prec_hr,& + nanx,nlat,nlon,ic,nlatt,nlont,id,slat,& + slon,rlat,rlon,slatt,slont,nd,um,vm,gu92,gv92,& + gu52,gv52,nger,Vdmin,Vref,ipos2,new_mi,new_ccm,& + new_kvars,new_corrpar) + +! !!!!!!!!!!!!!!! +!* 0. DECLARATIONS +! ------------ +! MODULES with constants and functions +use mod_csts +use mod_funcs + +implicit none + +!!!!!!!!!!!!!!!!!!!!!!!! +! INPUT ARGUMENTS +!!!!!!!!!!!!!!!!!!!!!!!! +integer, intent(in) :: nanx +!*********************************************** +! DOMAIN variables +!************************************************ +! sinoptic grid +integer, intent(in) :: nlat +integer, intent(in) :: nlon +integer, intent(in) :: ic +!----------------------------------------------- +! low resolution grid +integer, intent(in) :: nlatt +integer, intent(in) :: nlont +integer, intent(in) :: id +!------------------------------------------------ +double precision, intent(in) :: slat +double precision, intent(in) :: slon +double precision, intent(in) :: rlat +double precision, intent(in) :: rlon +!------------------------------------------------ +double precision, intent(in) :: slatt +double precision, intent(in) :: slont +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!! +! TIME variables +integer, intent(in) :: nd +!------------------------------------------------ +! Reanlysis fields +double precision, intent(in) :: u500(nd,ic) +double precision, intent(in) :: v500(nd,ic) +double precision, intent(in) :: t500(nd,ic) +double precision, intent(in) :: t850(nd,ic) +double precision, intent(in) :: msl_si(nd,ic) +double precision, intent(in) :: q700(nd,ic) +!------------------------------------------------ +! AEMET high resolution observational dat +double precision, intent(in) :: lon_hr(nptos) +double precision, intent(in) :: lat_hr(nptos) +double precision, intent(in) :: prec_hr(nd,nptos) +!------------------------------------------------ +double precision, intent(in) :: um(nd,ic) +double precision, intent(in) :: vm(nd,ic) +!------------------------------------------------ +integer, intent(in) :: nger +!------------------------------------------------ +double precision, intent(in) :: gu92(nger,ic) +double precision, intent(in) :: gv92(nger,ic) +double precision, intent(in) :: gu52(nger,ic) +double precision, intent(in) :: gv52(nger,ic) +!------------------------------------------------ +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!! +! OUTPUT ARGUMENTS +!!!!!!!!!!!!!!!!!!!!!!!! +integer, intent(out) :: ipos2 +!!!!!!!!!!!!!!!!!!!!!!!! +double precision, intent(out) :: Vdmin(nptos,4) +integer, intent(out) :: Vref(nptos,4) +!------------------------------------------------ +integer, intent(out) :: new_mi(nger,nptos) +double precision, intent(out) :: new_ccm(nger,nptos) +integer, intent(out) :: new_kvars(nger,nptos,npx) +double precision, intent(out) :: new_corrpar(nger,nptos,npx) +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! INNER FORTRAN VARIABLES +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +integer :: ipos +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +real :: x(ic) +real :: y(ic) +!------------------------------------------------ +real :: xcand(nptos) +real :: ycand(nptos) +!------------------------------------------------ +integer :: iri(nptos) +!------------------------------------------------ + +call utm_ERA(ic,nlat,nlon,slat,slon,rlat,rlon,x,y) + +call utm_obs(lon_hr,lat_hr,xcand,ycand) + +call ptos_ref_4(ic,x,y,xcand,ycand,Vdmin,Vref,ipos2) + +call ptos_ref(ic,x,y,xcand,ycand,iri,ipos) + +call sig_predic(nlat,nlon,nlatt,nlont,slat,slon,rlat,rlon,slatt,& + slont,nd,ic,id,prec_hr,nger,um,vm,gu92,gv92,gu52,& + gv52,iri,u500,v500,msl_si,q700,t500,t850,nanx,& + ipos,new_mi,new_ccm,new_kvars,new_corrpar) + +end subroutine + diff --git a/src/training_temp.f90 b/src/training_temp.f90 new file mode 100755 index 0000000000000000000000000000000000000000..009a5f366e5a9241d71360c758a5deb68b2a032c --- /dev/null +++ b/src/training_temp.f90 @@ -0,0 +1,103 @@ +! !!!!!!!!!!!!!!! + subroutine training_temp(t1000,msl_si,msl_lr,lon_hr,lat_hr,& + ngridd,nlat,nlon,ic,nlatt,nlont,id,slat,& + slon,rlat,rlon,slatt,slont,nd,day,month,& + um,vm,insol,Vdmin,Vref,ipos) + +! !!!!!!!!!!!!!!! +!* 0. DECLARATIONS +! ------------ +! MODULES with constants and functions +use mod_csts +use mod_funcs + +implicit none + +!!!!!!!!!!!!!!!!!!!!!!!! +! INPUT ARGUMENTS +!!!!!!!!!!!!!!!!!!!!!!!! +integer, intent(in) :: ngridd +!*********************************************** +! DOMAIN variables +!************************************************ +! sinoptic grid +integer, intent(in) :: nlat +integer, intent(in) :: nlon +integer, intent(in) :: ic +!----------------------------------------------- +! low resolution grid +integer, intent(in) :: nlatt +integer, intent(in) :: nlont +integer, intent(in) :: id +!------------------------------------------------ +double precision, intent(in) :: slat +double precision, intent(in) :: slon +double precision, intent(in) :: rlat +double precision, intent(in) :: rlon +!------------------------------------------------ +double precision, intent(in) :: slatt +double precision, intent(in) :: slont +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!! +! TIME variables +integer, intent(in) :: nd +integer, intent(in) :: day(nd) +integer, intent(in) :: month(nd) +!------------------------------------------------ +! Reanlysis fields +double precision, intent(in) :: t1000(nd,ic) +double precision, intent(in) :: msl_si(nd,ic) +double precision, intent(in) :: msl_lr(nd,id) +!------------------------------------------------ +! AEMET high resolution observational dat +double precision, intent(in) :: lon_hr(nptos) +double precision, intent(in) :: lat_hr(nptos) +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!! +! OUTPUT ARGUMENTS +!!!!!!!!!!!!!!!!!!!!!!!! +double precision, intent(out) :: um(nd,ic) +double precision, intent(out) :: vm(nd,ic) +double precision, intent(out) :: insol(nd) +!!!!!!!!!!!!!!!!!!!!!!!! +double precision, intent(out) :: Vdmin(nptos,4) +integer, intent(out) :: Vref(nptos,4) +!!!!!!!!!!!!!!!!!!!!!!!! +integer, intent(out) :: ipos + +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! INNER FORTRAN VARIABLES +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!------------------------------------------------ +real :: den(nd,ic) +!------------------------------------------------ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +real :: x(ic) +real :: y(ic) +!------------------------------------------------ +real :: xcand(nptos) +real :: ycand(nptos) +!------------------------------------------------ + + +!print*,"" +!print*,"*** TRAINING PROCESS ***" +!print*,"" + +call calc_tempes_densi_sealev(ic,nd,msl_si,t1000,den) + +call geos(ic,nd,id,slatt,slont,slat,slon,rlat,& + rlon,rlat,rlon,nlatt,nlont,nlat,nlon,den,msl_lr,ngridd,& + um,vm) + +call insolation(nd,day,month,insol) + +call utm_ERA(ic,nlat,nlon,slat,slon,rlat,rlon,x,y) + +call utm_obs(lon_hr,lat_hr,xcand,ycand) + +call ptos_ref_4(ic,x,y,xcand,ycand,Vdmin,Vref,ipos) + +end subroutine +