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

case = "b40.1955-2005.2deg.wcm.002"
year = "2005"

fin  = addfile(case+".cice.r."+year+"-01-01-00000.nc","r")

; These fields go straight through
aicen = fin->aicen
vicen = fin->vicen
vsnon = fin->vsnon
Tsfcn = fin->Tsfcn
uvel = fin->uvel
vvel = fin->vvel
scale_factor = fin->scale_factor
coszen = fin->coszen
swvdr = fin->swvdr
swvdf = fin->swvdf
swidr = fin->swidr
swidf = fin->swidf
strocnxT = fin->strocnxT
strocnyT = fin->strocnyT
stressp_1 = fin->stressp_1
stressp_2 = fin->stressp_2
stressp_3 = fin->stressp_3
stressp_4 = fin->stressp_4
stressm_1 = fin->stressm_1
stressm_2 = fin->stressm_2
stressm_3 = fin->stressm_3
stressm_4 = fin->stressm_4
stress12_1 = fin->stress12_1
stress12_2 = fin->stress12_2
stress12_3 = fin->stress12_3
stress12_4 = fin->stress12_4
iceumask = fin->iceumask

; Tracers
aerosnossl1 = fin->aerosnossl1
aerosnoint1 = fin->aerosnoint1
aeroicessl1 = fin->aeroicessl1
aeroiceint1 = fin->aeroiceint1
aerosnossl2 = fin->aerosnossl2
aerosnoint2 = fin->aerosnoint2
aeroicessl2 = fin->aeroicessl2
aeroiceint2 = fin->aeroiceint2
aerosnossl3 = fin->aerosnossl3
aerosnoint3 = fin->aerosnoint3
aeroicessl3 = fin->aeroicessl3
aeroiceint3 = fin->aeroiceint3
iage = fin->iage
FY = fin->FY

; These need to be converted
eicen = fin->eicen
esnon = fin->esnon
apnd = fin->apondn
hpnd = fin->hpondn

; These need to be created

; We don't have frz_onst from LE runs
frz_onset = uvel*0.

ndims1 = dimsizes(aicen)
ndims2 = dimsizes(eicen)
ndims3 = dimsizes(esnon)
print(ndims1)
print(ndims2)
ncat = ndims1(0)
nilyr = ndims2(0)/ncat
nslyr = ndims3(0)/ncat
;nilyr2 = nilyr*2
nj = ndims1(1)
ni = ndims1(2)

qice = new((/nilyr,ncat,nj,ni/),double)
sice = new((/nilyr,ncat,nj,ni/),double)
qsno = new((/nslyr,ncat,nj,ni/),double)

; Convert energy to enthalpy
do n=0,ncat-1
do k=0,nilyr-1
   nt = n*nilyr+k
   tmp = where(vicen(n,:,:).ne.0.,vicen(n,:,:)/int2dble(nilyr),default_fillvalue("double"))
   tmp@_FillValue = default_fillvalue("double")
;  For 8-layer case
;  qice(2*k,n,:,:) = eicen(nt,:,:)*0.5/tmp
;  qice(2*k,n,:,:) = where(ismissing(qice(2*k,n,:,:)),0.,qice(2*k,n,:,:))
;  qice(2*k+1,n,:,:) = eicen(nt,:,:)*0.5/tmp
;  qice(2*k+1,n,:,:) = where(ismissing(qice(2*k+1,n,:,:)),0.,qice(2*k+1,n,:,:))

;  For 4-layer case
   qice(k,n,:,:) = eicen(nt,:,:)/tmp
   qice(k,n,:,:) = where(ismissing(qice(k,n,:,:)),0.,qice(k,n,:,:))
   print(n+1)
   print(k+1)
   tmp1 = ndtooned(qice(k,n,:,:))
   ii = minind(tmp1)
   tmp2 = ndtooned(eicen(nt,:,:))
   tmp3 = ndtooned(vicen(n,:,:))
   print(tmp1(ii))
   print(tmp2(ii))
   print(tmp3(ii))
end do
end do
do n=0,ncat-1
do k=0,nslyr-1
   nt = n*nslyr+k
   tmp = where(vsnon(n,:,:).ne.0.,vsnon(n,:,:)/int2dble(nslyr),default_fillvalue("double"))
   tmp@_FillValue = default_fillvalue("double")
   qsno(k,n,:,:) = esnon(nt,:,:)/tmp
   qsno(k,n,:,:) = where(ismissing(qsno(k,n,:,:)),0.,qsno(k,n,:,:))
end do
end do


; Salinity
saltmax = 3.2d0
nsal = 0.407d0
msal = 0.573d0
pi = atan(1.0d0)*4.0d0

salinz = new((/nilyr/),double)
do k=0,nilyr-1
   zn = (int2dble(k+1)-0.5d0)/int2dble(nilyr)
   salinz(k) = (saltmax/2.d0)*(1.d0-cos(pi*zn^(nsal/(msal+zn))))
   sice(k,:,:,:) = salinz(k)
end do

delete(qice@_FillValue)
delete(sice@_FillValue)
delete(qsno@_FillValue)

file_atts = getvaratts(fin)
print(file_atts)
natts = dimsizes(file_atts)

fout = addfile(case+"_4lyr.cice5.r."+year+"-01-01-00000.nc","c")

setfileoption(fout,"DefineMode",True)

do iatt=0,natts-1
   fout@$file_atts(iatt)$ = fin@$file_atts(iatt)$
end do

dimNames = (/"ncat","nj","ni"/)
dimSizes = (/ncat,nj,ni/)
dimUnlim = (/False,False,False/)
filedimdef(fout,dimNames,dimSizes,dimUnlim)

filevardef(fout,"aicen",typeof(aicen),(/"ncat","nj","ni"/))
filevardef(fout,"vicen",typeof(vicen),(/"ncat","nj","ni"/))
filevardef(fout,"vsnon",typeof(vsnon),(/"ncat","nj","ni"/))
filevardef(fout,"Tsfcn",typeof(Tsfcn),(/"ncat","nj","ni"/))
filevardef(fout,"aerosnossl001",typeof(aerosnossl1),(/"ncat","nj","ni"/))
filevardef(fout,"aerosnoint001",typeof(aerosnoint1),(/"ncat","nj","ni"/))
filevardef(fout,"aeroicessl001",typeof(aeroicessl1),(/"ncat","nj","ni"/))
filevardef(fout,"aeroiceint001",typeof(aeroiceint1),(/"ncat","nj","ni"/))
filevardef(fout,"aerosnossl002",typeof(aerosnossl2),(/"ncat","nj","ni"/))
filevardef(fout,"aerosnoint002",typeof(aerosnoint2),(/"ncat","nj","ni"/))
filevardef(fout,"aeroicessl002",typeof(aeroicessl2),(/"ncat","nj","ni"/))
filevardef(fout,"aeroiceint002",typeof(aeroiceint2),(/"ncat","nj","ni"/))
filevardef(fout,"aerosnossl003",typeof(aerosnossl3),(/"ncat","nj","ni"/))
filevardef(fout,"aerosnoint003",typeof(aerosnoint3),(/"ncat","nj","ni"/))
filevardef(fout,"aeroicessl003",typeof(aeroicessl3),(/"ncat","nj","ni"/))
filevardef(fout,"aeroiceint003",typeof(aeroiceint3),(/"ncat","nj","ni"/))
filevardef(fout,"iage",typeof(iage),(/"ncat","nj","ni"/))
filevardef(fout,"FY",typeof(FY),(/"ncat","nj","ni"/))

filevardef(fout,"uvel",typeof(uvel),(/"nj","ni"/))
filevardef(fout,"vvel",typeof(vvel),(/"nj","ni"/))
filevardef(fout,"scale_factor",typeof(scale_factor),(/"nj","ni"/))
filevardef(fout,"coszen",typeof(coszen),(/"nj","ni"/))
filevardef(fout,"swvdr",typeof(swvdr),(/"nj","ni"/))
filevardef(fout,"swvdf",typeof(swvdf),(/"nj","ni"/))
filevardef(fout,"swidr",typeof(swidr),(/"nj","ni"/))
filevardef(fout,"swidf",typeof(swidf),(/"nj","ni"/))
filevardef(fout,"strocnxT",typeof(strocnxT),(/"nj","ni"/))
filevardef(fout,"strocnyT",typeof(strocnyT),(/"nj","ni"/))
filevardef(fout,"stressp_1",typeof(stressp_1),(/"nj","ni"/))
filevardef(fout,"stressp_2",typeof(stressp_2),(/"nj","ni"/))
filevardef(fout,"stressp_3",typeof(stressp_3),(/"nj","ni"/))
filevardef(fout,"stressp_4",typeof(stressp_4),(/"nj","ni"/))
filevardef(fout,"stressm_1",typeof(stressm_1),(/"nj","ni"/))
filevardef(fout,"stressm_2",typeof(stressm_2),(/"nj","ni"/))
filevardef(fout,"stressm_3",typeof(stressm_3),(/"nj","ni"/))
filevardef(fout,"stressm_4",typeof(stressm_4),(/"nj","ni"/))
filevardef(fout,"stress12_1",typeof(stress12_1),(/"nj","ni"/))
filevardef(fout,"stress12_2",typeof(stress12_2),(/"nj","ni"/))
filevardef(fout,"stress12_3",typeof(stress12_3),(/"nj","ni"/))
filevardef(fout,"stress12_4",typeof(stress12_4),(/"nj","ni"/))
filevardef(fout,"iceumask",typeof(iceumask),(/"nj","ni"/))

filevardef(fout,"qice001",typeof(qice),(/"ncat","nj","ni"/))
filevardef(fout,"qice002",typeof(qice),(/"ncat","nj","ni"/))
filevardef(fout,"qice003",typeof(qice),(/"ncat","nj","ni"/))
filevardef(fout,"qice004",typeof(qice),(/"ncat","nj","ni"/))
;filevardef(fout,"qice005",typeof(qice),(/"ncat","nj","ni"/))
;filevardef(fout,"qice006",typeof(qice),(/"ncat","nj","ni"/))
;filevardef(fout,"qice007",typeof(qice),(/"ncat","nj","ni"/))
;filevardef(fout,"qice008",typeof(qice),(/"ncat","nj","ni"/))
filevardef(fout,"sice001",typeof(sice),(/"ncat","nj","ni"/))
filevardef(fout,"sice002",typeof(sice),(/"ncat","nj","ni"/))
filevardef(fout,"sice003",typeof(sice),(/"ncat","nj","ni"/))
filevardef(fout,"sice004",typeof(sice),(/"ncat","nj","ni"/))
;filevardef(fout,"sice005",typeof(sice),(/"ncat","nj","ni"/))
;filevardef(fout,"sice006",typeof(sice),(/"ncat","nj","ni"/))
;filevardef(fout,"sice007",typeof(sice),(/"ncat","nj","ni"/))
;filevardef(fout,"sice008",typeof(sice),(/"ncat","nj","ni"/))
filevardef(fout,"qsno001",typeof(qsno),(/"ncat","nj","ni"/))
filevardef(fout,"apnd",typeof(apnd),(/"ncat","nj","ni"/))
filevardef(fout,"hpnd",typeof(hpnd),(/"ncat","nj","ni"/))
filevardef(fout,"frz_onset",typeof(frz_onset),(/"nj","ni"/))

setfileoption(fout,"DefineMode",False)

fout->aicen = aicen
fout->vicen = vicen
fout->vsnon = vsnon
fout->Tsfcn = Tsfcn
fout->aerosnossl001 = (/aerosnossl1/)
fout->aerosnoint001 = (/aerosnoint1/)
fout->aeroicessl001 = (/aeroicessl1/)
fout->aeroiceint001 = (/aeroiceint1/)
fout->aerosnossl002 = (/aerosnossl2/)
fout->aerosnoint002 = (/aerosnoint2/)
fout->aeroicessl002 = (/aeroicessl2/)
fout->aeroiceint002 = (/aeroiceint2/)
fout->aerosnossl003 = (/aerosnossl3/)
fout->aerosnoint003 = (/aerosnoint3/)
fout->aeroicessl003 = (/aeroicessl3/)
fout->aeroiceint003 = (/aeroiceint3/)
fout->iage = iage
fout->FY = FY
fout->uvel = uvel
fout->vvel = vvel
fout->scale_factor = scale_factor
fout->coszen = coszen
fout->swvdr = swvdr
fout->swvdf = swvdf
fout->swidr = swidr
fout->swidf = swidf
fout->strocnxT = strocnxT
fout->strocnyT = strocnyT
fout->stressp_1 = stressp_1
fout->stressp_2 = stressp_2
fout->stressp_3 = stressp_3
fout->stressp_4 = stressp_4
fout->stressm_1 = stressm_1
fout->stressm_2 = stressm_2
fout->stressm_3 = stressm_3
fout->stressm_4 = stressm_4
fout->stress12_4 = stress12_4
fout->stress12_1 = stress12_1
fout->stress12_2 = stress12_2
fout->stress12_3 = stress12_3
fout->iceumask = iceumask

fout->frz_onset = (/frz_onset/)
fout->apnd = (/apnd/)
fout->hpnd = (/hpnd/)
fout->qice001 = (/qice(0,:,:,:)/)
fout->qice002 = (/qice(1,:,:,:)/)
fout->qice003 = (/qice(2,:,:,:)/)
fout->qice004 = (/qice(3,:,:,:)/)
;fout->qice005 = (/qice(4,:,:,:)/)
;fout->qice006 = (/qice(5,:,:,:)/)
;fout->qice007 = (/qice(6,:,:,:)/)
;fout->qice008 = (/qice(7,:,:,:)/)
fout->sice001 = (/sice(0,:,:,:)/)
fout->sice002 = (/sice(1,:,:,:)/)
fout->sice003 = (/sice(2,:,:,:)/)
fout->sice004 = (/sice(3,:,:,:)/)
;fout->sice005 = (/sice(4,:,:,:)/)
;fout->sice006 = (/sice(5,:,:,:)/)
;fout->sice007 = (/sice(6,:,:,:)/)
;fout->sice008 = (/sice(7,:,:,:)/)
fout->qsno001 = (/qsno(0,:,:,:)/)

end