Contact person:  Michael Steele, University of Washington, mas@apl.washington.edu, 
phone:  (206) 543-6586.  October, 1999.

This folder contains 2 data files (nfr97_data and nfr98_data) and 2 postscript 
files (nfr97.ps and nfr98.ps) that show the monthly net freezing rate (nfr) over 
the entire Arctic Ocean, Barents Sea, and Greenland-Iceland-Norwegian Seas for the 
years 1997 and 1998.  The data have been generated using a numerical sea ice - 
ocean model described by Zhang et al. (Geophys. Res. Lett., 25, 1745-1748, 1998).  
Forcing was provided by the NCEP/NCAR reanalysis.  Net freezing rate is defined as 
the net growth minus melt over all ice thickness categories, where positive values 
represent net growth.

The postscript files are (hopefully) self-explanatory.  The data are provided as
unformatted binary files and are arranged on a cartesian (x,y) grid based on a 
polar stereographic projection.  The data may be read using the following IDL code 
snippets.

-------------------------

nx=130
ny=102

month_data(nx,ny)
all_data(12,nx,ny)

for each month, read in a nx*ny grid;
for i=0,11 do begin
readu,"filename",Month_data
all_dat(i,*,*)=month_data(*,*)
endfor

Each grid cell is 40km x 40km.  To convert from grid points (i,j) to kilometers 
(x,y): 

for i=0,nx-1 do begin
   x(i)=(i-62.5)*40.0
endfor
for j=0,ny-1 do begin
   y(j)=(j-54.5)*40.0
endfor

The model domain is drawn on SSM/I polar stereographic projection, rotated 55 
degrees, so that the positive X axis is along 55E longitude.  To convert from 
(x,y) to (lat,lon), use the IDL program XYtoLL.  The inverse is provided by 
LLtoXY.  

-------------------------------------

;  *** XXtoLL ***
;
;  Polar stereo projection of lat,lon onto an x,y grid in kilometers
;  (From Fortran subroutine rmaps.f)
;
;  Here we assume northern hemisphere!
;  This version allows
;  vector arguments for x and y, with vectors alat,alon returned.
; 
;  We use rot=55, which rotates our maps so that North America is at the bottom.
;
;  !DTOR = "degrees to radians", ie pi/180
;  !RADEG = "radians to degrees" ie 180/pi
;
PRO XYtoLL, x, y, alat, alon, rot

   nel = n_elements(x)
   alat = fltarr(nel)
   alon = fltarr(nel)
   tt = fltarr(nel)
   chi = fltarr(nel)

;  Definition of constants

        re=6378.273		; Radius of earth (km)  -Hughes Ellipsoid
        e2=DOUBLE(0.006693883)	; Eccentricity of earth -Hughes Ellipsoid
        e=sqrt(e2)
        slat=70.		; Standard parallel

;   rho is a vector
        rho = sqrt(x*x + y*y)
;   cm and t are scalars
        cm=COS(!DTOR*slat)/SQRT(1.-e2*(SIN(!DTOR*slat)^2))
        t=TAN((!DPI/4.)-(slat/(2.*!RADEG)))/         $
          ((1.-e*SIN(!DTOR*slat))/                   $
           (1.+e*SIN(!DTOR*slat)))^(e/2)
;   tt and chi are vectors
        tt=rho*t/(re*cm)
        chi=(!DPI/2.0)-2.0*ATAN(tt)

        alat = chi + ((e2/2)+(5*(e2^2)/24)+((e2^3)/12)) * SIN(2*chi) $
                      + ((7*(e2^2)/48)+(29*(e2^3)/240)) * SIN(4*chi) $
                                       + (7*(e2^3)/120) * SIN(6*chi)
        alat = alat * !RADEG
        alon = ATAN(x,-y)*!RADEG - rot
        neg = where(alon LE -180.)
        if min(neg) NE -1 then alon(neg) = alon(neg) + 360.

	RETURN
	END

-------------------------------------------

;  *** LLtoXY ***
;
;  Polar stereo projection of lat,lon onto an x,y grid in kilometers
;  (From Fortran subroutine rmaps.f)
;
;  Here we assume northern hemisphere!
;
;  see other notes above for XXtoLL
;
PRO LLtoXY, alat, alon, x, y, rot

   nel = n_elements(alat)
   x = fltarr(nel)
   y = fltarr(nel)
   t0 = fltarr(nel)
   rho = fltarr(nel)
   NP = lonarr(nel)
   XP = lonarr(nel)

;  Definition of constants

        re=6378.273		; Radius of earth (km)  -Hughes Ellipsoid
        e2=DOUBLE(0.006693883)	; Eccentricity of earth -Hughes Ellipsoid
        e=sqrt(e2)
        slat=70.		; Standard parallel

        NP = where(alat GE 89.995)
        XP = where(alat LT 89.995)

        if min(XP) NE -1 then begin   ; some points are below 89.995

;          t0 is a vector because alat is a vector
           t0(XP)=TAN((!DPI/4.)-(alat(XP)/(2.*!RADEG)))/     $
                  ((1.-e*SIN(!DTOR*alat(XP)))/               $
                   (1.+e*SIN(!DTOR*alat(XP))))^(e/2)
;          t1 is a scalar because slat is a scalar
           t1=    TAN((!DPI/4.)-(slat/(2.*!RADEG)))/         $
                  ((1.-e*SIN(!DTOR*slat))/                   $
                   (1.+e*SIN(!DTOR*slat)))^(e/2)
           cm=COS(!DTOR*slat)/SQRT(1.-e2*(SIN(!DTOR*slat)^2))
           rho(XP)=re*cm*t0(XP)/t1

           x(XP) =  rho(XP)*SIN(!DTOR*(alon(XP)+rot))
           y(XP) = -rho(XP)*COS(!DTOR*(alon(XP)+rot))

        endif

        if min(NP) NE -1 then begin   ; some points are above 89.995
           x(NP) = replicate(0.0, n_elements(NP))
           y(NP) = replicate(0.0, n_elements(NP))
        endif

	RETURN
	END