;	PRIMARY PRODUCTION MODEL TESTING
;	 USING MIKE BEHRENFELD'S MODEL

pro model,Chl_file,Temp_file,Irr_file,Julian_date,Chl_sat

win_xsize=2048
win_ysize=1024
win_xsize1=win_xsize-1
win_ysize1=win_ysize-1

;	Input the Chlorophyll data
;
Chl_sat=fltarr(win_xsize,win_ysize)
openr,1,Chl_file
readu,1,Chl_sat
close,1
Chl_sat=rotate(Chl_sat,7)
;tvscl,rebin(Chl_sat,win_xsize/4,win_ysize/4,/sample)

;	Input the Sea Surface Temperature data
;
Temp_surf=fltarr(win_xsize,win_ysize)
openr,1,Temp_file
readu,1,Temp_surf
close,1
;tvscl,rebin(Temp_surf,win_xsize/4,win_ysize/4,/sample)

;	Input the Day Length data
;
Daytime=fltarr(win_ysize,366)
Daylength=fltarr(win_ysize)
Day_file='day_len_1024.dat'
;read,'File name for Daylength data:',Day_file
openr,1,Day_file
readu,1,Daytime
close,1
Daylength(*)=Daytime(*,Julian_date)
;tvscl,rebin(Daytime,win_ysize/2,366/2,/sample)
Daytime=0



;	Input the Surface Irradiance data
;
;Irr_surf=fltarr(win_ysize,366)
Irr_surf=fltarr(win_xsize,win_ysize)
;Irr_file='irradiance_1024.dat'
;read,'File name for Surface Irradiance data:',Irr_file
openr,1,Irr_file
readu,1,Irr_surf
close,1
;tvscl,rebin(Irr_surf,win_xsize/4,win_ysize/4,/sample)

;	Input the Julian day of the year
;
;Julian_date=15
;read,'Julian date of the year:',Julian_date

;	Read in the name of the Primary Production Output file
;
;;Prim_Prod=fltarr(win_xsize,win_ysize)
;PP_file='PP_Jan.tst'
;read,'Enter the Output file name for Prim. Prod.:',PP_file

;	Definition of other variables
;
Chl_tot=0.0
Z_eu=0.0
Pb_opt=4.54
Temp_surf2=0.
sqr2=sqrt(2.)
sqr2n1=sqr2+1.
;;Pb_opt_total=0.
;;Pb_opt_count=0.

c0=1.2956
c1=2.75 ;c1=0.2749
c2=6.17 ;c2=6.17e-2
c3=-20.5 ;c3=-2.05e-2
c4=24.62 ;c4=2.46e-3
c5=-13.48 ;c5=-1.35e-4
c6=3.4132 ;c6=3.42e-6
c7=-0.327 ;c7=-3.28e-8

;	Looping over the world
;
;tvscl,rebin(Chl_sat,win_xsize/4,win_ysize/4,/sample)
Z_eu=fltarr(win_xsize,win_ysize)
Pb_opt=fltarr(win_xsize,win_ysize)
logic_arr=intarr(win_xsize,win_ysize)

;   Calculate Chl_tot from Satellite Surface Chlorophyll Data
Z_eu=38.0*Float(Chl_sat lt 1.)*(Chl_sat > 0.)^0.425 $
     +40.2*Float(Chl_sat ge 1.)*(Chl_sat > 0.)^0.507			;Chl_tot

;	Calculate Euphotic Depth with Morel's Case I model
Z_eu=(Z_eu > 0.)-float(Z_eu le 0.)
Z_eu=1/Z_eu
Z_eu=(Z_eu > 0.)
;Z_eu=(568.2*(Z_eu le 0.1)*Z_eu^(.746)+ 200.*(Z_eu gt 0.1)*Z_eu^(.293))	;Z_eu
Pb_opt=568.2*(Z_eu^(.746))
Z_eu=200.*(Z_eu^(.293))
Z_eu=Z_eu*float(Z_eu gt 102.)+Pb_opt*float(Z_eu le 102.)

;	Calculate Euphotic Depth with Mike's Equation
;Z_eu=exp(-0.3554*ALOG(Chl_sat > 0.)+3.5805)				;Z_eu
;Z_eu=35.8915*((Chl_sat > 0.)^(-0.3554))				;Z_eu
;;Z_eu=35.8915*exp(-0.3554*ALOG(Chl_sat > 0.))				;Z_eu

;    Calculate the Pb_opt from Satellite Surface Temperature Data
Pb_opt=4.*(Temp_surf gt 28.5)+1.13*(Temp_surf lt -1.)*(Temp_surf gt -10.) ;Pb_opt
logic_arr=(Temp_surf le 28.5)*(Temp_surf ge -1.)
Temp_surf=(0.1*Temp_surf > 1.e-5)
Pb_opt=float(logic_arr eq 0)*Pb_opt+float(logic_arr)* $
	(c7*Temp_surf^7.+c6*Temp_surf^6.+c5*Temp_surf^5. $
	+c4*Temp_surf^4.+c3*Temp_surf^3.+c2*Temp_surf^2. $
	+c1*Temp_surf+c0)					;Pb_opt

;	Calculate the Primary Production
;
for j=0,win_ysize1 do Chl_sat(*,j)=Daylength(j)*Chl_sat(*,j)
Chl_sat = 0.66125*Pb_opt $
;	*Irr_surf(j,Julian_date)/(Irr_surf(j,Julian_date)+4.1) $
	*Irr_surf $
;;	*(0.9888*Irr_surf+0.00349*Irr_surf^2-0.000013*Irr_surf^3) $
	/(Irr_surf+4.1) $
	*(Chl_sat > 0.)*Z_eu + (Chl_sat < 0.)
;tvscl,rebin(Chl_sat,win_xsize/4,win_ysize/4,/sample)

;;print,Temp_file,' Ave_Pb_opt=',Pb_opt_total/Pb_opt_count
Daylength=0
Irr_surf=0

;	To save the new temperature data
;openw,2,Temp_file
;writeu,2,Temp_surf
;close,2
;Temp_surf=0

;	To write the Primary Production data
;		into the output file
;openw,2,PP_file
;writeu,2,Chl_sat
;close,2

return
end