Update calcrn_vect.m

calcrn_vect.m

 function [albs,arns,brns,crns,arnv,brnv,crnv,cras,crav]=calcrn_vect(rg,ratm,emiss,emisv,albe,albv,sigmav,sigmava)
 
 
-    albs=(albe-sigmav*albv)./((1-sigmav).^2+sigmav*albv.*albe-(sigmav*albv).^2);
+   albs=(albe-sigmav*albv)./((1-sigmav).^2+sigmav*albv.*albe-(sigmav*albv).^2);
+    i=or(albs>0.4,albs<0.05);
+
+    if i==1
+    albs=max(min(albs,0.05),0.4);
+    a=albs*sigmav^2;
+    b=-sigmav*(1+albs*albe);
+    c=albe-albs*(1-sigmav)^2;
+
+    f = @(x) a*x^2+b*x+c;
+
+    lbv=fzero(f,albe);
+    albv=min(max(albv,0.05),0.3);
+    end
+
+    albe=sigmav*albv+(albs*(1-sigmav)^2)/(1-sigmav*albs*albv);
 
 	v1=1-albv*albs.*sigmav;
     v1a=1-albv*albs.*sigmava;
@@ -22,4 +37,4 @@ function [albs,arns,brns,crns,arnv,brnv,crnv,cras,crav]=calcrn_vect(rg,ratm,emis
     crnv=rg*(1-albv)*sigmav.*(1+(albs.*v4)./v1)+sigmava*emisv*ratm.*(1+(v4a.*v3)./v5a);
     crav=sigmava*emisv*ratm.*(1+(v4a.*v3)./v5a);
 
-end
\ No newline at end of file
+end