!for get_vgrid_*

      inquire(file=vgrid,exist=lexist)
      if(lexist) then
        open(19,file=vgrid,status='old')
      else
        inquire(file='../'//vgrid,exist=lexist)
        if(lexist) then
          open(19,file='../'//vgrid,status='old')
        else
          write(*,*)'Unable to find vgrid.in'
          stop
        endif
      endif

      read(19,*)ivcor
      select case(ivcor)
        case(2) !SZ
          read(19,*) nvrt,kz,h_s !kz>=1
          if(nvrt<2) stop 'nvrt<2'
          if(kz<1) then !.or.kz>nvrt-2) then
            write(*,*)'Wrong kz:',kz
            stop
          endif
          if(h_s<6) then
            write(*,*)'h_s needs to be larger:',h_s
            stop
          endif

          ! # of z-levels excluding "bottom" at h_s
          read(19,*) !for adding comment "Z levels"
          do k=1,kz-1
            read(19,*)j,ztot(k)
            if(ztot(k)>=-h_s) then
              print*, 'Illegal Z level:',k
              stop
            endif
            if(k>1) then; if(ztot(k)<=ztot(k-1)) then
              print*, 'z-level inverted:',k
              stop
            endif; endif
          enddo !k
          read(19,*) !level kz       
          ! In case kz=1, there is only 1 ztot(1)=-h_s
          ztot(kz)=-h_s

          nsig=nvrt-kz+1 !# of S levels (including "bottom" & f.s.)
          read(19,*) !for adding comment "S levels"
          read(19,*)h_c,theta_b,theta_f
          if(h_c<5.or.h_c>=h_s) then !large h_c to avoid 2nd type abnormality
            print*, 'h_c needs to be larger avoid 2nd type abnormality; &
     &do u want to continue? Enter 1 to continue, or ctrl-C to abort:'
            read*, itmp
          endif
          if(theta_b<0.or.theta_b>1) then
            write(*,*)'Wrong theta_b:',theta_b
            stop
          endif
          if(theta_f<=0) then
            write(*,*)'Wrong theta_f:',theta_f
            stop
          endif

          sigma(1)=-1 !bottom
          sigma(nsig)=0 !surface
          read(19,*) !level kz
          do k=kz+1,nvrt-1
            kin=k-kz+1
            read(19,*) j,sigma(kin)
            if(sigma(kin)<=sigma(kin-1).or.sigma(kin)>=0) then
              write(*,*)'Check sigma levels at:',k,sigma(kin),sigma(kin-1)
              stop
            endif
          enddo !k
          read(19,*) !level nvrt

        case(1) !localized sigma
          read(19,*)nvrt
          read(19,*)kbp(1:np)
          do i=1,nvrt
            read(19,*)j,sigma_lcl(i,1:np)
          enddo !i
        case default
          write(*,*)'Unknown ivcor:',ivcor
          stop
      end select
      close(19)