;	Annual Primary Production (Pb_opt is a function of T)
;		where Pb_opt = 1.54*10.^(0.0275*T - 0.07) 
;	the result is converted into logarithm scale:
;	0<count<246  =  (log10(Prim_Prod)+1.3367)/0.0357
;	for the case  0.005<Prim_Prod<2.83291e+07
;
@model_expT_cld_inh_ful.pro
img_name='Primary Production Map'

win_xsize=2048
win_ysize=1024
win_xsize1=win_xsize-1
win_ysize1=win_ysize-1
;;window,0,color=256,xsize=win_xsize/4,ysize=win_ysize/4,title='Global.PP'
;;restore,'annual.rgb'
;;tvlct,r,g,b

;	Input the Monthly Average Primary Production
;
Chl_file=strarr(14)
Temp_file=strarr(14)
Irr_file=strarr(14)
Julian_date=intarr(14)
day_of_month=intarr(14)
Temperary_data=fltarr(win_xsize,win_ysize)
;;Annual_PP=fltarr(win_xsize,win_ysize)
;;Season_PP=fltarr(win_xsize,win_ysize)
;;Season_PP_lat=fltarr(win_ysize)

;Chl_file(0)='/user5/czcs/full_fil_12.chl'
Chl_file(0)='/user5/newchl/new3_mean12.flt'
Temp_file(0)='/user5/sst/ful_Dec.sst'
Irr_file(0)='/user5/seawifs/ful_8312.irr'
day_of_month(0)=31
Julian_date(0)=349

;Chl_file(1)='/user5/czcs/full_fil_01.chl'
Chl_file(1)='/user5/newchl/new3_mean01.flt'
Temp_file(1)='/user5/sst/ful_Jan.sst'
Irr_file(1)='/user5/seawifs/ful_8401.irr'
day_of_month(1)=31
Julian_date(1)=15

;Chl_file(2)='/user5/czcs/full_fil_02.chl'
Chl_file(2)='/user5/newchl/new3_mean02.flt'
Temp_file(2)='/user5/sst/ful_Feb.sst'
Irr_file(2)='/user5/seawifs/ful_8402.irr'
day_of_month(2)=28
Julian_date(2)=Julian_date(1)+day_of_month(1)

;Chl_file(3)='/user5/czcs/full_fil_03.chl'
Chl_file(3)='/user5/newchl/new3_mean03.flt'
Temp_file(3)='/user5/sst/ful_Mar.sst'
Irr_file(3)='/user5/seawifs/ful_8403.irr'
day_of_month(3)=31
Julian_date(3)=Julian_date(2)+day_of_month(2)

;Chl_file(4)='/user5/czcs/full_fil_04.chl'
Chl_file(4)='/user5/newchl/new3_mean04.flt'
Temp_file(4)='/user5/sst/ful_Apr.sst'
Irr_file(4)='/user5/seawifs/ful_8404.irr'
day_of_month(4)=30
Julian_date(4)=Julian_date(3)+day_of_month(3)

;Chl_file(5)='/user5/czcs/full_fil_05.chl'
Chl_file(5)='/user5/newchl/new3_mean05.flt'
Temp_file(5)='/user5/sst/ful_May.sst'
Irr_file(5)='/user5/seawifs/ful_8405.irr'
day_of_month(5)=31
Julian_date(5)=Julian_date(4)+day_of_month(4)

;Chl_file(6)='/user5/czcs/full_fil_06.chl'
Chl_file(6)='/user5/newchl/new3_mean06.flt'
Temp_file(6)='/user5/sst/ful_Jun.sst'
Irr_file(6)='/user5/seawifs/ful_8406.irr'
day_of_month(6)=30
Julian_date(6)=Julian_date(5)+day_of_month(5)

;Chl_file(7)='/user5/czcs/full_fil_07.chl'
Chl_file(7)='/user5/newchl/new3_mean07.flt'
Temp_file(7)='/user5/sst/ful_Jul.sst'
Irr_file(7)='/user5/seawifs/ful_8307.irr'
day_of_month(7)=31
Julian_date(7)=Julian_date(6)+day_of_month(6)

;Chl_file(8)='/user5/czcs/full_fil_08.chl'
Chl_file(8)='/user5/newchl/new3_mean08.flt'
Temp_file(8)='/user5/sst/ful_Aug.sst'
Irr_file(8)='/user5/seawifs/ful_8308.irr'
day_of_month(8)=31
Julian_date(8)=Julian_date(7)+day_of_month(7)

;Chl_file(9)='/user5/czcs/full_fil_09.chl'
Chl_file(9)='/user5/newchl/new3_mean09.flt'
Temp_file(9)='/user5/sst/ful_Sep.sst'
Irr_file(9)='/user5/seawifs/ful_8309.irr'
day_of_month(9)=30
Julian_date(9)=Julian_date(8)+day_of_month(8)

;Chl_file(10)='/user5/czcs/full_fil_10.chl'
Chl_file(10)='/user5/newchl/new3_mean10.flt'
Temp_file(10)='/user5/sst/ful_Oct.sst'
Irr_file(10)='/user5/seawifs/ful_8310.irr'
day_of_month(10)=31
Julian_date(10)=Julian_date(9)+day_of_month(9)

;Chl_file(11)='/user5/czcs/full_fil_11.chl'
Chl_file(11)='/user5/newchl/new3_mean11.flt'
Temp_file(11)='/user5/sst/ful_Nov.sst'
Irr_file(11)='/user5/seawifs/ful_8311.irr'
day_of_month(11)=30
Julian_date(11)=Julian_date(10)+day_of_month(10)

;Chl_file(12)='/user5/czcs/full_fil_12.chl'
Chl_file(12)='/user5/newchl/new3_mean12.flt'
Temp_file(12)='/user5/sst/ful_Dec.sst'
Irr_file(12)='/user5/seawifs/ful_8312.irr'
day_of_month(12)=31
Julian_date(12)=Julian_date(11)+day_of_month(11)

;Chl_file(13)='/user5/czcs/full_fil_01.chl'
Chl_file(13)='/user5/newchl/new3_mean01.flt'
Temp_file(13)='/user5/sst/ful_Jan.sst'
Irr_file(13)='/user5/seawifs/ful_8401.irr'
day_of_month(13)=31
Julian_date(13)=15

;	Calculate the Annual Primary Productivity
;		from the average monthly

cell_area=19.*19.
pi=3.141592654
r_constant=pi/180.
;cell_area=cell_area*(2048./360.)*(2048./360.)

;;		Calculate Annual_PP
;	Spring:k_season=0, Summer:k_season=1,
;	Fall:k_season=2,   Winter:k_season=3
season_PP_global=fltarr(4)
season_start=intarr(4)
season_start(0)=3	; Spring
season_start(1)=6	; Summer
season_start(2)=9	; Fall
season_start(3)=0	; Winter
season_name=strarr(4)
;season_name(0)='/user5/PP/Spring_ful.dat'
;season_name(1)='/user5/PP/Summer_ful.dat'
;season_name(2)='/user5/PP/Fall_ful.dat'
;season_name(3)='/user5/PP/Winter_ful.dat'
Annual_PP=0
COS_r_lat=fltarr(win_ysize)
temp_lat=fltarr(win_ysize)
for j=0,win_ysize1 do begin
  r_lat=(float(j)/float(win_ysize1)-0.5)*pi
  COS_r_lat(j)=COS(r_lat)
endfor

for k_season=0,3 do begin
  Season_PP=0
  Season_PP_lat=0
  season_end=season_start(k_season)+2
  for k=season_start(k_season),season_end do begin
    model,Chl_file(k),Temp_file(k),Irr_file(k),Julian_date(k),Temperary_data
    Season_PP=TEMPORARY(Season_PP)+(Temperary_data > 0.)* $
	float(day_of_month(k))+(Temperary_data < 0.)
  endfor
  temp_lat=TOTAL((Season_PP > 0.),1)*COS_r_lat
  Season_PP_global(k_season)=TOTAL(temp_lat)*cell_area
  Annual_PP=TEMPORARY(Annual_PP)+Season_PP
;;  openw,2,season_name(k_season)
;;  writeu,2,Season_PP
;;  close,2
endfor
Annual_PP=(Annual_PP > 0.)-float(Annual_PP lt 0.)
PP_global=TOTAL(Season_PP_global)
Season_PP=0
Season_PP_lat=0

;	Display the Annual_PP
;;window,0,color=256,xsize=win_xsize/4,ysize=win_ysize/4,title='Global.PP'
;;restore,'annual.rgb'
;;tvlct,r,g,b
;;tvscl,rebin(Annual_PP,win_xsize/4,win_ysize/4,/sample)

;	Save the Annual Primary Production
;
openw,2,'/user3/PP/Annu_Glb_expT_cld_inh_ful.dat'
writeu,2,Annual_PP
close,2

;	Calculate the Average Chlorophyll
Average_chl=0
for k=1,12 do begin
    openr,1,chl_file(k)
    readu,1,Temperary_data
    close,1
    Average_chl=TEMPORARY(Average_chl)+Temperary_data*day_of_month(k)
endfor
Average_chl=rotate(Average_chl,7)/365.
Temperary_data=0
Chl_global=0.
Area_Ocean=0.
Area_Earth=0.
temp_lat=TOTAL((Average_chl > 0.),1)
Chl_global=TOTAL(temp_lat*COS_r_lat)
temp_lat=TOTAL((Average_chl ge 0.),1)
Area_Ocean=TOTAL(temp_lat*COS_r_lat)
Area_Earth=TOTAL(COS_r_lat)
Chl_global=Chl_global/Area_Ocean
Area_Ocean=cell_area*Area_Ocean
Area_Earth=cell_area*Area_Earth*float(win_xsize)

;	Calculate the Integrated Global Productivity

lat_10S=fix((90.-10.)/180.*float(win_ysize)+0.5)
lat_10N=fix((90.+10.)/180.*float(win_ysize)-0.5)
lat_23S=fix((90.-23.)/180.*float(win_ysize)+0.5)
lat_23N=fix((90.+23.)/180.*float(win_ysize)-0.5)
lat_50S=fix((90.-50.)/180.*float(win_ysize)+0.5)
lat_50N=fix((90.+50.)/180.*float(win_ysize)-0.5)

;;PP_global=0.
PP_Eq=0.
PP_Trop=0.
temp_lat=TOTAL((Annual_PP > 0.),1)*COS_r_lat
PP_Eq=TOTAL(temp_lat(lat_10S:lat_10N))*cell_area
PP_Trop=TOTAL(temp_lat(lat_23S:lat_23N))*cell_area
;;PP_global=TOTAL(temp_lat)*cell_area

PP_Oligotrophic50=0.
PP_Mesotrophic50=0.
PP_Eutrophic50=0.
temp_lat=TOTAL((Average_chl le 0.1)*(Annual_PP > 0.),1)*COS_r_lat
PP_Oligotrophic50=TOTAL(temp_lat(lat_50S:lat_50N))*cell_area
PP_Oligotrophic=TOTAL(temp_lat)*cell_area
temp_lat=TOTAL((Average_chl gt 0.1)*(Average_chl le 1.)*(Annual_PP > 0.),1) $
	*COS_r_lat
PP_Mesotrophic50=TOTAL(temp_lat(lat_50S:lat_50N))*cell_area
PP_Mesotrophic=TOTAL(temp_lat)*cell_area
temp_lat=TOTAL((Average_chl gt 1.)*(Annual_PP > 0.),1)*COS_r_lat
PP_Eutrophic50=TOTAL(temp_lat(lat_50S:lat_50N))*cell_area
PP_Eutrophic=TOTAL(temp_lat)*cell_area

Average_chl=0

MASK=bytarr(win_xsize,win_ysize)
openr,1,'/user5/mask/MASK_2048x1024.dat'
readu,1,MASK
close,1
MASK=rotate(MASK,7)
PP_Pacific=0.
PP_Atlantic=0.
PP_Indian=0.
PP_Antarctic=0.
PP_Arctic=0.
PP_Mediterranean=0.

temp_lat=TOTAL((MASK eq 1B)*(Annual_PP > 0.),1)
PP_Pacific=TOTAL(temp_lat*COS_r_lat)*cell_area
temp_lat=TOTAL((MASK eq 2B)*(Annual_PP > 0.),1)
PP_Atlantic=TOTAL(temp_lat*COS_r_lat)*cell_area
temp_lat=TOTAL((MASK eq 3B)*(Annual_PP > 0.),1)
PP_Indian=TOTAL(temp_lat*COS_r_lat)*cell_area
temp_lat=TOTAL((MASK eq 4B)*(Annual_PP > 0.),1)
PP_Antarctic=TOTAL(temp_lat*COS_r_lat)*cell_area
temp_lat=TOTAL((MASK eq 5B)*(Annual_PP > 0.),1)
PP_Arctic=TOTAL(temp_lat*COS_r_lat)*cell_area
temp_lat=TOTAL((MASK eq 16B)*(Annual_PP > 0.),1)
PP_Mediterranean=TOTAL(temp_lat*COS_r_lat)*cell_area
MASK=0
ANNUAL_PP=0

Sum_of_Oceans=PP_Pacific+PP_Atlantic+PP_Indian+PP_Antarctic $
		+PP_Arctic+PP_Mediterranean

print, 'Global_Annual_PP=',PP_global
print, 'Global_Spring_PP=',Season_PP_global(0)
print, 'Global_Summer_PP=',Season_PP_global(1)
print, 'Global_Fall_PP=',Season_PP_global(2)
print, 'Global_Winter_PP=',Season_PP_global(3)
print, 'Average_chl=',Chl_global
print, 'Global_Ocean_Area=',Area_Ocean
print, 'Global_Earth_Area=',Area_Earth
print, 'Pacific_Annual_PP=',PP_Pacific,',%=',PP_Pacific/Sum_of_Oceans*100.
print, 'Atlantic_Annual_PP',PP_Atlantic,',%=',PP_Atlantic/Sum_of_Oceans*100.
print, 'Indian_Annual_PP=',PP_Indian,',%=',PP_Indian/Sum_of_Oceans*100.
print, 'Antarctic_Annual_PP=',PP_Antarctic,',%=',PP_Antarctic/Sum_of_Oceans*100.
print, 'Arctic_Annual_PP=',PP_Arctic,',%=',PP_Arctic/Sum_of_Oceans*100.
print, 'Mediterranean_Annual_PP=',PP_Mediterranean,',%=',PP_Mediterranean/Sum_of_Oceans*100.
print, 'Sum_of_Oceans=',Sum_of_Oceans
print, 'Eq_Annual_PP=',PP_Eq
print, 'Trop_Annual_PP=',PP_Trop
print, 'Global_Annual_Oligotrophic(50S-50N)=',PP_Oligotrophic50
print, 'Global_Annual_Mesotrophic(50S-50N)=',PP_Mesotrophic50
print, 'Global_Annual_Eutrophic(50S-50N)=',PP_Eutrophic50
print, 'Global_Annual_Oligotrophic(90S-90N)=',PP_Oligotrophic
print, 'Global_Annual_Mesotrophic(90S-90N)=',PP_Mesotrophic
print, 'Global_Annual_Eutrophic(90S-90N)=',PP_Eutrophic
;;print, 'Quantum Yield=',PP_global/irr_global

;	Save the data to file 'annuNEW_T_cld_inh_ful.prn'
openw, 2, '/user3/PP/annuPP_expT_cld_inh_ful.prn'
printf,2, 'Annual_PP with T, with clouds, with inhibition:'
printf,2, 'Global_Annual_PP=',PP_global
printf,2, 'Global_Spring_PP=',Season_PP_global(0)
printf,2, 'Global_Summer_PP=',Season_PP_global(1)
printf,2, 'Global_Fall_PP=',Season_PP_global(2)
printf,2, 'Global_Winter_PP=',Season_PP_global(3)
printf,2, 'Average_chl=',Chl_global
printf,2, 'Global_Ocean_Area=',Area_Ocean
printf,2, 'Global_Earth_Area=',Area_Earth
printf,2, 'Pacific_Annual_PP=',PP_Pacific,',%=',PP_Pacific/Sum_of_Oceans*100.
printf,2, 'Atlantic_Annual_PP',PP_Atlantic,',%=',PP_Atlantic/Sum_of_Oceans*100.
printf,2, 'Indian_Annual_PP=',PP_Indian,',%=',PP_Indian/Sum_of_Oceans*100.
printf,2, 'Antarctic_Annual_PP=',PP_Antarctic,',%=',PP_Antarctic/Sum_of_Oceans*100.
printf,2, 'Arctic_Annual_PP=',PP_Arctic,',%=',PP_Arctic/Sum_of_Oceans*100.
printf,2, 'Mediterranean_Annual_PP=',PP_Mediterranean,',%=',PP_Mediterranean/Sum_of_Oceans*100.
printf,2, 'Sum_of_Oceans=',Sum_of_Oceans
printf,2, 'Eq_Annual_PP=',PP_Eq
printf,2, 'Trop_Annual_PP=',PP_Trop
printf,2, 'Global_Annual_Oligotrophic(50S-50N)=',PP_Oligotrophic50
printf,2, 'Global_Annual_Mesotrophic(50S-50N)=',PP_Mesotrophic50
printf,2, 'Global_Annual_Eutrophic(50S-50N)=',PP_Eutrophic50
printf,2, 'Global_Annual_Oligotrophic(90S-90N)=',PP_Oligotrophic
printf,2, 'Global_Annual_Mesotrophic(90S-90N)=',PP_Mesotrophic
printf,2, 'Global_Annual_Eutrophic(90S-90N)=',PP_Eutrophic
;;printf,2, 'Quantum Yield=',PP_global/irr_global
close,2

end