c======================================================================|
      subroutine model(idf,ro,pr,vx,vy,vz,bx,by,bz,gm
     &           ,gx,gxm,gy,gym,margin,x,ix,y,jx
     &   ,mf_params)
c======================================================================|
      implicit double precision (a-h,o-z)
c----------------------------------------------------------------------|
      dimension dxm(ix),x(ix)
      dimension dym(jx),y(jx)
      dimension ro(ix,jx),pr(ix,jx),vx(ix,jx)
      dimension vy(ix,jx),bx(ix,jx),by(ix,jx)
      dimension gx(ix,jx),gxm(ix,jx)
      dimension gy(ix,jx),gym(ix,jx)

      dimension bz(ix,jx),vz(ix,jx)

      dimension tem0(jx),pre0(jx),den0(jx),gym0(jx)
      dimension bmx0(jx),rbeta0(jx)

c----------------------------------------------------------------------|
c   parameters
c----------------------------------------------------------------------|
c   gm is gamma( ratio of specific heat )
      gm=1.05d0

c   ztr is height of Transition Region(Corona)
      ztr=18.d0
c   top of calculation box
      ymax=40.d0
c   bottom of calculation box
      ymin=0.d0
c   pi
      pi = acos(-1.0d0)
c----------------------------------------------------------------------|
c   grid
c   xmax is right side of calculation box
c   dx0 is x-grid size
c   dy0 is y-grid size
c----------------------------------------------------------------------|
      xmax=40.d0
      dx0=xmax/real(ix-margin*2-1)
      do i=1,ix
         dxm(i)=dx0
      enddo

      izero=margin+1
      x(izero)=0.
      do i=izero+1,ix
         x(i) = x(i-1)+dxm(i-1)
      enddo
      do i=izero-1,1,-1
         x(i) = x(i+1)-dxm(i)
      enddo

      dy0=(ymax-ymin)/real(jx-margin*2-1)
      do j=1,jx
         dym(j)=dy0
      enddo

      jzero=int(-ymin/dy0)+margin+2
      y(jzero)=0.
      do j=jzero+1,jx
        y(j) = y(j-1)+dym(j-1)
      enddo
      do j=jzero-1,1,-1
        y(j) = y(j+1)-dym(j)
      enddo

c----------------------------------------------------------------------|
c   gravitation
c----------------------------------------------------------------------|
      g0= -1.0d0/gm

      do j=1,jx
        gym0(j)=g0
      enddo

      do j=1,jx
      do i=1,ix
         gx(i,j)=0.0
         gxm(i,j)=0.0
         gy(i,j)=g0
         gym(i,j)=g0
      enddo
      enddo

      call bdspny(0,margin,gy,ix,jx)
      call bdsmny(0,margin,gym,ix,jx)
      call bdcnsy(1,margin,gy,0.d0,ix,jx)
      call bdcnsy(1,margin+1,gym,0.d0,ix,jx)

c-----------------------------------------------------------------------|
c   initial temperature, density, pressure distributions
c-----------------------------------------------------------------------|

c   temperature

c   width of transition region
      wtr=0.6
c   coronal temperature
      tcor=25.0
c   factor of adiabatic temperature gradient
      tadg=0.0
c   height of photoshere or bottom of simulation box
      zcnv= 0.0

      do j=1,jx
            tem0(j)=1.
     &        -tadg*(gm-1.)/gm*(y(j)-zcnv)*0.5*(1-tanh(y(j)-zcnv))
     &        +0.5*(tcor-1.)*(tanh((y(j)-ztr)/wtr)+1.)
      enddo

c   plasma beta

c   bottom of magnetic sheet
      zf1=4.0
c   top of magnetic sheet
      zf2=8.0

      rbetaf=1/1.
      do j=1,jx
          af1 = (  tanh((y(j)-zf1)/0.5) + 1.)/2
          af2 = ( -tanh((y(j)-zf2)/0.5) + 1.)/2
          rbeta0(j) = rbetaf*af1*af2
      enddo

c   density, pressure

      den0(jzero) = 1.
      pre0(jzero) = 1./gm * tem0(jzero)
      do j=jzero+1,jx
         den0(j) = den0(j-1)
     &            *((1+rbeta0(j-1))*tem0(j-1)+0.5*gm*gym0(j-1)*dym(j-1))
     &            /((1+rbeta0(j)  )*tem0(j)  -0.5*gm*gym0(j-1)*dym(j-1))
         pre0(j) = pre0(jzero)
     &        *(den0(j)/den0(jzero))*(tem0(j)/ tem0(jzero))
      enddo
      do j=jzero-1,1,-1
         den0(j) = den0(j+1)
     &             *((1+rbeta0(j+1))*tem0(j+1)-0.5*gm*gym0(j)*dym(j))
     &             /((1+rbeta0(j)  )*tem0(j)  +0.5*gm*gym0(j)*dym(j))
         pre0(j) = pre0(jzero)
     &        *(den0(j)/den0(jzero))*(tem0(j)/ tem0(jzero))
      enddo

c   magnetic field

      do j=1,jx
          bmx0(j)  = sqrt(8*pi*pre0(j)*rbeta0(j))
      enddo

c  set all of aboves 

c  magnetic direction angle( 0:only Parker mode, 90:only Interchange mode)
      angle=45.d0
      theta=angle*pi/180.d0

      do j=1,jx
      do i=1,ix
         ro(i,j)=den0(j)
         pr(i,j)=pre0(j)
         vx(i,j)=0.0
         vy(i,j)=0.0
         vz(i,j)=0.0
c  parker mode
         bx(i,j)=bmx0(j)
         bz(i,j)=0.0
c  mix mode( parker and interchange )
c         bx(i,j)=bmx0(j)*cos(theta)
c         bz(i,j)=bmx0(j)*sin(theta)
         by(i,j)=0.0
      enddo
      enddo

c----------------------------------------------------------------------|
c     write parameters to file
c----------------------------------------------------------------------|
c                                                     netcdf - start >>>
c     call ncputos(idf,'gm',gm)
c     call ncputos(idf,'ztr',ztr)
c     call ncputos(idf,'rbetaf',rbetaf)
c                                                     netcdf - end   <<<
      call dacputparamd(mf_params,'gm',gm)
      call dacputparamd(mf_params,'ztr',ztr)
      call dacputparamd(mf_params,'rbetaf',rbetaf)

      return
      end
