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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin

  FIRST     = True
  OVERWRITE = True
  

  ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  ;%%% input files
  ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ; the file that needs to be gridded
    filesIn = "OH_3Dglobal.geos5.72L.2x25.hemco.nc"
    f=addfile(filesIn, "r")
  	fsplit=str_split(filesIn,"/")
    fnameout=fsplit(dimsizes(fsplit)-1)
    fnameout = str_sub_str(fnameout, "geos5.72L.2x25","giss.F40")
    if ( .not.isfilepresent(fnameout) .or. OVERWRITE .eq. "True" ) then

    if ( FIRST ) then

    srcGridName = "GEOS.2x25.nc"
  	dstGridName = "GISS.F40.nc"
  	wgtFileName = "F40.wghts.nc"
    
    ;###############
    ; Source Grid
    ;###############
    Opt                = True
    
  	mf = addfile ( "./pressure.MERRA.nc","r"); this is 72 levels
    pIn                = mf->pmid
    Opt@ForceOverwrite = True
    Opt@PrintTimings   = True
    Opt@Title          = "GEOS-Chem 2x2.5 Grid"
    rectilinear_to_SCRIP(srcGridName,pIn&lat,pIn&lon,Opt)
    delete(Opt)
    ;printVarSummary(pIn)
    ;print(max(pIn(:,0,:,:)))
    ;print(min(pIn(:,0,:,:)))
    ;print(avg(pIn(:,0,:,:)))

    ;###############
    ; Output Grid
    ;###############
    Opt                = True
  	gf = addfile("./pressure.GCAP_2000s.F40.nc","r"); this is 40 levels
    pOut = gf->pmid
    Opt@ForceOverwrite = True
    Opt@PrintTimings   = True
  	Opt@Title          = "GISS ModelE F40 Grid"
    rectilinear_to_SCRIP(dstGridName,pOut&lat,pOut&lon,Opt)
    delete(Opt)
    delete(gf)
    
    ;##################
    ; Generate Weights
    ;##################
    Opt                      = True
    Opt@InterpMethod         = "bilinear" ; default
    Opt@ForceOverwrite       = True
    Opt@PrintTimings         = True
    ESMF_regrid_gen_weights(srcGridName,dstGridName,wgtFileName,Opt)
    delete(Opt)
    Opt                      = True
    Opt@PrintTimings         = True
    pMid = ESMF_regrid_with_weights(pIn,wgtFileName,Opt)    

    FIRST = False

    end if
    ;----------------------------
    ;print out the pressure info
    ;------------------------------
    ;printVarSummary(pOut)
    ;print("pout is")
    ;print(max(pOut(:,0,:,:)))
    ;print(min(pOut(:,0,:,:)))  ; 
    ;print(avg(pOut(:,0,:,:)))
    ;printVarSummary(pMid)
    ;print("pmid is")
    ;print(max(pMid(:,0,:,:)))
    ;print(min(pMid(:,0,:,:)))
    ;print(avg(pMid(:,0,:,:)))
    ;compare pmid and pOut surface
    ;pdiff = pMid(:,0,:,:) - pOut(:,0,:,:)
    ;print("pdiff is")
    ;print(max(pdiff(:,:,:)))
    ;print(min(pdiff(:,:,:)))
    ;print(avg(pdiff(:,:,:)))
    ;######################
    ; Regrid Horizontally
    ;######################
  	OH_in                   = f->OH
  	OH_mid                = ESMF_regrid_with_weights(OH_in,wgtFileName,Opt)
    ;########################
    ; Interpolate Vertically
    ;########################
    ; Linear interpolation, extrapolate near surface
    OH_out               = int2p_n_Wrap(pMid, OH_mid, pOut, -1,1)
    ;########################
    ; Output New File
    ;########################
    system( "/bin/rm -f "+fnameout )
    fout = addfile( fnameout, "c" )
  	fout->OH     =OH_out
    end if
    delete(f) ; Close input file
    system( "/bin/rm -f *.Log" )

  ;print out more information
  print(max(OH_in))
  print(min(OH_in))
  print(avg(OH_in))
  print(max(OH_out))
  print(min(OH_out)); this shows some negative values
  print(avg(OH_out))
  print(max(OH_mid))
  print(min(OH_mid))
  print(avg(OH_mid))
  ;test to see if there is any missing value
  is_msg1 = any(ismissing(OH_in))
  if (is_msg1) then
      print ("'OH_in' contains missing values.")
  else 
      print ("'OH_in' doesn't contain any missing values.")
  end if 

  is_msg2 = any(ismissing(OH_out))
  if (is_msg2) then
      print ("'OH_out' contains missing values.")
  else 
      print ("'OH_out' doesn't contain any missing values.")
  end if 

  is_msg3 = any(ismissing(OH_mid))
  if (is_msg3) then
      print ("'OH_mid' contains missing values.")
  else 
      print ("'OH_mid' doesn't contain any missing values.")
  end if 


  ;----------------------------------------------------------------------
  ;plot zonal mean on the orginal and regridded grid-these can be ignored
  ;----------------------------------------------------------------------
  ;Notes
  ; the zonal mean plot is little bit werid, it seems the zonalmean plot on the original code displays upside Down
  ; I think it maybe because the coordinates leve is in pressure in the original grid while on the regridded grid
  ;the leve is in actual levels
  ; I have used IDL to display the zonalmean plot for both grids and they seem to be fine
  plot = False
if (plot .eq. "True") then
  zave_in   = dim_avg_Wrap(OH_in(5,:,:,:))  
  zave_out   = dim_avg_Wrap(OH_out(5,:,:,:)) 
  ; Notes: the zonalmean plot is kind of weird 
  wks = gsn_open_wks("pdf", "oh_zonalmean_regrid_june")
  gsn_define_colormap(wks, "gui_default")

  res               = True
  
  res@gsnDraw         = False
  res@gsnFrame        = False

  res@cnFillOn        = True
  res@cnLinesOn       = False
  res@cnLineLabelsOn  = False

  ; plot data on original grid
  res@gsnAddCyclic    = False
  dims=tostring(dimsizes(zave_in))
  res@tiMainString = "GEOS5 2X25 OH in june"
  plot_orig = gsn_csm_contour(wks,zave_in(:,:), res)

  ;plot data on regridded data
  res@gsnAddCyclic    = False
  dims=tostring(dimsizes(zave_out))
  res@tiMainString = "ModelE F40 OH in june"
  plot_regrid = gsn_csm_contour(wks,zave_out(:,:), res)

  ; resouces for panelling
  pres        = True
  pres@gsnMaximize    = True
  pres@gsnPanelLabelBar = True

  gsn_panel(wks,(/plot_orig,plot_regrid/),(/2,1/),pres)
  delete(wks)
end if
  
  return
  
end