undef("coordcheck")
function coordcheck(x[*][*][*][*]:numeric)
begin 
; some versions of the model just have "degrees" for the units attribute
; of the lat/lon coordinate variable. This causes the gsn templates
; to produce a warning message. 

; this function adds the correct attribute to the lat/lon coordinates to
; remove the message. 

; fix attributes
  if(x&lat@units.ne."degrees_north")then
    x&lat@units="degrees_north"
  end if
  if(x&lon@units.ne."degrees_east")then
    x&lon@units="degrees_east"
  end if
  return(x)
end
;*****************************************************
undef("isDensity")
function isDensity(name:string)
local areDensity
begin 
  areDensity = (/"O2","OX","N4S","NOZ","CO","CO2","H2O","H2","HOX","OP",\
  "CH4","AR","HE","NAT","O21D","NO2","NO","O3","O1","OH","HO2","H","N2D",\
  "NE","O2P","H3P","HP","H2V1","H2V2","H2V3","H2V4","H2P"/)
  if(any(name.eq.areDensity))then
    return(True)
  else
    return(False)
  end if
end
;*********************************************************
undef("interp_middle")
; internal to halfzp
function interp_middle(x[*][*][*][*]:numeric,top:integer)
local k
begin
  do k = 1,top   ; starts at 4.5 vice 5.0 = x(2)
    x(:,k,:,:) = (/0.5*(x(:,k,:,:)+x(:,k+1,:,:))/)
  end do
  return(x)
end
;********************************************************
undef("extrap_top")
function extrap_top(name:string,x[*][*][*][*]:numeric,top:integer)
local xtop1,indices,xtop2,dims
begin
  x(:,top-1,:,:) = 2*x(:,top-1,:,:)-x(:,top-2,:,:)
  if(isDensity(name))then
    if (any(x(:,top-1,:,:).lt.0.))then       ; only execute if necessary 
      xtop1   = ndtooned(x(:,top-1,:,:))   
      indices = ind(xtop1.lt.0)      
      if(.not.any(ismissing(indices))) then 
        xtop2 = ndtooned(x(:,top-2,:,:))          
        xtop1(indices) = xtop2(indices) 
        dims = dimsizes(x)                      
        x(:,top-1,:,:) = (/onedtond( xtop1, (/dims(0),dims(2),dims(3)/))/)
      end if    
    end if 
  end if
  return(x)
end 
;*****************************************************
undef("halfzp")
function halfzp(name:string,x[*][*][*][*]:numeric)
local noproc,dims,nlev
begin 
; the model data is listed as being on the zp levels, when in 
; reality it is output on halfzp levels. This function puts the
; data on the full zp level by averaging the level above and below
; the requested level.

; note we use (/.../) syntax when copying one value to another in order
; not to effect the coordinate variables.
  

; variables we do not convert to full zp. we only check the units on the
; corresponding lat/lon coordinate arrays.   
  noproc = (/"Z","W","NE","POTEN"/)
  if(any(name.eq.noproc))then
    if(name.eq."Z")then
      x = (/x/100000./)                    ; unit conversion
      x@units = "km"
    end if
    x   = coordcheck(x)
    return(x)
  end if


; variables in which the top boundary is processed
  top    = (/"O2","OX","N4S","NOZ","CO","CO2","H2O","H2","HOX","CH4","AR",\
  "HE","NAT"/)

; variables in which the top and bottom boundary are processed.
  topbot  = (/"O21D","NO2","NO","O3","O1","OH","HO2","H","N2D","TI",\
  "TE"/)
  
  dims = dimsizes(x)
  nlev = dims(1)
;*****************************************************  
; we must process this from top to bottom, so we will reorder the array
; bottom [-17,-16.5,-16,-15.5.......4.0,4.5,5.0] top  ===>
; bottom [5.0,4.5,4.0,......-15.5,16,-16.5,-17] top
  x = x(:,::-1,:,:)
  print("converting "+name+" to full zp")
;*************************************
; TN UN VN
;*************************************
  if(name.eq."TN".or.name.eq.("UN").or.name.eq."VN")then
; the top and bott boundaries must be handled first and in that order!! They
; are involved in the middle calculation. The top most value is taken from the 
; bottom. This is a historical artifact that stored the array this way.

    x(:,nlev-1,:,:)=(/x(:,0,:,:)/)       ; top boundary

; The bottom boundary is taken from the first value up from the boundary
    x(:,0,:,:) = (/x(:,1,:,:)/)

    x = interp_middle(x,nlev-2)          ; interpolate the middle

;  back to original ordering
    x = x(:,::-1,:,:)
    if(name.eq."UN".or.name.eq."VN")then
      x = (/ x/100. /)                   ; unit conversion
      x@units = "m/s"
    end if
    x   = coordcheck(x)                  ; adjust unit atts on lat/lon arrays
    return(x)
  end if
;******************************************
; OP  O2P   : log interpolation
;******************************************
  if(name.eq."OP".or.name.eq."O2P")then
    if(name.eq."O2P")then
      x=(/mask(x,x.lt.0,1.e-12)/)
    end if
; log interpolate body 
; bottom boundary
    x(:,0,:,:)=(/sqrt(x(:,1,:,:)^3/x(:,2,:,:))/)
; middle
    do k=1,nlev-2
      x(:,k,:,:)=(/sqrt(x(:,k,:,:)*x(:,k+1,:,:))/)
    end do
; top boundary
    x(:,nlev-1,:,:) = (/x(:,nlev-1,:,:)^2/x(:,nlev-2,:,:)/)
;  back to original ordering
    x=x(:,::-1,:,:)
    x=coordcheck(x)
    return(x)
  end if
;******************************************
; top/bottom extrapolation
; "O21D","NO2","NO","O3","O1","OH","HO2","H","N2D","TI","TE"
;******************************************
  if(any(name.eq.topbot))then
    x = interp_middle(x,nlev-2)
    x = extrap_top(name,x,nlev-1)
    
; bottom boundary extrapolation
    x(:,0,:,:)=(/1.5*x(:,1,:,:)-0.5*x(:,2,:,:)/)
;  back to original ordering
    x=x(:,::-1,:,:) 
    x = coordcheck(x)
    return(x)
  end if
;*******************************************
; top extrapolation
; "O2","OX","N4S","NOZ","CO","CO2","H2O","H2","HOX","CH4","AR","HE","NAT"/)
;*******************************************
  if(any(name.eq.top))then
; bottom and middle x(:,1:nlev-2,:,:)
    x = interp_middle(x,nlev-2)
; top boundary extrapolation (recall that array is reversed)
    x(:,nlev-1,:,:)=(/2*x(:,nlev-1,:,:)-x(:,nlev-2,:,:)/)
;  back to original ordering
    x=x(:,::-1,:,:) 
    x=coordcheck(x)
    return(x)
  end if
;*******************************************
; unknown fields: may be user defined etc. we
; don't process these
;*******************************************
  x = coordcheck(x)                      
  return(x)
end