;***************************************;
; Convert ECHAM grib files to netcdf    *;
;***************************************;

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

undef("Pro2D")
function Pro2D(f:file,fout:file,name:string,oldname:string,itime:integer,\
               LAT1x1:numeric,LON1x1:numeric,LAT2x2:numeric,LON2x2:numeric,n:integer)
begin
if(oldname.eq."var179") then
  tmp1 = f->$oldname$(time|itime,lat|:,lon|:)*-1.0
else
  tmp1 = f->$oldname$(time|itime,lat|:,lon|:)
end if
tmp = tmp1(::-1,:) 
tmp@_FillValue = 9.96921e+36
tmp = where(tmp.lt.0.,0.0,tmp)
tmp!0 = "lat"
tmp!1 = "lon"
tmp&lat  = f->lat(::-1)
tmp&lon  = f->lon

x1 = name+"1x1"
x2 = name+"2x2"
;print("Regridding :"+x1+" "+x2)

if(oldname.eq."TSIR" .or. oldname.eq."SI") then
  opt = True
  opt@critpc = 50                ; require only 50% of the values to be precent 
  out1x1 = area_hi2lores_Wrap (tmp&lon,tmp&lat,tmp,True,1,LON1x1,LAT1x1,opt)
  out2x2 = area_hi2lores_Wrap (tmp&lon,tmp&lat,tmp,True,1,LON2x2,LAT2x2,opt)
else
  out1x1 = area_conserve_remap_Wrap (tmp&lon,tmp&lat,tmp,LON1x1,LAT1x1,False)
  out2x2 = area_conserve_remap_Wrap (tmp&lon,tmp&lat,tmp,LON2x2,LAT2x2,False)
end if

if(itime.eq.0) then 
  setfileoption(fout,"DefineMode",True)
  filevardef(fout,name,typeof(tmp),(/"time","lat","lon"/))
  filevarattdef(fout,name,tmp) 
  filevardef(fout,x1,typeof(tmp),(/"time","lat1","lon1"/))
  filevarattdef(fout,x1,tmp) 
  filevardef(fout,x2,typeof(tmp),(/"time","lat2","lon2"/))
  filevarattdef(fout,x2,tmp) 
  setfileoption(fout,"DefineMode",False)
end if
if(itime.eq.0 .and. n.eq.13) then  
  setfileoption(fout,"DefineMode",True)
  filevardef(fout,"TP","float",(/"time","lat","lon"/))
  filevardef(fout,"TP1x1","float",(/"time","lat1","lon1"/))
  filevardef(fout,"TP2x2","float",(/"time","lat2","lon2"/))
  filevardef(fout,"ALBEDO","float",(/"time","lat","lon"/))
  filevardef(fout,"ALBEDO1x1","float",(/"time","lat1","lon1"/))
  filevardef(fout,"ALBEDO2x2","float",(/"time","lat2","lon2"/))
  setfileoption(fout,"DefineMode",False)
end if

;print("Writing 2d data to file...")
fout->$name$(itime,:,:) = (/tmp/)
fout->$x1$(itime,:,:)   = (/out1x1/)
fout->$x2$(itime,:,:)   = (/out2x2/)
      
if(n.eq.13) then
      ;print("Merging the precip variables...")
      cp  = fout->CP(itime,:,:)
      cp@_FillValue = 9.96921e+36
      cp = where(cp.lt.0.0,0.0,cp)
      lsp = fout->LSP(itime,:,:)
      lsp@_FillValue = 9.96921e+36
      lsp = where(lsp.lt.0.0,0.0,lsp)
      tp  = (cp + lsp)*3600.0 
      fout->TP(itime,:,:) = (/tp/)

      cp1  = fout->CP1x1(itime,:,:)
      cp1@_FillValue = 9.96921e+36
      cp1 = where(cp1.lt.0.0,0.0,cp1)
      lsp1 = fout->LSP1x1(itime,:,:)
      lsp1@_FillValue = 9.96921e+36
      lsp1 = where(lsp1.lt.0.0,0.0,lsp1)
      tp1x1  = (cp1 + lsp1)*3600.0 
      fout->TP1x1(itime,:,:) = (/tp1x1/)

      cp2  = fout->CP2x2(itime,:,:)
      cp2@_FillValue = 9.96921e+36
      cp2 = where(cp2.lt.0.0,0.0,cp2)
      lsp2 = fout->LSP2x2(itime,:,:)
      lsp2@_FillValue = 9.96921e+36
      lsp2 = where(lsp2.lt.0.0,0.0,lsp2)
      tp2x2  = (cp2 + lsp2)*3600.0 
      fout->TP2x2(itime,:,:) = (/tp2x2/)

      print("Calculate the ALBEDO...")
      tsir  = fout->TSIR(itime,:,:)
      tsir@_FillValue = 9.96921e+36
      si = fout->SI(itime,:,:)
      si@_FillValue = 9.96921e+36
      tmp3 = mask(si,si.le.0.0,False)
      albedo = where(.not.ismissing(tmp3),(1 - (tsir/tmp3)),si@_FillValue)  
      fout->ALBEDO(itime,:,:) = (/albedo/)

      tsir1 = fout->TSIR1x1(itime,:,:)
      tsir1@_FillValue = 9.96921e+36
      si1 = fout->SI1x1(itime,:,:)
      si1@_FillValue = 9.96921e+36
      tmp4 = mask(si1,si1.le.0.0,False)
      albedo1x1 = where(.not.ismissing(tmp4),(1 - (tsir1/tmp4)),si1@_FillValue)  
      fout->ALBEDO1x1(itime,:,:) = (/albedo1x1/)

      tsir2  = fout->TSIR2x2(itime,:,:)
      tsir2@_FillValue = 9.96921e+36
      si2 = fout->SI2x2(itime,:,:)
      si2@_FillValue = 9.96921e+36
      tmp5 = mask(si2,si2.le.0.0,False)
      albedo2x2 = where(.not.ismissing(tmp5),(1 - (tsir2/tmp5)),si2@_FillValue)  
      fout->ALBEDO2x2(itime,:,:) = (/albedo2x2/)
end if
return(fout)
end

undef("Pro3D")
function Pro3D(f:file,fout:file,name:string,oldname:string,itime:integer,level:integer,plev:double,\
               LAT1x1:numeric,LON1x1:numeric,LAT2x2:numeric,LON2x2:numeric,n:integer)
begin
if(oldname.eq."var21") then
  tmp1 = f->$oldname$(time|itime,lev_2|level,lat|:,lon|:)*100.
else
  tmp1 = f->$oldname$(time|itime,lev_2|level,lat|:,lon|:)
end if

tmp = tmp1(::-1,:) 
tmp@_FillValue = 9.96921e+36
tmp = where(tmp.lt.0.,0.0,tmp)
tmp!0 = "lat"
tmp!1 = "lon"
tmp&lat  = f->lat(::-1)
tmp&lon  = f->lon

x1 = name+"1x1"
x2 = name+"2x2"
;print("Regridding :"+x1+" "+x2)

if(oldname.eq."TSIR" .or. oldname.eq."SI") then
  opt = True
  opt@critpc = 50                ; require only 50% of the values to be precent 
  out1x1 = area_hi2lores_Wrap (tmp&lon,tmp&lat,tmp,True,1,LON1x1,LAT1x1,opt)
  out2x2 = area_hi2lores_Wrap (tmp&lon,tmp&lat,tmp,True,1,LON2x2,LAT2x2,opt)
else
  out1x1 = area_conserve_remap_Wrap (tmp&lon,tmp&lat,tmp,LON1x1,LAT1x1,False)
  out2x2 = area_conserve_remap_Wrap (tmp&lon,tmp&lat,tmp,LON2x2,LAT2x2,False)
end if

if(itime.eq.0 .and. level.eq.0) then 
  setfileoption(fout,"DefineMode",True)
  filevardef(fout,name,typeof(tmp),(/"time","lev","lat","lon"/))
  filevarattdef(fout,name,tmp) 
  filevardef(fout,x1,typeof(tmp),(/"time","lev","lat1","lon1"/))
  filevarattdef(fout,x1,tmp) 
  filevardef(fout,x2,typeof(tmp),(/"time","lev","lat2","lon2"/))
  filevarattdef(fout,x2,tmp) 
  setfileoption(fout,"DefineMode",False)
end if
if(itime.eq.0 .and. level.eq.0 .and. n.eq.0) then
  setfileoption(fout,"DefineMode",True)
  filevardef(fout,"IWC1x1","float",(/"time","lev","lat1","lon1"/))
  filevardef(fout,"IWC2x2","float",(/"time","lev","lat2","lon2"/))
  filevardef(fout,"IWC","float",(/"time","lev","lat","lon"/))
  setfileoption(fout,"DefineMode",False)
end if

;print("Writing 3d data to file...")
fout->$name$(itime,level,:,:) = (/tmp/)
fout->$x1$(itime,level,:,:)   = (/out1x1/)
fout->$x2$(itime,level,:,:)   = (/out2x2/)

;print("Writing 3d iwc data to file...")
if(n.eq.21) then 
  q = fout->Q(itime,level,:,:)
  q@_FillValue = 9.96921e+36
  t = fout->T(itime,level,:,:)
  q = where(q.le.0.0,q@_FillValue,q) 
  ciwc = fout->CIWC(itime,level,:,:)
  rhom = ((lev(level)/(8.314*t)) * (1.0/((1-q)/(q*0.02896)+(1.0/0.01802))))
  rhod = where(q.gt.0.0,rhom * ((1-q)/q),q@_FillValue)
  rho = rhod + rhom
  rho = where(rho.le.0.0,q@_FillValue,rho)
  iwc = where(rho.gt.0.0,(ciwc/rho)*1.0e6,q@_FillValue); iwc in mg/m3
  iwc@_FillValue = 9.96921e+36
  iwc = where(iwc.lt.0.,iwc@_FillValue,iwc)
  fout->IWC(itime,level,:,:) = (/iwc/)

  q1 = fout->Q1x1(itime,level,:,:)
  q1@_FillValue = 9.96921e+36
  t1 = fout->T1x1(ts,level,:,:)
  ciwc1 = fout->CIWC1x1(itime,level,:,:)
  q1 = where(q1.le.0.0,q1@_FillValue,q1)
  rhom1 = ((lev(level)/(8.314*t1)) * (1.0/((1-q1)/(q1*0.02896)+(1.0/0.01802))))
  rhod1 = where(q1.gt.0.0,rhom1 * ((1-q1)/q1),q1@_FillValue)
  rho1 = rhod1 + rhom1
  rho1 = where(rho1.le.0.0,q1@_FillValue,rho1)
  iwc1x1 = where(rho1.gt.0.0,(ciwc1/rho1)*1.0e6,q1@_FillValue); iwc in mg/m3
  iwc1x1@_FillValue = 9.96921e+36
  iwc1x1 = where(iwc1x1.lt.0.,iwc1x1@_FillValue,iwc1x1)
  fout->IWC1x1(itime,level,:,:) = (/iwc1x1/)

  q2 = fout->Q2x2(itime,level,:,:)
  q2@_FillValue = 9.96921e+36
  t2 = fout->T2x2(itime,level,:,:)
  ciwc2 = fout->CIWC2x2(itime,level,:,:)
  q2 = where(q2.le.0.0,q2@_FillValue,q2)
  rhom2 = ((lev(level)/(8.314*t2)) * (1.0/((1-q2)/(q2*0.02896)+(1.0/0.01802))))
  rhod2 = where(q2.gt.0.0,rhom2 * ((1-q2)/q2),q2@_FillValue)
  rho2 = rhod2 + rhom2
  rho2 = where(rho2.le.0.0,q2@_FillValue,rho2)
  iwc2x2 = where(rho.gt.0.0,(ciwc/rho)*1.0e6,q@_FillValue); iwc in mg/m3
  iwc2x2@_FillValue = 9.96921e+36
  iwc2x2 = where(iwc2x2.lt.0.,iwc2x2@_FillValue,iwc2x2)
  fout->IWC2x2(itime,level,:,:) = (/iwc2x2/)
end if
  return(fout)
end

;*************
;Main program*
;*************
begin

name =    (/"var26",\
            "var142",\
            "var143",\
            "var178",\
            "var179",\
	    "var184",\
            "var245",\
            "var246",\
            "var247",\
            "var248",\
	    "var249",\
            "var251",\
            "var252",\
            "var253",\
            "var21",\
            "var22",\
            "var23",\
            "var24",\
            "var25",\
            "var223",\
            "var243",\
            "var244"/)

Lname =(/"surface pressure (in Pa)",\
     	 "large-scale precipitation rate (in kg m-2 s-1)",\
 	 "convective precipitation rate (in kg m-2 s-1)",\
   	 "net solar radiation at the TOA (W m-2)",\
	 "outgoing LW radiation, multiplied by -1 (Wm-2)",\
	 "incoming solar radiation at TOA (Wm-2)",\
	 "sunlit fraction for ISCCP simulator",\
	 "sunlit fraction x total cloud fraction from ISCCP simulator",\ 
	 "sunlit fraction x total cloud fraction x cloud top pressure from ISCCP simulator",\ 
	 "sunlit fraction x total cloud fraction x cloud albedo from ISCCP simulator",\
	 "sunlit fraction x total cloud fraction x cloud optical depth from ISCCP simulator",\
	 "frequency of deep convection",\
	 "frequency of shallow convection",\
	 "frequency of midlevel convection",\
         "relative humidity (0...1)",\
	 "specific humidity (kg/kg)",\
	 "cloud liquid water in (kg/kg) (grid-mean values, not in-cloud)",\
	 "cloud ice in (kg/kg) (grid-mean values, not in-cloud)",\
	 "temperature (K)",\
	 "cloud fraction (ECHAM6s ordinary cloud fraction output)",\
         "cloud fraction input to ISCCP simulator",\
         "cloud fraction x in-cloud optical depth input to ISCCP simulator"/)
Sname  = (/"SP",\
           "LSP",\
           "CP",\
           "TSIR",\
           "TTR",\
           "SI",\
           "sunlit_ISCCP","TCC_ISCCP",\ 
	   "CTP_ISCCP","CALB_ISCCP","TAU_ISCCP",\
           "DCFLAG","SCFLAG","MCFLAG",\
           "R","Q",\
           "CLWC",\
           "CIWC","T","CC",\
           "CC_ISCCP_3D","TAU_ISCCP_3D"/)

; Constants
NLAT  = 75	
NLON  = 360  
NLAT2 = 38
NLON2 = 180	
nlat  = 86
nlon  = 384 

LAT1x1 = fspan(-37.0,37.0,NLAT)
LAT1x1@units = "Degrees north"   
LAT1x1@axis = "Y"   
NLAT         = dimsizes(LAT1x1)
LON1x1       = fspan(0.0,359.0,NLON)
LON1x1!0     = "lon1"
LON1x1@units = "Degrees east"   
LON1x1@axis = "X"   

LAT2x2       = fspan(-37.0,37.0,NLAT2)
NLAT2        = dimsizes(LAT2x2)
LON2x2       = fspan(0.0,359.0,NLON2)
LON2x2!0     = "lon2"
LAT2x2@axis = "Y"
LON2x2@units = "Degrees east"
LON2x2@axis = "X"      

mlev = 36
nhym = 95
nhyi = 96
lev = (/10000.,15000.,20000.,25000.,30000.,35000.,40000.,45000.,50000./)
m_lev = (/60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,\ 
          77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95/)
;***********************************************
; get variable names from grib file
;***********************************************
diri = "/nobackup/rossby16/sm_marst/zproject/echam/"
diro = "/nobackup/rossby15/sm_marst/zproject/echam/"
outfile = diro+"ECHAM6_regridded_conserv_2007-2008.nc"
system ("/bin/rm -f " + outfile)    ; remove any pre-exist file    

f = addfile(diri+"echam_indata_2007-2008.nc","r")
   
names = getfilevarnames(f); extract all variable names

plev = dimsizes(f->lev_2)
nlev = dimsizes(f->lev)
msize = dimsizes(f->hyam)
isize = dimsizes(f->hyai)

print("Names in nc file: "+name) 

print("Create a large file to hold the data")
setfileoption("nc","Format","LargeFile")
fout = addfile(outfile,"c")
  
setfileoption(fout,"DefineMode",True)
dimNames = (/"time","lat","lon","lat1","lon1","lat2","lon2","lev","mlev","nhym","nhyi"/)
dimSizes = (/-1    , nlat, nlon, NLAT , NLON , NLAT2, NLON2, plev, nlev,  msize, isize/) 
dimUnlim = (/True  , False,False,False, False, False, False, False,False, False, False/) 
filedimdef(fout,dimNames,dimSizes,dimUnlim)
   
filevardef(fout, "time","double","time") 
filevardef(fout, "lev", "double","lev")                          
filevardef(fout, "mlev","double","mlev")                          
filevardef(fout, "lat", "double","lat")                          
filevardef(fout, "lon", "double","lon")                          
filevardef(fout, "lat1","double","lat1")                          
filevardef(fout, "lon1","double","lon1")                          
filevardef(fout, "lat2","double","lat2")                          
filevardef(fout, "lon2","double","lon2")                         
filevardef(fout, "hyam","double","nhym")    
filevardef(fout, "hybm","double","nhym")    
filevardef(fout, "hyai","double","nhyi")          
filevardef(fout, "hybi","double","nhyi")          
   
filevarattdef(fout,"time",f->time)             
fout->time@units = "hours since 2007-01-01 01:32:00" ; change the units to hours
filevarattdef(fout,"lat",f->lat)               
filevarattdef(fout,"lon",f->lon)               
filevarattdef(fout,"lat1",f->lat)              
filevarattdef(fout,"lon1",f->lon)              
filevarattdef(fout,"lat2",f->lat)              
filevarattdef(fout,"lon2",f->lon)              
filevarattdef(fout,"lev",f->lev_2)
filevarattdef(fout,"mlev",f->lev)	
filevarattdef(fout,"hyai",f->hyai)		
filevarattdef(fout,"hybi",f->hybi)		
filevarattdef(fout,"hyam",f->hyam)		
filevarattdef(fout,"hybm",f->hybm)		
setfileoption(fout,"DefineMode",False)

fout->lev  = (/f->lev_2/)
fout->mlev = (/f->lev/)
fout->lat  = (/f->lat(::-1)/)
fout->lon  = (/f->lon/) 
fout->lat1 = (/LAT1x1/)
fout->lon1 = (/LON1x1/) 
fout->lat2 = (/LAT2x2/)
fout->lon2 = (/LON2x2/) 
  
fout->hyai = (/f->hyai/)
fout->hybi = (/f->hybi/)
fout->hyam = (/f->hyam/)
fout->hybm = (/f->hybm/)
fout->time = (/f->time/)/60.0 ; change the time to hours

endtime = dimsizes(f->time) - 1  
nsize = dimsizes(name) - 1
nlevs = dimsizes(f->lev_2) - 1  

print("Regridding the variables...")       
do itime = 0, 3;endtime  
  print("Processing time step: "+itime)
  do n = 0, nsize ; loop over the different variables
   print("Processing variable:-> "+name(n)+" mapping to:->"+Sname(n))
   if(name(n).eq."var26" .or.\
      name(n).eq."var142" .or.\
      name(n).eq."var143" .or.\
      name(n).eq."var178" .or.\
      name(n).eq."var179" .or.\
      name(n).eq."var184" .or.\
      name(n).eq."var245" .or.\
      name(n).eq."var246" .or.\
      name(n).eq."var247" .or.\
      name(n).eq."var248" .or.\
      name(n).eq."var249" .or.\
      name(n).eq."var251" .or.\
      name(n).eq."var252" .or.\
      name(n).eq."var253") then
      fout = Pro2D(f,fout,Sname(n),name(n),itime,LAT1x1,LON1x1,LAT2x2,LON2x2,n)
    end if
  
   if(name(n).eq."var21" .or.\
      name(n).eq."var22" .or.\
      name(n).eq."var23" .or.\
      name(n).eq."var24" .or.\
      name(n).eq."var25" .or.\
      name(n).eq."var223") then
      do level = 0, nlevs
        fout = Pro3D(f,fout,Sname(n),name(n),itime,level,f->lev_2(level),LAT1x1,LON1x1,LAT2x2,LON2x2,n)
      end do
    end if
 end do

 ;print("Time in file: "+fout->time(itime))
end do

end