scalar scalar_lambda=0.1; scalar beta=0.0065;
cout<<"lambda= "; cin>>scalar_lambda;  cout<<"beta= "; cin>>beta;

volScalarField zeta
(
     IOobject
         (
	     "zeta",
	     runTime.timeName(),
	     mesh,
	     IOobject::MUST_READ,
	     IOobject::AUTO_WRITE
	 ),
	     mesh
);

volScalarField lambda
(
     IOobject
         (
	     "lambda",
	     runTime.timeName(),
	     mesh,
	     IOobject::MUST_READ,
	     IOobject::AUTO_WRITE
	 ),
	     mesh
);

volScalarField gg
(
     IOobject
         (
	     "gg",
	     runTime.timeName(),
	     mesh,
	     IOobject::MUST_READ,
	     IOobject::AUTO_WRITE
	 ),
	     mesh
);

volScalarField gg_int
(
     IOobject
         (
	     "gg_int",
	     runTime.timeName(),
	     mesh,
	     IOobject::MUST_READ,
	     IOobject::AUTO_WRITE
	 ),
	     mesh
);

ii=0;
forAll(lambda,ii){
	lambda[ii]=scalar_lambda;
}
Info<<"Lambda "<<endl;

// ALFA1
scalar jj=0; scalar sum=0; scalar kk=0; ii=0; scalar ww=0;
forAll(alfa,cellI){
         alfa[cellI]=0;
	    }
 while(jj<(nx*ny*nz)){ 
	 sum=sum+fi[jj]*fi[jj]*(1/vv)*dx*dz;
    
	 
	 if (kk==(nx-1))
	 {   
         alfa[ii]=alfa[ii]+sum;
		 kk=0; 
		 ii=ii+1;
		 sum=0;
	 } else
	 {
	      kk=kk+1;
	 }
	    if(ww==(nx*ny-1)){
		 ww=0; 
		 ii=0;
	 }else{
		 ww=ww+1;		 
	 } 
	 jj=jj+1;
	 
    }
	Info<<"alfa "<<endl;
	alfa.write();
	
// DELTA
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
         delta[cellI]=0;
	    }
 while(jj<(nx*ny*nz)){ 
     sum=sum+fi[jj]*(DD[jj]+dDD[jj])*fi[jj]*dx*dz;
     
	 
	 if (kk==(nx-1))
	 {   
         delta[ii]=delta[ii]+sum;
		 kk=0; 
		 ii=ii+1;
		 sum=0;
	 } else
	 {
	     kk=kk+1;
	 }
	 if(ww==(nx*ny-1)){
		 ww=0; 
		 ii=0;
	 }else{
		 ww=ww+1;		 
	 }
	 jj=jj+1;
	
    }
delta.write();
Info<<"delta "<<endl;
	

// GAMMA
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
         gamma[cellI]=0;
	    }
 while(jj<(nx*ny*nz)){ 
	 sum=sum+fi[jj]*((DD[jj]+dDD[jj])*grad_fi_y[jj]+grad_DD_fi_y[jj])*dx*dz;
	 
	 
	 if (kk==(nx-1))
	 {   
         gamma[ii]=gamma[ii]+sum;
		 kk=0; 
		 ii=ii+1;
		 sum=0;
	 } else
	 {
	     kk=kk+1;
	 }
	 
	 if(ww==(nx*ny-1)){
		 ww=0; 
		 ii=0;
	 }else{
		 ww=ww+1;		 
	 }

	 jj=jj+1;
	
    } 
gamma.write();
	Info<<"gamma "<<endl;
		
// ETA
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
         eta[cellI]=0;
	    }
 while(jj<(nx*ny*nz)){ 
       sum=sum+fi[jj]*(lap_dDD_fi[jj]+((1-beta)*dnisf[jj]-ds[jj])*fi[jj]-beta*nisf[jj]*fi[jj])*dx*dz;
	 
	 if (kk==(nx-1))
	 {   
         eta[ii]=eta[ii]+sum;
		 kk=0; 
		 ii=ii+1;
		 sum=0;
	 } else
	 {
	     kk=kk+1;
	 }
    if(ww==(nx*ny-1)){
		 ww=0; 
		 ii=0;
	 }else{
		 ww=ww+1;		 
	 }
	 jj=jj+1;
	
 }       
eta.write();
	Info<<"eta "<<endl;
		
			
	// ZETA
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
         zeta[cellI]=0;
	    }
 while(jj<(nx*ny*nz)){ 
	 sum=sum+fi[jj]*fi[jj]*beta*(nisf[jj]+dnisf[jj])*dx*dz;
	 
	 
	 if (kk==(nx-1))
	 {   
         zeta[ii]=zeta[ii]+sum;
		 kk=0; 
		 ii=ii+1;
		 sum=0;
	 } else
	 {
	     kk=kk+1;
	 }

if(ww==(nx*ny-1)){
		 ww=0; 
		 ii=0;
	 }else{
		 ww=ww+1;		 
	 }
	 jj=jj+1;
	 
    } 
zeta.write();
	Info<<"Zeta "<<endl;
	
	 // GG
jj=0; sum=0; kk=0; ii=0; ww=0; 
forAll(alfa,cellI){
         gg[cellI]=0;
	    }
 while(jj<(nx*ny*nz)){ 
	 sum=sum+fi[jj]*(beta/scalar_lambda)*(nisf[jj]*fi[jj])*dx*dz;
	 	 
	 if (kk==(nx-1))
	 {   
         gg[ii]=gg[ii]+sum;
		 kk=0; 
		 ii=ii+1;
		 sum=0;
	 } else
	 {
	     kk=kk+1;
	 }
if(ww==(nx*ny-1)){
		 ww=0; 
		 ii=0;
	 }else{
		 ww=ww+1;		 
	 }
	 jj=jj+1;
	 
    }  
gg.write();
	scalar gg_0;
	  forAll(gg,cellI)
        {
	     gg_0=gg_0+gg[cellI]*dy/HH;
	    }
	Info<<"GG "<<gg_0<<endl;
	
forAll(gg_int,cellI)
          {
	   gg_int[cellI]=gg_0;
          }

    gg_int.write();

