Test page. Here is a simple model, in Fortran 90, to compute radiative fluxes (PAR, NIR, TIR) in a plant canopy). If one is interested, I can explain it (it was written in haste, undocumented). It was also created in Compaq Visual Fortran 6.0, which accounts for the weird indenting when viewed in text mode:

The "\" before the comment delimiter, "!", is to prevent pmwiki from making this a major heading

\! canopy_abs+scatt.f90 ~/fortran/radtpt 21-22 Feb. 09

      real K,ks,I,I0,Kdiff,Kdir
      write(6,'("aleaf,sleaf,xf,D0,I0,fleaf,Kdiff,Kdir,r")')
	  read(5,*)aleaf,sleaf,xf,D0,I0,fleaf,Kdiff,Kdir,r
	  write(6,'("Recap: aleaf, sleaf, xf, D0, I0, fleaf, Kdiff, Kdir, r=",/,2x,&
	 & 9f8.3)')aleaf,sleaf,xf,D0,I0,fleaf,Kdiff,Kdir,r
	  a=aleaf*Kdiff
	  s=sleaf*Kdiff
	  adir=aleaf*Kdir
	  sdir=sleaf*Kdir
	  fdir=fleaf*Kdir
	  part1=-(1.+a/s)
	  K=sqrt(2.*a*s+a*a)
	  bplus=part1+K/s
	  bminus=part1-K/s
	  ks=adir+sdir+fdir
	  Eplus=(exp(K*xf)-exp(-ks*xf))/(K+ks)
	  Eminus=(exp(-ks*xf)-exp(-K*xf))/(K-ks)
	  fac=exp(-K*xf)
	  P1=(fdir-bminus*sdir)*I0*Eplus*fac
	  ratio=(1.+r*bminus)/(1.+r*bplus)
	  P2=(fdir-bplus*sdir)*I0*Eminus*ratio*fac
	  P3=(bminus-bplus)*r*I0*exp(-ks*xf)*fac/(1.+r*bplus)
	  denom=bminus-bplus*exp(-2.*K*xf)*(1.+r*bminus)/(1.+r*bplus)
	  write(6,'("K,bplus,bminus,ks,Eplus,Eminus=",/,2x,6f9.3)')&
	 &  K,bplus,bminus,ks,Eplus,Eminus
	  write(6,'("fac,P1,ratio,P2,P3,denom=",/,2x,6f10.3)')&
	 & fac,P1,ratio,P2,P3,denom
	  Aplus0=((bminus-bplus)*D0+bplus*(P2-P1+P3))/(bminus-bplus*fac*fac*ratio)
	  Aminus0=(Aplus0*bminus-(bminus-bplus)*D0)/bplus
	  U0=(Aminus0-Aplus0)/(bminus-bplus)
	  write(6,'("Aplus0=",f8.3,", Aminus0=",f8.3,", U0=",f8.3)')&
	 & Aplus0,Aminus0,U0
	  dx=0.01*xf
	  abs=0.
	  do  x=0.5*dx,xf-0.5*dx,dx
	    Aplus=Aplus0*exp(-K*x)+(fdir-bplus*sdir)*I0*(exp(-ks*x)-exp(-K*x))/(K-ks)
		Aminus=Aminus0*exp(K*x)+(fdir-bminus*sdir)*I0*(exp(K*x)-exp(-ks*x))/(K+ks)
	    D=(bminus*Aplus-bplus*Aminus)/(bminus-bplus)
		U=(Aminus-Aplus)/(bminus-bplus)
		I=I0*exp(-ks*x)
		absx=a*(D+U)+adir*I
		abs=abs+absx
		write(6,'("  x=",f6.3,": Aplus=",f8.3,", Aminus=",f8.3,/,&
	   & 4x,"D=",f8.3,", U=",f8.3,", I=",f8.3,", absrate=",f8.3)')&
	   & x,Aplus,Aminus,D,U,I,absx
	  enddo

\! Accounting

      abs_can=abs*dx
      Refl_toc=U0
	  Aplusxf=Aplus0*fac+(fdir-bplus*sdir)*I0*Eminus
	  Aminusxf=Aminus0*exp(K*xf)+(fdir-bminus*sdir)*I0*Eplus
	  Dxf=(bminus*Aplusxf-bplus*Aminusxf)/(bminus-bplus)
	  write(6,'("Aplusxf=",f8.3,", Aminusxf=",f8.3,", Dxf=",f8.3)')&
	 & Aplusxf,Aminusxf,Dxf
	  abssoil=(Dxf+I0*exp(-ks*xf))*(1.-r)
	  Totalout=Refl_toc+abs+abssoil
	  Totalin=I0+D0
	  write(6,'("Input: I0+D0 = ",f8.3)')Totalin
	  write(6,'("Outputs:",/,3x,"Refl_toc=",f8.3,/,3x,"abs_can=",f8.3,&
	 & /,3x,"abs_soil=",f8.3,/,"Total outputs=",f8.3)')Refl_toc,abs_can,abssoil,&
	 & Refl_toc+abs_can+abssoil
	  stop
	  end