pro mkdiffs
;
; Read 2 history files, calculate differences of selected 4d vars, 
; and write new output file containing the diffs.
;
 file1 = 'timegcm_edyn.nc'            ; edynamo
 file0 = 'timegcm_pdyn.s_apex.nc'     ; pdynamo
 filediffs = 'edyn-pdyn_diffs.nc'     ; raw diffs pdyn apex wts - pdyn esmf wts
 percent_diffs = 0
;
;file1 = 'edynamo_apex_time.nc'       ; edynamo with apex_subs from timegcm
;file0 = 'edynamo_apex_tie.nc'        ; edynamo with apex_subs from tiegcm
;filediffs = 'apex_time-tie_diffs.nc' ; raw diffs 
;percent_diffs = 0
;
; Fields to read and calculate file1-file0 differences 
; (vars must be on both file0 and file1):
;
;fields = ['PHIM3D','ED13D','ED23D','EPHI3D','ELAM3D','EMZ3D']

;fields = ['PHIM2D','ED1_a','ED2_a','EPHI_a','ELAM_a','ED1','ED2','EPHI','ELAM',$
;          'PHIM3D','ED13D','ED23D','EPHI3D','ELAM3D','EMZ3D']

;fields = ['PHIM2D','ED1_a','ED2_a','EPHI_a','ELAM_a','ED1','ED2','EPHI','ELAM',$
;          'THETAM','ZPOTENM3D','PHIM3D','ED13D','ED23D','EPHI3D','EMZ3D']

;fields = ['PHIHM','PFRAC','PHIM2D','ZIGM11_GLB','ZIGM22_GLB','ZIGMC_GLB','ZIGM2_GLB','RHS_GLB',$
;          'C0_01','C0_10','COFUM_01','COFUM_09','RIM1_PRESOLV','RIM2_PRESOLV',$
;          'RIM1_SOLV','RIM2_SOLV']

;fields = ['C0_01','C0_02','C0_03','C0_04','C0_05','C0_06','C0_07','C0_08','C0_09','C0_10',$
;          'COFUM_01','COFUM_02','COFUM_03','COFUM_04','COFUM_05','COFUM_06','COFUM_07',$
;          'COFUM_08','COFUM_09']

;fields = ['ZIGM11_GLB','ZIGM22_GLB','ZIGMC_GLB','ZIGM2_GLB','RHS_GLB']

;fields = ['RIM1_rhs','RIM2_rhs','RHS_N','RHS_glb']

;fields = ['ZIGM11_b','ZIGM22_b','ZIGMC_b','ZIGM2_b','RIM1_b','RIM2_b',$
;          'ZIGM11_c','ZIGM22_c','ZIGMC_c','ZIGM2_c','RIM1_c','RIM2_c',$
;          'ZIGM11_d','ZIGM22_d','ZIGMC_d','ZIGM2_d','RIM1_d','RIM2_d',$
;          'ZIGM11_e','ZIGM22_e','ZIGMC_e','ZIGM2_e','RIM1_e','RIM2_e']

;fields = ['PED_MAG','HAL_MAG','ZPOT_MAG','ADOTV1_MAG','ADOTV2_MAG',$
;          'SINI_MAG','ADOTA1_MAG','ADOTA2_MAG','BE3_MAG',$
;          'ZIGM11_b','ZIGM22_b','ZIGMC_b','ZIGM2_b','RIM1_b','RIM2_b']

;fields = ['ADOTV1','ADOTV2','ADOTA1','ADOTA2','A1DTA2','BE3','SINI']
;          'ZB_DIAG','BMOD_DIAG'] 

;fields = ['TN','UN','VN','OMEGA','Z','BARMASS','PED','HALL',$
;          'PED_MAG','HAL_MAG','ZPOT_MAG','ADOTV1_MAG','ADOTV2_MAG',$
;          'ZIGM11_b','ZIGM22_b','ZIGMC_b','ZIGM2_b','RIM1_b','RIM2_b']

 fields =  ['PED_MAG','HAL_MAG','ZPOT_MAG','ADOTV1_MAG','ADOTV2_MAG',$
            'ZIGM11_b','ZIGM22_b','ZIGMC_b','ZIGM2_b','RIM1_b','RIM2_b']

;fields = [ 'PED_GEO','HAL_GEO','ZPOT_GEO','SCHT_GEO','ADOTV1_GEO','ADOTV2_GEO',$
;           'SINI_GEO','ADOTA1_GEO','ADOTA2_GEO','A1DTA2_GEO','BE3_GEO',$
;           'PED_MAG','HAL_MAG','ZPOT_MAG','ADOTV1_MAG','ADOTV2_MAG',$
;           'SINI_MAG','ADOTA1_MAG','ADOTA2_MAG','A1DTA2_MAG','BE3_MAG']

;fields = ['TN','OMEGA','Z','BARMASS','SCHEIGHT','WN',$
;          'PED_GEO','HAL_GEO','ZPOT_GEO','SCHT_GEO','ADOTV1_GEO','ADOTV2_GEO',$
;          'SINI_GEO','ADOTA1_GEO','ADOTA2_GEO','A1DTA2_GEO','BE3_GEO']

;fields = ['TN_DYN','W_DYN','Z_DYN','BARM_DYN','PED_DYN','HAL_DYN','SCHEIGHT','WN',$
;          'ADOTV1','ADOTV2','SINI','ADOTA1','ADOTA2','A1DTA2',$
;          'PED_GEO','HAL_GEO','ZPOT_GEO','SCHT_GEO','ADOTV1_GEO','ADOTV2_GEO',$
;          'SINI_GEO','ADOTA1_GEO','ADOTA2_GEO','A1DTA2_GEO','BE3_GEO']
;
nflds = n_elements(fields)
;
; Open perturbed and control files read-only:
;
ncid1 = ncdf_open(file1) ; default is nowrite
print,'Opened perturbed file (file1) ',file1
ncid0 = ncdf_open(file0) ; default is nowrite
print,'Opened control file (file0) ',file0
file1struct = ncdf_inquire(ncid1) ; get info for perturbed file
;
; Create new diffs file and define dimensions and coord variables from
; perturbed file:
;
nciddiffs = ncdf_create(filediffs,/clobber) ; overwrite if pre-existing
for i=0,file1struct.ndims-1 do begin
  ncdf_diminq,ncid1,i,name,size
  id = ncdf_dimdef(nciddiffs,name,size)
;
; Define coordinate variable if it exists (assumed to be double):
;
  id1 = ncdf_varid(ncid1,name)
  if (id1 ne -1) then begin
    vinfo = ncdf_varinq(ncid1,name)
    id = ncdf_vardef(nciddiffs,vinfo.name,vinfo.dim,/double)
; Copy coord var attributes:
    for ii=0,vinfo.natts-1 do begin
      attname = ncdf_attname(ncid1,id1,ii)
      iatt = ncdf_attcopy(ncid1,id1,attname,nciddiffs,id)
    endfor
  endif
endfor
;
; Define vars mtime, time, and some global attributes from
; perturbed file1:
;
ncdf_varget,ncid1,"mtime",mtime
ncdf_varget,ncid1,"year",year
ncdf_varget,ncid1,"day",day
ncdf_varget,ncid1,"time",time
ncdf_varget,ncid1,"ut",ut
ntime = n_elements(time)

info = ncdf_varinq(ncid1,'mtime')
id   = ncdf_vardef(nciddiffs,'mtime',info.dim,/long)
info = ncdf_varinq(ncid1,'year')
id   = ncdf_vardef(nciddiffs,'year',info.dim,/long)
info = ncdf_varinq(ncid1,'day')
id   = ncdf_vardef(nciddiffs,'day',info.dim,/long)
info = ncdf_varinq(ncid1,'ut')
id   = ncdf_vardef(nciddiffs,'ut',info.dim,/long)
id   = ncdf_vardef(nciddiffs,'time',/double)
print,'ntime=',ntime,' mtime=' & print,mtime
;
; Global attribute model_name:
ncdf_attget,ncid1,'model_name',model_name,/global
ncdf_attput,nciddiffs,'model_name',model_name,/global,/char
;
; Global attribute model_version:
ncdf_attget,ncid1,'model_version',model_version,/global
ncdf_attput,nciddiffs,'model_version',model_version,/global,/char
;
; Global attribute missing_value:
ncdf_attget,ncid1,'missing_value',missing_value,/global
ncdf_attput,nciddiffs,'missing_value',missing_value,/global,/double
;
; Add global attribute with names of file0,1:
;
attribute = "Difference fields of perturbed file "+file1+ $
  " minus control file "+file0
ncdf_attput,nciddiffs,'description',attribute,/global,/char
;
; Define requested fields on diff file, with attributes:
; (assumed to be doubles)
;
for i=0,nflds-1 do begin
  id1 = ncdf_varid(ncid1,fields[i])
  if (id1 eq -1) then begin
    print,'>>> WARNING: Could not find var ',fields[i],$
    ' on file1 ',file1
  endif else begin
    vinfo = ncdf_varinq(ncid1,fields[i])
;
; If diffs of Z, name var on diffs file ZDIFFS, to avoid naming
; conflict w/ perturbed Z: 
    if (vinfo.name eq 'Z') then begin
      iddiff = ncdf_vardef(nciddiffs,"ZDIFFS",vinfo.dim,/double)
    endif else begin
      iddiff = ncdf_vardef(nciddiffs,vinfo.name,vinfo.dim,/double)
    endelse

    for ii=0,vinfo.natts-1 do begin
      attname = ncdf_attname(ncid1,id1,ii)
      ncdf_attget,ncid1,id1,attname,attval
;
; Add diffs to long_name attribute:
      if (attname eq 'long_name') then begin
        if (strlen(attval) eq 0) then attval = vinfo.name

        if (percent_diffs le 0) then longname = 'Diffs: '+string(attval) $
        else longname = '%Diffs: '+string(attval)

        longname = byte(longname)
        ncdf_attput,nciddiffs,iddiff,attname,longname,/char
      endif else begin
        iatt = ncdf_attcopy(ncid1,id1,attname,nciddiffs,iddiff)
      endelse
    endfor

  endelse
endfor
;
; Define geopotential Z from file1 (copy attributes):
;
idz1 = ncdf_varid(ncid1,'Z')
if (idz1 eq -1) then begin
  print,'>>> WARNING: Could not find var Z on file1 ',file1
endif else begin
  zinfo = ncdf_varinq(ncid1,'Z')
  idz = ncdf_vardef(nciddiffs,zinfo.name,zinfo.dim,/double)
  for ii=0,zinfo.natts-1 do begin
    attname = ncdf_attname(ncid1,idz1,ii)
    if (attname eq 'long_name') then begin
      zname = "Perturbed Geopotential Height"
      ncdf_attput,nciddiffs,idz,attname,zname,/char
    endif else begin
      iatt = ncdf_attcopy(ncid1,idz1,attname,nciddiffs,idz)
    endelse
  endfor
endelse
;
; Take new diffs file out of define mode, into data mode:
;
ncdf_control,nciddiffs,/endef
;
; Provide values of coord vars from file1:
;
for i=0,file1struct.ndims-1 do begin
  ncdf_diminq,ncid1,i,name,size    ; dim on file1
  idv = ncdf_varid(nciddiffs,name) ; is there a coord var w/ this dim name?
  if (idv ne -1) then begin        ; copy it from file1 to filediffs
    ncdf_varget,ncid1,name,cvar
    ncdf_varput,nciddiffs,idv,cvar
  endif
endfor
;
; Put mtime var on filediffs (was read from file1 above):
;
id = ncdf_varid(nciddiffs,'mtime')
ncdf_varput,nciddiffs,id,mtime
id = ncdf_varid(nciddiffs,'year')
ncdf_varput,nciddiffs,id,year
id = ncdf_varid(nciddiffs,'day')
ncdf_varput,nciddiffs,id,day
id = ncdf_varid(nciddiffs,'ut')
ncdf_varput,nciddiffs,id,ut
;
; Put full field Z on filediffs from file1 (ids were set above):
;
ncdf_varget,ncid1,idz1,z
ncdf_varput,nciddiffs,idz,z
;
; Loop over fields, calculate diffs, and put diffs on output file:
;
for i=0,nflds-1 do begin
  ncdf_varget,ncid1,fields[i],f1
  ncdf_varget,ncid0,fields[i],f0
  diffs = f1-f0
  percent = (f1-f0)/f0*100.
  minpercent = min(percent)
  maxpercent = max(percent)
  if (percent_diffs le 0) then begin
    print,format="('Field ',a16,' f0=',2e12.4,' f1=',2e12.4,' raw diffs=',2e12.4)",$
      fields[i],min(f0),max(f0),min(f1),max(f1),min(diffs),max(diffs)
  endif else begin
    print,format="('Field ',a16,' f0=',2e12.4,' f1=',2e12.4,' %diffs=',2e12.4)",$
      fields[i],min(f0),max(f0),min(f1),max(f1),minpercent,maxpercent
  endelse
;
; Put diffs on diffs file:
;
  if (fields[i] eq 'Z') then begin
    id = ncdf_varid(nciddiffs,'ZDIFFS')
  endif else begin
    id = ncdf_varid(nciddiffs,fields[i])
  endelse
  if (id ne -1) then begin
    if (percent_diffs le 0) then ncdf_varput,nciddiffs,id,diffs $
    else ncdf_varput,nciddiffs,id,percent
  endif
endfor
;
; Close out files:
;
ncdf_close,ncid0
ncdf_close,ncid1
ncdf_close,nciddiffs
if (percent_diffs le 0) then $
     print,'Wrote output file with raw differences to ',filediffs $
else print,'Wrote output file with percent differences to ',filediffs
end