!m_vqs: # of master grids
      m_vqs=8
      dz_bot_min=4 !min. bottom layer thickness [m]
      allocate(hsm(m_vqs),nv_vqs(m_vqs),a_vqs(m_vqs))
      hsm=(/5,10,50,100,350,1050,2000,5001/)
      nv_vqs(1:5)=(/3,4,9,12,21/) !# of levels for each master grid (increasing with depth)
      nv_vqs(6)=nv_vqs(5)+9
      nv_vqs(7)=nv_vqs(6)+5
      nv_vqs(8)=nv_vqs(7)+5

      if(m_vqs<2) then 
        write(*,*)'Check vgrid.in:',m_vqs
        stop
      endif
      if(hsm(1)<0) stop 'hsm(1)<0'
      do m=2,m_vqs
        if(hsm(m)<=hsm(m-1)) then
          write(*,*)'Check hsm:',m,hsm(m),hsm(m-1)
          stop
        endif
      enddo !m

!     Other consts.
!     Stretching const. for the 1st master grid and also for depth <= hsm(1)
!     |a_vqs0|<=1 (1: skew toward bottom; -1: toward surface; 0: no bias)
      a_vqs0=0

!     Generate a master vgrid (z_mas)
      etal=0 !used in master grid only; elev.
      if(etal<=-hsm(1)) then
        write(*,*)'elev<hsm:',etal
        stop
      endif

      nvrt_m=nv_vqs(m_vqs)
      print*, 'nvrt in master vgrid=',nvrt_m
      allocate(z_mas(nvrt_m,m_vqs))
      z_mas=-1.e5
      do m=1,5 !m_vqs
        do k=1,nv_vqs(m)
          sigma=(k-1.0)/(1.0-nv_vqs(m))

!         Alternative transformations below
!         Option 1: quadratic 
!          a_vqs(m)=max(-1.d0,a_vqs0-(m-1)*0.03)
!          tmp=a_vqs(m)*sigma*sigma+(1+a_vqs(m))*sigma !transformed sigma
!          z_mas(k,m)=tmp*(etal+hsm(m))+etal

!          Option 2: S
          theta_b=0
          if(m<=4) then
            theta_f=2.5
          else if(m==5) then
            theta_f=2.65
          endif
          cs=(1-theta_b)*sinh(theta_f*sigma)/sinh(theta_f)+ &
     &theta_b*(tanh(theta_f*(sigma+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5)
          z_mas(k,m)=etal*(1+sigma)+hsm(1)*sigma+(hsm(m)-hsm(1))*cs
        enddo !k
      enddo !m

      !Last 2 master grids
      z_mas(1:nv_vqs(5),6)=z_mas(1:nv_vqs(5),5)
      z_mas(1:nv_vqs(6),7)=z_mas(1:nv_vqs(5),5)
      z_mas(1:nv_vqs(6),8)=z_mas(1:nv_vqs(5),5)
      z_mas(1+nv_vqs(5):nv_vqs(6),6)=(/-400,-460,-520,-590,-660,-740,-830,-930,-1050/)
      z_mas(1+nv_vqs(5):nv_vqs(7),7)=(/-400,-460,-520,-590,-660,-740,-830,-930,-1050,-1200,-1400,-1600,-1800,-2000/)
      z_mas(1+nv_vqs(5):nv_vqs(8),8)=(/-400,-460,-520,-590,-660,-740,-830,-930,-1050,-1200,-1400,-1600,-1800,-2000,-2400,-2900,-3500,-4200,-5001./)