C     gateway subroutine
      
      subroutine mexfunction(nlhs, plhs, nrhs, prhs)
      implicit none
C-----------------------------------------------------------------------
C     (integer) Replace integer by integer*8 on 64-bit platforms
      integer*8 plhs(*), prhs(*)
      integer*8 i_ut, i_alt, i_theta, i_phi, i_dst, i_f107
      integer*8 i_dnhmf, i_dnlmf, i_dcord, i_dpred, i_dcurr
      integer*8 i_bmdl_x, i_bmdl_y, i_bmdl_z 
      integer*8 i_jmdl_x, i_jmdl_y
      integer*8 mxGetPr, MXCREATEDOUBLEMATRIX, mxIsChar, mxGetString
C-----------------------------------------------------------------------
C
      integer*8 nlhs, nrhs, mxGetM, mxGetN
      integer*8 m_in, n_in, N_data
      character*80 path
      integer*8 status, strlen
         
      open(-1, file='STD_OUT', status='unknown')

      m_in = mxGetM(prhs(1))
      n_in = mxGetN(prhs(1))
      N_data = m_in * n_in

c      print*,m_in    

      i_ut = mxGetPr(prhs(1))
      i_alt = mxGetPr(prhs(2))
      i_theta = mxGetPr(prhs(3))
      i_phi = mxGetPr(prhs(4))
      i_dst = mxGetPr(prhs(5))
      i_f107 = mxGetPr(prhs(6))

C     The input  must be a string.
      if(mxIsChar(prhs(7)) .ne. 1) then
         call mexErrMsgTxt('Input must be a string.')
C     The input must be a row vector.
      elseif (mxGetM(prhs(7)) .ne. 1) then
         call mexErrMsgTxt('Input must be a row vector.')
      endif 
C     Get the length of the input string.
      strlen = mxGetM(prhs(7))*mxGetN(prhs(7))
C     Get the string contents (dereference the input integer).
      status = mxGetString(prhs(7), path, 80)
C     Check if mxGetString is successful.
      if (status .ne. 0) then 
         call mexErrMsgTxt('Filename length must be less than 80')
      endif
      
      i_dnhmf = mxGetPr(prhs(8))
      i_dnlmf = mxGetPr(prhs(9))
      i_dcord = mxGetPr(prhs(10))
      i_dpred = mxGetPr(prhs(11))
      i_dcurr = mxGetPr(prhs(12))
      
      
      plhs(1) = MXCREATEDOUBLEMATRIX(m_in, 7, 0)
      i_bmdl_x = mxGetPr(plhs(1))
      plhs(2) = MXCREATEDOUBLEMATRIX(m_in, 7, 0)
      i_bmdl_y = mxGetPr(plhs(2))
      plhs(3) = MXCREATEDOUBLEMATRIX(m_in, 7, 0)
      i_bmdl_z = mxGetPr(plhs(3))
      plhs(4) = MXCREATEDOUBLEMATRIX(m_in, 4, 0)
      i_jmdl_x = mxGetPr(plhs(4))
      plhs(5) = MXCREATEDOUBLEMATRIX(m_in, 4, 0)
      i_jmdl_y = mxGetPr(plhs(5))

      
C     call the computational routine
      call make_cm4(%val(i_bmdl_x), %val(i_bmdl_y), %val(i_bmdl_z), 
     *     %val(i_jmdl_x), %val(i_jmdl_y), 
     *     %val(i_ut), %val(i_alt), %val(i_theta), 
     *     %val(i_phi), %val(i_dst), %val(i_f107), 
     *     N_data, path, 
     *     %val(i_dnhmf), %val(i_dnlmf), %val(i_dcord), %val(i_dpred), 
     *     %val(i_dcurr))
      
      close(-1)
      return
      end
      
C     computational subroutine
      subroutine make_cm4(bmdl_x, bmdl_y, bmdl_z, 
     *     jmdl_x, jmdl_y, 
     *     ut, alt, theta, 
     *     phi, dst, f107, 
     *     N_data, path, 
     *     dnhmf, dnlmf, dcord, dpred, dcurr)
      implicit real*8 (a-h,o-z)
      real*8  ut(*), alt(*), theta(*), phi(*)
      real*8 dst(*), f107(*)
      real*8 bmdl_x(N_data, 7), bmdl_y(N_data, 7), bmdl_z(N_data, 7)
      real*8 jmdl_x(N_data, 4), jmdl_y(N_data, 4) 
      real*8 bmdl(3,7), jmdl(3,4)
      
      character*80 path(1)
      
c     begin CM4 declarations ------------------------------------------------------
      integer*8 unit(3)
      logical load(3),indx(2),gmut,cord,pred(6),curr,coef
      integer*8 nhmf(2),nlmf(2),perr,oerr,cerr
      real*8 dnhmf(2), dnlmf(2), dpred(6)
      real*8 mut      
      real*8 gmdl(1)      ! note, an array of length one      
c     end CM4 declarations --------------------------------------------------------
      
c     begin CM4 assignments -------------------------------------------------------      
      unit(1)=11                ! logical unit for model coefficients
      load(1)=.true.            ! set to initially read model coefficients
      indx(1)=.false.           ! get Dst from argument
      indx(2)=.false.           ! get solar flux from argument
      gmut   =.true.            ! compute MUT from UT
      coef   =.false.           ! no coefficients generated
      perr   = 1                ! error print flag
      oerr   = 50               ! error logical unit
      cerr   = 0                ! error flag
c     end CM4 assignments ---------------------------------------------------------
      
      open(oerr, file='cm4.msg', status = 'unknown')
      write(oerr, *)
      
      do i=1,6
         if(dpred(i).ge.0.5) then
            pred(i) = .true.
         else
            pred(i) = .false.
         endif
      enddo
      
      if(dcord.ge. 0.5) then
		cord = .true.
	else
		cord = .false.
	endif

      if(dcurr.ge. 0.5) then
		curr = .true. ! current density J fields computed
	else
		curr = .false.
	endif

      nhmf(1) = dnhmf(1)
      nhmf(2) = dnhmf(2)
      nlmf(1) = dnlmf(1)
      nlmf(2) = dnlmf(2)

      
      do i=1,N_data
c	print*,'cm4_m.f:', ut(i)
         call CM4FIELD(PATH,UNIT,LOAD,INDX,GMUT,CORD,PRED,CURR,
     1        COEF,NHMF,NLMF,UT(i),MUT,THETA(i),PHI(i),ALT(i),DST(i),
     2        F107(i),BMDL,JMDL,GMDL,PERR,OERR,CERR)
         do is = 1,7
            bmdl_x(i,is) = bmdl(1,is)
            bmdl_y(i,is) = bmdl(2,is)
            bmdl_z(i,is) = bmdl(3,is)
         enddo
c sheet current density
         do is = 1,4
            jmdl_x(i,is) = jmdl(1,is)
            jmdl_y(i,is) = jmdl(2,is)
         enddo
      enddo
 999  close(oerr)
      end
      
      include 'cm4field.f'