/*====================================================================================
! Purpose(s): 
!
! Contain(s):  c2ndstep
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!====================================================================================*/

#include "CArrayMasks.h"
#include "CVars.h"
#include "CMisc.h"
#include "CErrors.h"
#include <stdio.h>

#ifdef __INTEL_COMPILER
#define PYM c2ndstep
#else 
#define PYM c2ndstep_
#endif

//extern "C" void _FORTRAN c2ndstep_ (_DOUBLE_ARRAY AAR,     _DOUBLE_ARRAY AAT,
extern "C" void _FORTRAN PYM (_DOUBLE_ARRAY AAR,     _DOUBLE_ARRAY AAT,    
                                    _DOUBLE_ARRAY AXCoefs, _DOUBLE_ARRAY AYCoefs, 
                                    _DOUBLE_ARRAY AXi,     _DOUBLE_ARRAY AYj, 
                                    _DOUBLE_ARRAY ADelX,   _DOUBLE_ARRAY ADelY,
                                    _DOUBLE_ARRAY AU,      _DOUBLE_ARRAY AV, 
                                    _DOUBLE_ARRAY AF, 
                                    _DOUBLE_ARRAY AP, 
                                    _INT_ARRAY ANf, 
                                    _INT_ARRAY Axnpc,             
                                    _INT APorosityModel, 
                                    _DOUBLE_ARRAY AXPorosity, _DOUBLE_ARRAY AGc, 
                                    _DOUBLE ADelT, _DOUBLE Afrsurf, _DOUBLE ARhoF, 
                                    _INT AKb, _INT AKt, _INT AKl, _INT AKr, _INT AIMax, _INT AJMax)
{
/*!=======================================================================
! Purpose(s): 
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!=======================================================================*/
  
  int PIdxL, PIdxR, PIdxT, PIdxB;

  // Creates the mesh object.
  TMesh Mesh(AIMax, AJMax, AXCoefs, AYCoefs, AXi, AYj, ADelX, ADelY);

  // Creates the masks.
  TRMask      wAR        (Mesh, AAR);
  TTMask      wAT        (Mesh, AAT);
  TCMask      wF         (Mesh, AF);
  TCMask      wP         (Mesh, AP);
  TRMask      wU         (Mesh, AU);
  TTMask      wV         (Mesh, AV);
  
  TIntMask    wNf        (Mesh, ANf);
  TIntMask    wXPorIdx   (Mesh, Axnpc);

  Mesh.DefineArea(1, AIMax, 1, AJMax);

  Mesh.First2();
  do
  {
    if (Mesh.IsInDefinedArea())
    // If the cell is inside the area of calculus...
    {

      // Gets the porosity index at every side of the cell.
      IdxUnpack(wXPorIdx.c(), PIdxL, PIdxR, PIdxT, PIdxB);

      double rhoR = ARhoF*wF.IR();
      double rhoT = ARhoF*wF.IT();

      // UPDATES U COMPONENT.
      if (rhoR>=Afrsurf && wAR.c()>=OBSTACLE) 
      // If it is a valid cell...
      {

        double DelU = ADelT*R_Dx(wP)/rhoR;

        if (PIdxR>0) 
        // If right side has porous material...
        {
          if (APorosityModel==pmNonlinearFriction) DelU = DelU/AGc[PIdxR];
          else if (APorosityModel==pmLinearFriction) printf("Non linear frictional porosity not implemented.", __FILE__, __LINE__);
        }
        
        *wU.C=wU.c()-DelU;

      }
      else
      // If it is not a valid cell...
      {
        *wU.C=0.0;
      }

      // UPDATES V COMPONENT.
      if (Mesh.I!=1 || AKl!=bcSpecInflowOutFlow)
      {
        if (rhoT>=Afrsurf && wAT.c()>=OBSTACLE) 
        // If it is a valid cell...
        {

          double DelV = ADelT*T_Dy(wP)/rhoT;

          if (PIdxT>0)
          // If top side has porous material...
          {
            if (APorosityModel==pmNonlinearFriction) DelV = DelV/AGc[PIdxT];
            else if (APorosityModel==pmLinearFriction) printf("Non linear frictional porosity not implemented.", __FILE__, __LINE__);
          }
          
          if (Mesh.I>1) *wV.C=wV.c()-DelV;

        }
        else
        // If it is not a valid cell...
        {
          *wV.C=0.0;
        }

      } //if


    } //if


  } while(Mesh.Next2());



}