From e43dcd09fe99208c4e497fe9e8d7c32f42e7b45b Mon Sep 17 00:00:00 2001 From: mrke Date: Sat, 14 May 2022 10:31:20 +1000 Subject: [PATCH] fixing nasty issue emerging with R version 4.2 for endotherm model where failing to put a decimal point after some values was causing them to be treated as integers --- R/endoR_devel.R | 1 + man/ectotherm.Rd | 2 +- man/micro_aust.Rd | 10 +- man/micro_era5.Rd | 10 +- man/micro_global.Rd | 20 ++-- man/micro_ncep.Rd | 10 +- man/micro_nz.Rd | 20 ++-- man/micro_terra.Rd | 20 ++-- man/micro_uk.Rd | 10 +- man/micro_usa.Rd | 10 +- man/onelump.Rd | 14 +-- man/onelump_var.Rd | 2 +- man/twolump.Rd | 2 +- src/CONV_ENDO.f | 22 ++-- src/GEOM_ENDO.f | 126 ++++++++++++----------- src/GETKFUR.f | 2 +- src/IRPROP.f | 22 ++-- src/RESPFUN.f | 10 +- src/SEVAP_ENDO.f | 4 +- src/SIMULSOL.f | 246 ++++++++++++++++++++++---------------------- src/SOLVENDO.f | 69 +++++++++---- src/WETAIR.f | 2 +- 22 files changed, 333 insertions(+), 301 deletions(-) diff --git a/R/endoR_devel.R b/R/endoR_devel.R index 5b0d2fd0..504e22b2 100644 --- a/R/endoR_devel.R +++ b/R/endoR_devel.R @@ -404,6 +404,7 @@ endoR_devel <- function( #PANTSTEP <- 0 # check if heat stressed already (to save computation) QGEN <- 0 + QRESP <- 0 TC_REF <- TC QBASREF <- QBASAL diff --git a/man/ectotherm.Rd b/man/ectotherm.Rd index 3a03bca0..24e3f704 100644 --- a/man/ectotherm.Rd +++ b/man/ectotherm.Rd @@ -552,7 +552,7 @@ ytick<-seq(-20, -10, by=2) axis(side=2, at=ytick, labels = FALSE) mtext(text = rev(seq(0, 100, 20)), side = 2, line = 1, at = seq(-20, -10, 2), las = 2) abline(h = -10, lty = 2, col = 'grey') -mtext(text = c('body temperature (°C)', 'activity', 'shade (\%)', 'depth (cm)'), side = 2, line = 2.5, at = c(30, 9, 0, -15)) +mtext(text = c('body temperature (°C)', 'activity', 'shade (\%)', 'depth (cm)'), side = 2, line = 2.5, at = c(30, 9, 0, -15)) text(-0.2, c(ecto$T_F_max + 1, ecto$T_F_min + 1), c('T_F_max', 'T_F_min'), col = c('red', 'blue'), cex = 0.75) # seasonal activity plot (dark blue = night, light blue = basking, orange = foraging) diff --git a/man/micro_aust.Rd b/man/micro_aust.Rd index 8fc420ab..ca2ff832 100644 --- a/man/micro_aust.Rd +++ b/man/micro_aust.Rd @@ -309,11 +309,11 @@ soilmoist <- cbind(dates, soilmoist) minshade <- micro$minshade[1] # plotting above-ground conditions in minimum shade -with(metout, {plot(TALOC ~ dates, xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout, {plot(TALOC ~ dates, xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l" ,main = paste("air and sky temperature, ", minshade, "\% shade", sep = ""), ylim = c(-20, 60))}) -with(metout, {points(TAREF ~ dates, xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout, {points(TAREF ~ dates, xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l", lty = 2, col = 'blue')}) -with(metout,{points(TSKYC ~ dates, xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{points(TSKYC ~ dates, xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l", col = 'light blue', main = paste("sky temperature, ", minshade, "\% shade", sep = ""))}) with(metout, {plot(RHLOC ~ dates, xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l", ylim = c(0, 100), main = paste("humidity, ", minshade, "\% shade", sep = ""))}) @@ -331,11 +331,11 @@ with(metout, {plot(SNOWDEP ~ dates, xlab = "Date and Time", ylab = "Snow Depth ( # plotting soil temperature for(i in 1:10){ if(i == 1){ - plot(soil[, i + 3] ~ soil[, 1], xlab = "Date and Time", ylab = "Soil Temperature (°C)" + plot(soil[, i + 3] ~ soil[, 1], xlab = "Date and Time", ylab = "Soil Temperature (°C)" , col = i, type = "l", main = paste("soil temperature ", minshade, "\% shade",sep="")) }else{ points(soil[, i + 3] ~ soil[, 1], xlab = "Date and Time", ylab = "Soil Temperature - (°C)", col = i, type = "l") + (°C)", col = i, type = "l") } } diff --git a/man/micro_era5.Rd b/man/micro_era5.Rd index 9f8afd9d..dd32a8a2 100644 --- a/man/micro_era5.Rd +++ b/man/micro_era5.Rd @@ -332,11 +332,11 @@ soil <- cbind(dates,soil) soilmoist <- cbind(dates, soilmoist) # plotting above-ground conditions in minimum shade -with(metout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",main=paste("air and sky temperature",sep=""), ylim = c(-20, 60))}) -with(metout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",lty=2,col='blue')}) -with(metout,{points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",col='light blue',main=paste("sky temperature",sep=""))}) with(metout,{plot(RHLOC ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",ylim=c(0,100),main=paste("humidity",sep=""))}) @@ -354,11 +354,11 @@ with(metout,{plot(SNOWDEP ~ dates,xlab = "Date and Time", ylab = "Snow Depth (cm # plotting soil temperature for(i in 1:10){ if(i==1){ - plot(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)" + plot(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)" ,col=i,type = "l",main=paste("soil temperature",sep="")) }else{ points(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature - (°C)",col=i,type = "l") + (°C)",col=i,type = "l") } } diff --git a/man/micro_global.Rd b/man/micro_global.Rd index bbacfab4..d8f79470 100644 --- a/man/micro_global.Rd +++ b/man/micro_global.Rd @@ -309,15 +309,15 @@ minshade <- micro$minshade[1] maxshade <- micro$maxshade[1] # plotting above-ground conditions in minimum shade -with(metout, {plot(TALOC ~ micro$dates, xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(metout, {plot(TALOC ~ micro$dates, xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l", main = paste("air temperature, ", minshade, "\% shade",sep = ""))}) -with(metout, {points(TAREF ~ micro$dates, xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(metout, {points(TAREF ~ micro$dates, xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l",lty = 2, col = 'blue')}) with(metout, {plot(RHLOC ~ micro$dates, xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",ylim = c(0, 100),main = paste("humidity, ",minshade, "\% shade",sep=""))}) with(metout, {points(RH ~ micro$dates, xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",col = 'blue',lty = 2, ylim = c(0, 100))}) -with(metout, {plot(TSKYC ~ micro$dates, xlab = "Date and Time", ylab = "Sky Temperature (°C)" +with(metout, {plot(TSKYC ~ micro$dates, xlab = "Date and Time", ylab = "Sky Temperature (°C)" , type = "l", main = paste("sky temperature, ", minshade, "\% shade", sep=""))}) with(metout, {plot(VREF ~ micro$dates, xlab = "Date and Time", ylab = "Wind Speed (m/s)" , type = "l", main = "wind speed", col = 'blue',ylim = c(0, 15))}) @@ -331,34 +331,34 @@ with(metout, {plot(SOLR ~ micro$dates,xlab = "Date and Time", ylab = "Solar Radi # plotting soil temperature for minimum shade for(i in 1:10){ if(i==1){ - plot(soil[,i + 2] ~ micro$dates, xlab = "Date and Time", ylab = "Soil Temperature (°C)" + plot(soil[,i + 2] ~ micro$dates, xlab = "Date and Time", ylab = "Soil Temperature (°C)" ,col = i, type = "l", main = paste("soil temperature ", minshade, "\% shade", sep="")) }else{ points(soil[,i + 2] ~ micro$dates, xlab = "Date and Time", ylab = "Soil Temperature - (°C)", col = i, type = "l") + (°C)", col = i, type = "l") } } # plotting above-ground conditions in maximum shade -with(shadmet,{plot(TALOC ~ micro$dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(shadmet,{plot(TALOC ~ micro$dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l", main = "air temperature, sun")}) -with(shadmet,{points(TAREF ~ micro$dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(shadmet,{points(TAREF ~ micro$dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l", lty = 2, col = 'blue')}) with(shadmet,{plot(RHLOC ~ micro$dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l", ylim = c(0, 100),main = "humidity, shade")}) with(shadmet,{points(RH ~ micro$dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l", col = 'blue',lty = 2, ylim = c(0, 100))}) -with(shadmet,{plot(TSKYC ~ micro$dates,xlab = "Date and Time", ylab = "Sky Temperature (°C)", +with(shadmet,{plot(TSKYC ~ micro$dates,xlab = "Date and Time", ylab = "Sky Temperature (°C)", type = "l", main = "sky temperature, shade")}) # plotting soil temperature for maximum shade for(i in 1:10){ if(i==1){ plot(shadsoil[,i + 2] ~ micro$dates, xlab = "Date and Time", ylab = "Soil Temperature - (°C)", col = i, type = "l", main = paste("soil temperature ", maxshade, "\% shade", sep="")) + (°C)", col = i, type = "l", main = paste("soil temperature ", maxshade, "\% shade", sep="")) }else{ points(shadsoil[,i + 2] ~ micro$dates, xlab = "Date and Time", ylab = "Soil Temperature - (°C)", col = i, type = "l") + (°C)", col = i, type = "l") } } } diff --git a/man/micro_ncep.Rd b/man/micro_ncep.Rd index ada67ab7..67299d11 100644 --- a/man/micro_ncep.Rd +++ b/man/micro_ncep.Rd @@ -324,11 +324,11 @@ soil <- cbind(dates,soil) soilmoist <- cbind(dates, soilmoist) # plotting above-ground conditions in minimum shade -with(metout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",main=paste("air and sky temperature",sep=""), ylim = c(-20, 60))}) -with(metout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",lty=2,col='blue')}) -with(metout,{points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",col='light blue',main=paste("sky temperature",sep=""))}) with(metout,{plot(RHLOC ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",ylim=c(0,100),main=paste("humidity",sep=""))}) @@ -346,11 +346,11 @@ with(metout,{plot(SNOWDEP ~ dates,xlab = "Date and Time", ylab = "Snow Depth (cm # plotting soil temperature for(i in 1:10){ if(i==1){ - plot(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)" + plot(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)" ,col=i,type = "l",main=paste("soil temperature",sep="")) }else{ points(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature - (°C)",col=i,type = "l") + (°C)",col=i,type = "l") } } diff --git a/man/micro_nz.Rd b/man/micro_nz.Rd index 21920c56..e0d75b75 100644 --- a/man/micro_nz.Rd +++ b/man/micro_nz.Rd @@ -314,15 +314,15 @@ minshade<-micro$minshade[1] maxshade<-micro$maxshade[1] # plotting above-ground conditions in minimum shade -with(plotmetout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(plotmetout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l",main=paste("air temperature, ",minshade,"\% shade",sep=""))}) -with(plotmetout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(plotmetout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l",lty=2,col='blue')}) with(plotmetout,{plot(RHLOC ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",ylim=c(0,100),main=paste("humidity, ",minshade,"\% shade",sep=""))}) with(plotmetout,{points(RH ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",col='blue',lty=2,ylim=c(0,100))}) -with(plotmetout,{plot(TSKYC ~ dates,xlab = "Date and Time", ylab = "Sky Temperature (°C)" +with(plotmetout,{plot(TSKYC ~ dates,xlab = "Date and Time", ylab = "Sky Temperature (°C)" , type = "l",main=paste("sky temperature, ",minshade,"\% shade",sep=""))}) with(plotmetout,{plot(VREF ~ dates,xlab = "Date and Time", ylab = "Wind Speed (m/s)" , type = "l",main="wind speed",ylim = c(0, 15))}) @@ -336,34 +336,34 @@ with(plotmetout,{plot(SOLR ~ dates,xlab = "Date and Time", ylab = "Solar Radiati # plotting soil temperature for minimum shade for(i in 1:10){ if(i==1){ - plot(plotsoil[,i+3]~plotsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)" + plot(plotsoil[,i+3]~plotsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)" ,col=i,type = "l",main=paste("soil temperature ",minshade,"\% shade",sep="")) }else{ points(plotsoil[,i+3]~plotsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature - (°C)",col=i,type = "l") + (°C)",col=i,type = "l") } } # plotting above-ground conditions in maximum shade -with(plotshadmet,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(plotshadmet,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l",main="air temperature, sun")}) -with(plotshadmet,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(plotshadmet,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l",lty=2,col='blue')}) with(plotshadmet,{plot(RHLOC ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",ylim=c(0,100),main="humidity, shade")}) with(plotshadmet,{points(RH ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",col='blue',lty=2,ylim=c(0,100))}) -with(plotshadmet,{plot(TSKYC ~ dates,xlab = "Date and Time", ylab = "Sky Temperature (°C)", +with(plotshadmet,{plot(TSKYC ~ dates,xlab = "Date and Time", ylab = "Sky Temperature (°C)", type = "l",main="sky temperature, shade")}) # plotting soil temperature for maximum shade for(i in 1:10){ if(i==1){ plot(plotshadsoil[,i+3]~plotshadsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature - (°C)",col=i,type = "l",main=paste("soil temperature ",maxshade,"\% shade",sep="")) + (°C)",col=i,type = "l",main=paste("soil temperature ",maxshade,"\% shade",sep="")) }else{ points(plotshadsoil[,i+3]~plotshadsoil[,1],xlab = "Date and Time", ylab = "Soil Temperature - (°C)",col=i,type = "l") + (°C)",col=i,type = "l") } } } diff --git a/man/micro_terra.Rd b/man/micro_terra.Rd index cd108104..2089df1b 100644 --- a/man/micro_terra.Rd +++ b/man/micro_terra.Rd @@ -311,15 +311,15 @@ minshade <- micro$minshade[1] maxshade <- micro$maxshade[1] # plotting above-ground conditions in minimum shade -with(metout, {plot(TALOC ~ micro$dates3, xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(metout, {plot(TALOC ~ micro$dates3, xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l", main = paste("air temperature, ", minshade, "\% shade",sep = ""))}) -with(metout, {points(TAREF ~ micro$dates3, xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(metout, {points(TAREF ~ micro$dates3, xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l",lty = 2, col = 'blue')}) with(metout, {plot(RHLOC ~ micro$dates3, xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",ylim = c(0, 100),main = paste("humidity, ",minshade, "\% shade",sep=""))}) with(metout, {points(RH ~ micro$dates3, xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",col = 'blue',lty = 2, ylim = c(0, 100))}) -with(metout, {plot(TSKYC ~ micro$date3s, xlab = "Date and Time", ylab = "Sky Temperature (°C)" +with(metout, {plot(TSKYC ~ micro$date3s, xlab = "Date and Time", ylab = "Sky Temperature (°C)" , type = "l", main = paste("sky temperature, ", minshade, "\% shade", sep=""))}) with(metout, {plot(VREF ~ micro$dates3, xlab = "Date and Time", ylab = "Wind Speed (m/s)" , type = "l", main = "wind speed", col = 'blue',ylim = c(0, 15))}) @@ -333,34 +333,34 @@ with(metout, {plot(SOLR ~ micro$dates3,xlab = "Date and Time", ylab = "Solar Rad # plotting soil temperature for minimum shade for(i in 1:10){ if(i==1){ - plot(soil[,i + 2] ~ micro$dates3, xlab = "Date and Time", ylab = "Soil Temperature (°C)" + plot(soil[,i + 2] ~ micro$dates3, xlab = "Date and Time", ylab = "Soil Temperature (°C)" ,col = i, type = "l", main = paste("soil temperature ", minshade, "\% shade", sep="")) }else{ points(soil[,i + 2] ~ micro$dates3, xlab = "Date and Time", ylab = "Soil Temperature - (°C)", col = i, type = "l") + (°C)", col = i, type = "l") } } points(metout$SNOWDEP ~ micro$dates, type = 'h', col = 'light blue') # plotting above-ground conditions in maximum shade -with(shadmet,{plot(TALOC ~ micro$dates3,xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(shadmet,{plot(TALOC ~ micro$dates3,xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l", main = "air temperature, sun")}) -with(shadmet,{points(TAREF ~ micro$dates3,xlab = "Date and Time", ylab = "Air Temperature (°C)" +with(shadmet,{points(TAREF ~ micro$dates3,xlab = "Date and Time", ylab = "Air Temperature (°C)" , type = "l", lty = 2, col = 'blue')}) with(shadmet,{plot(RHLOC ~ micro$dates3,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l", ylim = c(0, 100),main = "humidity, shade")}) with(shadmet,{points(RH ~ micro$dates3,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l", col = 'blue',lty = 2, ylim = c(0, 100))}) -with(shadmet,{plot(TSKYC ~ micro$dates3,xlab = "Date and Time", ylab = "Sky Temperature (°C)", +with(shadmet,{plot(TSKYC ~ micro$dates3,xlab = "Date and Time", ylab = "Sky Temperature (°C)", type = "l", main = "sky temperature, shade")}) # plotting soil temperature for maximum shade for(i in 1:10){ if(i==1){ plot(shadsoil[,i + 2] ~ micro$dates3, xlab = "Date and Time", ylab = "Soil Temperature - (°C)", col = i, type = "l", main = paste("soil temperature ", maxshade, "\% shade", sep="")) + (°C)", col = i, type = "l", main = paste("soil temperature ", maxshade, "\% shade", sep="")) }else{ points(shadsoil[,i + 2] ~ micro$dates3, xlab = "Date and Time", ylab = "Soil Temperature - (°C)", col = i, type = "l") + (°C)", col = i, type = "l") } } points(shadmet$SNOWDEP ~ micro$dates, type = 'h', col = 'light blue') diff --git a/man/micro_uk.Rd b/man/micro_uk.Rd index 55114186..d06d75a1 100644 --- a/man/micro_uk.Rd +++ b/man/micro_uk.Rd @@ -315,11 +315,11 @@ soilmoist <- cbind(dates, soilmoist) minshade<-micro$minshade # plotting above-ground conditions in minimum shade -with(metout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",main=paste("air and sky temperature, ",minshade,"\% shade",sep=""), ylim = c(-20, 60))}) -with(metout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",lty=2,col='blue')}) -with(metout,{points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout,{points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l",col='light blue',main=paste("sky temperature, ",minshade,"\% shade",sep=""))}) with(metout,{plot(RHLOC ~ dates,xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l",ylim=c(0,100),main=paste("humidity, ",minshade,"\% shade",sep=""))}) @@ -337,11 +337,11 @@ with(metout,{plot(SNOWDEP ~ dates,xlab = "Date and Time", ylab = "Snow Depth (cm # plotting soil temperature for(i in 1:10){ if(i==1){ - plot(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)" + plot(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature (°C)" ,col=i,type = "l",main=paste("soil temperature ",minshade,"\% shade",sep="")) }else{ points(soil[,i+3]~soil[,1],xlab = "Date and Time", ylab = "Soil Temperature - (°C)",col=i,type = "l") + (°C)",col=i,type = "l") } } diff --git a/man/micro_usa.Rd b/man/micro_usa.Rd index 6833d8a1..3f408cd3 100644 --- a/man/micro_usa.Rd +++ b/man/micro_usa.Rd @@ -307,11 +307,11 @@ soilmoist <- cbind(dates, soilmoist) minshade<-micro$minshade # plotting above-ground conditions in minimum shade -with(metout, {plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout, {plot(TALOC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l", main = paste("air and sky temperature, ", minshade, "\% shade", sep = ""), ylim = c(-20, 60))}) -with(metout, {points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout, {points(TAREF ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l", lty = 2,col = 'blue')}) -with(metout, {points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" +with(metout, {points(TSKYC ~ dates,xlab = "Date and Time", ylab = "Temperature (°C)" , type = "l", col = 'light blue', main = paste("sky temperature, ", minshade, "\% shade", sep = ""))}) with(metout, {plot(RHLOC ~ dates, xlab = "Date and Time", ylab = "Relative Humidity (\%)" , type = "l", ylim = c(0, 100), main = paste("humidity, ", minshade, "\% shade", sep = ""))}) @@ -329,11 +329,11 @@ with(metout, {plot(SNOWDEP ~ dates, xlab = "Date and Time", ylab = "Snow Depth ( # plotting soil temperature for(i in 1:10){ if(i==1){ - plot(soil[,i+3] ~ soil[,1], xlab = "Date and Time", ylab = "Soil Temperature (°C)" + plot(soil[,i+3] ~ soil[,1], xlab = "Date and Time", ylab = "Soil Temperature (°C)" ,col = i,type = "l", main = paste("soil temperature ", minshade, "\% shade", sep = "")) }else{ points(soil[, i + 3] ~ soil[, 1], xlab = "Date and Time", ylab = "Soil Temperature - (°C)", col = i, type = "l") + (°C)", col = i, type = "l") } } diff --git a/man/onelump.Rd b/man/onelump.Rd index 7c23fffb..13485f22 100644 --- a/man/onelump.Rd +++ b/man/onelump.Rd @@ -81,10 +81,10 @@ tmins <- t/60 par(mfrow = c(1,2)) Ww_g <- 5 # body weight, g -Tc_init <- 20 # initial body temperature, °C +Tc_init <- 20 # initial body temperature, °C geom <- 2 # shape (2 = ellipsoid) -Tair <- 20 # air temperature, °C -Trad <- Tair # radiant temperature, °C +Tair <- 20 # air temperature, °C +Trad <- Tair # radiant temperature, °C vel <- 1 # wind speed, m/s Qsol <- 500 # horizontal plane solar radiation, W/m2 Zen <- 20 # zenith angle of sun, degrees @@ -92,7 +92,7 @@ alpha <- 0.85 # solar absorptivity, - Tbs<-onelump(t=t, alpha = alpha, Tc_init = Tc_init, Ww_g = Ww_g, geom = geom, Tair = Tair, Trad = Trad, vel = vel, Qsol = Qsol, Zen = Zen) -plot(Tbs$Tc ~ tmins, type= 'l' ,col = 1, ylim = c(20, 30), ylab = 'Temperature, °C',xlab='time, min', las = 1) +plot(Tbs$Tc ~ tmins, type= 'l' ,col = 1, ylim = c(20, 30), ylab = 'Temperature, °C',xlab='time, min', las = 1) text(80, 27, " 500 g") text(80, 24.5, "5 g") text(90, 20.5, "Tair for both sizes", col = "blue") @@ -107,16 +107,16 @@ abline(Tair,0, lty = 1, col = 'blue') abline(h = Tair + .1, lty = 2, col = 'blue') Ww_g <- 5 # body weight, g -Tair <- 25 # air temperature, °C +Tair <- 25 # air temperature, °C vel <-0.5 # wind speed, m/s Tbs<-onelump(t=t, alpha = alpha, Tc_init = Tc_init, Ww_g = Ww_g, geom = geom, Tair = Tair, Trad = Trad, vel = vel, Qsol = Qsol, Zen = Zen) -plot(Tbs$Tc~tmins,type='l',col=1,ylim=c(20,30),ylab='Temperature, °C',xlab='time, min', las = 1) +plot(Tbs$Tc~tmins,type='l',col=1,ylim=c(20,30),ylab='Temperature, °C',xlab='time, min', las = 1) abline(h = Tair, lty = 1, col = 'blue') Ww_g <- 500 # body weight, g -Tair <- 20 # air temperature, °C +Tair <- 20 # air temperature, °C vel <-1 # wind speed, m/s Tbs<-onelump(t=t, alpha = alpha, Tc_init = Tc_init, Ww_g = Ww_g, diff --git a/man/onelump_var.Rd b/man/onelump_var.Rd index f10065b4..b203e54c 100644 --- a/man/onelump_var.Rd +++ b/man/onelump_var.Rd @@ -144,7 +144,7 @@ for (i in 1:length(DOYs)) { colnames(Tbs_ode) <- c('time', 'Tc', 'Tcf', 'tau', 'dTdt') Tbs_ode$time <- Tbs_ode$time / 3600 # convert to hours - with(Tbs_ode, plot(Tc ~ time, type = 'l', col = '1', ylim = c(-10, 80), xlim = c(0, 23), ylab='Temperature, °C',xlab = 'hour of day', main = paste0(round(loc[1], 1), round(loc[2], 1), ", ", mons[i], ", ", Ww_g," g"))) + with(Tbs_ode, plot(Tc ~ time, type = 'l', col = '1', ylim = c(-10, 80), xlim = c(0, 23), ylab='Temperature, °C',xlab = 'hour of day', main = paste0(round(loc[1], 1), round(loc[2], 1), ", ", mons[i], ", ", Ww_g," g"))) with(Tbs_ode, points(Tcf ~ time, type = 'l', col = '2')) points(Tairf(time) ~ hours, type = 'l', col = 'blue', lty = 2) legend(0,70, c("Tc", "Tcf", "Tair"), lty = c(1, 1, 2), lwd = c(2.5, 2.5, 2.5), col = c("black", "red", "blue")) diff --git a/man/twolump.Rd b/man/twolump.Rd index b71cb2ab..798ae9a7 100644 --- a/man/twolump.Rd +++ b/man/twolump.Rd @@ -157,7 +157,7 @@ for (i in 1:length(DOYs)) { colnames(Tbs_ode) = c("time", "Tc", "Tsk", "Ts", "Tcf") Tbs_ode$time <- Tbs_ode$time / 3600 # convert to hours - with(Tbs_ode, plot(Tc ~ time, type = 'l', col = '1', ylim = c(-10, 80), xlim = c(0, 23), ylab='Temperature, °C',xlab = 'hour of day', main = paste0(loc, ", ", mons[i], ", ", Ww_g," g"))) + with(Tbs_ode, plot(Tc ~ time, type = 'l', col = '1', ylim = c(-10, 80), xlim = c(0, 23), ylab='Temperature, °C',xlab = 'hour of day', main = paste0(loc, ", ", mons[i], ", ", Ww_g," g"))) with(Tbs_ode, points(Ts ~ time, lty = 2, type = 'l', col = '1')) with(Tbs_ode, points(Tcf ~ time, type = 'l', col = '2')) points(Tairf(time) ~ hours, type = 'l', col = 'blue', lty = 2) diff --git a/src/CONV_ENDO.f b/src/CONV_ENDO.f index 1f0a9971..11354897 100644 --- a/src/CONV_ENDO.f +++ b/src/CONV_ENDO.f @@ -124,7 +124,7 @@ SUBROUTINE CONV_ENDO(TS,TENV,NGEOM,SURFAR,FLTYPE,FURTST,D,TFA, PR=CP*VISDYN/THCOND C COMPUTING SCHMIDT NUMBER - IF (FLTYPE .EQ. 0.00) THEN + IF (INT(FLTYPE) .EQ. 0) THEN C AIR SC=VISDYN/(DENSTY*DIFVPR) ELSE @@ -143,11 +143,11 @@ SUBROUTINE CONV_ENDO(TS,TENV,NGEOM,SURFAR,FLTYPE,FURTST,D,TFA, ENDIF C STABILITY CHECK - IF(DELTAT .EQ. 0.0000000E+00)THEN + IF(DELTAT .LE. 0.0000000E+00)THEN DELTAT = DELTAT + 0.00001 ENDIF - GR=((DENSTY**2)*BETA*GRAV*(D**3)*DELTAT)/(VISDYN**2) + GR=((DENSTY**2.)*BETA*GRAV*(D**3.)*DELTAT)/(VISDYN**2.) C CORRECTING IF NEGATIVE DELTAT GR = ABS(GR) C RAYLEIGH NUMBER @@ -161,13 +161,13 @@ SUBROUTINE CONV_ENDO(TS,TENV,NGEOM,SURFAR,FLTYPE,FURTST,D,TFA, C IF A TRUNCATED CONE (5) C WE WILL USE THE CYLINDER EQUATIONS (NGEOM=1). - IF(NGEOM.EQ.5)THEN + IF(INT(NGEOM).EQ.5)THEN NGEOM = 1 ENDIF C ********************* FREE CONVECTION ******************** - IF(NGEOM.EQ.1)THEN + IF(INT(NGEOM).EQ.1)THEN C FREE CONVECTION IN CYLINDER C FROM P.334 KREITH (1965): MC ADAM'S 1954 RECOMMENDED COORDINATES C BUT CHECK OUT P. 443-445 IN BIRD, STEWART & LIGHTFOOT, 2002 @@ -183,7 +183,7 @@ SUBROUTINE CONV_ENDO(TS,TENV,NGEOM,SURFAR,FLTYPE,FURTST,D,TFA, ENDIF C FREE CONVECTION IN SPHERE OR ELLIPSOID - IF((NGEOM.EQ.2).OR.(NGEOM.EQ.4))THEN + IF((INT(NGEOM).EQ.2).OR.(INT(NGEOM).EQ.4))THEN C FREE CONVECTION IN SPHERE C FROM P.413 BIRD ET AL (1960) TRANSPORT PHENOMENA) ANU=2.+0.60* ((GR**.25)*(PR**.333)) @@ -200,7 +200,7 @@ SUBROUTINE CONV_ENDO(TS,TENV,NGEOM,SURFAR,FLTYPE,FURTST,D,TFA, ENDIF ENDIF - IF(NGEOM.EQ.3)THEN + IF(INT(NGEOM).EQ.3)THEN C FREE CONVECTION FOR A PLATE (AS USED BY PHILLIPS & HEATH 1992, C GATES EQN. 9.77, ASSUMES TURBULENT FLOW, NOTE PHILLIPS & HEATH C HAD 2/3 EXPONENT NOT 1/3 AS IN GATES) @@ -232,7 +232,7 @@ SUBROUTINE CONV_ENDO(TS,TENV,NGEOM,SURFAR,FLTYPE,FURTST,D,TFA, C COMPUTING NUSSELT NUMBER C CYLINDER - IF(NGEOM.EQ.1)THEN + IF(INT(NGEOM).EQ.1)THEN IF(RE.LE.4.)THEN ANU=.891*RE**.33 ELSE IF((RE.GT.4.).AND.(RE.LE.40))THEN @@ -247,10 +247,10 @@ SUBROUTINE CONV_ENDO(TS,TENV,NGEOM,SURFAR,FLTYPE,FURTST,D,TFA, ENDIF ENDIF C SPHERE OR ELLIPSOID - IF((NGEOM.EQ.2).OR.(NGEOM.EQ.4))THEN + IF((INT(NGEOM).EQ.2).OR.(INT(NGEOM).EQ.4))THEN ANU=0.37*RE**0.6 ENDIF - IF(NGEOM.EQ.3)THEN + IF(INT(NGEOM).EQ.3)THEN C FORCED CONVECTION FOR A PLATE (AS USED BY PHILLIPS & HEATH 1992, C GATES EQN. 9.49, ASSUMES TURBULENT FLOW) ANU = 0.032*RE**0.8 @@ -274,7 +274,7 @@ SUBROUTINE CONV_ENDO(TS,TENV,NGEOM,SURFAR,FLTYPE,FURTST,D,TFA, C ENDIF C USING BIRD, STEWART, & LIGHTFOOT'S MIXED CONVECTION FORMULA (P. 445, TRANSPORT PHENOMENA, 2002) - NUTOTAL = (NUFREE**3 + NUFORCED**3)**(1./3.) + NUTOTAL = (NUFREE**3. + NUFORCED**3.)**(1./3.) HCCOMB = NUTOTAL*(THCOND/D) C DIFFERENT QCONV DEPENDING ON WHETHER THERE IS FUR OR NOT diff --git a/src/GEOM_ENDO.f b/src/GEOM_ENDO.f index 8f51de12..1dfa2726 100644 --- a/src/GEOM_ENDO.f +++ b/src/GEOM_ENDO.f @@ -71,13 +71,13 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT VOLFAT = MASFAT/FATDEN C CYLINDER - IF (SHAPE .EQ. 1) THEN - R1 = (VOL/(SHAPEB*PI*2))**(1./3.) !R_SKIN + IF (INT(SHAPE) .EQ. 1) THEN + R1 = (VOL/(SHAPEB*PI*2.))**(1./3.) !R_SKIN R2 = R1 + ZFUR ! R_FA - ALENTH=SHAPEB*R1*2 + ALENTH=SHAPEB*R1*2. C FUR-AIR AREA = END CIRCLE AREA + CIRCLE PERIMETER * LENGTH (LENGTH = 2*D2) - AREASKIN=2.*PI*R1**2+2.*PI*R1*ALENTH - ATOT=2.*PI*R2**2+2.*PI*R2*ALENTH + AREASKIN=2.*PI*R1**2.+2.*PI*R1*ALENTH + ATOT=2.*PI*R2**2.+2.*PI*R2*ALENTH AWIDTH=2.*R2 AHEIT=AWIDTH AREA=ATOT @@ -92,24 +92,24 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT FATTHK = 0. ENDIF ASILN = AWIDTH * ALENTH - ASILP = PI * R2**2 - IF(ORIENT.EQ.0)THEN ! AVERAGE POSTURE + ASILP = PI * R2**2. + IF(INT(ORIENT).EQ.0)THEN ! AVERAGE POSTURE ASIL = (ASILN + ASILP)/2. ENDIF - IF(ORIENT.EQ.1)THEN ! NORMAL TO SUN'S RAYS (HEAT MAXIMISING) + IF(INT(ORIENT).EQ.1)THEN ! NORMAL TO SUN'S RAYS (HEAT MAXIMISING) ASIL = ASILN ENDIF - IF(ORIENT.EQ.2)THEN ! PARALLEL TO SUN'S RAYS (HEAT MINIMISING) + IF(INT(ORIENT).EQ.2)THEN ! PARALLEL TO SUN'S RAYS (HEAT MINIMISING) ASIL = ASILP ENDIF - IF(ORIENT.EQ.3)THEN - ASIL=2*R2*ALENTH*sin(ZEN*pi/180.)+pi*R2**2*cos(ZEN*pi/180.) + IF(INT(ORIENT).EQ.3)THEN + ASIL=2.*R2*ALENTH*sin(ZEN*pi/180.)+pi*R2**2.*cos(ZEN*pi/180.) ENDIF D = VOL**(1./3.) GO TO 999 ENDIF - IF (SHAPE .EQ. 2) THEN + IF (INT(SHAPE) .EQ. 2) THEN C SPHERE C BODY, FUR DIMENSIONS FOR CONDUCTION-RADIATION C VOL = 4./3. * PI * R**3 @@ -132,10 +132,10 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT RFLESH = R1 ENDIF C SKIN SURFACE AREA - AREASKIN=4.*PI*R1**2 ! SPHERE WITHOUT FUR - AREA = 4.*PI*R2**2 ! SPHERE WITH FUR + AREASKIN=4.*PI*R1**2. ! SPHERE WITHOUT FUR + AREA = 4.*PI*R2**2. ! SPHERE WITH FUR ATOT = AREA - ASIL = PI * R2**2 + ASIL = PI * R2**2. ASILN = ASIL ASILP = ASIL ALENTH=R2*2. @@ -145,20 +145,20 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT ENDIF C FLAT PLATE - IF (SHAPE .EQ. 3) THEN + IF (INT(SHAPE) .EQ. 3) THEN AHEIT=(VOL/(SHAPEB*SHAPEC))**(1./3.) ALENTH=AHEIT*SHAPEB AWIDTH=AHEIT*SHAPEC AREASKIN=ALENTH*AWIDTH*2.+ALENTH*AHEIT*2.+AWIDTH*AHEIT*2. ASILN = (AWIDTH+ZFUR) * (ALENTH+ZFUR) ASILP = (AWIDTH+ZFUR) * (AHEIT+ZFUR) - IF((ORIENT.EQ.0).OR.(ORIENT.EQ.3))THEN ! AVERAGE POSTURE + IF((INT(ORIENT).EQ.0).OR.(INT(ORIENT).EQ.3))THEN ! AVERAGE POSTURE ASIL = (ASILN + ASILP)/2. ENDIF - IF(ORIENT.EQ.1)THEN ! NORMAL TO SUN'S RAYS (HEAT MAXIMISING) + IF(INT(ORIENT).EQ.1)THEN ! NORMAL TO SUN'S RAYS (HEAT MAXIMISING) ASIL = ASILN ENDIF - IF(ORIENT.EQ.2)THEN ! PARALLEL TO SUN'S RAYS (HEAT MINIMISING) + IF(INT(ORIENT).EQ.2)THEN ! PARALLEL TO SUN'S RAYS (HEAT MINIMISING) ASIL = ASILP ENDIF R1 = AWIDTH/2. @@ -181,7 +181,7 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT ENDIF C ELLIPSOID -100 IF (SHAPE .EQ. 4) THEN +100 IF (INT(SHAPE) .EQ. 4) THEN IF(INT(SUBQFAT).EQ.0.OR.(MASFAT.LE.0.00))THEN ! NO SUBCUTANEOUS FAT C ASSUMING A PROLATE SPHEROID C VOL = 4/3* PI*A*B*C @@ -193,27 +193,28 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT ALENTH=ASEMAJ*2. AHEIT=BSEMIN*2. AWIDTH=CSEMIN*2. - E = ((ASEMAJ**2 - CSEMIN**2)**0.5 )/ASEMAJ ! ECCENTRICITY + E = ((ASEMAJ**2. - CSEMIN**2)**0.5 )/ASEMAJ ! ECCENTRICITY ASILN = PI * (ASEMAJ + ZFUR) * (BSEMIN + ZFUR) ASILP = PI * (BSEMIN + ZFUR) * (CSEMIN + ZFUR) - IF(ORIENT.EQ.0)THEN ! AVERAGE POSTURE + IF(INT(ORIENT).EQ.0)THEN ! AVERAGE POSTURE ASIL = (ASILN + ASILP)/2. ENDIF - IF(ORIENT.EQ.1)THEN ! NORMAL TO SUN'S RAYS (HEAT MAXIMISING) + IF(INT(ORIENT).EQ.1)THEN ! NORMAL TO SUN'S RAYS (HEAT MAXIMISING) ASIL = ASILN ENDIF - IF(ORIENT.EQ.2)THEN ! PARALLEL TO SUN'S RAYS (HEAT MINIMISING) + IF(INT(ORIENT).EQ.2)THEN ! PARALLEL TO SUN'S RAYS (HEAT MINIMISING) ASIL = ASILP ENDIF - IF(ORIENT.EQ.3)THEN ! COMPUTE AS A FUNCTION OF ZENITH ANGLE + IF(INT(ORIENT).EQ.3)THEN ! COMPUTE AS A FUNCTION OF ZENITH ANGLE STH=sin(90.*pi/180.) CTH=cos(90.*pi/180.) SPH=sin(ZEN*pi/180.) CPH=cos(ZEN*pi/180.) - AA=CTH**2*(CPH**2/(ASEMAJ+ZFUR)**2+SPH**2/(BSEMIN+ZFUR)**2) - & +STH**2/(CSEMIN+ZFUR)**2 - TWOHH=2.*CTH*STH*CPH*(1./(BSEMIN+ZFUR)**2-1./(ASEMAJ+ZFUR)**2) - BB=SPH**2/(ASEMAJ+ZFUR)**2+CPH**2/(BSEMIN+ZFUR)**2 + AA=CTH**2.*(CPH**2./(ASEMAJ+ZFUR)**2.+SPH**2./ + & (BSEMIN+ZFUR)**2.)+STH**2./(CSEMIN+ZFUR)**2. + TWOHH=2.*CTH*STH*CPH*(1./(BSEMIN+ZFUR)**2.-1./ + & (ASEMAJ+ZFUR)**2.) + BB=SPH**2./(ASEMAJ+ZFUR)**2.+CPH**2/(BSEMIN+ZFUR)**2. PS=0.5*atan2(TWOHH,AA-BB) SPS=sin(PS) CPS=cos(PS) @@ -224,10 +225,10 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT ASIL=pi*SEMAX1*SEMAX2 ENDIF C SKIN SURFACE AREA - IF(SAMODE.EQ.1.)THEN + IF(INT(SAMODE).EQ.1)THEN AREASKIN=10.*(AMASS*1000.)**0.667/10000. ! WALSBERG G. E. & KING J. E. (1978) THE RELATIONSHIP OF THE EXTERNAL SURFACE AREA OF BIRDS TO SKIN SURFACE AREA AND BODY MASS. JOURNAL OF EXPERIMENTAL BIOLOGY 76 , 185-189. ELSE - IF(SAMODE.EQ.2)THEN + IF(INT(SAMODE).EQ.2)THEN AREASKIN=1110.*AMASS**0.65/10000. ! STAHL W. R. (1967) SCALING OF RESPIRATORY VARIABLES IN MAMMALS. JOURNAL OF APPLIED PHYSIOLOGY 22 , 453-460. ELSE AREASKIN=2.*PI*BSEMIN**2.+2.*PI*(ASEMAJ*BSEMIN/E)*ASIN(E) ! ELLIPSOID WITHOUT FUR @@ -237,19 +238,19 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT R2 = R1 + ZFUR A2 = ASEMAJ + ZFUR C2 = R1 + ZFUR - E2 = ((A2**2 - C2**2)**0.5 )/A2 ! ECCENTRICITY FOR TOTAL AREA + E2 = ((A2**2. - C2**2.)**0.5 )/A2 ! ECCENTRICITY FOR TOTAL AREA C AREA AT THE FUR TIPS - IF(SAMODE.EQ.1.)THEN + IF(INT(SAMODE).EQ.1)THEN AREA=8.11*(AMASS*1000.)**0.667/10000. ! WALSBERG G. E. & KING J. E. (1978) THE RELATIONSHIP OF THE EXTERNAL SURFACE AREA OF BIRDS TO SKIN SURFACE AREA AND BODY MASS. JOURNAL OF EXPERIMENTAL BIOLOGY 76 , 185-189. ELSE - IF(SAMODE.EQ.2)THEN + IF(INT(SAMODE).EQ.2)THEN AREA=1110.*AMASS**0.65* ! STAHL W. R. (1967) SCALING OF RESPIRATORY VARIABLES IN MAMMALS. JOURNAL OF APPLIED PHYSIOLOGY 22 , 453-460. - & (2.*PI*R2**2 + 2.*PI*(A2*R2/E2)*ASIN(E2)/ ! OUTER SURFACE AREA IF AN ELLIPSOID WITH FUR - & (2.*PI*BSEMIN**2 + 2.*PI*(ASEMAJ*BSEMIN/E)*ASIN(E))) ! OUTER SURFACE AREA OF ELLIPSOID WITHOUT FUR + & (2.*PI*R2**2. + 2.*PI*(A2*R2/E2)*ASIN(E2)/ ! OUTER SURFACE AREA IF AN ELLIPSOID WITH FUR + & (2.*PI*BSEMIN**2. + 2.*PI*(ASEMAJ*BSEMIN/E)*ASIN(E))) ! OUTER SURFACE AREA OF ELLIPSOID WITHOUT FUR & /10000. ! INCREASE IN AREA DUE TO FUR ELSE - AREA = 2.*PI*R2**2 + 2.*PI*(A2*R2/E2)*ASIN(E2) ! ELLIPSE WITH FUR + AREA = 2.*PI*R2**2. + 2.*PI*(A2*R2/E2)*ASIN(E2) ! ELLIPSE WITH FUR ENDIF ENDIF ATOT = AREA @@ -275,14 +276,14 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT CFA=1.0 CFB=(SHAPEB*FLSHBSEMIN)+(2*FLSHBSEMIN) CFC=(2.*SHAPEB*(FLSHBSEMIN**2.))+(FLSHBSEMIN**2.) - CFD=(SHAPEB*(FLSHBSEMIN**3.))-(((VOLFAT+FLSHVL)*3)/(4*PI)) + CFD=(SHAPEB*(FLSHBSEMIN**3.))-(((VOLFAT+FLSHVL)*3.)/(4.*PI)) C BREAKING UP THE FORMULA INTO COMPONENT PARTS FOR CLEARER CODING. CFT1A=((-1.*CFB)**3.)/(27.*(CFA**3.)) CFT1B=(CFB*CFC)/(6*(CFA**2.)) - CFT1C=(CFD/(2*CFA)) + CFT1C=(CFD/(2.*CFA)) CFT1= CFT1A+CFT1B-CFT1C CFT2A=(CFT1**2.) - CFT2B=((CFC/(3*CFA))-((CFB**2.)/(9*(CFA**2.))))**3. + CFT2B=((CFC/(3.*CFA))-((CFB**2.)/(9*(CFA**2.))))**3. C CORRECTING FOR WHEN CFT2A+CFT2B IS LESS THAN 0.0, WHICH CREATES A PROBLEM IN THE C FAT THICKNESS CALCULATION. THIS OCCURS WHEN FAT IS GETTING CLOSE TO NOT BEING ENOUGH TO COVER THE ELLIPSOID. IF(CFT2A+CFT2B.GE.0.0)THEN @@ -290,17 +291,17 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT ELSE CFT2=0.0 ENDIF - CFT3=(CFB/(3*CFA)) + CFT3=(CFB/(3.*CFA)) C FATTHK=((CFT1+CFT2)**0.333)+((CFT1-CFT2)**0.333)-CFT3 C FORTRAN CAN'T SEEM TO HANDLE THE CUBE ROOTS OF NEGATIVE NUMBERS SO WE NEED TO C DO SOME MODIFCATIONS IF((CFT1+CFT2).LT.0.0)THEN - CFT4= -1.*((-1*(CFT1+CFT2))**(1./3.)) + CFT4= -1.*((-1.*(CFT1+CFT2))**(1./3.)) ELSE CFT4= (CFT1+CFT2)**(1./3.) ENDIF IF((CFT1-CFT2).LT.0.0)THEN - CFT5= -1.*((-1*(CFT1-CFT2))**(1./3.)) + CFT5= -1.*((-1.*(CFT1-CFT2))**(1./3.)) ELSE CFT5= (CFT1-CFT2)**(1./3.) ENDIF @@ -310,7 +311,7 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT C SETTING SUBQFAT TO 0 C SETTING FATTHK TO 0.0 FATTHK=0.0 - SUBQFAT=0 + SUBQFAT=0. GO TO 100 ENDIF ASEMAJ=FLSHASEMAJ+FATTHK @@ -320,27 +321,28 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT ALENTH=ASEMAJ*2. AHEIT=BSEMIN*2. AWIDTH=CSEMIN*2. - E = ((ASEMAJ**2 - CSEMIN**2)**0.5 )/ASEMAJ ! ECCENTRICITY + E = ((ASEMAJ**2. - CSEMIN**2.)**0.5 )/ASEMAJ ! ECCENTRICITY ASILN = PI * (ASEMAJ + ZFUR) * (BSEMIN + ZFUR) ASILP = PI * (BSEMIN + ZFUR) * (CSEMIN + ZFUR) - IF(ORIENT.EQ.0)THEN ! AVERAGE POSTURE + IF(INT(ORIENT).EQ.0)THEN ! AVERAGE POSTURE ASIL = (ASILN + ASILP)/2. ENDIF - IF(ORIENT.EQ.1)THEN ! NORMAL TO SUN'S RAYS (HEAT MAXIMISING) + IF(INT(ORIENT).EQ.1)THEN ! NORMAL TO SUN'S RAYS (HEAT MAXIMISING) ASIL = ASILN ENDIF - IF(ORIENT.EQ.2)THEN ! PARALLEL TO SUN'S RAYS (HEAT MINIMISING) + IF(INT(ORIENT).EQ.2)THEN ! PARALLEL TO SUN'S RAYS (HEAT MINIMISING) ASIL = ASILP ENDIF - IF(ORIENT.EQ.3)THEN ! COMPUTE AS A FUNCTION OF ZENITH ANGLE + IF(INT(ORIENT).EQ.3)THEN ! COMPUTE AS A FUNCTION OF ZENITH ANGLE STH=sin(90.*pi/180.) CTH=cos(90.*pi/180.) SPH=sin(ZEN*pi/180.) CPH=cos(ZEN*pi/180.) - AA=CTH**2*(CPH**2/(ASEMAJ+ZFUR)**2+SPH**2/(BSEMIN+ZFUR)**2) - & +STH**2/(CSEMIN+ZFUR)**2 - TWOHH=2.*CTH*STH*CPH*(1./(BSEMIN+ZFUR)**2-1./(ASEMAJ+ZFUR)**2) - BB=SPH**2/(ASEMAJ+ZFUR)**2+CPH**2/(BSEMIN+ZFUR)**2 + AA=CTH**2.*(CPH**2./(ASEMAJ+ZFUR)**2.+SPH**2./ + & (BSEMIN+ZFUR)**2.)+STH**2./(CSEMIN+ZFUR)**2. + TWOHH=2.*CTH*STH*CPH*(1./(BSEMIN+ZFUR)**2.-1./ + & (ASEMAJ+ZFUR)**2.) + BB=SPH**2./(ASEMAJ+ZFUR)**2.+CPH**2./(BSEMIN+ZFUR)**2. PS=0.5*atan2(TWOHH,AA-BB) SPS=sin(PS) CPS=cos(PS) @@ -351,10 +353,10 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT ASIL=pi*SEMAX1*SEMAX2 ENDIF C SKIN SURFACE AREA - IF(SAMODE.EQ.1.)THEN + IF(INT(SAMODE).EQ.1)THEN AREASKIN=10.*(AMASS*1000.)**0.667/10000. ! WALSBERG G. E. & KING J. E. (1978) THE RELATIONSHIP OF THE EXTERNAL SURFACE AREA OF BIRDS TO SKIN SURFACE AREA AND BODY MASS. JOURNAL OF EXPERIMENTAL BIOLOGY 76 , 185-189. ELSE - IF(SAMODE.EQ.2)THEN + IF(INT(SAMODE).EQ.2)THEN AREASKIN=1110.*AMASS**0.65/10000. ! STAHL W. R. (1967) SCALING OF RESPIRATORY VARIABLES IN MAMMALS. JOURNAL OF APPLIED PHYSIOLOGY 22 , 453-460. ELSE AREASKIN=2.*PI*BSEMIN**2.+2.*PI*(ASEMAJ*BSEMIN/E)*ASIN(E) ! ELLIPSOID WITHOUT FUR @@ -364,18 +366,18 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT R2 = R1 + ZFUR A2 = ASEMAJ + ZFUR C2 = R1 + ZFUR - E2 = ((A2**2 - C2**2)**0.5 )/A2 ! ECCENTRICITY FOR TOTAL AREA + E2 = ((A2**2. - C2**2.)**0.5 )/A2 ! ECCENTRICITY FOR TOTAL AREA C AREA AT THE FUR TIPS - IF(SAMODE.EQ.1.)THEN + IF(INT(SAMODE).EQ.1)THEN AREA=8.11*(AMASS*1000.)**0.667/10000. ! WALSBERG G. E. & KING J. E. (1978) THE RELATIONSHIP OF THE EXTERNAL SURFACE AREA OF BIRDS TO SKIN SURFACE AREA AND BODY MASS. JOURNAL OF EXPERIMENTAL BIOLOGY 76 , 185-189. ELSE - IF(SAMODE.EQ.2)THEN + IF(INT(SAMODE).EQ.2)THEN AREA=1110.*AMASS**0.65* ! STAHL W. R. (1967) SCALING OF RESPIRATORY VARIABLES IN MAMMALS. JOURNAL OF APPLIED PHYSIOLOGY 22 , 453-460. - & (2.*PI*R2**2 + 2.*PI*(A2*R2/E2)*ASIN(E2)/ ! OUTER SURFACE AREA IF AN ELLIPSOID WITH FUR - & (2.*PI*BSEMIN**2 + 2.*PI*(ASEMAJ*BSEMIN/E)*ASIN(E))) ! OUTER SURFACE AREA OF ELLIPSOID WITHOUT FUR + & (2.*PI*R2**2. + 2.*PI*(A2*R2/E2)*ASIN(E2)/ ! OUTER SURFACE AREA IF AN ELLIPSOID WITH FUR + & (2.*PI*BSEMIN**2. + 2.*PI*(ASEMAJ*BSEMIN/E)*ASIN(E))) ! OUTER SURFACE AREA OF ELLIPSOID WITHOUT FUR & /10000. ! INCREASE IN AREA DUE TO FUR ELSE - AREA = 2.*PI*R2**2 + 2.*PI*(A2*R2/E2)*ASIN(E2) ! ELLIPSE WITH FUR + AREA = 2.*PI*R2**2. + 2.*PI*(A2*R2/E2)*ASIN(E2) ! ELLIPSE WITH FUR ENDIF ENDIF ATOT = AREA @@ -394,7 +396,7 @@ SUBROUTINE GEOM_ENDO(AMASS,ANDENS,FATPCT,SHAPE,ZFUR,SUBQFAT C (AREA/HAIR)*(HAIRS/M2)*M2 DHAIR = DHARA RHO = RHOARA - HAIRAR = (PI * (DHAIR/2.)**2) * (RHO * AREASKIN) + HAIRAR = (PI * (DHAIR/2.)**2.) * (RHO * AREASKIN) C SKIN AREA IS TOTAL AREA - HAIR AREA CONVSK = AREASKIN - HAIRAR ELSE diff --git a/src/GETKFUR.f b/src/GETKFUR.f index d78dbfb0..735c809a 100644 --- a/src/GETKFUR.f +++ b/src/GETKFUR.f @@ -63,7 +63,7 @@ SUBROUTINE GETKFUR(RHOAR,LHAR,ZZFUR,DHAR,KAIR,KHAIR,RESULTS) IF (HAIRSP .LT. 0.0000000) THEN C RECALCULATING DENSITY BASED ON SPECIFIED HAIR DIAMETER LUNIT = DHAR - RHOEFF = 1./LUNIT**2 + RHOEFF = 1./LUNIT**2. RHOAR = (RHOEFF*ZZFUR)/LHAR RHOCM2 = RHOAR/10000. C RESETTING HAIR DENSITY TO MAX. POSBL VALUE.' diff --git a/src/IRPROP.f b/src/IRPROP.f index f1692ee1..4c53908b 100644 --- a/src/IRPROP.f +++ b/src/IRPROP.f @@ -61,11 +61,11 @@ SUBROUTINE IRPROP(TA,DHAIRD,DHAIRV,LHAIRD,LHAIRV,ZFURD,ZFURV, FURTST = RHOD*DHAIRD*LHAIRD*ZFURD C COMPUTING AVERAGE VALUES OF PARAMETERS (ALREADY IN SI UNITS) - RHO = RHOD*(1-PVEN) + (RHOV*PVEN) - DHAIR = (DHAIRD*(1-PVEN)) + DHAIRV*PVEN - LHAIR = (LHAIRD*(1-PVEN)) + LHAIRV*PVEN - ZFUR = (ZFURD*(1-PVEN)) + ZFURV*PVEN - REFL = (REFLD*(1-PVEN)) + REFLV*PVEN + RHO = RHOD*(1.-PVEN) + (RHOV*PVEN) + DHAIR = (DHAIRD*(1.-PVEN)) + DHAIRV*PVEN + LHAIR = (LHAIRD*(1.-PVEN)) + LHAIRV*PVEN + ZFUR = (ZFURD*(1.-PVEN)) + ZFURV*PVEN + REFL = (REFLD*(1.-PVEN)) + REFLV*PVEN C NOW PUT THE DATA FROM THE 3 PARTS OF THE ANIMAL INTO THE APPROPRIATE ARRAY C INDEX IS AVERAGE (1), DORSAL(2), VENTRAL(3) @@ -85,16 +85,16 @@ SUBROUTINE IRPROP(TA,DHAIRD,DHAIRV,LHAIRD,LHAIRV,ZFURD,ZFURV, C VENTRAL FUR VALUES OF BODY PART C ACCOUNTING FOR PTVEN IN VENTRAL FUR PROPERTIES - PVENV = PVEN*2 ! BECAUSE ANIMAL IS CONSIDERED AS THE AVERAGE OF A WHOLLY DORSAL AND 'WHOLLY' VENTRAL CALC, BUT VENTRAL FUR MAY NOT COVER ALL OF VENTRAL HALF + PVENV = PVEN*2. ! BECAUSE ANIMAL IS CONSIDERED AS THE AVERAGE OF A WHOLLY DORSAL AND 'WHOLLY' VENTRAL CALC, BUT VENTRAL FUR MAY NOT COVER ALL OF VENTRAL HALF IF(PVENV.GT.1.0)THEN PVENV=1.0 ENDIF - DHAR(3) = DHAIRD*(1-PVENV) + DHAIRV*PVENV - LHAR(3) = LHAIRD*(1-PVENV) + LHAIRV*PVENV - RHOAR(3) = RHOD*(1-PVENV) + RHOV*PVENV - ZZFUR(3) = ZFURD*(1-PVENV) + ZFURV*PVENV - REFLFR(3) = REFLD*(1-PVENV) + REFLV*PVENV + DHAR(3) = DHAIRD*(1.-PVENV) + DHAIRV*PVENV + LHAR(3) = LHAIRD*(1.-PVENV) + LHAIRV*PVENV + RHOAR(3) = RHOD*(1.-PVENV) + RHOV*PVENV + ZZFUR(3) = ZFURD*(1.-PVENV) + ZFURV*PVENV + REFLFR(3) = REFLD*(1.-PVENV) + REFLV*PVENV DO 9, L=1,3 C INDEX,L,IS THE AVERAGE(1),FRONT(2), BACK(3) OR DORSAL(2), VENTRAL(3) OF THE BODY PART diff --git a/src/RESPFUN.f b/src/RESPFUN.f index 8b45c7cc..6e20d187 100644 --- a/src/RESPFUN.f +++ b/src/RESPFUN.f @@ -71,17 +71,17 @@ SUBROUTINE RESPFUN(TAIREF,X,O2GAS,N2GAS,CO2GAS,BARPRS,QMIN,RQ, P_N2 = RP_N2 P_CO2 = RP_CO2 C ALLOWING USER TO MODIFY GAS VALUES FOR BURROW, ETC. CONDITIONS - IF (P_O2 .NE. (O2GAS/100.))THEN + IF((P_O2.GT.(O2GAS/100.)).OR.(P_O2.LT.(O2GAS/100.)))THEN P_O2 = O2GAS/100. ELSE P_O2 = RP_O2 ENDIF - IF (P_N2 .NE. (N2GAS/100.))THEN + IF((P_N2.GT.(N2GAS/100.)).OR.(P_N2.LT.(N2GAS/100.)))THEN P_N2 = N2GAS/100. ELSE P_N2 = RP_N2 ENDIF - IF (P_CO2 .NE. (CO2GAS/100.))THEN + IF((P_CO2.GT.(CO2GAS/100.)).OR.(P_CO2.LT.(CO2GAS/100.)))THEN P_CO2 = CO2GAS/100. ELSE P_CO2 = RP_CO2 @@ -113,7 +113,7 @@ SUBROUTINE RESPFUN(TAIREF,X,O2GAS,N2GAS,CO2GAS,BARPRS,QMIN,RQ, C OXYGEN CONSUMPTION BASED ON HEAT GENERATION ESTIMATE TO MAINTAIN C BODY TEMPERATURE, CORRECTED FOR SUBSTRATE UTILIZED. C LITERS OF O2/S @ STP: (DATA FOR EQUIVALENCIES FROM KLEIBER, 1961) - IF (RQ .EQ. 1.0) THEN + IF(RQ .GE. 1.0) THEN C CARBOHYDRATE: (ASIDE) CARBOHYDRATES WORTH 4193 CAL/G C LITER(STP)/S = J/S*(CAL/J)*(KCAL/CAL)*(LITER O2/KCAL)) O2STP = RESPGEN*(1./4.185)*(1./1000.)*(1.0/5.057) @@ -209,7 +209,7 @@ SUBROUTINE RESPFUN(TAIREF,X,O2GAS,N2GAS,CO2GAS,BARPRS,QMIN,RQ, C 2.22E-03* C (EDWARDS & HAINES.1978. J. COMP. PHYSIOL. 128: 177-184 IN WELCH, 1980) C FOR A 0.01 KG ANIMAL, THE MAX. RATE WOULD BE 1.67**10^-6 G/S - TESVAL = 2.22E-03*GMASS/1000.*15 + TESVAL = 2.22E-03*GMASS/1000.*15. IF (GEVAP .GT. TESVAL) THEN GEVAP = TESVAL ELSE diff --git a/src/SEVAP_ENDO.f b/src/SEVAP_ENDO.f index f6ee6d2b..a916e7b1 100644 --- a/src/SEVAP_ENDO.f +++ b/src/SEVAP_ENDO.f @@ -75,7 +75,7 @@ SUBROUTINE SEVAP_ENDO(BARPRS,TA,RELHUM,VEL,TC,TSKIN,ALT,SKINW, AIRVD = VD C EYES OPEN - WEYES = HD * (PCTEYES / 100) * CONVSK * (SURFVD - AIRVD) ! KG/S + WEYES = HD * (PCTEYES / 100.) * CONVSK * (SURFVD - AIRVD) ! KG/S C COMPUTE CUTANEOUS WATER LOSS BASED ON AEFF @@ -88,7 +88,7 @@ SUBROUTINE SEVAP_ENDO(BARPRS,TA,RELHUM,VEL,TC,TSKIN,ALT,SKINW, C DIVIDE INTO LOSS FROM FURRED SURFACES AND BARE SURFACES C BARE = FORCED + FREE C PROPORTION WETTED AREA THAT IS FUR - SKINRAF = SKINW * (1 - PCTBAREVAP * 0.01) * 0.01 + SKINRAF = SKINW * (1. - PCTBAREVAP * 0.01) * 0.01 C PROPORTION WETTED AREA THAT IS BARE SKINRAHF = (SKINW * PCTBAREVAP * 0.01) * 0.01 diff --git a/src/SIMULSOL.f b/src/SIMULSOL.f index 8c04646a..867340b8 100644 --- a/src/SIMULSOL.f +++ b/src/SIMULSOL.f @@ -37,7 +37,7 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, DOUBLE PRECISION AK1,AK2,ALT,ASEMAJ,ASQG,BETARA,BG,BL,BP,BR,BS DOUBLE PRECISION BSEMIN,BSQG,CD,CF,CONVAR,CONVRES,CONVSK,CSEMIN - DOUBLE PRECISION CSQG,D,DIFTOL,DV1,DV2,DV3,DV3T1,DV3T2,DV3T3,DV4 + DOUBLE PRECISION CSQG,D,DIFTOL,DV1,DV2,DV3,DV4 DOUBLE PRECISION DV5,EMIS,ENVVARS,FABUSH,FAGRD,FASKY,FATTHK,FAVEG DOUBLE PRECISION FLSHASEMAJ,FLSHBSEMIN,FLSHCSEMIN,FLTYPE,FLYHR DOUBLE PRECISION FURTHRMK,FURTST,FURVARS,FURWET,GEOMVARS,HC,HD @@ -48,7 +48,7 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, DOUBLE PRECISION RSKIN,SEVAPRES,SHAPE,SIG,SKINW,SOLCT,SOLPRO,SSQG DOUBLE PRECISION SUBQFAT,SUCCESS,SURFAR,TA,TBUSH,TC,TCONDSB,TFA DOUBLE PRECISION TFACALC,TFADIFF,TFAT1,TFAT2,TFAT3,TFAT4,TLOWER,TR - DOUBLE PRECISION TRAITS,TRAPPX,TS,TS2T1,TS2T2,TSKCALC,TSKCALC1 + DOUBLE PRECISION TRAITS,TRAPPX,TS,TSKCALC,TSKCALC1 DOUBLE PRECISION TSKCALC2,TSKCALCAV,TSKDIFF,TSKIN,TSKT1,TSKT2 DOUBLE PRECISION TSKT3,TSKY,TVEG,VEL,VOL,XR,ZFUR,ZL DOUBLE PRECISION LHAIR,DHAIR,RHO,REFL,KHAIR @@ -149,7 +149,7 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, C BEGIN CALCULATING TSKIN AND TFA VALUES FOR FURRED BODY PARTS **************************************************************************************************************** 5 CONTINUE - NTRY = NTRY + 1 + NTRY = NTRY + 1. DO 105, I=1,20 11 CONTINUE CALL CONV_ENDO(TS,TA,SHAPE,SURFAR,FLTYPE,FURTST,D,TFA,VEL,ZFUR, @@ -170,76 +170,76 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, KFUR = FURTHRMK ELSE C NEED A TRAD APPROXIMATION FOR CALCULATING KRAD. - TRAPPX = (TSKIN*(1-XR))+(TFA*XR) + TRAPPX = (TSKIN*(1.-XR))+(TFA*XR) KRAD = (16.0*SIG*(TRAPPX+273.15)**3.)/(3.*BETARA(1)) ! EQ7 IN CONLEY AND PORTER 1986 KFUR = KEFF+KRAD ENDIF - IF(IPT.EQ.1)THEN ! CYLINDER GEOMETRY - CF=(PCOND*2*PI*KFURCMPRS*LEN)/(DLOG(RFURCMP/RSKIN)) + IF(INT(IPT).EQ.1)THEN ! CYLINDER GEOMETRY + CF=(PCOND*2.*PI*KFURCMPRS*LEN)/(DLOG(RFURCMP/RSKIN)) IF(PCOND.GT.0.0)THEN TFACMP=(CF*TSKIN+CD*TCONDSB)/(CD+CF) ELSE TFACMP=0.0 ENDIF CD1= (KFURCMPRS/DLOG(RFURCMP/RSKIN))*PCOND+ - & (KFUR/DLOG(RFUR/RSKIN))*(1-PCOND) + & (KFUR/DLOG(RFUR/RSKIN))*(1.-PCOND) CD2= (KFURCMPRS/DLOG(RFURCMP/RSKIN))*PCOND - CD3= (KFUR/DLOG(RFUR/RSKIN))*(1-PCOND) + CD3= (KFUR/DLOG(RFUR/RSKIN))*(1.-PCOND) - DV1=1+((2*PI*LEN*RFLESH**2.*CD1)/(4*AK1*VOL))+ - & ((2*PI*LEN*RFLESH**2.*CD1)/(2*AK2*VOL))* + DV1=1.+((2.*PI*LEN*RFLESH**2.*CD1)/(4.*AK1*VOL))+ + & ((2.*PI*LEN*RFLESH**2.*CD1)/(2.*AK2*VOL))* & DLOG(RSKIN/RFLESH) - DV2=((QSEVAP)*((RFLESH**2.*CD1)/(4*AK1*VOL)))+ - & (QSEVAP)*((RFLESH**2.*CD1)/(2*AK2*VOL))* + DV2=((QSEVAP)*((RFLESH**2.*CD1)/(4.*AK1*VOL)))+ + & (QSEVAP)*((RFLESH**2.*CD1)/(2.*AK2*VOL))* & DLOG(RSKIN/RFLESH) - DV3=(((2*PI*LEN)/DV1)* - & (TC*CD1-DV2-TFACMP*CD2-TFA*CD3)*RFLESH**2.)/(2*VOL) + DV3=(((2.*PI*LEN)/DV1)* + & (TC*CD1-DV2-TFACMP*CD2-TFA*CD3)*RFLESH**2.)/(2.*VOL) IF(XR.LT.1.0)THEN - DV4= CD2+(KFUR/DLOG(RFUR/RRAD))*(1-PCOND) + DV4= CD2+(KFUR/DLOG(RFUR/RRAD))*(1.-PCOND) TR = DV3/DV4+((TFACMP*CD2)/DV4)+ - & (TFA*((KFUR/DLOG(RFUR/RRAD))*(1-PCOND)))/DV4 + & (TFA*((KFUR/DLOG(RFUR/RRAD))*(1.-PCOND)))/DV4 ELSE TR = TFA ENDIF ENDIF - IF(IPT.EQ.2)THEN ! SPHERE GEOMETRY - CF=(PCOND*4*PI*KFURCMPRS*RFURCMP*RSKIN)/(RFURCMP-RSKIN) + IF(INT(IPT).EQ.2)THEN ! SPHERE GEOMETRY + CF=(PCOND*4.*PI*KFURCMPRS*RFURCMP*RSKIN)/(RFURCMP-RSKIN) IF(PCOND.GT.0.0)THEN TFACMP=(CF*TSKIN+CD*TCONDSB)/(CD+CF) ELSE TFACMP=0.0 ENDIF CD1= ((KFURCMPRS*RFURCMP)/(RFURCMP-RSKIN))*PCOND+ - & ((KFUR*RFUR)/(RFUR-RSKIN))*(1-PCOND) + & ((KFUR*RFUR)/(RFUR-RSKIN))*(1.-PCOND) CD2= ((KFURCMPRS*RFURCMP)/(RFURCMP-RSKIN))*PCOND - CD3= ((KFUR*RFUR)/(RFUR-RSKIN))*(1-PCOND) - DV1=1+((4*PI*RSKIN*RFLESH**2.*CD1)/(6*AK1*VOL))+ - & ((4*PI*RSKIN*RFLESH**3.*CD1)/(3*AK2*VOL))* + CD3= ((KFUR*RFUR)/(RFUR-RSKIN))*(1.-PCOND) + DV1=1.+((4.*PI*RSKIN*RFLESH**2.*CD1)/(6.*AK1*VOL))+ + & ((4.*PI*RSKIN*RFLESH**3.*CD1)/(3.*AK2*VOL))* & ((RSKIN-RFLESH)/(RFLESH*RSKIN)) - DV2=((QSEVAP)*((RFLESH**2.*CD1)/(6*AK1*VOL)))+ - & (QSEVAP)*((RFLESH**3.*CD1)/(3*AK2*VOL))* + DV2=((QSEVAP)*((RFLESH**2.*CD1)/(6.*AK1*VOL)))+ + & (QSEVAP)*((RFLESH**3.*CD1)/(3.*AK2*VOL))* & ((RSKIN-RFLESH)/(RFLESH*RSKIN)) - DV3=(((4*PI*RSKIN)/DV1)* - & (TC*CD1-DV2-TFACMP*CD2-TFA*CD3)*RFLESH**3.)/(3*VOL*RRAD) + DV3=(((4.*PI*RSKIN)/DV1)* + & (TC*CD1-DV2-TFACMP*CD2-TFA*CD3)*RFLESH**3.)/(3.*VOL*RRAD) IF(XR.LT.1.0)THEN - DV4= CD2+((KFUR*RFUR)/(RFUR-RRAD))*(1-PCOND) + DV4= CD2+((KFUR*RFUR)/(RFUR-RRAD))*(1.-PCOND) TR = DV3/DV4+((TFACMP*CD2)/DV4)+ - & (TFA*(((KFUR*RFUR)/(RFUR-RRAD))*(1-PCOND)))/DV4 + & (TFA*(((KFUR*RFUR)/(RFUR-RRAD))*(1.-PCOND)))/DV4 ELSE TR = TFA ENDIF ENDIF - IF(IPT.GE.3)THEN ! ELLIPSOID GEOMETRY + IF(IPT.GE.3.)THEN ! ELLIPSOID GEOMETRY FLSHASEMAJ=ASEMAJ-FATTHK FLSHBSEMIN=BSEMIN-FATTHK FLSHCSEMIN=CSEMIN-FATTHK - IF((SUBQFAT.EQ.1.).AND.(FATTHK.GT.0.00))THEN + IF((INT(SUBQFAT).EQ.1).AND.(FATTHK.GT.0.00))THEN ASQG = FLSHASEMAJ**2. BSQG = FLSHBSEMIN**2. CSQG = FLSHCSEMIN**2. @@ -250,7 +250,7 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, ENDIF SSQG = (ASQG*BSQG*CSQG)/(ASQG*BSQG+ASQG*CSQG+BSQG*CSQG) C GETTING THE RADIUS IN THE "B" DIRECTION AT THE FLESH: - IF((SUBQFAT.EQ.1.).AND.(FATTHK.GT.0.00))THEN + IF((INT(SUBQFAT).EQ.1).AND.(FATTHK.GT.0.00))THEN BG=FLSHBSEMIN ELSE BG=BSEMIN @@ -260,8 +260,8 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, C GETTING THE RADIUS IN THE "B" DIRECTION AT THE FUR-AIR INTERFACE: BL=BSEMIN+ZL BR=BS+(XR*ZL) - CF=(PCOND*3*KFURCMPRS*VOL*BLCMP*BS)/ - & ((((3*SSQG)**0.5)**3.)*(BLCMP-BS)) + CF=(PCOND*3.*KFURCMPRS*VOL*BLCMP*BS)/ + & ((((3.*SSQG)**0.5)**3.)*(BLCMP-BS)) IF(PCOND.GT.0.0)THEN TFACMP=(CF*Tskin+CD*TCONDSB)/(CD+CF) ELSE @@ -269,23 +269,23 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, ENDIF CD1= ((KFURCMPRS*BLCMP)/(BLCMP-BS))*PCOND+ - & ((KFUR*BL)/(BL-BS))*(1-PCOND) + & ((KFUR*BL)/(BL-BS))*(1.-PCOND) CD2= ((KFURCMPRS*BLCMP)/(BLCMP-BS))*PCOND - CD3= ((KFUR*BL)/(BL-BS))*(1-PCOND) + CD3= ((KFUR*BL)/(BL-BS))*(1.-PCOND) - DV1=1+((3*BS*SSQG*CD1)/(2*AK1*(((3*SSQG)**0.5)**3.)))+ + DV1=1.+((3*BS*SSQG*CD1)/(2.*AK1*(((3*SSQG)**0.5)**3.)))+ & ((BS*CD1)/(AK2))* & ((BS-BG)/(BS*BG)) - DV2=((QSEVAP)*((SSQG*CD1)/(2*AK1*VOL)))+ - & (QSEVAP)*(((((3*SSQG)**0.5)**3.)*CD1)/(3*AK2*VOL))* + DV2=((QSEVAP)*((SSQG*CD1)/(2.*AK1*VOL)))+ + & (QSEVAP)*(((((3.*SSQG)**0.5)**3.)*CD1)/(3.*AK2*VOL))* & ((BS-BG)/(BS*BG)) DV3=((BS/DV1)* & (TC*CD1-DV2-TFACMP*CD2-TFA*CD3))/(BR) IF(XR.LT.1.0)THEN - DV4= CD2+((KFUR*BL)/(BL-BR))*(1-PCOND) + DV4= CD2+((KFUR*BL)/(BL-BR))*(1.-PCOND) TR = DV3/DV4+((TFACMP*CD2)/DV4)+ - & (TFA*(((KFUR*BL)/(BL-BR))*(1-PCOND)))/DV4 + & (TFA*(((KFUR*BL)/(BL-BR))*(1.-PCOND)))/DV4 ELSE TR = TFA ENDIF @@ -294,75 +294,75 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, TR = TR+273.15 C THESE QR VARIABLES INCORPORATE THE VARIOUS HR VALUES FOR RADIANT EXCHANGE WITH SKY, GROUND, ETC. - QR1=CONVAR*(FASKY*4.*EMIS*SIG*((TR+(TSKY+273.15))/2)**3) - QR2=CONVAR*(FABUSH*4.*EMIS*SIG*((TR+(TBUSH+273.15))/2)**3) - QR3=CONVAR*(FAVEG*4.*EMIS*SIG*((TR+(TVEG+273.15))/2)**3) - QR4=CONVAR*(FAGRD*4.*EMIS*SIG*((TR+(TLOWER+273.15))/2)**3) + QR1=CONVAR*(FASKY*4.*EMIS*SIG*((TR+(TSKY+273.15))/2.)**3) + QR2=CONVAR*(FABUSH*4.*EMIS*SIG*((TR+(TBUSH+273.15))/2.)**3) + QR3=CONVAR*(FAVEG*4.*EMIS*SIG*((TR+(TVEG+273.15))/2.)**3) + QR4=CONVAR*(FAGRD*4.*EMIS*SIG*((TR+(TLOWER+273.15))/2.)**3) IF(PCOND.LT.1)THEN !FOLLOWING CALCULATIONS ARE FOR WHEN THERE IS LESS THAN 100% CONDUCTION C INCLUDES TERM (QFSEVAP) FOR HEAT LOST DUE TO EVAPORATION FROM THE FUR SURFACE TO C CALCULATIONS OF TFA (E.G. WET FUR FROM RAIN) - IF(IPT.EQ.1.)THEN ! CYLINDER GEOMETRY + IF(INT(IPT).EQ.1)THEN ! CYLINDER GEOMETRY IF(XR.LT.1)THEN TFAT1=QR1*TSKY+QR2*TBUSH+QR3*TVEG+QR4*TLOWER- & (QR1+QR2+QR3+QR4)*((DV3/DV4)+((TFACMP*CD2)/DV4)) - TFAT2=((2*PI*LEN)/DV1)*(TC*CD1-DV2-TFACMP*CD2) + TFAT2=((2.*PI*LEN)/DV1)*(TC*CD1-DV2-TFACMP*CD2) TFAT3=HC*CONVAR*TA-CD*TFACMP+CD*TCONDSB-QFSEVAP+QSLR - TFAT4=(2*PI*LEN*CD3)/DV1+(QR1+QR2+QR3+QR4)* - & (((KFUR/DLOG(RFUR/RRAD))*(1-PCOND))/DV4)+HC*CONVAR + TFAT4=(2.*PI*LEN*CD3)/DV1+(QR1+QR2+QR3+QR4)* + & (((KFUR/DLOG(RFUR/RRAD))*(1.-PCOND))/DV4)+HC*CONVAR TFACALC=(TFAT1+TFAT2+TFAT3)/TFAT4 TR2 = DV3/DV4+((TFACMP*CD2)/DV4)+ - & (TFACALC*((KFUR/DLOG(RFUR/RRAD))*(1-PCOND)))/DV4 + & (TFACALC*((KFUR/DLOG(RFUR/RRAD))*(1.-PCOND)))/DV4 ELSE TFAT1=QR1*TSKY+QR2*TBUSH+QR3*TVEG+QR4*TLOWER - TFAT2=((2*PI*LEN)/DV1)*(TC*CD1-DV2-TFACMP*CD2) + TFAT2=((2.*PI*LEN)/DV1)*(TC*CD1-DV2-TFACMP*CD2) TFAT3=HC*CONVAR*TA-CD*TFACMP+CD*TCONDSB-QFSEVAP+QSLR - TFAT4=(2*PI*LEN*CD3)/DV1+(QR1+QR2+QR3+QR4)+HC*CONVAR + TFAT4=(2.*PI*LEN*CD3)/DV1+(QR1+QR2+QR3+QR4)+HC*CONVAR TFACALC=(TFAT1+TFAT2+TFAT3)/TFAT4 TR2 = TFACALC ENDIF ENDIF - IF(IPT.EQ.2.)THEN ! SPHERE GEOMETRY - IF(XR.LT.1)THEN - TFAT1=((4*PI*RSKIN)/DV1)*(TC*CD1-DV2-TFACMP*CD2) + IF(INT(IPT).EQ.2)THEN ! SPHERE GEOMETRY + IF(XR.LT.1.)THEN + TFAT1=((4.*PI*RSKIN)/DV1)*(TC*CD1-DV2-TFACMP*CD2) TFAT2=QR1*TSKY+QR2*TBUSH+QR3*TVEG+QR4*TLOWER- & (QR1+QR2+QR3+QR4)*((DV3/DV4)+((TFACMP*CD2)/DV4)) TFAT3=HC*CONVAR*TA-CD*TFACMP+CD*TCONDSB-QFSEVAP+QSLR - TFAT4=(4*PI*RSKIN*CD3)/DV1+(QR1+QR2+QR3+QR4)* - & ((((KFUR*RFUR)/(RFUR-RRAD))*(1-PCOND))/DV4)+HC*CONVAR + TFAT4=(4.*PI*RSKIN*CD3)/DV1+(QR1+QR2+QR3+QR4)* + & ((((KFUR*RFUR)/(RFUR-RRAD))*(1.-PCOND))/DV4)+HC*CONVAR TFACALC=(TFAT1+TFAT2+TFAT3)/TFAT4 TR2 = DV3/DV4+((TFACMP*CD2)/DV4)+ - & (TFACALC*(((KFUR*RFUR)/(RFUR-RRAD))*(1-PCOND)))/DV4 + & (TFACALC*(((KFUR*RFUR)/(RFUR-RRAD))*(1.-PCOND)))/DV4 ELSE - TFAT1=((4*PI*RSKIN)/DV1)*(TC*CD1-DV2-TFACMP*CD2) + TFAT1=((4.*PI*RSKIN)/DV1)*(TC*CD1-DV2-TFACMP*CD2) TFAT2=QR1*TSKY+QR2*TBUSH+QR3*TVEG+QR4*TLOWER TFAT3=HC*CONVAR*TA-CD*TFACMP+CD*TCONDSB-QFSEVAP+QSLR - TFAT4=(4*PI*RSKIN*CD3)/DV1+(QR1+QR2+QR3+QR4)+HC*CONVAR + TFAT4=(4.*PI*RSKIN*CD3)/DV1+(QR1+QR2+QR3+QR4)+HC*CONVAR TFACALC=(TFAT1+TFAT2+TFAT3)/TFAT4 TR2 = TFACALC ENDIF ENDIF IF(IPT.GE.3.)THEN ! ELLIPSOID GEOMETRY - IF(XR.LT.1)THEN + IF(XR.LT.1.)THEN TFAT1=QR1*TSKY+QR2*TBUSH+QR3*TVEG+QR4*TLOWER- & (QR1+QR2+QR3+QR4)*((DV3/DV4)+((TFACMP*CD2)/DV4)) - TFAT2=((3*VOL*BS)/((((3*SSQG)**0.5)**3.)*DV1))* + TFAT2=((3.*VOL*BS)/((((3.*SSQG)**0.5)**3.)*DV1))* & (TC*CD1-DV2-TFACMP*CD2) TFAT3=HC*CONVAR*TA-CD*TFACMP+CD*TCONDSB-QFSEVAP+QSLR - TFAT4=(3*VOL*BS*CD3)/((((3*SSQG)**0.5)**3.)*DV1)+ + TFAT4=(3.*VOL*BS*CD3)/((((3.*SSQG)**0.5)**3.)*DV1)+ & (QR1+QR2+QR3+QR4)* & ((((KFUR*BL)/(BL-BR))*(1-PCOND))/DV4)+HC*CONVAR TFACALC=(TFAT1+TFAT2+TFAT3)/TFAT4 TR2 = DV3/DV4+((TFACMP*CD2)/DV4)+ - & (TFACALC*(((KFUR*BL)/(BL-BR))*(1-PCOND)))/DV4 + & (TFACALC*(((KFUR*BL)/(BL-BR))*(1.-PCOND)))/DV4 ELSE - TFAT1=((3*VOL*BS)/((((3*SSQG)**0.5)**3.)*DV1))* + TFAT1=((3.*VOL*BS)/((((3.*SSQG)**0.5)**3.)*DV1))* & (TC*CD1-DV2-TFACMP*CD2) TFAT2=QR1*TSKY+QR2*TBUSH+QR3*TVEG+QR4*TLOWER TFAT3=HC*CONVAR*TA-CD*TFACMP+CD*TCONDSB-QFSEVAP+QSLR - TFAT4=(3*VOL*BS*CD3)/((((3*SSQG)**0.5)**3.)*DV1)+ + TFAT4=(3.*VOL*BS*CD3)/((((3.*SSQG)**0.5)**3.)*DV1)+ & (QR1+QR2+QR3+QR4)+HC*CONVAR TFACALC=(TFAT1+TFAT2+TFAT3)/TFAT4 TR2 = TFACALC @@ -377,80 +377,80 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, QCONV = HC*CONVAR*(TFACALC-TA) QCOND = CD*(TFACMP-TCONDSB) ELSE !BEGIN CALCULATIONS FOR 100% CONDUCTION. NEED A DIFFERENT SOLUTION APPROACH SINCE THERE IS NO UNCOMPRESSED TFA TO SOLVE FOR - IF(IPT.EQ.1.)THEN ! CYLINDER GEOMETRY - CF1=(2*PI*KFURCMPRS*LEN)/(DLOG(RFURCMP/RSKIN)) - DV5=1+((CF1*RFLESH**2.)/(4*AK1*VOL))+ - & ((CF1*RFLESH**2.)/(2*AK2*VOL))*DLOG(RSKIN/RFLESH) + IF(INT(IPT).EQ.1)THEN ! CYLINDER GEOMETRY + CF1=(2.*PI*KFURCMPRS*LEN)/(DLOG(RFURCMP/RSKIN)) + DV5=1.+((CF1*RFLESH**2.)/(4.*AK1*VOL))+ + & ((CF1*RFLESH**2.)/(2.*AK2*VOL))*DLOG(RSKIN/RFLESH) ENDIF - IF(IPT.EQ.2.)THEN ! SPHERE GEOMETRY - CF1=(4*PI*KFURCMPRS*RFURCMP)/(RFURCMP-RSKIN) - DV5=1+((CF1*RFLESH**2.)/(6*AK1*VOL))+ - & ((CF1*RFLESH**3.)/(3*AK2*VOL))*((RSKIN-RFLESH)/ + IF(INT(IPT).EQ.2)THEN ! SPHERE GEOMETRY + CF1=(4.*PI*KFURCMPRS*RFURCMP)/(RFURCMP-RSKIN) + DV5=1.+((CF1*RFLESH**2.)/(6.*AK1*VOL))+ + & ((CF1*RFLESH**3.)/(3.*AK2*VOL))*((RSKIN-RFLESH)/ & (RSKIN-RFLESH)) ENDIF - IF(IPT.EQ.3.)THEN ! ELLIPSOID GEOMETRY - CF1=(3*KFURCMPRS*VOL*BLCMP*BS)/ + IF(INT(IPT).EQ.3)THEN ! ELLIPSOID GEOMETRY + CF1=(3.*KFURCMPRS*VOL*BLCMP*BS)/ & ((((3*SSQG)**0.5)**3.)*(BL-BS)) DV5=1+((CF1*SSQG)/(2*AK1*VOL))+ - & ((CF1*(((3*SSQG)**0.5)**3.))/(3*AK2*VOL))*((BS-BG)/ + & ((CF1*(((3*SSQG)**0.5)**3.))/(3.*AK2*VOL))*((BS-BG)/ & (BS*BG)) ENDIF TFACMPT1=(CF1/DV5)*TC+CD*TCONDSB TFACMPT2=CD+CF1/DV5 TFACMP=TFACMPT1/TFACMPT2 - QRAD=0 - QCONV=0 - QFSEVAP = 0 - QSLR = 0 + QRAD=0. + QCONV=0. + QFSEVAP = 0. + QSLR = 0. QCOND=CD*(TFACMP-TCONDSB) ENDIF QENV = QRAD+QCONV+QCOND+QFSEVAP-QSLR - IF(IPT.EQ.1.)THEN ! CYLINDER GEOMETRY - IF(PCOND.LT.1)THEN - TSKCALC1 = TC-(((QENV+QSEVAP)*RFLESH**2.)/(4*AK1*VOL))- - & (((QENV+QSEVAP)*RFLESH**2.)/(2*AK2*VOL))* + IF(INT(IPT).EQ.1)THEN ! CYLINDER GEOMETRY + IF(PCOND.LT.1.)THEN + TSKCALC1 = TC-(((QENV+QSEVAP)*RFLESH**2.)/(4.*AK1*VOL))- + & (((QENV+QSEVAP)*RFLESH**2.)/(2.*AK2*VOL))* & DLOG(RSKIN/RFLESH) - TSKCALC2 = ((QENV*RFLESH**2.)/(2*CD1*VOL))+((TFACMP*CD2)/CD1)+ + TSKCALC2 =((QENV*RFLESH**2.)/(2.*CD1*VOL))+((TFACMP*CD2)/CD1)+ & ((TFACALC*CD3)/CD1) ELSE - TSKCALC1=TC-((QENV*RFLESH**2.)/(4*AK1*VOL))- - & ((QENV*RFLESH**2.)/(2*AK2*VOL))* + TSKCALC1=TC-((QENV*RFLESH**2.)/(4.*AK1*VOL))- + & ((QENV*RFLESH**2.)/(2.*AK2*VOL))* & DLOG(RSKIN/RFLESH) - TSKCALC2=(((QENV*RFLESH**2.)/(2*KFURCMPRS*VOL))* + TSKCALC2=(((QENV*RFLESH**2.)/(2.*KFURCMPRS*VOL))* & DLOG(RFURCMP/RSKIN))+TFACMP ENDIF ENDIF - IF(IPT.EQ.2.)THEN ! SPHERE GEOMETRY - IF(PCOND.LT.1)THEN - TSKCALC1 = TC-(((QENV+QSEVAP)*RFLESH**2.)/(6*AK1*VOL))- - & (((QENV+QSEVAP)*RFLESH**3.)/(3*AK2*VOL))* + IF(INT(IPT).EQ.2)THEN ! SPHERE GEOMETRY + IF(PCOND.LT.1.)THEN + TSKCALC1 = TC-(((QENV+QSEVAP)*RFLESH**2.)/(6.*AK1*VOL))- + & (((QENV+QSEVAP)*RFLESH**3.)/(3.*AK2*VOL))* & ((RSKIN-RFLESH)/(RSKIN*RFLESH)) - TSKCALC2 = ((QENV*RFLESH**3.)/(3*CD1*VOL*RSKIN))+ + TSKCALC2 = ((QENV*RFLESH**3.)/(3.*CD1*VOL*RSKIN))+ & ((TFACMP*CD2)/CD1)+((TFACALC*CD3)/CD1) ELSE - TSKCALC1=TC-((QENV*RFLESH**2.)/(6*AK1*VOL))- - & ((QENV*RFLESH**3.)/(3*AK2*VOL))* + TSKCALC1=TC-((QENV*RFLESH**2.)/(6.*AK1*VOL))- + & ((QENV*RFLESH**3.)/(3.*AK2*VOL))* & ((RSKIN-RFLESH)/(RSKIN*RFLESH)) - TSKCALC2=(((QENV*RFLESH**3.)/(3*KFURCMPRS*VOL))* + TSKCALC2=(((QENV*RFLESH**3.)/(3.*KFURCMPRS*VOL))* & ((RFURCMP-RSKIN)/(RFURCMP*RSKIN)))+TFACMP ENDIF ENDIF - IF(IPT.GE.3.)THEN ! ELLIPSOID GEOMETRY - IF(PCOND.LT.1)THEN - TSKCALC1 = TC-(((QENV+QSEVAP)*SSQG)/(2*AK1*VOL))- - & (((QENV+QSEVAP)*(((3*SSQG)**0.5)**3.))/(3*AK2*VOL))* + IF(INT(IPT).GE.3)THEN ! ELLIPSOID GEOMETRY + IF(PCOND.LT.1.)THEN + TSKCALC1 = TC-(((QENV+QSEVAP)*SSQG)/(2.*AK1*VOL))- + & (((QENV+QSEVAP)*(((3.*SSQG)**0.5)**3.))/(3.*AK2*VOL))* & ((BS-BG)/(BS*BG)) - TSKCALC2 = ((QENV*(((3*SSQG)**0.5)**3.))/(3*CD1*VOL*BS))+ + TSKCALC2 = ((QENV*(((3.*SSQG)**0.5)**3.))/(3.*CD1*VOL*BS))+ & ((TFACMP*CD2)/CD1)+((TFACALC*CD3)/CD1) ELSE - TSKCALC1=TC-((QENV*SSQG)/(2*AK1*VOL))- - & ((QENV*(((3*SSQG)**0.5)**3.))/(3*AK2*VOL))* + TSKCALC1=TC-((QENV*SSQG)/(2.*AK1*VOL))- + & ((QENV*(((3.*SSQG)**0.5)**3.))/(3.*AK2*VOL))* & ((BS-BG)/(BS*BG)) - TSKCALC2=(((QENV*(((3*SSQG)**0.5)**3.))/(3*KFURCMPRS*VOL))* + TSKCALC2=(((QENV*(((3.*SSQG)**0.5)**3.))/(3.*KFURCMPRS*VOL))* & ((BLCMP-BS)/(BLCMP*BS)))+TFACMP ENDIF ENDIF @@ -466,11 +466,11 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, GO TO 16 ENDIF C IF NO, TRY ANOTHER INITIAL TFA GUESS - IF(SOLPRO.EQ.1.)THEN + IF(INT(SOLPRO).EQ.1)THEN C FIRST SOLUTION PROCEDURE IS TO SET TFA GUESS TO THE CALCULATED TFA TFA=TFACALC ELSE - IF(SOLPRO.EQ.2.)THEN + IF(INT(SOLPRO).EQ.2)THEN C SECOND SOLUTION PROCEDURE IS TO SET TFA GUESS TO AVERAGE OF PREVIOUS GUESS C AND CALCULATED TFA TFA=(TFACALC+TFA)/2. @@ -522,12 +522,12 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, SOLCT=SOLCT+1. IF(SOLCT.GE.100.)THEN - IF(SOLPRO.NE.3.0)THEN + IF(INT(SOLPRO).NE.3)THEN SOLCT=0. - SOLPRO=SOLPRO+1 + SOLPRO=SOLPRO+1. ELSE C EVEN THE SECOND WAY OF SOLVING FOR BALANCE DOESN'T WORK, INCREASE TOLERANCE - IF(DIFTOL.EQ.0.001)THEN + IF(DIFTOL.LE.0.001)THEN DIFTOL=0.01 SOLCT=0. SOLPRO=1. @@ -544,17 +544,17 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, 16 IF(TSKDIFF.LT.DIFTOL)THEN C IF YES, BOTH TFA AND TSK GUESSES ARE SIMILAR TO THE CALCULATED VALUES - IF(IPT.EQ.1.)THEN ! CYLINDER GEOMETRY + IF(INT(IPT).EQ.1)THEN ! CYLINDER GEOMETRY QGENNET = (TC-TSKCALCAV)/((RFLESH**2./(4.*AK1*VOL))+ & ((RFLESH**2./(2.*AK2*VOL))*LOG(RSKIN/RFLESH))) ENDIF - IF(IPT.EQ.2.)THEN ! SPHERE GEOMETRY + IF(INT(IPT).EQ.2)THEN ! SPHERE GEOMETRY QGENNET = (TC-TSKCALCAV)/((RFLESH**2./(6.*AK1*VOL))+ & ((RFLESH**3./(3.*AK2*VOL))*((RSKIN-RFLESH)/(RFLESH*RSKIN)))) ENDIF - IF(IPT.GE.3.)THEN ! ELLIPSOID GEOMETRY + IF(INT(IPT).GE.3)THEN ! ELLIPSOID GEOMETRY QGENNET = (TC-TSKCALCAV)/((SSQG/(2.*AK1*VOL))+ & (((((3.*SSQG)**0.5)**3.)/(3.*AK2*VOL))* & ((BS-BG)/(BG*BS)))) @@ -567,17 +567,17 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, GO TO 5 ELSE SUCCESS=0. - IF(IPT.EQ.1.)THEN ! CYLINDER GEOMETRY + IF(INT(IPT).EQ.1)THEN ! CYLINDER GEOMETRY QGENNET = (TC-TSKCALCAV)/((RFLESH**2./(4.*AK1*VOL))+ & ((RFLESH**2./(2.*AK2*VOL))*LOG(RSKIN/RFLESH))) ENDIF - IF(IPT.EQ.2.)THEN ! SPHERE GEOMETRY + IF(INT(IPT).EQ.2)THEN ! SPHERE GEOMETRY QGENNET = (TC-TSKCALCAV)/((RFLESH**2./(6.*AK1*VOL))+ & ((RFLESH**3./(3.*AK2*VOL))*((RSKIN-RFLESH)/(RFLESH*RSKIN)))) ENDIF - IF(IPT.GE.3.)THEN ! ELLIPSOID GEOMETRY + IF(INT(IPT).GE.3)THEN ! ELLIPSOID GEOMETRY QGENNET = (TC-TSKCALCAV)/((SSQG/(2.*AK1*VOL))+ & (((((3*SSQG)**0.5)**3.)/(3.*AK2*VOL))* & ((BS-BG)/(BG*BS)))) @@ -590,7 +590,7 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, C COMMENT: THIS LOOP IS FOR WHEN THERE IS NO FUR. 120 CONTINUE - NTRY = NTRY + 1 + NTRY = NTRY + 1. DO 140, I=1,20 125 CONTINUE CALL CONV_ENDO(TS,TA,SHAPE,CONVSK,FLTYPE,FURTST,D,TFA,VEL,ZFUR, @@ -605,13 +605,13 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, C THESE QR VARIABLES INCORPORATE THE VARIOUS HR VALUES FOR RADIANT EXCHANGE WITH SKY, GROUND, ETC. QR1=CONVSK*(FASKY*4.*EMIS*SIG* - & (((TSKIN+273.15)+(TSKY+273.15))/2)**3) + & (((TSKIN+273.15)+(TSKY+273.15))/2.)**3.) QR2=CONVSK*(FABUSH*4.*EMIS*SIG* - & (((TSKIN+273.15)+(TBUSH+273.15))/2)**3) + & (((TSKIN+273.15)+(TBUSH+273.15))/2.)**3.) QR3=CONVSK*(FAVEG*4.*EMIS*SIG* - & (((TSKIN+273.15)+(TVEG+273.15))/2)**3) + & (((TSKIN+273.15)+(TVEG+273.15))/2.)**3.) QR4=CONVSK*(FAGRD*4.*EMIS*SIG* - & (((TSKIN+273.15)+(TLOWER+273.15))/2)**3) + & (((TSKIN+273.15)+(TLOWER+273.15))/2.)**3.) TSKT1= ((4.*AK1*VOL)/(RSKIN**2.)*TC)-QSEVAP+HC*CONVSK* & TA+QSLR @@ -639,9 +639,9 @@ SUBROUTINE SIMULSOL(DIFTOL,IPT,FURVARS,GEOMVARS,ENVVARS,TRAITS, TSKIN=TSKCALC TSKCALCAV=TSKCALC TFA=TSKCALC - NTRY=NTRY+1 - IF(NTRY.EQ.101.)THEN - IF(DIFTOL.EQ.0.001)THEN + NTRY=NTRY+1. + IF(INT(NTRY).EQ.101)THEN + IF(DIFTOL.LE.0.001)THEN DIFTOL=0.01 NTRY=0. ELSE diff --git a/src/SOLVENDO.f b/src/SOLVENDO.f index 98988318..85179b88 100644 --- a/src/SOLVENDO.f +++ b/src/SOLVENDO.f @@ -57,6 +57,7 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) PI = ACOS(-1.0d0) ZBRENTout=0. + QRESP=0. QGEN=input(1) QBASAL=input(2) TA=input(3) @@ -150,6 +151,34 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) TSKINMAX=TC ! initialise Q10mult=1. ! initialise PANT_COST = 0.D0 ! initialise + D=0. + AHEIT=0. + AK1_LAST=0. + ALENTH=0. + AREACND=0. + VOL=0. + VMULT=0. + TLUNG=0. + TC_LAST=0. + SHAPEB_LAST=0. + R2=0. + R1=0. + PCTWET_LAST=0. + PANT_LAST=0. + MASFAT=0. + KFURCMPRS=0. + FLSHVL=0. + FATTHK=0. + FAGRD=0. + DMULT=0. + CONVSK=0. + CONVAR=0. + AWIDTH=0. + ASILP=0. + ASILN=0. + AREASKIN=0. + AREACND=0. + KEFARA = (/0.,0.,0./) do while(QGEN < QBASAL) @@ -217,8 +246,8 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) QNORM = QSOLR endif - ABSAND = 1 - REFLFR(2) !# solar absorptivity of dorsal fur (fractional, 0-1) - ABSANV = 1 - REFLFR(3) !# solar absorptivity of ventral fur (fractional, 0-1) + ABSAND = 1. - REFLFR(2) !# solar absorptivity of dorsal fur (fractional, 0-1) + ABSANV = 1. - REFLFR(3) !# solar absorptivity of ventral fur (fractional, 0-1) C CORRECT FASKY FOR % VEGETATION SHADE OVERHEAD, ASHADE FAVEG = FSKREF*(SHADE/100.) @@ -310,7 +339,7 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) !# set fur depth and conductivity !# index for KEFARA, the conductivity, is the average (1), front/dorsal (2), back/ventral(3) of the body part - if((QSOLR.GT.0).OR.(ZZFUR(2).NE.ZZFUR(3)))THEN + if((QSOLR.GT.0).OR.(INT(ZZFUR(2)).NE.INT(ZZFUR(3))))THEN if(S == 1)THEN ZL = ZZFUR(2) KEFF = KEFARA(2) @@ -328,13 +357,13 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) RFUR = R1 + ZL !# body radius including fur, m D = 2. * RFUR !# diameter, m RRAD = RSKIN + (XR * ZL) !# effective radiation radius, m - IF(SHAPE.NE.4)THEN ! For cylinder and sphere geometries + IF(INT(SHAPE).NE.4)THEN ! For cylinder and sphere geometries RFURCMP=RSKIN+ZFURCOMP ELSE RFURCMP=RFUR ! Note that this value is never used if conduction not being modeled, but need to have a value for the calculations ENDIF - IF(SHAPE.EQ.4)THEN ! For ellipsoid geometry + IF(INT(SHAPE).EQ.4)THEN ! For ellipsoid geometry BLCMP=BSEMIN+ZFURCOMP ELSE BLCMP=RFUR ! Note that this value is never used if conduction not being modeled, but need to have a value for the calculations @@ -343,13 +372,13 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) LEN = ALENTH !# length, m !# Correcting volume to account for subcutaneous fat - if((SUBQFAT.EQ.1.).AND.(FATTHK.GT.0.0))THEN + if((INT(SUBQFAT).EQ.1).AND.(FATTHK.GT.0.0))THEN VOL = FLSHVL ENDIF !# Calculating the "Cd" variable: Qcond = Cd(Tskin-Tsub), where Cd = Conduction area*ksub/subdepth IF(S==2)THEN ! doing ventral side, add conduction - AREACND = ATOT * (PCOND * 2) + AREACND = ATOT * (PCOND * 2.) CD = (AREACND * KSUB) / 0.025 !# assume conduction happens from 2.5 cm depth ELSE !# doing dorsal side, no conduction. No need to adjust areas used for convection. AREACND = 0. @@ -368,13 +397,13 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) & PCTEYES/) !# set IPT, the geometry assumed in SIMULSOL: 1 = cylinder, 2 = sphere, 3 = ellipsoid - if((SHAPE.eq.1.).or.(SHAPE.eq.3.))THEN + if((INT(SHAPE).eq.1).or.(INT(SHAPE).eq.3))THEN IPT = 1. ENDIF - if(SHAPE.eq.2.)THEN + if(INT(SHAPE).eq.2)THEN IPT = 2. ENDIF - If(SHAPE.eq.4.)THEN + If(INT(SHAPE).eq.4)THEN IPT = 3. ENDIF @@ -408,7 +437,7 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) TLUNG =(TC + TS) * 0.5 !# average of skin and core TAEXIT = min(TA + DELTAR, TLUNG) !# temperature of exhaled air, deg C - IF(RESPIRE.EQ.1.)THEN + IF(INT(RESPIRE).EQ.1)THEN !# now guess for metabolic rate that balances the heat budget while allowing metabolic rate !# to remain at or above QBASAL, via 'shooting method' ZBRENT QMIN = QBASAL @@ -436,7 +465,7 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) TC_LAST = TC PANT_LAST = PANT PCTWET_LAST = PCTWET - IF(THERMOREG.EQ.1)THEN + IF(INT(THERMOREG).EQ.1)THEN if(SHAPEB.lt.SHAPEB_MAX)THEN SHAPEB = SHAPEB + UNCURL else @@ -454,14 +483,14 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) Q10mult = Q10**((TC - TC_REF)/10.) if(PANT.lt.PANT_MAX)THEN PANT = PANT + PANT_INC - PANT_COST=((PANT-1D0)/(PANT_MAX-1D0)*(PANT_MULT-1)*QBASREF) + PANT_COST=((PANT-1.)/(PANT_MAX-1.)*(PANT_MULT-1.)*QBASREF) QBASAL = QBASREF * Q10mult + PANT_COST else PANT = PANT_MAX - PANT_COST=((PANT-1D0)/(PANT_MAX-1D0)*(PANT_MULT-1)*QBASREF) + PANT_COST=((PANT-1.)/(PANT_MAX-1.)*(PANT_MULT-1.)*QBASREF) QBASAL = QBASREF * Q10mult + PANT_COST PCTWET = PCTWET + PCTWET_INC - if((PCTWET.GT.PCTWET_MAX).OR.(PCTWET_INC.eq.0.))THEN + if((PCTWET.GT.PCTWET_MAX).OR.(INT(PCTWET_INC).eq.0))THEN PCTWET = PCTWET_MAX RETURN ENDIF @@ -526,17 +555,17 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) HTOVPR = 2.5012E+06 - 2.3787E+03 * TA SWEATGS = (SIMULSOLout(1,6) + SIMULSOLout(2,6)) * 0.5 - & / HTOVPR * 1000 + & / HTOVPR * 1000. EVAPGS = ZBRENTout(3) + SWEATGS HTOVPR=2.5012E+06 - 2.3787E+03 * TA ! latent heat of vapourisation, W/kg/C - SWEATGS=(QSEVAPD + QSEVAPV) * 0.5 / HTOVPR * 1000 ! water lost from skin, g/s + SWEATGS=(QSEVAPD + QSEVAPV) * 0.5 / HTOVPR * 1000. ! water lost from skin, g/s EVAPGS=GEVAP + SWEATGS ! total evaporative water loss, g/s sigma=5.6697E-8 QIR=QIRIN - QIROUT - QIROUTD=sigma * EMISAN * AREASKIN * (TSKCALCAVD + 273.15)**4 + QIROUTD=sigma * EMISAN * AREASKIN * (TSKCALCAVD + 273.15)**4. QIRIND=QRADD * (-1.) + QIROUTD - QIROUTV=sigma * EMISAN * AREASKIN * (TSKCALCAVD + 273.15)**4 + QIROUTV=sigma * EMISAN * AREASKIN * (TSKCALCAVD + 273.15)**4. QIRINV=QRADV * (-1.) + QIROUTV QSOL=QSLRD * DMULT + QSLRV * VMULT ! solar, W @@ -562,7 +591,7 @@ subroutine SOLVENDO(INPUT,TREG,MORPH,ENBAL,MASBAL) !names(enbal)=c("QSOL", "QIRIN", "QGEN", "QEVAP", "QIROUT", "QCONV", "QCOND", "ENB", "NTRY", "SUCCESS") MASBAL=(/AIRVOL, O2STP, GEVAP, SWEATGS, O2MOL1, - & O2MOL2, N2MOL1, N2MOL2, AIRML1, AIRML2/) * 3600 + & O2MOL2, N2MOL1, N2MOL2, AIRML1, AIRML2/) * 3600. !names(masbal)=c("AIR_L", "O2_L", "H2OResp_g", "H2OCut_g", "O2_mol_in", "O2_mol_out", "N2_mol_in", "N2_mol_out", "AIR_mol_in", "AIR_mol_out") RETURN diff --git a/src/WETAIR.f b/src/WETAIR.f index 80f282af..afaeb0dd 100644 --- a/src/WETAIR.f +++ b/src/WETAIR.f @@ -93,6 +93,6 @@ SUBROUTINE WETAIR(DB,WB,RH,DP,BP,E,ESAT,VD,RW,TVIR,TVINC,DENAIR, IF (RH .LE. 0.0) GO TO 500 WTRPOT = 4.615E+5 * TK * DLOG(RH / 100.0) GO TO 600 -500 WTRPOT = -999 +500 WTRPOT = -999. 600 RETURN END \ No newline at end of file