/***************************************************************************
                 module classe solsor_dynspg_flt.h  -  description
***************************************************************************/
// Mohamed Berrada
// locean-ipsl.upmc, Paris, April 23, 2009
//===========================================================================
//                          methode forward
forward ()
{
  if(Yt==TU)
    return;
  else{
    if(Yi==0){
      int t,jn,jj,ji, icount = 0,ncut=0,niter=0,ishift;
      niter=0;ncut=0;
      if(Yt==TU+1 && neuler==0)
	t=0;
      else
	t=1;
      //---------------------
      //initialisation de gcx
      for(jj=0;jj<NY;jj++)
	for(ji=0;ji<NX;ji++){
	  //	gcxb(ji,jj)=gcx(ji,jj);
	  gcx(ji,jj)=YS_gcx_dynspg_flt(0,ji,jj,Yt);
	}
      //if(Yt==6) printf("------------------gcx(ji,jj,0,t)  =%lf\n",gcx(28,2) );
      //---------------------
      //boundarie conditions -- ceci doit etre a l'interieur de loop nmax mais cas special  
      for(jj=0;jj<NY;jj++){// ! East-West boundaries
	gcx(0,jj)=0.;gcx(NX-1,jj)=0.;
      }
      for(ji=0;ji<NX;ji++){// ! North-South boundaries
	gcx(ji,0)=0.;gcx(ji,NY-1)=0.;
      }
      
      double ztmp,zres,zres2,res=0.;
      //---------------------
      //boucle 
      //xsolmat_init();
      for(jn=1;jn<=nmax;jn++){
	// cas special        ! applied the lateral boundary conditions
	// CALL lbc_lnk_e( gcx, c_solver_pt, 1. ) ;  //c_solver_pt='T'    
	//         ! Residus    
	//         ! Guess black update
	
	for(jj=1;jj<NY-1;jj++){
	  ishift = jj%2;
	  for(ji=1+ishift;ji<NX-1;ji+=2){
	    ztmp = YS_gcb_dynspg_flt(0,ji,jj,Yt) 
	      - gcp(ji,jj,0,t) * gcx(ji  ,jj-1)   
	      - gcp(ji,jj,1,t) * gcx(ji-1,jj  )   
	      - gcp(ji,jj,2,t) * gcx(ji+1,jj  )   
	      - gcp(ji,jj,3,t) * gcx(ji  ,jj+1);
	    //     ! Estimate of the residual
	    zres = ztmp - gcx(ji,jj);
	    gcr(ji,jj) = zres * gcdmat(ji,jj,t) * zres;
	    //               ! Guess update
	    gcx(ji,jj) = sor * ztmp + (1.-sor) * gcx(ji,jj);
	  }
	}
	icount = icount + 1; 
	
	//cas special     ! applied the lateral boundary conditions
	//   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )  ;
	
	//   ! Guess red update
	for(jj=1;jj<NY-1;jj++){
	  ishift = (jj-1)%2;
	  for(ji=1+ishift;ji<NX-1;ji+=2){
	    ztmp = YS_gcb_dynspg_flt(0,ji,jj,Yt) 
	      - gcp(ji,jj,0,t) * gcx(ji  ,jj-1)   
	      - gcp(ji,jj,1,t) * gcx(ji-1,jj  )   
	      - gcp(ji,jj,2,t) * gcx(ji+1,jj  )   
	      - gcp(ji,jj,3,t) * gcx(ji  ,jj+1);
	    //     ! Estimate of the residual
	    zres = ztmp - gcx(ji,jj);
	    gcr(ji,jj) = zres * gcdmat(ji,jj,t) * zres;
	    //               ! Guess update
	    gcx(ji,jj) = sor * ztmp + (1.-sor) * gcx(ji,jj);
	  }
	}
	//	if(jn>=nmax-1) xtest("gcx__",gcx(28,3),28,3,0,Yt);
	icount = icount + 1;
	//     ! test of convergence
	if( jn > nmin && (jn-nmin)%nmod == 0 ) {
	  zres2=0.;
	  for(jj=1;jj<NY-1;jj++)
	    for(ji=1;ji<NX-1;ji++)
	      if(zres2<gcr(ji,jj)) zres2=gcr(ji,jj);
	  //      zres2 = MAXVAL( gcr(2:nlci-1,2:nlcj-1) );
	  //               ! test of convergence
	  if( zres2 < resmax || jn == nmax ){
	    res = sqrt( zres2 );
	    niter = jn; tniter[Yt]=jn;
	    ncut = 999;
	    //	    printf("FOR=======tniter(%d)=%d\n",Yt,tniter[Yt]);
	  }
	}
	//	if(jn>=nmax-1) xtest("gcx____",gcx(28,3),28,3,0,Yt);
	
	//         ! indicator of non-convergence or explosion
	//         IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
	if( ncut == 999 ) break;
      }
      //    cout<<"niter="<<niter<<endl;
      //    printf("res=%24.16e\n",res);
      // CALL lbc_lnk_e( gcx, c_solver_pt, 1. );//    ! boundary conditions
      /* FILE* fp;
	 fp=fopen("testh.pof","w");
	 for(int i=0;i<NX;i++)
	 for(int j=0;j<NY;j++){
	 fprintf(fp,"%4d  %4d %4d  %24.16e\n",Yt+1,i,j,gcr(i,j));
	 }
	 fclose(fp);*/
      /*      xtest("gcx",gcx(28,3),28,3,0,Yt);
	      xtest("gcb",YS_gcb_dynspg_flt(0,28,3,Yt),28,3,0,Yt); */
    }
  }
  //
}

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

backward ()
{

  if(Ycurward==BACKWARD){
    if(Yt==TU)
      return;
    else{
      if(Yi==0){
	double G_ztmp;
	int t,jn,jj,ji, ishift,istart;
	if(Yt==TU+1 && neuler==0)
	  t=0;
	else
	  t=1;
	//---------------------
	for(ji=0;ji<NX;ji++){
	  for(jj=0;jj<NY;jj++){
	    YG_gcx_dynspg_flt(0,ji,jj,Yt)=0.;
	    YG_gcb_dynspg_flt(0,ji,jj,Yt)=0.;
	  }
	}
	for(jj=0;jj<NY;jj++){
	  G_gcx(0,jj)=0.;G_gcx(NX-1,jj)=0.;
	  G_gcx0(0,jj)=0.;G_gcx0(NX-1,jj)=0.;
	}
	for(ji=0;ji<NX;ji++){
	  G_gcx(ji,0)=0.;G_gcx(ji,NY-1)=0.;
	  G_gcx0(ji,0)=0.;G_gcx0(ji,NY-1)=0.;
	}

	for(ji=1;ji<NX-1;ji++)
	  for(jj=1;jj<NY-1;jj++){
	    //	    YG_gcx2(0,ji,jj,Yt)=Ytb_gcx(ji,jj);
	    G_gcx(ji,jj)=YG_gcx2(0,ji,jj,Yt);
	    G_gcx0(ji,jj)=0.;
	  }

	for(jn=tniter[Yt];jn>=1;jn--){
	  //	printf("BACK1=======YG_gcx_dynspg_flt(%d)=%23.16e\n",Yt,G_gcx(9,9));
	  //	printf("BACK2=======YG_gcb_dynspg_flt(%d)=%23.16e\n",Yt,YG_gcb_dynspg_flt(0,9,9,Yt));

	  //  BACK   ! Guess red update
	  for(jj=NY-2;jj>=1;jj--){
	    ishift = (jj-1)%2;
	    istart=NX-2-(ishift-1)%2*(NX%2-1)%2-ishift*NX%2;
	    for(ji=istart;ji>=1;ji-=2){
	      G_ztmp=G_gcx(ji,jj)*sor;
	      G_gcx0(ji,jj)=G_gcx(ji,jj)*(1.-sor);
	      YG_gcb_dynspg_flt(0,ji,jj,Yt)+=G_ztmp;
	      G_gcx0(ji  ,jj-1)+=- G_ztmp*gcp(ji,jj,0,t) ;
	      G_gcx0(ji-1,jj  )+=- G_ztmp*gcp(ji,jj,1,t) ;
	      G_gcx0(ji+1,jj  )+=- G_ztmp*gcp(ji,jj,2,t) ;
	      G_gcx0(ji  ,jj+1)+=- G_ztmp*gcp(ji,jj,3,t) ;
	    }
	  }

	  //   BACK      ! Guess black update
	  
	  for(jj=NY-2;jj>=1;jj--){
	    ishift = jj%2;
	    istart=NX-2-(ishift-1)%2*(NX%2-1)%2-ishift*NX%2;
	    for(ji=istart;ji>=1;ji-=2){
	      G_ztmp=(G_gcx(ji,jj)+G_gcx0(ji,jj))*sor;
	      G_gcx0(ji,jj)=(G_gcx(ji,jj)+G_gcx0(ji,jj))*(1.-sor);
	      YG_gcb_dynspg_flt(0,ji,jj,Yt)+=G_ztmp;
	      G_gcx0(ji  ,jj-1)+=- G_ztmp*gcp(ji,jj,0,t) ;
	      G_gcx0(ji-1,jj  )+=- G_ztmp*gcp(ji,jj,1,t) ;
	      G_gcx0(ji+1,jj  )+=- G_ztmp*gcp(ji,jj,2,t) ;
	      G_gcx0(ji  ,jj+1)+=-G_ztmp* gcp(ji,jj,3,t) ;
	    }
	  }
	for(jj=0;jj<NY;jj++)
	  for(ji=0;ji<NX;ji++){
	    G_gcx(ji,jj)=G_gcx0(ji,jj);
	    G_gcx0(ji,jj)=0.;
	  }
  

	}
	for(jj=1;jj<NY-1;jj++)
	  for(ji=1;ji<NX-1;ji++){
	    YG_gcx_dynspg_flt(0,ji,jj,Yt)+=G_gcx(ji,jj);
	  }
	

	
      }
    }
    
  }
  else if(Ycurward==LINWARD){
    if(Yt==TU)
      return;
    else{
      if(Yi==0){
	int t,jn,jj,ji, ishift;
	if(Yt==TU+1 && neuler==0)
	  t=0;
	else
	  t=1;
	//---------------------
	//initialisation de gcx
	for(jj=0;jj<NY;jj++)
	  for(ji=0;ji<NX;ji++){
	    Ytb_gcx(ji,jj)=YG_gcx_dynspg_flt(0,ji,jj,Yt);
	  }
	//---------------------
	//boundarie conditions -- ceci doit etre a l'interieur de loop nmax mais cas special  
	for(jj=0;jj<NY;jj++){// ! East-West boundaries
	  Ytb_gcx(0,jj)=0.;Ytb_gcx(NX-1,jj)=0.;
	}
	for(ji=0;ji<NX;ji++){// ! North-South boundaries
	  Ytb_gcx(ji,0)=0.;Ytb_gcx(ji,NY-1)=0.;
	}
	
	double Ytb_ztmp;
	//---------------------
	//boucle 
	//xsolmat_init();
	//	printf("LIN=======tniter(%d)=%d\n",Yt,tniter[Yt]);
	for(jn=1;jn<=tniter[Yt];jn++){
	  // cas special        ! applied the lateral boundary conditions
	  // CALL lbc_lnk_e( gcx, c_solver_pt, 1. ) ;  //c_solver_pt='T'    
	  //         ! Residus    
	  //         ! Guess black update
	  
	  for(jj=1;jj<NY-1;jj++){
	    ishift = jj%2;//MOD( jj, 2 );
	    for(ji=1+ishift;ji<NX-1;ji+=2){
	      Ytb_ztmp = YG_gcb_dynspg_flt(0,ji,jj,Yt) 
		- gcp(ji,jj,0,t) * Ytb_gcx(ji  ,jj-1)   
		- gcp(ji,jj,1,t) * Ytb_gcx(ji-1,jj  )   
		- gcp(ji,jj,2,t) * Ytb_gcx(ji+1,jj  )   
		- gcp(ji,jj,3,t) * Ytb_gcx(ji  ,jj+1);
	      //               ! Guess update
	      Ytb_gcx(ji,jj) = sor * Ytb_ztmp + (1.-sor) * Ytb_gcx(ji,jj);
	    }
	  }
	  
	  //cas special     ! applied the lateral boundary conditions
	  //   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )  ;
	  
	  //   ! Guess red update
	  for(jj=1;jj<NY-1;jj++){
	    ishift = (jj-1)%2;
	    for(ji=1+ishift;ji<NX-1;ji+=2){
	      Ytb_ztmp = YG_gcb_dynspg_flt(0,ji,jj,Yt) 
		- gcp(ji,jj,0,t) * Ytb_gcx(ji  ,jj-1)   
		- gcp(ji,jj,1,t) * Ytb_gcx(ji-1,jj  )   
		- gcp(ji,jj,2,t) * Ytb_gcx(ji+1,jj  )   
		- gcp(ji,jj,3,t) * Ytb_gcx(ji  ,jj+1);
	      //               ! Guess update
	      Ytb_gcx(ji,jj) = sor * Ytb_ztmp + (1.-sor) * Ytb_gcx(ji,jj);
	    }
	  }
	  
	}
	for(ji=1;ji<NX-1;ji++)
	  for(jj=1;jj<NY-1;jj++)
	    YG_gcx2(0,ji,jj,Yt)=Ytb_gcx(ji,jj);
      }
    }
    //
    return;
  }
  //
}

//===========================================================================
//**********************  FIN DU MODULE solsor_dynspg_flt   *****************