program reading
       implicit none
       integer idim,jdim,kdim,tdim
       parameter(idim=512,jdim=386,kdim=30,tdim=12)
       real height(40)
       data height /-5.01, -15.07, -25.28, -35.76, -46.61,
     &              -57.98, -70.02, -82.92, -96.92, -112.32,
     &              -129.49, -148.96, -171.4, -197.79, -229.48,
     &              -268.46, -317.65, -381.39, -465.91, -579.3,
     &              -729.35, -918.37, -1139.15, -1378.57, -1625.7,
     &              -1875.11, -2125.01, -2375., -2625., -2875.,
     &              -3125., -3375., -3625., -3875., -4125.,
     &              -4375., -4625., -4875., -5125., -5375./
       real*8 sst8(idim,jdim,kdim)
       real*8 h8(513,386),mask8(idim,jdim)
       real   h4(idim,jdim),srho(kdim),mask(idim,jdim)
       real   sst4(idim,jdim,kdim)
       real   sh(idim,jdim,kdim)
       real   sst1000(idim,jdim)
       integer iyear,imth,lev,idend,idy
       character*4 year
       character*2 month,day

       integer i,j,k,t

       print*,height

       open(11,file='../../h-rho.dat',form='unformatted')
       read(11)h8
       close(11)
       do j=1,jdim
        do i=1,idim
           h4(i,j)=h8(i,j)
        enddo
       enddo

       open(11,file='../../masku.dat',form='unformatted')
       read(11)mask8
       close(11)
       do j=1,jdim
        do i=1,idim
           mask(i,j)=mask8(i,j)
        enddo
       enddo
       call zs(h4,sh)
        do iyear=1995,1995
         write(year,'(i4)')iyear
           do imth=1,tdim
C ---          do imth=3,tdim
          write(month,'(i2)')imth
        if((imth.eq.1) .OR. (imth.eq.3) .OR. (imth.eq.5).OR.(imth.eq.7)
     1    .OR.(imth.eq.8) .OR. (imth.eq.10) .OR. (imth.eq.12)) idend=31
        if((imth.eq.4) .OR. (imth.eq.6) .OR. (imth.eq.9).OR.
     1     (imth.eq.11)) idend=30
          IF (imth.eq.2)then
           IF(MOD(iyear,100).NE.0.AND.MOD(iyear,4).EQ.0) THEN
             WRITE(*,*) 'LEAP YEAR'
             idend=29
            ELSEIF(MOD(iyear,400).EQ.0) THEN
             WRITE(*,*) 'LEAP YEAR'
             idend=29
            ELSE
             WRITE(*,*) 'NOT A LEAP YEAR'
             idend=28
           ENDIF
           ENDIF
           print*,idend
          if(imth.lt.10) month(1:1)='0'
          do idy=1,idend
          write(day,'(i2)')idy
          if(idy.lt.10) day(1:1)='0'
          open(11,file='../uraw/ucur-'//year//month//day//'-s.dat',
     1         form='unformatted')
          open(51,file='u-'//year//month//day//'.dat',
     1         form='unformatted')
          read(11)sst8
          do k=1,kdim
           do j=1,jdim
            do i=1,idim
             sst4(i,j,k)=sst8(i,j,k)
            enddo
           enddo
          enddo
          do lev=1,40
            call vinterp(height(lev),mask,sh,sst4,sst1000)
            write(51)sst1000
          enddo ! lev
          close(11)
          close(51)
         enddo  ! idy
        enddo  ! tdim
       enddo    ! iyear
      end
      subroutine zs(h, z_r)
       implicit none
       integer Lp, Mp, NL
       parameter(Lp=512,Mp=386,NL=30)
       real    ds, hc, hmin, Tcline, theta_s, theta_b
       parameter(hmin=30, Tcline=75, theta_s=5.0, theta_b=0.4)
       real    hinv, z_r0, z_r(Lp*Mp*NL), h(Lp*Mp), zeta(Lp*Mp)
       real    cff1, cff2
       real    sc_w(nl), Cs_w(nl), sc_r(nl), Cs_r(nl)
       integer i, j, k, ij, ip
!-----------------------------------------------------------------------
!  Depths from S-coordinates transformation.
!-----------------------------------------------------------------------
!  Define S-Curves in domain [-1 < sc < 0] at vertical W- and
!  RHO-points.
        ds=1.0/real(NL)
        hc=MIN(hmin,Tcline)
            print*, 'hmin, hc=', hmin, hc
        cff1=1.0/SINH(theta_s)
        cff2=0.5/TANH(0.5*theta_s)
        sc_w(0)=-1.0
        Cs_w(0)=-1.0
        do k=1,NL
          sc_w(k)=ds*real(k-NL)
          Cs_w(k)=(1.0-theta_b)*cff1*SINH(theta_s*sc_w(k))+
     &            theta_b*(cff2*TANH(theta_s*(sc_w(k)+0.5))-0.5)
          sc_r(k)=ds*(real(k-NL)-0.5)
          Cs_r(k)=(1.0-theta_b)*cff1*SINH(theta_s*sc_r(k))+
     &            theta_b*(cff2*TANH(theta_s*(sc_r(k)+0.5))-0.5)
        enddo
        do j=1,Mp
         do i=1,Lp
           ij=i+(j-1)*Lp
           zeta(ij)=0.
         enddo
        enddo
!
!  Compute model grid depths (meters) at RHO-points.
!
        do k=1,NL
          cff1=hc*(sc_r(k)-Cs_r(k))
          cff2=Cs_r(k)
          do j=1,Mp
            do i=1,Lp
              ij=i+(j-1)*Lp
              ip=ij+(k-1)*Lp*Mp
              hinv=1.0/h(ij)
              z_r0=cff1+cff2*h(ij)
              z_r(ip)=z_r0+zeta(ij)*(1.0+z_r0*hinv)
            enddo
          enddo
        enddo
        print*,'lihq end of z_s'
c           print*,z_r
!
      endsubroutine zs

      subroutine vinterp(height,mask,sh,input,output)
       implicit none
       integer idim,jdim,kdim,level
       parameter(idim=512,jdim=386,kdim=30)
       real    height,sh(idim,jdim,kdim),mask(idim,jdim)
       real    input(idim,jdim,kdim), output(idim,jdim)
       real    d1,d2
       integer i,j,k

       do j=1,jdim
        do i=1,idim
         if(mask(i,j).eq.1)then
          if(sh(i,j,1).ge.height) then
            output(i,j)=-9999.
          elseif(sh(i,j,kdim).le.height) then
            output(i,j)=input(i,j,kdim)
          else
           do k=2,kdim
            if((sh(i,j,k).ge.height).and.(sh(i,j,k-1).lt.height))then
              level=k
              goto 1000
            endif
           enddo
1000       continue
           d1=sh(i,j,level-1)-height
           d2=height-sh(i,j,level)
           output(i,j)=input(i,j,level-1)*d2/(d1+d2) +
     1                 input(i,j,level  )*d1/(d1+d2)
          endif  ! height vs. sh
         else    ! mask
            output(i,j)=-9999.
         endif   ! mask
        enddo
       enddo
      endsubroutine vinterp