c mex subroutine to calculate fluxes and adjust met variables 
c using matlab data in the fortran program bulk_mat_v25b.f
c
c  To compile:
c  
c     fmex bulk_mat_v25b.f bulk_mat_v25b_function.f 
c  Output:
c     bulk_mat_v25b.mex4  or bulk_mat_v25b.solmex  **
c                ** exact form depends on compiler
c
c  To run, in matlab, type:
c
c [QH,QE,RF,TAU,TAUr,QE_webb,ws_h,Ta_h,qq_h,T0, ...
c      tstar,qstar,ustar,Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter,wg]= ...
c       bulk_mat_v25b(hUm,hTm,hUs,hTs,ts_depth,ws,sst,atb,rh,pr,Zi,...
c                  rs,rl,rain,jcool,jwarm,time,lat,lon)
c
c  or use bulk_v25.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
      integer dt_wrm,dter,Cd,Ch,Ce,wg,TAUr,QE_webb,lat,lon
c
      m=mxGetM(prhs(1))
      n=mxGetN(prhs(1))
c
c create matrices for the return arguments
      if (nlhs .eq. 26) 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)
        plhs(24)=mxCreateFull(m,n,0)
        plhs(25)=mxCreateFull(m,n,0)
        plhs(26)=mxCreateFull(m,n,0)
      elseif (nlhs .eq. 10) 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)
      else
        call mexErrMsgTxt('either 26 or 10 output arguments required')
      endif
c
c reference arguments to get array pointers
c   1  2 3   4    5    6      7    8    9   10  11    12    13
c [QH,QE,RF,TAU,TAUr,QE_webb,ws_h,Ta_h,qq_h,T0,tstar,qstar,ustar, ...
c       14 15 16 17 18 19 20 21 22  23   24    25   26
c       Cd,Ch,Ce,RR,RT,RQ,Zl,Zo,Zot,Zoq,dt_wrm,dter,wg]= ...
c        bulk_mat_v25b(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,lat,lon)
c                       14    15   16   17   18  19
c
c Set pointers to the output variables
      if (nlhs .eq. 26) then
        QH=mxGetPr(plhs(1))
        QE=mxGetPr(plhs(2))
        RF=mxGetPr(plhs(3))
        TAU=mxGetPr(plhs(4))
        TAUr=mxGetPr(plhs(5))
        QE_webb=mxGetPr(plhs(6))
        ws_h=mxGetPr(plhs(7))
        Ta_h=mxGetPr(plhs(8))
        qq_h=mxGetPr(plhs(9))
        T0=mxGetPr(plhs(10))
        tstar=mxGetPr(plhs(11))
        qstar=mxGetPr(plhs(12))
        ustar=mxGetPr(plhs(13))
        Cd=mxGetPr(plhs(14))
        Ch=mxGetPr(plhs(15))
        Ce=mxGetPr(plhs(16))
        RR=mxGetPr(plhs(17))
        RT=mxGetPr(plhs(18))
        RQ=mxGetPr(plhs(19))
        Zl=mxGetPr(plhs(20))
        Zo=mxGetPr(plhs(21))
        Zot=mxGetPr(plhs(22))
        Zoq=mxGetPr(plhs(23))
        dt_wrm=mxGetPr(plhs(24))
        dter=mxGetPr(plhs(25))
        wg=mxGetPr(plhs(26))
      else
        QH=mxGetPr(plhs(1))
        QE=mxGetPr(plhs(2))
        RF=mxGetPr(plhs(3))
        TAU=mxGetPr(plhs(4))
        TAUr=mxGetPr(plhs(5))
        QE_webb=mxGetPr(plhs(6))        
        ws_h=mxGetPr(plhs(7))
        Ta_h=mxGetPr(plhs(8))
        qq_h=mxGetPr(plhs(9))
        T0=mxGetPr(plhs(10))
      endif
c
c Set pointers to the input variables
      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))
      lat=mxGetPr(prhs(18))
      lon=mxGetPr(prhs(19))
c 
c  MAKE SURE THAT BULK_V25 is called correctly. This is from bulk_mat_v2_5b.f:
c      subroutine bulk_mat_v25b(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,lat,lon
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,wg,TAUr,QE_webb)
c
      call bulk_mat_v25b(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(lat),%VAL(lon),
     & %VAL(QH),%VAL(QE),
     & %VAL(RF),%VAL(TAU),%VAL(ustar),%VAL(tstar),%VAL(qstar),
     & %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),%VAL(wg),%VAL(TAUr),%VAL(QE_webb))
c
      return
      end