c Code to show how to read PIOMAS files in binary format and model grid
c information files
c
c advect.H: monthly change in ice thickness due to convergence (meter of ice/s)
c iceprod.H: monthly net ice production (m of ice/s)
c oflux.H: monthly ocean heat flux (m of ice/s) **do not use**
c tice0.H: temperature at the surface of ice, snow, or open ocean (K)
c gice.H: monthly 12-category ice thickness distribution (area fraction)
c heff.H:    monthly ice thickness (m)
c area.H:    month ice concentration (fraction)
c icevel.H:  monthly ice velocity (m/s)
c ssh.H:     monthly SSH (cm)
c uo1_10.H:     upper 10 levels of monthly ocean velocity (cm/s)
c otemp1_10.H:  upper 10 levels of monthly ocean temperature (C)
c osali1_10.H:  upper 10 levels of monthly ocean salinity (multiplied by 1000 to get psu)

      parameter(nx1=360,ny1=120,nx=nx1,ny=ny1,imt=nx1,jmt=ny1,km=30)
      parameter(m1=12)
      parameter(nyear1=1980,nyear2=2023,nday=365,nmon=12)

      dimension advect(imt,jmt),prod(imt,jmt),oflux(imt,jmt)
      dimension tice0(imt,jmt)
      dimension heff(imt,jmt),area(imt,jmt),gice(imt,jmt,m1)
      dimension uice(imt,jmt),vice(imt,jmt),ssh(imt,jmt)
      dimension uo(imt,jmt,km),vo(imt,jmt,km)
      dimension uo1(imt,jmt,km),vo1(imt,jmt,km)
      dimension to(imt,jmt,km),so(imt,jmt,km)
      dimension clon(imt,jmt),clat(imt,jmt),kmt(imt,jmt)
      dimension kmu(imt,jmt),iwork(imt,jmt)
      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)
      dimension dz(km),ZDZ (KM),ZDZZ(KM+1)
      character *80 fopen(5),f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12
      character *4 cyear(1900:2100),cyear1(1900:2100)
      character *12 char
      integer SLEN

c advect.H: monthly change in ice thickness due to convergence (meter of ice/s)
c iceprod.H: monthly net ice production (m of ice/s)
c oflux.H: monthly ocean heat flux (m of ice/s)
c tice0.H: temperature at the surface of ice, snow, or open ocean (K)
c gice.H: monthly 12-category ice thickness distribution (area fraction)
c heff.H:    monthly ice thickness (m)
c area.H:    month ice concentration (fraction)
c icevel.H:  monthly ice velocity (m/s)
c ssh.H:     monthly SSH (cm)
c uo1_10.H:     first 10 levels of monthly ocean velocity (cm/s)
c otemp1_10.H:  first 10 levels of monthly ocean temperature (C)
c osali1_10.H:  first 10 levels of monthly ocean salinity (multiplied by 1000 to get psu)

      f1='advect.H'
      f2='iceprod.H'
      f3='oflux.H'
      f4='tice0.H'
      f5='gice.H'
      f6='heff.H'
      f7='area.H'
      f8='icevel.H'
      f9='ssh.H'
      f10='uo1_10.H'
      f11='otemp1_10.H'
      f12='osali1_10.H'

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)

c read lon and lat for vector fields (like sea ice and ocean veclocities)
      open(24,file='grid.dat.pop')
        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 do rotation like this for vector uo and vo
c        uo1(:,:,k)=uo(:,:,k)*cosd(-angle)+vo(:,:,k)*sind(-angle)
c        vo1(:,:,k)=vo(:,:,k)*cosd(-angle)-uo(:,:,k)*sind(-angle)
        read(24,'(10f8.2)') ((angle(i,j),i=1,nx),j=1,ny)
      close(24)

c read model scaler grid mask; ocean levels > 0, land = 0
      open(20,file='io.dat_360_120.output')
      read(20,'(360i2)') kmt
      close(20)
c get model vector grid mask; ocean levels > 0, land = 0
      IWORK=eoshift(KMT,dim=1,shift=1)
      KMU = min(KMT,IWORK)
      IWORK=eoshift(KMU,dim=2,shift=1)
      KMU = min(KMU,IWORK)  ! KMU = minimum of surrounding KMTs

c read ocean level thickness dz(k) in cm
c zdz(kmt(i,j)) is the ocean depth at scaler grid cell (i,j)
c zdz(kmu(i,j)) is the ocean depth at vector cell (i,j)
      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)

      do 999 iyear=nyear1,nyear2

      write(unit=cyear(iyear),fmt='(i4)') iyear

      i=slen(f1)
      open(1,file=f1(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')

      i=slen(f2)
      open(2,file=f2(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')

      i=slen(f3)
      open(3,file=f3(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')

      i=slen(f4)
      open(4,file=f4(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
      i=slen(f5)
      open(5,file=f5(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
      i=slen(f6)
      open(16,file=f6(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
      i=slen(f7)
      open(7,file=f7(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
      i=slen(f8)
      open(8,file=f8(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
      i=slen(f9)
      open(9,file=f9(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
      i=slen(f10)
      open(10,file=f10(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
      i=slen(f11)
      open(11,file=f11(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
      i=slen(f12)
      open(12,file=f12(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')
 
c Note: some machines need a record number to read data like the following
c      irec = 0
c      do imon=1,nmon
c        do k=1,m1 
c          irec = irec+1
c          read(5,rec=irec)((gice(i,j,k),i=1,nx1),j=1,ny1)
c        end do
c      end do

c The following is for machines that do not need a record number to read data
      do imon=1,nmon
        read(1)((advect(i,j),i=1,nx1),j=1,ny1)
        read(2)((prod(i,j),i=1,nx1),j=1,ny1)
        read(3)((oflux(i,j),i=1,nx1),j=1,ny1)
        read(4)((tice0(i,j),i=1,nx1),j=1,ny1)
        do k=1,m1 
          read(5)((gice(i,j,k),i=1,nx1),j=1,ny1)
        end do
        read(16)((heff(i,j),i=1,nx1),j=1,ny1)
        read(7)((area(i,j),i=1,nx1),j=1,ny1)
        read(8)((uice(i,j),i=1,nx1),j=1,ny1)
        read(8)((vice(i,j),i=1,nx1),j=1,ny1)
        read(9)((ssh(i,j),i=1,nx1),j=1,ny1)
c        do k=1,km
        do k=1,10
          read(10)((uo(i,j,k),i=1,nx1),j=1,ny1)
        end do
c        do k=1,km
        do k=1,10
          read(10)((vo(i,j,k),i=1,nx1),j=1,ny1)
        end do
c        do k=1,km
        do k=1,10
          read(11)((to(i,j,k),i=1,nx1),j=1,ny1)
        end do
c        do k=1,km
        do k=1,10
          read(12)((so(i,j,k),i=1,nx1),j=1,ny1)
        end do
      end do

      i=nx/2
      j=ny/2
      write(*,'(2i6, 2f12.3, 2f8.2)') iyear, iday
     &, clon(i,j), clat(i,j), uo(i,j,1), heff(i,j)

      close(2)
      close(3)
      close(4)
      close(5)
      close(16)
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
