%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%%%								                                   %%%%
      %%%%                    Author = Dr Andrea Boghi	         		   %%%%
      %%%%								                                   %%%%
      %%%%                           LEGEND	   		    	               %%%%
      %%%%								                                   %%%%
      %%%%                        C(1) = CO_2		     		           %%%%
      %%%%                        C(2) = CH_4		    		           %%%%
      %%%%                        C(3) = N_2				               %%%%
      %%%%								                                   %%%%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INPUT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      function [thetaG,CL,CG,r,z,nr,nz] = modelThetaG(x,y)

      nr = 60;
      nz = 200;
      nit = 64000; %64000 
      kE = x(1)	     %%% 
      kT = x(2)	     %%%  (1-100)
      kV = x(3)	     %%% mol/(dm s) (1-100) 
      Ktheta = x(4)
      % % % SCO20 =x(4) 	 %%% mol/l s-1 (1-100)
      % % % Ktheta = x(5)	 %%% mol/l s-1 (1-100)
      % % % Ktheta = 1.608245131690465d0*SCO20/kE
      % % % SCO20 = 6.217957575589712d-01*Ktheta*kE %%% alphaCH4 = 2.d0
      SCO20 = 4.292683555238209d-01*Ktheta*kE %%% alphaCH4 = 1.d0
	  PH2O = 5.d0;     %%% kPa
      Lz = 1.7d0;   %%% dm
      Lr = 0.5d0;   %%% dm
      TC = 23.d0;   %%% �C
      theta =0.66d0;
      fL = 0.55d0;
      pH = 7.d0;
      CFL = 4.5d-01;
     
      %%%%% GRID
      dr = Lr/(nr-1);
      for i=1:nr;
         r(i) = dr*(i-1);
      end;
      dr2 = dr*dr;

      dz = Lz/(nz-1);
      for j=1:nz;
         z(j) = dz*(j-1) -Lz;
      end;
      dz2 = dz*dz;

      alpha = dr2*dz2/(dr2 +dz2);
      beta = dr2/(dr2 +dz2);

      %%%% Calculate velocities
      [ur,uz] = calculateVelocities(r,z,nr,nz);
	  
      %%%% Parameters
      [CL0,thetaG0,thetaL0,Lv,CH2O,CD,Sm,Dm,KH,C0,DL] = Parameters(theta,TC,kE,kT,kV,SCO20,Ktheta,PH2O,pH,y,r,z,nr,nz);

      for j=1:nz;
	 for i=1:nr;
	    for k=1:3;
	       CL(i,j,k) = CL0(k);
	    end
	    [thetaL(i,j),thetaG(i,j)] = calculateThetas(CL(i,j,:),KH,CD,CH2O,theta);
	 end;
      end;

      for i=1:nr;
	     for k=1:3;
			CL(i,nz,k) = C0(k);
		 end;
         [thetaL(i,nz),thetaG(i,nz)] = calculateThetas(CL(i,nz,:),KH,CD,CH2O,theta);
      end;

      %%%%% iterative cycle
      for iit=1:nit;

         %%% Under-Relaxation
	 for k=1:3;
	    for j=1:nz;
	       for i=1:nr;
                  CL_old(i,j,k) = CL(i,j,k);
	       end;
        end;
     end;

  	 %%% R Derivative
     for j=1:nz;
	    dthetaLdr(:,j) = firstOrdDer(thetaL(:,j),dr,nr);
	 end;

         %%% Z Derivative
	 for i=1:nr;
	    dthetaLdz(i,:) = firstOrdDer(thetaL(i,:),dz,nz);
	 end;

	 for j=1:nz;
	    for i=1:nr;
	       [thetaL(i,j),thetaG(i,j)] = calculateThetas(CL(i,j,:),KH,CD,CH2O,theta);
	    end;
	 end;

	 %%% Solve
	 for k=1:3;
	    for j=2:nz-1;
	       for i=2:nr-1;
                  [ap,Sp] = coefficientsPP(Sm(i,j,k),Dm(i,j,k),KH(k),DL(k),kE,theta,fL,thetaL(i,j),alpha);
                  [as,an] = coefficientsNS(DL(k),uz(i,j),fL,thetaL(i,j),dthetaLdz(i,j),dz,beta);
                  [ae,aw] = coefficientsWE(DL(k),r(i),ur(i,j),fL,thetaL(i,j),dthetaLdr(i,j),dr,beta);
          	  CL(i,j,k) = (ae*CL(i-1,j,k) +aw*CL(i+1,j,k) +as*CL(i,j-1,k) +an*CL(i,j+1,k) +Sp)/ap;
	       end; 
        end;

	    for j=2:nz-1;
               [ap,Sp] = coefficientsPP(Sm(1,j,k),Dm(1,j,k),KH(k),DL(k),kE,theta,fL,thetaL(1,j),alpha);
               [as,an] = coefficientsNS(DL(k),uz(1,j),fL,thetaL(1,j),dthetaLdz(1,j),dz,beta);
	       CL(1,j,k) = (4.d0*(1.d0 -beta)*CL(2,j,k) +as*CL(1,j-1,k) +an*CL(1,j+1,k) +Sp)/(ap+2.d0*(1.d0 -beta));
	    end;

	    for j=2:nz-1;
               [ap,Sp] = coefficientsPP(Sm(nr,j,k),Dm(nr,j,k),KH(k),DL(k),kE,theta,fL,thetaL(nr,j),alpha);
               [as,an] = coefficientsNS(DL(k),uz(nr,j),fL,thetaL(nr,j),dthetaLdz(nr,j),dz,beta);
	       CL(nr,j,k) = (2.d0*(1.d0 -beta)*CL(nr-1,j,k) +as*CL(nr,j-1,k) +an*CL(nr,j+1,k) +Sp)/ap;
	    end;

	    for i=2:nr-1;
               [ap,Sp] = coefficientsPP(Sm(i,1,k),Dm(i,1,k),KH(k),DL(k),kE,theta,fL,thetaL(i,1),alpha);
               [ae,aw] = coefficientsWE(DL(k),r(i),ur(i,1),fL,thetaL(i,1),dthetaLdr(i,1),dr,beta);
	       CL(i,1,k) = (ae*CL(i-1,1,k) +aw*CL(i+1,1,k) +2.d0*beta*CL(i,2,k) +Sp)/ap;
	    end;

            [ap,Sp] = coefficientsPP(Sm(1,1,k),Dm(1,1,k),KH(k),DL(k),kE,theta,fL,thetaL(1,1),alpha);
	    CL(1,1,k) = (4.d0*(1.d0 -beta)*CL(2,1,k) +2.d0*beta*CL(1,2,k) +Sp)/(ap+2.d0*(1.d0 -beta));
            [ap,Sp] = coefficientsPP(Sm(nr,1,k),Dm(nr,1,k),KH(k),DL(k),kE,theta,fL,thetaL(nr,1),alpha);
	    CL(nr,1,k) = (2.d0*(1.d0 -beta)*CL(nr-1,1,k) +2.d0*beta*CL(nr,2,k) +Sp)/ap;

	    %%% Fluxes
	    for i=1:nr;
	       Jn(i,k) = -thetaL(i,nz)*fL*DL(k)*(3.d0*CL(i,nz,k) -4.d0*CL(i,nz-1,k) ........
	                 +CL(i,nz-2,k))/(2.d0*dz) +uz(i,nz)*CL(i,nz,k);
	       Js(i,k) = -thetaL(i,1)*fL*DL(k)*(-3.d0*CL(i,1,k) +4.d0*CL(i,2,k) ........
		         -CL(i,3,k))/(2.d0*dz) +uz(i,1)*CL(i,1,k);
            end;

	    for j=1:nz;
	       Je(j,k) = -thetaL(nr,j)*fL*DL(k)*(3.d0*CL(nr,j,k) -4.d0*CL(nr-1,j,k) ........
		         +CL(nr-2,j,k))/(2.d0*dr) +ur(nr,j)*CL(nr,j,k);
            end;

	 end;
	
	 %%% sources
         Gm = sources(Sm,Dm,CL,KH,thetaG,kE,nr,nz);

         Gint = volumeInt(Gm,r,dr,dz,nr,nz);

	 %%% fluxes
         JintZ = surfaceZInt(Jn-Js,r,dr,nr);

	 JintR = surfaceRInt(Je,r,dz,nr,nz);

	 Jint = JintZ +JintR;

         %%% Under-Relaxation
	 for k=1:3;
	    for j=1:nz;
	       for i=1:nr;
                  CL(i,j,k) = CFL*CL(i,j,k) +(1.d0-CFL)*CL_old(i,j,k);
	       end;
            end;
         end;

      end;
      
      for k=1:3;
	    for j=1:nz;
	       for i=1:nr;
          	  CG(i,j,k) = CL(i,j,k)/KH(k);
	       end; 
        end;
      end;


      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						           			 * 
      %%%*                             Parameters 		 		     	 *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function [CL0,thetaG0,thetaL0,Lv,CH2O,CD,Sm,Dm,KH,C0,DL] = Parameters(theta,TC,kE,kT,kV,SCO20,Ktheta,PH2O,pH,y,r,z,nr,nz)

          for j=1:nz;
             for i=1:nr;
                Lv(i,j) = 1.d0*y(1)*exp( -( ........
                          ( ( (r(i) -y(2))*cos(y(6)) -(z(j) -y(4))*sin(y(6)) )^2.d0 )/(2.d0*y(3)^2.d0) ........
                       +  ( ( (r(i) -y(2))*sin(y(6)) +(z(j) -y(4))*cos(y(6)) )^2.d0 )/(2.d0*y(5)^2.d0) ) );
             end;
          end;

	  TK = TC +273.15;
	  RU = 8.314d0; %%% dm3 kPa K-1 mol-1
	  alphaCH4 = 1.d0; %%%5.d-02;
	  K1 = 4.45d-07; %%% mol/l
	  H2CO30 = 2.d-02; %%% mol/l
      H2CO3i = 8.d-02; %%% mol/l

	  DG(1) = 1.55d-03; %%% dm2 s-1
      DG(2) = 2.20d-03; %%% dm2 s-1
      %DG(3) = 2.04d-03; %%% dm2 s-1
	  %%%%% ATTENTION: THE ABSENCE OF N2 TRANSFER REQUIRES TO MODIFY DG(3)
	  DG(3) = 0.d0; %%% dm3 s-1

      DL(1) = (1.93d-07 +(1.18d-07)*K1*10.d0^pH); %%% dm2 s-1
	  %DL(1) = 1.18d-07; %%% dm2 s-1
      DL(2) = 1.73d-07; %%% dm2 s-1
      DL(3) = 2.02d-07; %%% dm2 s-1

      %%%% ATTENTION: ACIDS
      KH(1) = 8.30d-01;
      KH(2) = 3.15d-02;
      KH(3) = 1.62d-02;
	  %%%% ATTENTION: ACIDS
      KH(1) = KH(1)*(1.d0 +K1*10.d0^pH);

      % % % DON'T CHANGE
	  CG0(1) = 1.20d-05; %%% mol/l
      CG0(2) = 9.21d-08; %%% mol/l
      CG0(3) = 3.22d-02; %%% mol/l
      
      %%% Initialization
      thetaG0 = (Ktheta*PH2O/2.d0)*(1.d0 +sqrt(1.d0 +(4.d0*RU*TK*SCO20*(1.d0+alphaCH4))/(kE*Ktheta*PH2O*PH2O)));
      thetaL0 = theta -thetaG0;
      CL0(1) = KH(1)*SCO20/(kE*thetaG0);
      CL0(2) = KH(2)*alphaCH4*SCO20/(kE*thetaG0);
      CL0(3) = KH(3)*0.d0;

      for j=1:nz;
         for i=1:nr;
	        S(i,j,1) = SCO20 +kV*Lv(i,j);
            S(i,j,2) = alphaCH4*S(i,j,1);
	        S(i,j,3) = 0.d0;
         end;
      end;

	  %%% OUT
	  CH2O = PH2O/(RU*TK); %%% mol/l
	  CD = 1.d0/(Ktheta*RU*TK); %%% mol/l

	  for k=1:3;
         for j=1:nz;
            for i=1:nr;
	     	   Dm(i,j,k) = kT*Lv(i,j)*DG(k);
		       Sm(i,j,k) = S(i,j,k) +Dm(i,j,k)*CG0(k);
            end;
         end;
	  end;

      % % % DON'T CHANGE

      C0(1) = 2.d-03;  %%% mol/l 8.d-02  !!! mol/l	  2.d-03  !!! mol/l
	  C0(2) = 2.9d-09; %%% mol/l
 	  C0(3) = 5.0d-04; %%% mol/l

	  for k=1:3;
	     CG0(k) = C0(k)/KH(k);
	  end;


	  %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						                										 * 
      %%%*                          Calculate thetas        		     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

	  function [thetaL,thetaG] = calculateThetas(CL,KH,CD,CH2O,theta)

	  thetaG = ( CL(1)/KH(1) + CL(2)/KH(2) + CL(3)/KH(3) + CH2O )/CD;
	  thetaL = theta - thetaG;
 
	  %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						                										 * 
      %%%*                          Calculate velocities        		     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

	  function [ur,uz] = calculateVelocities(r,z,nr,nz)

	  for j=1:nz;
	     for i=1:nr;
		    ur(i,j) = 0.d0;
			uz(i,j) = 0.d0;
		 end;
	  end;

	  %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						                										 * 
      %%%*                          First Oder Derivative        		     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

	  function dads = firstOrdDer(a,ds,ns)

      dads(1) = (-3.d0*a(1) +4.d0*a(2) -a(3))/(2.d0*ds);
      for j=2:ns-1;
         dads(j) = (a(j+1)-a(j-1))/(2.d0*ds);
      end;
      dads(ns) = (3.d0*a(ns) -4.d0*a(ns-1) +a(ns-2))/(2.d0*ds);

      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						               								 * 
      %%%*                          coefficients        		     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function [ap,Sp] = coefficientsPP(Sm,Dm,KH,DL,kE,theta,fL,thetaL,alpha)
    
      Sp = alpha*Sm/(thetaL*fL*DL);
      ap = 2.d0 +alpha*(Dm +kE*(theta -thetaL))/(KH*thetaL*fL*DL);

      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*	 									 * 
      %%%*                          coefficients        		     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function [as,an] = coefficientsNS(DL,uz,fL,thetaL,dthetaLdz,dz,beta)

      an = beta*(1.d0 +(dthetaLdz/thetaL -uz/(thetaL*fL*DL))*dz/2.d0);
      as = beta*(1.d0 -(dthetaLdz/thetaL -uz/(thetaL*fL*DL))*dz/2.d0);

      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						              			 * 
      %%%*                          coefficients        		     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function [ae,aw] = coefficientsWE(DL,r,ur,fL,thetaL,dthetaLdr,dr,beta)

      aw = (1.d0 -beta)*(1.d0 +(dthetaLdr/thetaL +1.d0/r -ur/(thetaL*fL*DL))*dr/2.d0);
      ae = (1.d0 -beta)*(1.d0 -(dthetaLdr/thetaL +1.d0/r -ur/(thetaL*fL*DL))*dr/2.d0);

	  %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						              											 * 
      %%%*                          Sources						     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function Gm = sources(Sm,Dm,CL,KH,thetaG,kE,nr,nz)

 	  for k=1:3;
         for j=1:nz;
            for i=1:nr;
	     	   Gm(i,j,k) = Sm(i,j,k) -((Dm(i,j,k) +kE*thetaG(i,j))/KH(k))*CL(i,j,k);
            end;
         end;
	  end;

	  %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						              											 * 
      %%%*                          Volume Integral						     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function Gint = volumeInt(Gm,r,dr,dz,nr,nz)

 	  for k=1:3;
	     Gint(k) = 0.d0;
         for j=1:nz-1;
            for i=1:nr-1;
	     	   Gint(k) = Gint(k) +((Gm(i,j,k) + Gm(i+1,j,k) + Gm(i,j+1,k) ........
			           + Gm(i+1,j+1,k))/4.d0)*((r(i)+r(i+1))/2.d0)*dr*dz;
            end;
         end;
	  end;

	  %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						              											 * 
      %%%*                          SurfaceZInt						     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function Jint = surfaceZInt(Jw,r,dr,nr)

 	  for k=1:3;
	     Jint(k) = 0.d0;
         for i=1:nr-1;
	     	Jint(k) = Jint(k) +((Jw(i,k) + Jw(i+1,k))/2.d0)*((r(i)+r(i+1))/2.d0)*dr;
         end;
	  end;

	  %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						              											 * 
      %%%*                          SurfaceRInt						     	         *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function Jint = surfaceRInt(Jw,r,dz,nr,nz)

 	  for k=1:3;
	     Jint(k) = 0.d0;
         for j=1:nz;
	     	Jint(k) = Jint(k) +Jw(j,k)*r(nr)*dz;
         end;
	  end;

	  %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
      %%%*						               											 * 
      %%%*                          Calculate Source-Sink						     	 *
      %%%*                                                                               *
      %%%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*

      function [St,Rt,Et] = calcSoSi(Lv,thetaG,CL,kE,kT,kV,SCO20,pH,nr,nz)

	  alphaCH4 = 1.d0; %%%5.d-02;
	  K1 = 4.45d-07; %%% mol/l

	  DG(1) = 1.55d-03 %%% dm2 s-1
      DG(2) = 2.20d-03; %%% dm2 s-1
      %DG(3) = 2.04d-03; %%% dm2 s-1
	  %%%%% ATTENTION: THE ABSENCE OF N2 TRANSFER REQUIRES TO MODIFY DG(3)
	  DG(3) = 0.d0; %%% dm3 s-1

      %%%% ATTENTION: ACIDS
      KH(1) = 8.30d-01;
      KH(2) = 3.15d-02;
      KH(3) = 1.62d-02;
	  %%%% ATTENTION: ACIDS
      KH(1) = KH(1)*(1.d0 +K1*10.d0^pH);

      % % % DON'T CHANGE
	  CG0(1) = 1.20d-05; %%% mol/l
      CG0(2) = 9.21d-08; %%% mol/l
      CG0(3) = 3.22d-02; %%% mol/l

      for j=1:nz;
	     for i=1:nr;
	        St(i,j,1) = SCO20 +kV*Lv(i,j);
			%print*, ' SCO2=',St(i,j,1);
            St(i,j,2) = alphaCH4*St(i,j,1);
	        St(i,j,3) = 0.d0;
		 end;
      end;

      for k=1:3;
         for j=1:nz;
		    for i=1:nr;
	           Rt(i,j,k) = kT*Lv(i,j)*DG(k)*(CL(i,j,k)/KH(k) -CG0(k));
               Et(i,j,k) = kE*thetaG(i,j)*CL(i,j,k)/KH(k);
			end;
         end;
      end;