/***************************************************************************
                 module classe sbcgyre_utau_vtau_wndm.h  -  description
***************************************************************************/
// Mohamed Berrada
// locean-ipsl.upmc, Paris, Mars 18, 2009
//===========================================================================
//                          methode forward
forward (YREAL x1,YREAL x2,YREAL x3,YREAL x4)
{/*	1	from	sbcgyre_utau_vtau	1  i-1 	j         t
	2	from	sbcgyre_utau_vtau	2  i 	j-1       t
	3	from	sbcgyre_utau_vtau	1  i 	j         t
	4	from	sbcgyre_utau_vtau	2  i 	j         t*/
 /* verification norm(xY-xF),norm(xY-xF)/norm(xF)
   wndm   3.299904462350809e-09,     2.039381910331513e-11
*/
  if(Yt==TU){
    YS1=0.;
  }
  else{
    //      ! Estimation of wind speed as a function of wind stress ( |tau|=rhoa*Cd*|U|^2 )
    double zrhoa  = 1.22;//         ! Air density kg/m3
    double zcdrag = 1.5e-3;//       ! drag coefficient
    double zcoef = 0.5 / ( zrhoa * zcdrag ); 
    if(Yj==NY-1 || Yj==0 || Yi==NX-1 || Yi==0)
      YS1=0.;
    else{
      double  ztx = x1+x3;//utau(ji-1,jj  ) + utau(ji,jj) 
      double  zty = x2+x4;// vtau(ji  ,jj-1) + vtau(ji,jj) 
      double  ztau = sqrt( ztx * ztx + zty * zty );
      YS1=sqrt(ztau * zcoef ) * tmask(Yi,Yj,0);
    }
    // CALL lbc_lnk( wndm(:,:) , 'T', 1. ) pour l'instant "closed domain"
  }
  //
}

//===========================================================================
//                         methode  backward

backward (YREAL x1,YREAL x2,YREAL x3,YREAL x4)
{
  YJ1I1=0.;     YJ1I2=0.;     YJ1I3=0.;     YJ1I4=0.;
  if(Yt!=TU){
    //      ! Estimation of wind speed as a function of wind stress ( |tau|=rhoa*Cd*|U|^2 )
    double zrhoa  = 1.22;//         ! Air density kg/m3
    double zcdrag = 1.5e-3;//       ! drag coefficient
    double zcoef = 0.5 / ( zrhoa * zcdrag ); 
    if(Yj==NY-1 || Yj==0 || Yi==NX-1 || Yi==0){
      YJ1I1=0.;     YJ1I2=0.;     YJ1I3=0.;     YJ1I4=0.;
    }
    else{
      double  ztx = x1+x3,dztx_x1=1.,dztx_x3=1.;//utau(ji-1,jj  ) + utau(ji,jj) 
      double  zty = x2+x4,dzty_x2=1.,dzty_x4=1.;// vtau(ji  ,jj-1) + vtau(ji,jj) 
      double  ztau = sqrt( ztx * ztx + zty * zty ), 
	dztau_ztx  = ztx/ztau,
	dztau_zty  = zty/ztau,
	dztau_x1   = dztau_ztx*dztx_x1,
	dztau_x2   = dztau_zty*dzty_x2,
	dztau_x3   = dztau_ztx*dztx_x3,
	dztau_x4   = dztau_zty*dzty_x4;
      YJ1I1 = 0.5 / sqrt(ztau * zcoef ) * zcoef * dztau_x1 * tmask(Yi,Yj,0);
      YJ1I2 = 0.5 / sqrt(ztau * zcoef ) * zcoef * dztau_x2 * tmask(Yi,Yj,0);
      YJ1I3 = 0.5 / sqrt(ztau * zcoef ) * zcoef * dztau_x3 * tmask(Yi,Yj,0);
      YJ1I4 = 0.5 / sqrt(ztau * zcoef ) * zcoef * dztau_x4 * tmask(Yi,Yj,0);
    }
    
  }
}
//===========================================================================
//********************* FIN DU MODULE sbcgyre_utau_vtau_wndm **********************