load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

; This ncl script reads in aerosol column burdens (loads) from two versions 
; of NorESM/CAM-Oslo and makes global plots of the annually averaged loads,
; including global average as a number in the title line for each figure. 

; Model independent constants
 g=9.80665
 pi=3.1415926
 re=6378.39e3            ; earth radius in m
 coffa=pi*re^2./180.      
 area1=4.*pi*re^2          

; *************************************************************************
; **** To be edited by the user if the ncl script is run interactively ****
;
; Define plot type and plot output format
 if (.not. isvar("plot_type")) then      ; is plot_type on command line?
  plot_type = 6                      ; 0 => SO2 load
                                     ; 1 => SO4 load
                                     ; 2 => BC  load
                                     ; 3 => POM load
                                     ; 4 => SS  load
                                     ; 5 => DU  load
                                     ; 6 => H2O load
                                     ; 7 => Excess look-up table load
 end if 
 if (.not. isvar("format")) then        ; is format on command line?
  format = "ps"
;  format = "eps"
;  format = "png"
;  format = "pdf"
 end if
;
; *************************************************************************
; No changes by the user should be necessary below...
; *************************************************************************
;old  all_files_I  = systemfunc ("ls /media/BackupAK/aerocomA2r128-tester/CTRL2000/aerocomA2r128_2006.cam2.h0.0007-*.nc")
      all_files_I  = systemfunc ("ls " + filepath_I + filenamep_I + "*")
      all_files_II = systemfunc ("ls " + filepath_II + filenamep_II + "*")
      f0_I      = addfile (filepath_I+filename_I, "r")
      f0_II     = addfile (filepath_II+filename_II, "r")
      f1_I      = addfiles (all_files_I, "r")   ; note the "s" of addfile
      f1_II     = addfiles (all_files_II, "r")   ; note the "s" of addfile

; Reading Gaussian weights and other required model variables
 gw0_I=doubletofloat(f0_I->gw)
 gw0_II=doubletofloat(f0_II->gw)

 lon_I=f0_I->lon
 dlon_I=360./dimsizes(lon_I)
 lon_II=f0_II->lon
 dlon_II=360./dimsizes(lon_II)
  
; Initialization (and obtain correct variable dimensions)
 tmp_I=f1_I[:]->PS
 tmp_II=f1_II[:]->PS         
 load_I=tmp_I
 load_II=tmp_II

 if (plot_type.eq.0) then
 var="C_SO2"              ; name of input-variable and plot
 varname="SO2"            ; variable name used in text string:
 if(ModI.eq."CAM4-Oslo") then
   load_I=(/(f1_I[:]->C_SO2)/)*1.e6     ; variable to be plotted from I
 else
   load_I=(/(f1_I[:]->cb_SO2)/)*1.e6/1.998   ; variable to be plotted from I
 end if
 load_II=(/(f1_II[:]->cb_SO2)/)*1.e6/1.998   ; variable to be plotted from II
 else if (plot_type.eq.1) then
 var="C_SO4"              ; name of input-variable and plot
 varname="SO4"            ; variable name used in text string:
 if(ModI.eq."CAM4-Oslo") then
   load_I=(/(f1_I[:]->C_SO4)/)*1.e6     ; variable to be plotted from I
 else
   load_I=(/(f1_I[:]->cb_SO4_A1)/)/3.06 + (/(f1_I[:]->cb_SO4_A2)/)/3.59 + (/(f1_I[:]->cb_SO4_AC)/)/3.06 + (/(f1_I[:]->cb_SO4_NA)/)/3.06 + (/(f1_I[:]->cb_SO4_PR)/)/3.06 + (/(f1_I[:]->cb_SO4_A1_OCW)/)/3.06 + (/(f1_I[:]->cb_SO4_A2_OCW)/)/3.59 + (/(f1_I[:]->cb_SO4_AC_OCW)/)/3.06 + (/(f1_I[:]->cb_SO4_NA_OCW)/)/3.06 + (/(f1_I[:]->cb_SO4_PR_OCW)/)/3.06
   load_I = load_I*1.e6
 end if
 load_II=(/(f1_II[:]->cb_SO4_A1)/)/3.06 + (/(f1_II[:]->cb_SO4_A2)/)/3.59 + (/(f1_II[:]->cb_SO4_AC)/)/3.06 +  (/(f1_II[:]->cb_SO4_NA)/)/3.06 + (/(f1_II[:]->cb_SO4_PR)/)/3.06 + (/(f1_II[:]->cb_SO4_A1_OCW)/)/3.06 + (/(f1_II[:]->cb_SO4_A2_OCW)/)/3.59 + (/(f1_II[:]->cb_SO4_AC_OCW)/)/3.06  + (/(f1_II[:]->cb_SO4_NA_OCW)/)/3.06 + (/(f1_II[:]->cb_SO4_PR_OCW)/)/3.06
 load_II = load_II*1.e6
 else if (plot_type.eq.2) then
 var="C_BC"               ; name of input-variable and plot
 varname="BC"             ; variable name used in text string:
 if(ModI.eq."CAM4-Oslo") then
   load_I=(/(f1_I[:]->C_BC)/)*1.e6     ; variable to be plotted from I
 else
   load_I=(/(f1_I[:]->cb_BC)/) + (/(f1_I[:]->cb_BC_A_OCW)/) + (/(f1_I[:]->cb_BC_AC_OCW)/) + (/(f1_I[:]->cb_BC_AI_OCW)/) + (/(f1_I[:]->cb_BC_NI_OCW)/) + (/(f1_I[:]->cb_BC_N_OCW)/)
   load_I = load_I*1.e6
 end if
 load_II=(/(f1_II[:]->cb_BC)/) + (/(f1_II[:]->cb_BC_A_OCW)/) + (/(f1_II[:]->cb_BC_AC_OCW)/) + (/(f1_II[:]->cb_BC_AI_OCW)/) + (/(f1_II[:]->cb_BC_NI_OCW)/) + (/(f1_II[:]->cb_BC_N_OCW)/)
 load_II = load_II*1.e6
 else if (plot_type.eq.3) then
 var="C_POM"              ; name of input-variable and plot
 varname="POM"            ; variable name used in text string:
 if(ModI.eq."CAM4-Oslo") then
   load_I=(/(f1_I[:]->C_POM)/)*1.e6    ; variable to be plotted from I
 else
   load_I=(/(f1_I[:]->cb_OM)/) + (/(f1_I[:]->cb_OM_AI_OCW)/) + (/(f1_I[:]->cb_OM_AC_OCW)/) + (/(f1_I[:]->cb_OM_NI_OCW)/) \
          + (/f1_I[:]->cb_SOA_A1_OCW/) + (/f1_I[:]->cb_SOA_NA_OCW/)
   load_I = load_I*1.e6
 end if
 load_II=(/(f1_II[:]->cb_OM)/) + (/(f1_II[:]->cb_OM_AI_OCW)/) + (/(f1_II[:]->cb_OM_AC_OCW)/) + (/(f1_II[:]->cb_OM_NI_OCW)/) \
        + (/f1_II[:]->cb_SOA_A1_OCW/) + (/f1_II[:]->cb_SOA_NA_OCW/)
 load_II = load_II*1.e6
 else if (plot_type.eq.4) then
 var="C_SS"               ; name of input-variable and plot
 varname="Sea-salt"       ; variable name used in text string:
 if(ModI.eq."CAM4-Oslo") then
   load_I=(/(f1_I[:]->C_SS)/)*1.e6     ; variable to be plotted from I
 else
   load_I=(/(f1_I[:]->cb_SALT)/) + (/(f1_I[:]->cb_SS_A1_OCW)/) + (/(f1_I[:]->cb_SS_A2_OCW)/) + (/(f1_I[:]->cb_SS_A3_OCW)/) 
   load_I = load_I*1.e6
 end if
 load_II=(/(f1_II[:]->cb_SALT)/) + (/(f1_II[:]->cb_SS_A1_OCW)/) + (/(f1_II[:]->cb_SS_A2_OCW)/) + (/(f1_II[:]->cb_SS_A3_OCW)/) 
 load_II = load_II*1.e6
 else if (plot_type.eq.5) then
 var="C_DUST"             ; name of input-variable and plot
 varname="Dust"           ; variable name used in text string:
 if(ModI.eq."CAM4-Oslo") then
   load_I=(/(f1_I[:]->C_DUST)/)*1.e6   ; variable to be plotted from I
 else
   load_I=(/(f1_I[:]->cb_DUST)/) + (/(f1_I[:]->cb_DST_A2_OCW)/) + (/(f1_I[:]->cb_DST_A3_OCW)/)
   load_I = load_I*1.e6
 end if
 load_II=(/(f1_II[:]->cb_DUST)/) + (/(f1_II[:]->cb_DST_A2_OCW)/) + (/(f1_II[:]->cb_DST_A3_OCW)/)
 load_II = load_II*1.e6
 else if (plot_type.eq.6) then
 var="DAERH2O"             ; name of input-variable and plot
 varname="Condensed H2O"   ; variable name used in text string:
 load_I=(/(f1_I[:]->DAERH2O)/)       ; variable to be plotted from I
 load_II=(/(f1_II[:]->DAERH2O)/)     ; variable to be plotted from II 
 else if (plot_type.eq.7) then
 var="AKCXS"             ; name of input-variable and plot
 varname="AKCXS"   ; variable name used in text string:
 load_I=(/(f1_I[:]->AKCXS)/)       ; variable to be plotted from I
 load_II=(/(f1_II[:]->AKCXS)/)     ; variable to be plotted from II 
 end if
 end if
 end if
 end if
 end if
 end if
 end if
 end if

; Calculating area weighted loads

 load_Ia=load_I              ; initialization of global average variable
 load_IIa=load_II

 xdims_I = dimsizes(gw0_I)
 ;print(xdims_I)
 ydims_I = dimsizes(load_Ia)
 ;print(ydims_I)
 do i=0,dimsizes(gw0_I)-1
 load_Ia(:,i,:)=load_I(:,i,:)*coffa*dlon_I*gw0_I(i)
 end do

 xdims_II = dimsizes(gw0_II)
 ;print(xdims_I)
 ydims_II = dimsizes(load_IIa)
 ;print(ydims_II)
 do i=0,dimsizes(gw0_II)-1
 load_IIa(:,i,:)=load_II(:,i,:)*coffa*dlon_II*gw0_II(i)
 end do

; Defining color scales for each load variable
if (var .eq. "C_SO2" .or.  var .eq. "C_SO4") then
 digg=(/0.1,0.2,0.3,0.5,1,2,3,5,10,20/)                      ; C_SO2 & C_SO4
 else if (var .eq. "C_BC") then
  digg=(/0.03,0.05,0.1,0.2,0.3,0.5,1,2,3,5/)                 ; C_BC
  else if (var .eq. "C_POM" .or. var .eq. "C_SS") then
   digg=(/0.3,0.5,1,2,3,5,10,20,30,50/)                      ; C_POM & C_SS
   else if (var .eq. "C_DUST") then
    digg=(/0.5,1,2.5,5,10,25,50,100,250,500/)                ; C_DUST
    else if (var .eq. "DAERH2O") then
     digg=(/0.5,1,2.5,5,10,25,50,100,250,500/)               ; C_H2O
    else  ; AKCXS
     digg=(/0.01,0.02,0.03,0.05,.1,.2,.3,.5/)                         ; Replace with error message
    end if
   end if
  end if
 end if
end if

;;;;;;;;;;;;;;;;;;;;;;;;;
;
;    Make the Plot
;
;;;;;;;;;;;;;;;;;;;;;;;;;

  wks = gsn_open_wks(format,var)

  gsn_define_colormap(wks,"amwg_blueyellowred")
;  gsn_define_colormap(wks,"BlueDarkRed18")
;  gsn_define_colormap(wks,"precip2_15lev")
;  gsn_define_colormap(wks,"gui_default")
;  gsn_define_colormap(wks,"hotres")
  plot=new(2,graphic)
  
  res                   = True                  ; plot mods desired
  res@gsnSpreadColors   = False                 ; use full colormap
  res@mpFillOn          = False                 
  res@cnFillOn          = True                  ; color fill 
  res@cnLinesOn         = False                 ; no contour lines
  res@cnLineLabelsOn    = False                 
  res@gsnFrame          = False                 ; Do not draw plot
  res@gsnDraw           = False                 ; Do not advance frame
  res@lbLabelBarOn      = False 
  res@tmXBOn =False
  res@tmXTOn =False
  res@tmYLOn =False
  res@tmYROn =False
  res@cnMissingValFillPattern = 0
  res@cnMissingValFillColor = 16
  res@tiMainFontHeightF = 0.03
  res@tiMainFontThicknessF = 2
  res@txFontHeightF       = 0.02
  res@cnFillMode        = "RasterFill"  ; Turn on raster fill
  res@tiMainFont = "helvetica"
  res@tmYRMode             = "Automatic"
  res@cnInfoLabelOn        = False
  res@cnLevelSelectionMode      = "ExplicitLevels"       ; manual levels
;  res@cnFillColors   = (/3,4,5,6,7,8,9,0,10,11,12,13,14,15,16/)   ; gir hvitt midt i ?
;  res@cnFillColors   = (/2,3,4,5,6,7,8,9,10,11,12,13,14,15,16/) 
  res@cnFillColors   = (/3,5,6,8,9,10,11,12,13,14,15,16/) 
;  res@cnLevels       =   sprintf("%4.1f",digg)             ; min level
  res@cnLevels       =   sprintf("%5.2f",digg)             ; min level

;  res@tiMainString   = "CAM4-Oslo"
if (var .eq. "C_SO2" .or.  var .eq. "C_SO4") then
  res@gsnRightString   = "avg = "+sprintf("%4.2f",(sum(dim_avg_n(load_Ia,0))/area1))+" mg S m~S~-2~N~"
else if (var .eq. "AKCXS") then
  res@gsnRightString  = "avg = "+sprintf("%6.4f",(sum(dim_avg_n(load_Ia,0))/area1))+" mg m~S~-2~N~"
else
  res@gsnRightString   = "avg = "+sprintf("%4.2f",(sum(dim_avg_n(load_Ia,0))/area1))+" mg m~S~-2~N~"
end if
end if
  res@gsnLeftString   = varname + " column burden"
  plot(0) = gsn_csm_contour_map_ce(wks,dim_avg_n_Wrap(load_I,0),res)    ; create the plot

;  res@tiMainString   = "CAM5-Oslo"
if (var .eq. "C_SO2") then
  res@gsnRightString  = "avg = "+sprintf("%4.2f",(sum(dim_avg_n(load_IIa,0))/area1))+" mg S m~S~-2~N~"
else if (var .eq. "C_SO4") then
  res@gsnRightString  = "avg = "+sprintf("%4.2f",(sum(dim_avg_n(load_IIa,0))/area1))+" mg S m~S~-2~N~"
else if (var .eq. "AKCXS") then
  res@gsnRightString  = "avg = "+sprintf("%6.4f",(sum(dim_avg_n(load_IIa,0))/area1))+" mg m~S~-2~N~"
else
  res@gsnRightString  = "avg = "+sprintf("%4.2f",(sum(dim_avg_n(load_IIa,0))/area1))+" mg m~S~-2~N~"
end if
end if
end if
  res@gsnLeftString   = varname + " column burden"
  plot(1) = gsn_csm_contour_map_ce(wks,dim_avg_n_Wrap(load_II,0),res)   ; create the plot

  pres                            = True        ; panel plot mods desired
;  pres@gsnMaximize                = True        ; fill the page
;  pres@txString   =   var
  pres@txFontHeightF              =0.015
  pres@txFontThicknessF              =2
  pres@gsnPanelLabelBar = True
  pres@lbLabelFontHeightF = 0.01
;  pres@lbOrientation        ="Vertical"

  gsn_panel(wks,plot,(/1,2/),pres)              ; create panel plot

end