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

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

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


volScalarField zeta
(
     IOobject
         (
	     "zeta",
	     runTime.timeName(),
		  "/Weights",
	     mesh,
	     IOobject::NO_READ,
	     IOobject::AUTO_WRITE
	 ),
	     mesh
);

volScalarField lambda
(
     IOobject
         (
	     "lambda",
	     runTime.timeName(),
		  "/Weights",
	     mesh,
	     IOobject::NO_READ,
	     IOobject::AUTO_WRITE
	 ),
	     mesh
);

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

// ALFA1
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
	alfa[cellI]=0;
}
 while(jj<(nx*ny*nz)){ 
	 sum=sum+(psi1[jj]*fi1[jj]*(1/vv1)+psi2[jj]*fi2[jj]*(1/vv2))*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;
	 
    }
alfa.write();
Info<<"alfa "<<endl;
		
// DELTA
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
	delta[cellI]=0;
}
 while(jj<(nx*ny*nz)){ 
     sum=sum+(psi1[jj]*(DD1[jj]+dDD1[jj])*fi1[jj]+psi2[jj]*(DD2[jj]+dDD2[jj])*fi2[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;
		
// GAMMA1
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
	gamma[cellI]=0;
}
 while(jj<(nx*ny*nz)){ 
	 sum=sum+(psi1[jj]*((DD1[jj]+dDD1[jj])*grad_fi1_y[jj]+grad_DD1_fi1_y[jj])+psi2[jj]*((DD2[jj]+dDD2[jj])*grad_fi2_y[jj]+grad_DD2_fi2_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;
	
	
// ETA1
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
	eta[cellI]=0;
}
 while(jj<(nx*ny*nz)){ 
       sum=sum+(psi1[jj]*(lap_dDD1_fi1[jj]+((1-beta)*dnisf1[jj]-ds1[jj])*fi1[jj]+(1-beta)*dnisf2[jj]*fi2[jj]-beta*nisf1[jj]*fi1[jj]-beta*nisf2[jj]*fi2[jj])+psi2[jj]*(lap_dDD2_fi2[jj]+ds12[jj]*fi1[jj]-ds2[jj]*fi2[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;
	
		
	// ZETA1
jj=0; sum=0; kk=0; ii=0; ww=0;
forAll(alfa,cellI){
	zeta[cellI]=0;
}
 while(jj<(nx*ny*nz)){ 
	 sum=sum+(psi1[jj]*fi1[jj]*beta*(nisf1[jj]+dnisf1[jj])+psi2[jj]*fi2[jj]*beta*(nisf2[jj]+dnisf2[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;
scalar gg_0=0;
forAll(alfa,cellI){
	gg[cellI]=0;
}
 while(jj<(nx*ny*nz)){ 
	 sum=sum+psi1[jj]*(beta/scalar_lambda)*(nisf1[jj]*fi1[jj]+nisf2[jj]*fi2[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();
          forAll(gg,cellI)
          {
	     gg_0=gg_0+gg[cellI]*dy/HH;
	  }
	Info<<"GG "<<gg_0<<endl;
	
	          forAll(gg,ii)
	          {
              gg_int[ii]=gg_0;
                  }
       

scalar AA0=0;

	for(ii=0;ii<(nx*ny*nz);ii++) 
	{
	AA0+= (1-beta)*(nisf1[ii]*fi1[ii]+nisf2[ii]*fi2[ii])*dx*dy*dz;
	}
	cout<<"AA0="<<AA0;
