search

EDmodel / ED2

Ecosystem Demography Model
78 stars 112 forks source link

R utility to check met forcings #119

Open fabeit opened 9 years ago

fabeit commented 9 years ago

I have written a little R script to check if met variables in hdf5 files are within boundaries. This avoids crashing a run after some time because something is out of bounds. If it was me I would run this script in the initialization of every ED run. Would this be of interest to be added to the repository and in case where do I add it?

mdietze commented 9 years ago

Yes of interest. Presumably should be dropped into ED2/R-utils.

fabeit commented 9 years ago

@mdietze given that some of the met forcings are calculated at run time so a crash can't be prevented, wouldn't it make sense to add these checks before every ED run?

mdietze commented 9 years ago

There's nothing in ED2 that I'm aware of that's changing the contents of the hdf5 met files, so I'm not sure what you mean by "some of the met forcings are calculated at run time".

fabeit commented 9 years ago

@mdietze I was referring to the interpolation of the short-wave radiation, isn't that done in real time?

mdietze commented 9 years ago

All input met variables (not just SW) are interpolated in real time according to the rules you specify in the ED met header file, but that's done internally in ED and doesn't change the content of the hdf5 met files, so I'm still not sure I see why running an R script on the hdf5 met files every time at run time is preferable to running the check once when you build the met inputs.

Don't get me wrong, the script is still of interest since that check should be round at build time, so please do submit a pull request.

fabeit commented 9 years ago

@mdietze that makes perfect sense, I will put this script in R utils folder.

gklarenberg commented 7 years ago

@fabeit @mdietze What is this file called in R-utils...? I am having the strangest problem: ED2 ran fine yesterday evening (freshly downloaded and processed meteo data), and all of a sudden today it halts at 20 Feb 1996 with the message that my Total SW radiation does not make sense (reported to be 1.50798E+03, higher than the ok max). I haven't changed my inputs, but I want to check them anyway. The only thing I changed in ED2IN is to not write all the detailed budget and photosynthesis outputs etc: I changed IDETAILED=63 to IDETAILED=0....

fabeit commented 7 years ago

@gklarenberg I have not uploaded this file in fact, was your crash always happening or only after you have changed some run options? Here is the R script, it will read h5 files, but you will need to modify it a bit to load your files because your files might have different names.

sink("bad data.txt") library(rhdf5) l=list.files(pattern="OL1") l2=list.files(pattern="OL2") for (i in 1:length(l)) { cat("working on", l[i],"\n")

longwave rad

b=h5read(l[i],"dlwrf") if(min(b)< 40 || max(b)> 600) cat("DLWR out of range, min&max",l[i],min(b),max(b),"\n")

precip rate [kg/m2/s]

b=h5read(l[i],"prate") if(min(b)< 0 || max(b)> 0.1111) cat("prec rate out of range, min&max",l[i],min(b),max(b),"\n")

pressure [Pa]

b=h5read(l[i],"pres") if(min(b)< 45000 || max(b)> 110000) cat("pressure out of range, min&max",l[i],min(b),max(b),"\n")

spec humidity [Kg_h2o/Kg_air]

b=h5read(l[i],"sh") if(min(b)< 0.000001|| max(b)> 0.032) cat("spec humidity out of range, min&max",l[i],min(b),max(b),"\n")

temperature [K]

b=h5read(l[i],"tmp") if(min(b)< 260 || max(b)> 321) cat("pressure rate out of range, min&max",l[i],min(b),max(b),"\n")

ugrd&vgrd and wind speed

u=h5read(l[i],"ugrd") if(min(u)< -20 || max(u)> 20) cat("u wind component out of range, min&max,",l[i],min(u),max(u),"\n")

v=h5read(l[i],"vgrd") if(min(v)< -20 || max(v)> 20) cat("v wind component out of range, min&max,",l[i],min(v),max(v),"\n")

speed=sqrt(u^2+v^2) if(min(speed)< 0 || max(speed)> 50) cat("wind speed out of range, min&max,",l[i],min(speed),max(speed),"\n")

Shortwave Nera-IR beam radiation

nb=h5read(l2[i],"nbdsf") if(min(nb)< 0 || max(nb)> 1500) cat("Near-IR beam radiation out of range, min&max,",l[i],min(nb),max(nb),"\n")

Shortwave Near-IR diffuse radiation

nd=h5read(l2[i],"nddsf") if(min(nd)< 0 || max(nd)> 1500) cat("Near-IR beam radiation out of range, min&max,",l[i],min(nd),max(nd),"\n")

Shortwave Visible beam radiation

vb=h5read(l2[i],"vbdsf") if(min(vb)< 0 || max(vb)> 1500) cat("Visible diffuse radiation out of range, min&max,",l[i],min(vb),max(vb),"\n")

Shortwave Visible diffuse radiation

vd=h5read(l2[i],"vddsf") if(min(vd)< 0 || max(vd)> 1500) cat("Visible diffuse radiation out of range, min&max,",l[i],min(vd),max(vd),"\n")

Total shortwave

swt=vd+vb+nd+nb if(min(swt)< 0 || max(swt)> 1500) cat("Total shortwave radiation out of range, min&max,",l[i],min(swt),max(swt),"\n")

} sink()