c mex subroutine to calculate fluxes and adjust met variables 
c using matlab data in the fortran program bulk_v2a.f
c
c  To compile:
c  
c     fmex bulk_v2a.f bulk_v2ag.f 
c
c
c
c  To run, in matlab, type:
c
c [QH,QE,RF,TAU,ws_h,Ta_h,qq_h,T0, ...
c      tstar,qstar,ustar,Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter]= ...
c       bulk_v2a(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,Zi,...
c                  rs,rl,rain,jcool,jwarm,time)
c
c  or use bulk_v2a_drv.m
c
      subroutine mexFunction(nlhs,plhs,nrhs,prhs)
      integer plhs(*),prhs(*)
      integer nlhs, nrhs
      integer mxGetM, mxGetN, mxGetPr, mxCreateFull
c
      integer hUm,hTm,hUs,hTs,ws_h,Ta_h,qq_h,T0
      integer ws,sst,atb,rh,QH,QE,TAU,RF,m,n,pr,Zi
      integer rs,rl,rain,jcool,jwarm,ts_depth,time
      integer ustar,tstar,qstar,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter
      integer Cd,Ch,Ce
c
c check for proper number of input arguments
c      if (nrhs .ne. 16) then
c         call mexErrMsgTxt('Sixteen input arguments required')
c      elseif (nlhs .ne. 22) then
c         call mexErrMsgTxt('Twenty-two output arguments required')
c      endif
c
      m=mxGetM(prhs(1))
      n=mxGetN(prhs(1))
c
c create matrices for the return arguments
      if (nlhs .eq. 23) then
      plhs(1)=mxCreateFull(m,n,0)
      plhs(2)=mxCreateFull(m,n,0)
      plhs(3)=mxCreateFull(m,n,0)
      plhs(4)=mxCreateFull(m,n,0)
      plhs(5)=mxCreateFull(m,n,0)
      plhs(6)=mxCreateFull(m,n,0)
      plhs(7)=mxCreateFull(m,n,0)
      plhs(8)=mxCreateFull(m,n,0)
      plhs(9)=mxCreateFull(m,n,0)
      plhs(10)=mxCreateFull(m,n,0)
      plhs(11)=mxCreateFull(m,n,0)
      plhs(12)=mxCreateFull(m,n,0)
      plhs(13)=mxCreateFull(m,n,0)
      plhs(14)=mxCreateFull(m,n,0)
      plhs(15)=mxCreateFull(m,n,0)
      plhs(16)=mxCreateFull(m,n,0)
      plhs(17)=mxCreateFull(m,n,0)
      plhs(18)=mxCreateFull(m,n,0)
      plhs(19)=mxCreateFull(m,n,0)
      plhs(20)=mxCreateFull(m,n,0)
      plhs(21)=mxCreateFull(m,n,0)
      plhs(22)=mxCreateFull(m,n,0)
      plhs(23)=mxCreateFull(m,n,0)
      elseif (nlhs .eq. 8) then
      plhs(1)=mxCreateFull(m,n,0)
      plhs(2)=mxCreateFull(m,n,0)
      plhs(3)=mxCreateFull(m,n,0)
      plhs(4)=mxCreateFull(m,n,0)
      plhs(5)=mxCreateFull(m,n,0)
      plhs(6)=mxCreateFull(m,n,0)
      plhs(7)=mxCreateFull(m,n,0)
      plhs(8)=mxCreateFull(m,n,0)
      else
      call mexErrMsgTxt('either 23 or 8 output arguments required')
      endif
c
c Dereference arguments to get array pointers
c   1  2 3   4    5    6   7    8    9    10  11 
c [QH,QE,RF,TAU,ws_h,Ta_h,qq_h,T0,tstar,qstar,ustar, ...
c       12 13 14 15 16 17 18 19 20  21   22   23
c       Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter]= ...
c        bulk_v2a(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,Zi,rs,rl,
c                  1   2   3   4    5      6   7  8  9  10 11 12  13  
c                      rain,jcool,jwarm,time)
c                       14    15   16   17
c
      if (nlhs .eq. 23) then
      QH=mxGetPr(plhs(1))
      QE=mxGetPr(plhs(2))
      RF=mxGetPr(plhs(3))
      TAU=mxGetPr(plhs(4))
      ws_h=mxGetPr(plhs(5))
      Ta_h=mxGetPr(plhs(6))
      qq_h=mxGetPr(plhs(7))
      T0=mxGetPr(plhs(8))
      tstar=mxGetPr(plhs(9))
      qstar=mxGetPr(plhs(10))
      ustar=mxGetPr(plhs(11))
      Cd=mxGetPr(plhs(12))
      Ch=mxGetPr(plhs(13))
      Ce=mxGetPr(plhs(14))
      RR=mxGetPr(plhs(15))
      RT=mxGetPr(plhs(16))
      RQ=mxGetPr(plhs(17))
      Zl=mxGetPr(plhs(18))
      Zo=mxGetPr(plhs(19))
      Zot=mxGetPr(plhs(20))
      Zoq=mxGetPr(plhs(21))
      dt_wrm=mxGetPr(plhs(22))
      dter=mxGetPr(plhs(23))
      else
      QH=mxGetPr(plhs(1))
      QE=mxGetPr(plhs(2))
      RF=mxGetPr(plhs(3))
      TAU=mxGetPr(plhs(4))
      ws_h=mxGetPr(plhs(5))
      Ta_h=mxGetPr(plhs(6))
      qq_h=mxGetPr(plhs(7))
      T0=mxGetPr(plhs(8))
      endif
      hUm=mxGetPr(prhs(1))
      hTm=mxGetPr(prhs(2))
      hUs=mxGetPr(prhs(3))
      hTs=mxGetPr(prhs(4))
      ts_depth=mxGetPr(prhs(5))
      ws=mxGetPr(prhs(6))
      sst=mxGetPr(prhs(7))
      atb=mxGetPr(prhs(8))
      rh=mxGetPr(prhs(9))
      pr=mxGetPr(prhs(10))
      Zi=mxGetPr(prhs(11))
      rs=mxGetPr(prhs(12))
      rl=mxGetPr(prhs(13))
      rain=mxGetPr(prhs(14))
      jcool=mxGetPr(prhs(15))
      jwarm=mxGetPr(prhs(16))
      time=mxGetPr(prhs(17))
c 
c  MAKE SURE THAT CALL BULK_V2A correctly. This is from bulk_v2a.f:
c      subroutine bulk_v2a(m,time,
c     & hUm,hTm,hUs,hTs,ts_depth,              !    MC
c     & ws,sst,atb,qq,ws_h,Ta_h,qq_h,          !    MC
c     & rs,rl,rain,pp,zi,Jcool,jwarm,
c     & QH,QE,RF,TAU,
c     & Ustar,Tstar,Qstar,
c     & CD,CH,CE,
c     & RR,RT,RQ,ZL,ZO,zot,zoq,dt_wrm,dter,T0)
 
c

      call bulk_v2a(m,%VAL(time),%VAL(hUm),%VAL(hTm),%VAL(hUs),
     & %VAL(hTs),%VAL(ts_depth),%VAL(ws),%VAL(sst),%VAL(atb),%VAL(rh),
     & %VAL(ws_h),%VAL(Ta_h),%VAL(qq_h),%VAL(rs),%VAL(rl),%VAL(rain),
     & %VAL(pr),%VAL(Zi),%VAL(jcool),%VAL(jwarm),%VAL(QH),%VAL(QE),
     & %VAL(RF),%VAL(TAU),%VAL(ustar),%VAL(qstar),%VAL(tstar),
     & %VAL(Cd),%VAL(Ch),%VAL(Ce),%VAL(RR),%VAL(RT),%VAL(RQ),
     & %VAL(Zl),%VAL(Zo),%VAL(Zot),%VAL(Zoq),%VAL(dt_wrm),
     & %VAL(dter),%VAL(T0))
c
      return
      end