;=============================================================== ; ; Annual Primary Production with model_1: ; --- Pb_opt is a function of T ; the result is converted into logarithm scale: ; 0 0.)^0.425 $ +40.2*Float(Chl_sat ge 1.)*(Chl_sat > 0.)^0.507 ;============================================================ ; Calculate Euphotic Depth with Morel's Case I model ;============================================================ print,'Calculating Euphotic depth' Z_eu=(Z_eu > 0.)-float(Z_eu le 0.) Z_eu=1/Z_eu Z_eu=(Z_eu > 0.) 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=35.891tempdata5*((Chl_sat > 0.)^(-0.3554)) ;Z_eu=35.8915*exp(-0.3554*ALOG(Chl_sat > 0.)) ;================================================================ ; Calculate the Pb_opt from Satellite Surface Temperature Data ;================================================================ print,'Calculating Pb_opt' Pb_opt=4.*(Temp_surf gt 28.5)+1.13*(Temp_surf lt -1.)*(Temp_surf gt -10.) 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) ;================================================================ ; Calculate the Primary Production ;================================================================ print,'Calculating Primary Production' PP_data = 0.66125*Pb_opt*Irr_surf/(Irr_surf+4.1)*(Chl_sat > 0.)*Z_eu + (Chl_sat < 0.) PP_data = (Ice_data lt 25)*PP_data ;below 10% of ice coverage Chl_sat = (Ice_data lt 25)*Chl_sat ;below 10% of ice coverage Irr_surf=0 endif ;PP calculate end ;============================================================== ; Main program ;============================================================== ;============================ ; Color table PalFileIn = 'pal0.rgb' colors = bytarr(256,3) openr, 1, PalFileIn readu, 1, colors close,1 colors(250,*) = [200,200,200] tvlct,colors ;=========================== savePPmap = 1 ;flag for saving global PP map savePPres = 1 ;flag for saving global PP results savePPseason = 1 ;flag for saving seasonal PP map savePPmonth = 1 ;flag for saving monthly PP map PPCaclulate = 1 ;flag for calculating PP savemapname = '/usr1/PP/data_files/Annu_Glb_T_Inh_ful_Ice.bin' saveresname = '/usr1/PP/results/annuPP_Glb_T_Inh_ful_Ice.txt' img_name='Primary Production Map' win_xsize=2048 win_ysize=1024 win_xsize1=win_xsize-1 win_ysize1=win_ysize-1 Chl_file=strarr(13) Temp_file=strarr(13) Irr_file=strarr(13) Month_pp_name = strarr(32) Julian_date=intarr(13) Ice_mask=strarr(13) day_of_month=intarr(13) ; January: start of Winter ; Chl_file(1)= '/usr1/chlorophyll/c01_fill.bin' Temp_file(1)='/usr1/sst/ful_Jan.sst' Irr_file(1)= '/usr1/seawifs/ful_8401.irr' Ice_mask(1) = '/usr1/Mask/Ice_Mask_0' Month_pp_name(1) = '/usr1/PP/data_files/monthpp_01.bin' day_of_month(1)=31 Julian_date(1)=1 ; ; February ; Chl_file(2)='/usr1/chlorophyll/c02_fill.bin' Temp_file(2)='/usr1/sst/ful_Feb.sst' Irr_file(2)='/usr1/seawifs/ful_8402.irr' Ice_mask(2) = '/usr1/Mask/Ice_Mask_1' Month_pp_name(2) = '/usr1/PP/data_files/monthpp_02.bin' day_of_month(2)=28 Julian_date(2)=32 ;Julian_date(1)+day_of_month(1) ; ; March ; Chl_file(3)= '/usr1/chlorophyll/c03_fill.bin' Temp_file(3)='/usr1/sst/ful_Mar.sst' Irr_file(3)= '/usr1/seawifs/ful_8403.irr' Month_pp_name(3) = '/usr1/PP/data_files/monthpp_03.bin' Ice_mask(3) = '/usr1/Mask/Ice_Mask_2' day_of_month(3)=31 Julian_date(3)=60 ;Julian_date(2)+day_of_month(2) ; ; April:start of Spring ; Chl_file(4)= '/usr1/chlorophyll/c04_fill.bin' Temp_file(4)='/usr1/sst/ful_Apr.sst' Irr_file(4)= '/usr1/seawifs/ful_8404.irr' Month_pp_name(4) = '/usr1/PP/data_files/monthpp_04.bin' Ice_mask(4) = '/usr1/Mask/Ice_Mask_3' day_of_month(4)=30 Julian_date(4)=91 ;Julian_date(3)+day_of_month(3) ; ; May ; Chl_file(5)= '/usr1/chlorophyll/c05_fill.bin' Temp_file(5)='/usr1/sst/ful_May.sst' Irr_file(5)= '/usr1/seawifs/ful_8405.irr' Month_pp_name(5) = '/usr1/PP/data_files/monthpp_05.bin Ice_mask(5) = '/usr1/Mask/Ice_Mask_4' day_of_month(5)=31 Julian_date(5)=121 ;Julian_date(4)+day_of_month(4) ; ; June ; Chl_file(6)='/usr1/chlorophyll/c06_fill.bin' Temp_file(6)='/usr1/sst/ful_Jun.sst' Irr_file(6)='/usr1/seawifs/ful_8406.irr' Month_pp_name(6) = '/usr1/PP/data_files/monthpp_06.bin Ice_mask(6) = '/usr1/Mask/Ice_Mask_5' day_of_month(6)=30 Julian_date(6)=152 ;Julian_date(5)+day_of_month(5) ; ; July: start of Summer ; Chl_file(7)='/usr1/chlorophyll/c07_fill.bin' Temp_file(7)='/usr1/sst/ful_Jul.sst' Irr_file(7)='/usr1/seawifs/ful_8307.irr' Month_pp_name(7) = '/usr1/PP/data_files/monthpp_07.bin Ice_mask(7) = '/usr1/Mask/Ice_Mask_6' day_of_month(7)=31 Julian_date(7)= 182 ;Julian_date(6)+day_of_month(6) ; ; August ; Chl_file(8)='/usr1/chlorophyll/c08_fill.bin' Temp_file(8)='/usr1/sst/ful_Aug.sst' Irr_file(8)='/usr1/seawifs/ful_8308.irr' Month_pp_name(8) = '/usr1/PP/data_files/monthpp_8.bin Ice_mask(8) = '/usr1/Mask/Ice_Mask_7' day_of_month(8)=31 Julian_date(8)=213 ;Julian_date(7)+day_of_month(7) ; ; September ; Chl_file(9)='/usr1/chlorophyll/c09_fill.bin' Temp_file(9)='/usr1/sst/ful_Sep.sst' Irr_file(9)='/usr1/seawifs/ful_8309.irr' Month_pp_name(9) = '/usr1/PP/data_files/monthpp_09.bin' Ice_mask(9) = '/usr1/Mask/Ice_Mask_8' day_of_month(9)=30 Julian_date(9)=244 ;Julian_date(8)+day_of_month(8) ; ; October: start of Fall ; Chl_file(10)='/usr1/chlorophyll/c10_fill.bin' Temp_file(10)='/usr1/sst/ful_Oct.sst' Irr_file(10)='/usr1/seawifs/ful_8310.irr' Month_pp_name(10) = '/usr1/PP/data_files/monthpp_10.bin' Ice_mask(10) = '/usr1/Mask/Ice_Mask_9' day_of_month(10)=31 Julian_date(10)=274 ;Julian_date(9)+day_of_month(9) ; ; November ; Chl_file(11)='/usr1/chlorophyll/c11_fill.bin' Temp_file(11)='/usr1/sst/ful_Nov.sst' Irr_file(11)='/usr1/seawifs/ful_8311.irr' Month_pp_name(11) = '/usr1/PP/data_files/monthpp_11.bin' Ice_mask(11) = '/usr1/Mask/Ice_Mask_10' day_of_month(11)=30 Julian_date(11)=305 ; ; December ; Chl_file(12)='/usr1/chlorophyll/c12_fill.bin' Temp_file(12)='/usr1/sst/ful_Dec.sst' Irr_file(12)='/usr1/seawifs/ful_8312.irr' Month_pp_name(12) = '/usr1/PP/data_files/monthpp_12.bin' Ice_mask(12) = '/usr1/Mask/Ice_Mask_11' day_of_month(12)=31 Julian_date(12)=335 ;Julian_date(11)+day_of_month(11) ; Calculate the Annual Primary Productivity ; from the average monthly cell_area=19.546*19.546 pi=3.141592654 r_constant=pi/180. ;cell_area=cell_area*(2048./360.)*(1024./180.) ; 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)=4 ; Spring season_start(1)=7 ; Summer season_start(2)=10 ; Fall season_start(3)=1 ; Winter season_name=strarr(4) season_name(0)='/usr1/PP/data_files/Spring.bin' season_name(1)='/usr1/PP/data_files/Summer.bin' season_name(2)='/usr1/PP/data_files/Fall.bin' season_name(3)='/usr1/PP/data_files/Winter.bin' 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 ;================ Get daytime information Daytime=fltarr(win_ysize,366) Daylength=fltarr(win_ysize) Day_file='day_len_1024.dat' print, 'Getting daylength from ',Day_file openr,1,Day_file ;Read DayLength readu,1,Daytime close,1 window, 4, colors =255,xsize=366,ysize=256,xpos=0,ypos=750,title =Day_file tvscl,rotate(rebin(Daytime,win_ysize/4,366,/sample),1) Temp_data=fltarr(win_xsize,win_ysize) Chl_data =fltarr(win_xsize,win_ysize) Month_PP = fltarr(win_xsize, win_ysize) Season_PP = fltarr(win_xsize, win_ysize) Annual_PP = fltarr(win_xsize, win_ysize) Average_chl = fltarr(win_xsize, win_ysize) Annual_PP(*,*) = 0.0 Average_chl(*,*) = 0.0 No_Chl_Data = bytarr(win_xsize, win_ysize) No_Chl_Data(*,*) = 0 print, 'Model calculations for the first month only ' ; for k_season=0, 3 do begin ;four seasons ;for k_season = 0, 0 do begin ;one month Season_PP(*,*) = 0.0 season_end=season_start(k_season)+2 for k=season_start(k_season),season_end do begin ;four seasonss ;for k=season_start(k_season), season_start(k_season) do begin ;one month mmod,Chl_file(k),Temp_file(k),Irr_file(k),Ice_mask(k), Chl_data, Temp_data ; ==== Integrate over a month Month_PP(*,*) = 0.0 for day = Julian_date(k),Julian_date(k)+day_of_month(k)-1 do begin print,day, daytime(0,day),daytime(200,day),daytime(800,day),daytime(1023,day) Daylength(*)=Daytime(*,day) for j=0,win_ysize1 do Month_PP(*,j) = Month_PP(*,j) + Daylength(j) * Temp_data(*,j) endfor ;=============================================================== ; Saving monthly PP data ;=============================================================== if(SavePPMonth eq 1) then begin namePP = Month_pp_name(k) print, 'Name of the month file: ',namePP openw,2,namePP writeu,2,Month_PP close,2 endif Average_chl = Average_chl + Chl_data No_Chl_Data = No_Chl_Data + (Chl_data gt 0) Season_PP = TEMPORARY(Season_PP) + Month_PP endfor ;=============================================================== ; Saving Sasonal PP data ;=============================================================== if(SavePPSeason eq 1) then begin namePP = season_name(k_season) print,'Saving seasonal PP in ',namePP openw,2,namePP writeu,2,Season_PP close,2 endif 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) ;================================================= ; display the Annual_PP ;================================================= window,0,color=256,xsize=win_xsize/4,ysize=win_ysize/4,title='Global.PP' tvscl,rebin(Annual_PP,win_xsize/4,win_ysize/4,/sample) ;================================================= ; Save the Annual Primary Production ;================================================= if (savePPmap eq 1) then begin print, 'Saving global PP map into ',savemapname openw, 2, savemapname writeu,2, Annual_PP close, 2 endif ;====== Calculate the Average Chlorophyll ======== print, 'Calculating average chlorophyll' Average_chl = Average_chl/(No_Chl_Data + 0.000000001) ;================================================= ; Mask ;================================================= print, 'Applying mask for Ocean regions' MASK=bytarr(win_xsize,win_ysize) Mask_name = 'masknew.dat' openr,1,Mask_name readu,1,MASK close,1 MASK=rotate(MASK,7) Temp_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 gt 0.001),1) ;exclude ice cover Area_Ocean=TOTAL(temp_lat*COS_r_lat) ;from Ocean area ;for calculating ;grand Chl average 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 ====== print, 'Calculating Itegrated Global PP' 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_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_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 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 Sum_of_Oceans=PP_Pacific+PP_Atlantic+PP_Indian+PP_Antarctic $ +PP_Arctic+PP_Mediterranean print, '_________________________________________________________________' print, ' GLOBAL PRIMARY PRODUCTION RESULTS' print, '_________________________________________________________________' print, 'Annual Global Production = ',PP_global print, 'Seasonal Production: April-June = ',Season_PP_global(0) print, 'Seasonal Production: July-September = ',Season_PP_global(1) print, 'Seasonal Production: October-December = ',Season_PP_global(2) print, 'Seasonal Production: January-March = ',Season_PP_global(3) print, ' ' print, '_________________________________________________________________' print, ' OCEAN BASIN RESULTS' print, '_________________________________________________________________' print, 'Pacific = ',PP_Pacific,', %=',PP_Pacific/Sum_of_Oceans*100. print, 'Atlantic = ',PP_Atlantic,', %=',PP_Atlantic/Sum_of_Oceans*100. print, 'Indian = ',PP_Indian,', %=',PP_Indian/Sum_of_Oceans*100. print, 'Antarctic = ',PP_Antarctic,', %=',PP_Antarctic/Sum_of_Oceans*100. print, 'Arctic = ',PP_Arctic,', %=',PP_Arctic/Sum_of_Oceans*100. print, 'Mediterranean =',PP_Mediterranean,',%=',PP_Mediterranean/Sum_of_Oceans*100. print, 'Global = ',Sum_of_Oceans print, ' ' print, '________________________________________________________________' print, ' TROPHIC PRODUCTION' print, '________________________________________________________________' ;;print, 'Eq_Annual_PP=',PP_Eq ;;print, 'Trop_Annual_PP=',PP_Trop print, 'Oligotrophic (less than 0.1 mg Chl/m3) (50S-50N) =',PP_Oligotrophic50 print, 'Mesotrophic (between O.1 and 1 mg Chl/m3) (50S-50N)=',PP_Mesotrophic50 print, 'Eutrophic (more than 1mg Chl/m3) (50S-50N) =',PP_Eutrophic50 print, ' ' print, 'Oligotrophic(less than 0.1 mg Chl/m3) (90S-90N) =',PP_Oligotrophic print, 'Mesotrophic (between O.1 and 1 mg Chl/m3) (90S-90N)=',PP_Mesotrophic print, 'Eutrophic (between O.1 and 1 mg Chl/m3) (90S-90N) =',PP_Eutrophic print, ' ' print, '_________________________________________________________________' print, ' ACILILARY DATA' print, '_________________________________________________________________' print, 'Average Chlorophyll = ',Chl_global print, 'Global Ocean Area = ',Area_Ocean print, 'Global Earth Area = ',Area_Earth print, ' ' ;;printf,2, 'Quantum Yield=',PP_global/irr_global ;======================================================= ; Save the data to file 'annuPP_T_cld_inh_ful_Ice.txt' ;======================================================= if (savePPres eq 1) then begin print, 'Saving result data into ',saveresname openw, 2, saveresname printf,2, ' GLOBAL PRIMARY PRODUCTION RESULTS' printf,2, 'Annual Global Production=',PP_global printf,2, 'Seasonal Production: April-June=',Season_PP_global(0) printf,2, 'Seasonal Production: July-September=',Season_PP_global(1) printf,2, 'Seasonal Production: October-December =',Season_PP_global(2) printf,2, 'Seasonal Production: January-March=',Season_PP_global(3) printf,2, ' OCEAN BASIN RESULTS' printf,2, 'Pacific=',PP_Pacific,',%=',PP_Pacific/Sum_of_Oceans*100. printf,2, 'Atlantic',PP_Atlantic,',%=',PP_Atlantic/Sum_of_Oceans*100. printf,2, 'Indian=',PP_Indian,',%=',PP_Indian/Sum_of_Oceans*100. printf,2, 'Antarctic=',PP_Antarctic,',%=',PP_Antarctic/Sum_of_Oceans*100. printf,2, 'Arctic=',PP_Arctic,',%=',PP_Arctic/Sum_of_Oceans*100. printf,2, 'Mediterranean=',PP_Mediterranean,',%=',PP_Mediterranean/Sum_of_Oceans*100. printf,2, 'Global=',Sum_of_Oceans printf,2, ' TROPHIC PRODUCTION' printf,2, 'Eq_Annual_PP=',PP_Eq printf,2, 'Trop_Annual_PP=',PP_Trop printf,2, 'Oligotrophic (less than 0.1 mg Chl/m3 (50S-50N)=',PP_Oligotrophic50 printf,2, 'Mesotrophic (between O.1 and 1 mg Chl/m3 (50S-50N)=',PP_Mesotrophic50 printf,2, 'Eutrophic (more than 1mg Chl/m3) (50S-50N)=',PP_Eutrophic50 printf,2, 'Oligotrophic(less than 0.1 mg Chl/m3 (90S-90N)=',PP_Oligotrophic printf,2, 'Mesotrophic (between O.1 and 1 mg Chl/m3 (90S-90N)=',PP_Mesotrophic printf,2, 'Eutrophic (between O.1 and 1 mg Chl/m3 (90S-90N)=',PP_Eutrophic printf,2, ' ACILILARY DATA' printf,2, 'Average Chlorophyll=',Chl_global printf,2, 'Global Ocean Area=',Area_Ocean printf,2, 'Global Earth Area=',Area_Earth ;;printf,2, 'Quantum Yield=',PP_global/irr_global close,2 endif print, '' print, ' |||||' print, ' /=========\' print, ' / MPLET \' print, ' | CO ED |' print, ' | /\ |' print, ' | ** |' print, ' \ \____/ /' print, ' \=========/' end