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

pro model,Chl_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)

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

;	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)
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

;	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.

;	Looping over the world
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(0:win_xsize1,0:win_ysize1)=4.54

;	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.)

Daylength=0
Irr_surf=0

;	To write the Primary Production data
;		into the output file

;openw,2,PP_file
;writeu,2,Chl_sat
;close,2

return
end