/*Code by John Steffen, Austin Todd, and Jiangyi Hu 
 *Fall2012
 
 * Edited 2014 September 8 by Jason Keefer
 * - warnings will follow when compiling, but will give desired output
 
 *
  To compile c function:
  example #1 MATLAB 
 *
 *mex mexmft12.c (within Matlab interface)

*/

#include "mft12.c"
#include "mex.h"

void mexFunction( int nlhs, mxArray *argv2[],
                 int nrhs, const mxArray *argv1[])

{
 
  int    dyn_in_prm;
  float  dyn_in_val, dyn_in_val2, CONVECT, conv_crit, pressure;
  int    air_moist_prm;
  float  air_moist_val;
  int    sfc_moist_prm;
  float  sfc_moist_val, salinity;
  int    ss_prm;
  float  ss_val, t_air; 
  int    sst_prm;
  float  t_skin, ref_ht_wind, ref_ht_tq;
  int    astab; 
  float  Qnet;
  int    warn;
  int    mft_flag;
  float  dom_phs_spd, h_sig, lhf, q_star, shf, tau[2], t_star,
         u_star[2], wave_age, ww_stab, z_over_L, zo_m;
  int    eqv_neut_prm;
  float  z_wanted, u_at_z, t_at_z, q_at_z;
  int    flux_model, z0_mom_prm, z0_TQ_prm, stable_prm;
  
  float  phi_u;
  
  double *inMatrix;                /*1x27 input matrix*/
  double *outMatrix;              /* output matrix*/


  /* create a pointer to the real data in the input matrix */ 
  inMatrix = mxGetPr(argv1[0]);
  
  /* create the output matrix */
  argv2[0] = mxCreateDoubleMatrix(1,15,mxREAL);
 
  /* get a pointer to the real data in the output matrix */
  outMatrix = mxGetPr(argv2[0]);
  
  dyn_in_prm     = *(inMatrix);
  dyn_in_val     = *(inMatrix+1);
  dyn_in_val2    = *(inMatrix+2); 
  CONVECT        = *(inMatrix+3);
  conv_crit      = *(inMatrix+4); 
  pressure       = *(inMatrix+5);
  air_moist_prm  = *(inMatrix+6);
  air_moist_val  = *(inMatrix+7); 
  sfc_moist_prm  = *(inMatrix+8);
  sfc_moist_val  = *(inMatrix+9);
  salinity       = *(inMatrix+10);
  ss_prm         = *(inMatrix+11);
  ss_val         = *(inMatrix+12);
  t_air          = *(inMatrix+13);
  sst_prm        = *(inMatrix+14);  
  t_skin         = *(inMatrix+15);
  ref_ht_wind    = *(inMatrix+16);
  ref_ht_tq      = *(inMatrix+17);
  z_wanted       = *(inMatrix+18); 
  astab          = *(inMatrix+19);  
  eqv_neut_prm   = *(inMatrix+20);
  Qnet           = *(inMatrix+21);
  warn           = *(inMatrix+22);
  flux_model     = *(inMatrix+23);
  z0_mom_prm     = *(inMatrix+24);
  z0_TQ_prm      = *(inMatrix+25);
  stable_prm     = *(inMatrix+26);

  /*  
  conv_crit      = 0.00005; 
  sst_prm        = 1;  
  z_wanted       = 10; 
  eqv_neut_prm   = 0;
  Qnet           = 300;
  warn           = 0;
  flux_model     = 9;
  z0_mom_prm     = 0;
  z0_TQ_prm      = 0;
  stable_prm     = 0;
  */
  
    

  
  
/*  
mexPrintf( "test1\n" );
  
*/   
/*  
mexPrintf( "Input Variables: %i %f %f %f %f %f %i %f %i %f %f %i %f %f %i %f %f %f %i %f %i %i %f %i\n", 
   dyn_in_prm, dyn_in_val, dyn_in_val2,
   CON_P, conv_crit, press_sfc, air_moist_prm, air_moist_val, 
   sfc_moist_prm, sfc_moist_val, salinity, ss_prm, ss_val, ta, 
   sst_prm, tskin, z_wind_ref, z_temp_ref, astab, Qnet, warn, eqv_neut, z_wanted, flux_model); 

*/

    
  mft_flag = ht_adj_( dyn_in_prm, dyn_in_val, dyn_in_val2, CONVECT, 
              conv_crit, pressure,
              air_moist_prm, air_moist_val, sfc_moist_prm, sfc_moist_val,
              salinity, ss_prm, ss_val, t_air, sst_prm, t_skin, ref_ht_wind,
              ref_ht_tq, z_wanted, astab, eqv_neut_prm, Qnet, warn, flux_model,
              z0_mom_prm, z0_TQ_prm, stable_prm,
              &shf, &lhf, &tau, &u_star, &t_star,
              &q_star, &z_over_L, &wave_age, &dom_phs_spd, &h_sig, &ww_stab, &zo_m, 
              &u_at_z, &t_at_z, &q_at_z);
  
  /*
  mexPrintf( "Output variables: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", 
          shf, lhf, tau, u_star, t_star, q_star, z_over_L, wave_age, dom_phs_spd, h_sig, 
          ww_stab, zo_m, u_at_z, t_at_z, q_at_z);


 
 mexPrintf( "test2\n" );
 */

  /* should be the magnitude of u_star */ 
  outMatrix[0] = u_star[0]; /* should be first component of u_star */
  outMatrix[1] = u_star[1]; /* should be the second component of u_star */
  outMatrix[2] = t_star;
  outMatrix[3] = q_star;
  outMatrix[4] = z_over_L;
  outMatrix[5] = dom_phs_spd;
  outMatrix[6] = wave_age;
  outMatrix[7] = h_sig;
  outMatrix[8] = tau[0]; /* should be the first component of tau */
  outMatrix[9] = tau[1]; /* should be the second component of tau */
  outMatrix[10] = shf;
  outMatrix[11] = lhf;
  outMatrix[12] = u_at_z;
  outMatrix[13] = t_at_z;
  outMatrix[14] = q_at_z; 

/*  outMatrix[10] = ww_stab;
  outMatrix[11] = zo_m; */
  
/*     
  return( mft_flag );
*/  
  

}