c read wind (m/s), ice velocity(m/s), and surface ocean velocity (cm/s)
      parameter (im=360,jm=120,ix=im,iy=jm,KM=30,m1=12) 
      parameter (imt=im,jmt=jm,nx1=imt,ny1=jmt,nx=nx1,ny=ny1) 
      parameter (nyear1=2000,nyear2=2009,nday=365)
      common/bnd/ibnd(im,jm),kmt(im,jm),kmu(im,jm)
      dimension heff(im,jm)
      DIMENSION psurf(IM,JM),ssh(IM,JM),temp(IM,JM,km),sali(IM,JM,km)
      DIMENSION uo(IM,JM,km),vo(IM,JM,km),RHOP(IM,JM,km)
      DIMENSION uo1(IM,JM,km),vo1(IM,JM,km)
      DIMENSION ai(imt,jmt),hi(imt,jmt)
      DIMENSION gice(imt,jmt,m1)
      DIMENSION uice(imt,jmt),vice(imt,jmt)
      DIMENSION uo7(imt,jmt),vo7(imt,jmt)
      DIMENSION fiix(imt,jmt),fiiy(imt,jmt)
      DIMENSION ostrx(imt,jmt),ostry(imt,jmt)
      DIMENSION airstrxo(imt,jmt),airstryo(imt,jmt)
      DIMENSION s11(imt,jmt),s12(imt,jmt),s22(imt,jmt)
      DIMENSION uwind(imt,jmt),vwind(imt,jmt)
      DIMENSION clon(imt,jmt),clat(imt,jmt)
      DIMENSION sali1(IM,JM,km),salim(IM,JM,km,12)
      dimension ulat(imt,jmt),ulon(imt,jmt),HTN(imt,jmt),HTE(imt,jmt)
     &,HUS(imt,jmt),HUW(imt,jmt),angle(imt,jmt),dxt(imt,jmt)
     &,dyt(imt,jmt),alpha(imt,jmt)
      dimension dz(km),ZDZ (KM),ZDZZ(KM+1)
      dimension day_month(12), month_end(13)
      data day_month/31,28,31,30,31,30,31,31,30,31,30,31/
      data month_end/0, 31,59,90,120,151,181,212,243,273,304,334,365/
      character *80 fn(9)
      character *99 mask,fn1,fn2,fn3,fn4
      character *4 cyear(1900:2200),cyear1(1900:2200)
      character *2 cmon,cmon1
      character *12 char
      integer SLEN

      fn(1)='aiday.H'          ! daily ice concentration (frac)
      fn(2)='hiday.H'          ! daily mean ice thickness (m)
      fn(3)='giceday.H'        ! daily ice thickness distribution (frac) 
      fn(4)='uiday.H'          ! daily ice velocity (m/s)
      fn(5)='uo7day.H'         ! ocean velocity at 7th lev (30-36 m) (cm/s)
      fn(6)='fiiday.H'         ! daily ice interaction force(N/m*m)
      fn(7)='ostrday.H'        ! daily ocean surface stress (dyn/cm*cm)
      fn(8)='airstroday.H'     ! air stress on openwater (tauo*g(1)) (N/m*m)
      fn(9)='stressday.H'      ! Ice internal stress (N/m)

c read lon and lat for scalar fields (like sea ice thickness and concentration)
      open(20,file='grid.dat.rot')
      read(20,'(10f8.2)') ((clon(i,j),i=1,nx1),j=1,ny1)
      read(20,'(10f8.2)') ((clat(i,j),i=1,nx1),j=1,ny1)
      close(20)

      open(24,file='grid.dat.pop')
c ulat and ulon are lat and lon for vector fields (like sea ice and ocean veclocities)
        read(24,'(10f8.2)') ((ulat(i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((ulon(i,j),i=1,nx),j=1,ny)
c HTN, HTE are lengths of the northern and eastern sides of a scaler grid cell in km, HTN*HTE is the area of a scaler grid cell in km**2
c HUS, HUW are lengths of the southern and western sides of a vector grid cell in km
        read(24,'(10f8.2)') ((HTN  (i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((HTE  (i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((HUS  (i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((HUW  (i,j),i=1,nx),j=1,ny)
c angle is the angle between latitude line and  grid cell x-coordinate line, needed for vector rotation for plotting the vectors in spherical lat-lon coordinate system
c For example, to rotate surface ocean veclocity to a spherical lat-lon coordinate system, do the following:
c uo1(:,:,1)=uo(:,:,1)*cosd(-angle)+vo(:,:,1)*sind(-angle)
c vo1(:,:,1)=vo(:,:,1)*cosd(-angle)-uo(:,:,1)*sind(-angle)
c where uo and vo are based on model coordinate system, while uo1 and vo1 are based on spherical lat-lon coordinate syetem
        read(24,'(10f8.2)') ((angle(i,j),i=1,nx),j=1,ny)
      close(24)

c read model grid mask with ocean levels; ocean: levels>0; land: levels=0
c ocean levels can also be used to calculate ocean depth in the model
      open(20,file='io.dat_360_120.output')
      read(20,'(360i2)') kmt
      close(20)

c kmu is model grid mask for vector fields; ocean > 0, land = 0
c kmu has the number of ocean levels
c kmu can be used with dz to compute ocean depth for vector fields
      DO 730 J=1,ny-1
      DO 730 I=1,nx-1
      KMU(I,J)=MIN(kmt(I,J),kmt(I+1,J),kmt(I,J+1),kmt(I+1,J+1))
 730  CONTINUE
      DO J=1,ny-1
      DO I=nx,nx
      KMU(I,J)=MIN(kmt(I,J),kmt(1,J),kmt(I,J+1),kmt(1,J+1))
      end do
      end do


c read ocean level thickness dz(k) in cm
c dz can be used with kmt or kmu to calculate ocean depth
      open(20,file='dz.dta30')
      do k=1,km
        read(20,*)dz(k)
        dz(k)=dz(k)*0.01
      end do
      ZDZ(1)=DZ(1)
      do k=2,km
        ZDZ(K)=ZDZ(K-1)+DZ(K)
      end do
      do k=1,km
        zdzz(k)=zdz(k)-0.5*dz(k)
      end do
      close(20)

c  begin yearly loop
      do 999 iyear=nyear1,nyear2

      write(unit=cyear(iyear),fmt='(i4)') iyear
      write(unit=cyear1(iyear),fmt='(i4)') iyear+1

      do iii=1,9
        i=slen(fn(ii))
        open(ii+10,file=fn(ii)(1:i)//cyear(iyear)
     &  ,access='direct',form='unformatted',recl=nx1*ny1*4
     &  ,status='old')
      end do

      do 888 iday=1,nday

        read(11)((ai(i,j),i=1,nx1),j=1,ny1)
        read(12)((hi(i,j),i=1,nx1),j=1,ny1)
        do k=1,m1
          read(13)((gice(i,j,k),i=1,nx1),j=1,ny1)
        end do
        read(14)((uice(i,j),i=1,nx1),j=1,ny1)
        read(14)((vice(i,j),i=1,nx1),j=1,ny1)
        read(15)((uo7(i,j),i=1,nx1),j=1,ny1)
        read(15)((vo7(i,j),i=1,nx1),j=1,ny1)
        read(16)((fiix(i,j),i=1,nx1),j=1,ny1)
        read(16)((fiiy(i,j),i=1,nx1),j=1,ny1)
        read(17)((ostrx(i,j,1),i=1,nx1),j=1,ny1)
        read(17)((ostry(i,j,1),i=1,nx1),j=1,ny1)
        read(18)((airstrxo(i,j,1),i=1,nx1),j=1,ny1)
        read(18)((airstryo(i,j,1),i=1,nx1),j=1,ny1)
        read(19)((s11(i,j,1),i=1,nx1),j=1,ny1)
        read(19)((s12(i,j,1),i=1,nx1),j=1,ny1)
        read(19)((s22(i,j,1),i=1,nx1),j=1,ny1)

      i=imt/2
      j=jmt/2
      
      write(*,'(2i5,6f10.3)') iyear,iday
     &,ai(i,j),hi(i,j),uice(i,j),vice(i,j),fiix(i,j),fiiy(i,j)

888   continue
      do ii=1,9
        close(ii+10)
      end do
999   continue
      
      stop
      end

      INTEGER FUNCTION slen (string)
C ---
C --- this function computes the length of a character string less
C --- trailing blanks
C --- slen > 0, length of string less trailing blanks
C ---      = 0, character string is blank
C ---
      CHARACTER*(*) string
      CHARACTER*1 cblank
      INTEGER i
      DATA cblank/' '/
C ---
      DO 50 i = LEN(string), 1, -1
         IF (string(i:i) .NE. ' ')  GO TO 100
50    CONTINUE
      i = 0
100   CONTINUE
      slen = i
      RETURN
      END
