/* gptitan: photochimie */
/* GCCM */

/* tout est passe en simple precision */
/* sauf pour l'inversion de la matrice */

/* nitriles et hydrocarbures separes pour l'inversion */

#include "titan.h"

void gptitan_(const int *NLAT,
     float *RA, float *TEMP, float *NB, 
     char CORPS[][10], float Y[][NLEV], float *FTOP,
     float *YTOP, float *DECLIN, float *FIN, int *LAT, float *MASS, 
     float *botCH4, 
     float KRPD[][NLEV][RDISS+1][15], float KRATE[][NLEV],
     int reactif[][5], int *nom_prod, int *nom_perte, 
     int prod[][200], int perte[][200][2], int *aerprod, int *utilaer, 
     float MAER[][NLEV], float PRODAER[][NLEV], 
     float CSN[][NLEV], float CSH[][NLEV],
     int *htoh2, float *surfhaze)
{
   char   outlog[100],corps[100][10];
   int    dec,declinint,i,j,k,l;
   int    ireac,ncom1,ncom2;
   float  *fl,*fp,*mu,**jac,*ym1;
   float  *fluxtop,fluxCH4;
   float  cm,conv,cp,delta,deltamax,dv,dr,drp,drm;
   float  rr,np,nm,factdec,s,test,time,ts,v;
   double *fd,**jacd;
   char   str2[15];
   FILE   *out;

/* va avec htoh2 */
   float  dyh,dyh2;

/* va avec aer */
   float  dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
   float  **k_dep,**faer;
   float  *productaer,*csurn,*csurh,*mmolaer;

   if( (*aerprod) == 1 )
   {
    k_dep = rm2d( 1, 5, 1, 3 );     /* k en s-1, reactions d'initiation */
    faer  = rm2d( 1, 5, 1, 3 );     /* fraction de chaque compose */
    productaer  = rm1d( 0, 3 );     /* local production rate by pathways */
    mmolaer     = rm1d( 0, 3 );     /* local molar mass by pathways */
    csurn = rm1d( 0, 3 );           /* local C/N by pathways */
    csurh = rm1d( 0, 3 );           /* local C/H by pathways */
   }

/* DEBUG
      printf("CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
*/

   for( i = 0; i <= NC; i++)
   {
     strcpy( corps[i], CORPS[i] );
     corps[i][strcspn(CORPS[i], " ")] = '\0';
   }
   
   strcpy( outlog, "chimietitan" );
   strcat( outlog, ".log" );
   deltamax = 1.e5;

/*  DEBUG
            out = fopen( outlog, "a" );
   fprintf(out,"CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
            fclose( out );
*/
   ym1     = rm1d( 0,   NC-1 );
   fl      = rm1d( 0,   NC-1 );
   fp      = rm1d( 0,   NC-1 );
   fd      = dm1d( 0,   NC-1 );
   mu      = rm1d( 0, NLEV-1 );
   fluxtop = rm1d( 0,   NC-1 );
   jac     = rm2d( 0,   NC-1, 0,   NC-1 );
   jacd    = dm2d( 0,   NC-1, 0,   NC-1 );

/* DEBUG */
/*
            out = fopen( "err.log", "a" );
            fprintf( out,"%s\n", );
            fclose( out );
*/

/* initialisation krate pour dissociations */

   if( ( (*DECLIN) *10 + 267 ) < 14. ) 
   {
     declinint = 0;
     dec = 0;
   }
   else 
   {  
       if( ( (*DECLIN) * 10 + 267 ) > 520. ) 
       {
          declinint = 14;
          dec = 534;
       }
       else 
       {
          declinint = 1;
          dec       = 27;
          while( ( (*DECLIN) * 10 + 267 ) >= (float)(dec+20) ) 
          {
            dec += 40;
            declinint++;
          }
       }
   }
   if( ( (*DECLIN) >= -24. ) && ( (*DECLIN) <= 24. ) ) 
           factdec = ( (*DECLIN) - (dec-267)/10. ) / 4.;
   else
           factdec = ( (*DECLIN) - (dec-267)/10. ) / 2.7;

   for( i = 0; i <= RDISS; i++ )   /* RDISS correspond a la dissociation de N2 par les GCR */
      for( j = 0; j <= NLEV-1; j++ )
      { 
         if( factdec < 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint-1] * fabs(factdec) 
                                         + KRPD[*LAT][j][i][declinint] * ( 1 + factdec);
         if( factdec > 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint+1] * factdec 
                                         + KRPD[*LAT][j][i][declinint] * ( 1 - factdec);
         if( factdec == 0. ) KRATE[i][j] = KRPD[*LAT][j][i][declinint];
      }
        
/* initialisation mu, CH4 au sol */
     
   for( j = 0; j <= NLEV-1; j++ )
   {
      mu[j] = 0.0e0;
      for( i = 0; i <= ST-1; i++ ) 
      {
         if( ( strcmp(corps[i], "CH4") == 0 ) && ( Y[i][j] <= *botCH4 ) && ( j == 0 ) )
         {
             fluxCH4 = (*botCH4 - Y[i][j]);
             Y[i][j] = *botCH4;
         }
         mu[j] += ( MASS[i] * Y[i][j] );
      }
   }

/* ****************** */
/*  Main loop: level  */
/* ****************** */

   for( j = NLEV-1; j >= 0; j-- )
   {

/* DEBUG
            out = fopen( outlog, "a" );
        fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
            fclose( out );

            out = fopen( "profils.log", "w" );
    fprintf(out,"%d %e %e %e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
    for (i=0;i<=NREAC-1;i++) fprintf(out,"%d %e\n",i,KRATE[i][j]);
    for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
            fclose( out );
*/

      time     = ts = 0.0e0;
/*      delta    = (*FIN);  */
      delta    = 1.e-3;

      for( i = 0; i <= ST-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);

/* ++++++++++++ */
/*  time loop.  */
/* ++++++++++++ */

      while( time < (*FIN) )     
      {

/* Calcul de f et de la matrice jacobienne */
/* --------------------------------------- */

         for( i = 0; i <= ST-1; i++ )     /* productions et pertes chimiques */
         {
            Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */

            fp[i] = fl[i] = 0.0e0;                    /* init for next step */
            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;

            for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
            {
               ireac = prod[i][l];                  /* Number of the reaction involves. */
               ncom1 = reactif[ireac][0];           /* First compound which reacts. */
               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
               {
                  jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
                  fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
               }
               else                                 /* General case. */
               {
                  ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
                  jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
                  jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
                  fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
               }
            }
            
            for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
            {
               ireac = perte[i][l][0];              /* Reaction number. */
               ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
               {
                  jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
                  fl[i]     += ( KRATE[ireac][j] * NB[j] );
               }
               else                                 /* General case. */
               {
                  jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
                  jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
                  fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
               }
            }
         }

/* Aerosols */
/* -------- */
         if( (*aerprod) == 1 )
         {
             aer(corps,TEMP,NB,Y,&j,k_dep,faer,
              &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
              mmolaer,productaer,csurn,csurh);

             for( i = 0; i <= 3; i++ )
             {
               PRODAER[i][j] = productaer[i];
                  MAER[i][j] = mmolaer[i];
                   CSN[i][j] = csurn[i];
                   CSH[i][j] = csurh[i];
             }
/* DEBUG
printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
      printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
              fp[utilaer[2]],dyc2h2*NB[j]);
if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
      printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
              fp[utilaer[5]],dyhcn*NB[j]);
if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
      printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
              fp[utilaer[6]],dyhc3n*NB[j]);
if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
      printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
              fp[utilaer[13]],dynccn*NB[j]);
if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
      printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
              fp[utilaer[14]],dych3cn*NB[j]);
if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
      printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
              fp[utilaer[15]],dyc2h3cn*NB[j]);
*/

             fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
             fp[utilaer[5]] -= (    dyhcn * NB[j] );
             fp[utilaer[6]] -= (   dyhc3n * NB[j] );
             fp[utilaer[13]]-= (   dynccn * NB[j] );
             fp[utilaer[14]]-= (  dych3cn * NB[j] );
             fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
             if( Y[utilaer[2]][j]  != 0.0 )
       jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
             if( Y[utilaer[5]][j]  != 0.0 )
       jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
             if( Y[utilaer[6]][j]  != 0.0 )
       jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
             if( Y[utilaer[13]][j] != 0.0 )
     jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
             if( Y[utilaer[14]][j] != 0.0 )
     jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
             if( Y[utilaer[15]][j] != 0.0 )
     jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
         }
            
/* H -> H2 on haze particles */
/* ------------------------- */
         if( (*htoh2) == 1 )
         {
              heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
/* dyh <= 0 / 1.0 en adsor., 1 en reac. */

/* DEBUG
printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
*/

              fp[utilaer[0]] += ( dyh  * NB[j] );
              fp[utilaer[1]] += ( dyh2 * NB[j] );
              if( Y[utilaer[0]][j] != 0.0 )
       jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
         }

         for( i = 0; i <= ST-1; i++ )     /* finition pour inversion */
         {
            for( k = 0; k <= ST-1; k++ )
            {
               jac[i][k] *= ( -THETA * delta );     /* Correction time step. */
               if( k == i ) jac[k][k] += NB[j];     /* Correction diagonal. */
               jacd[i][k] = (double)jac[i][k];
            }
            
            fd[i] = (double)(delta * ( fp[i] - Y[i][j]*fl[i] ));
         }

/*         for( i = ST; i <= NC-1; i++ )    pas d'inversion (soot,prod): que faire?
         {
            Y[i][j] = ??? ;
         }
*/

/* Inversion of matrix cf method LU */
/* -------------------------------- */

/* inversion by blocs: */
/* Hydrocarbons */

         solve( jacd, fd, 0, NHC-1 );             

         for( i = 0; i <= NHC-1; i++ )
         {
            Y[i][j] += (float)fd[i];
            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
         }

/* Nitriles */

         solve( jacd, fd, NHC, ST-1 );             

         for( i = NHC+1; i <= ST-1; i++ )
         {
            Y[i][j] += (float)fd[i];
            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
         }

/* end inversion by blocs: */

         for( i = 0; i <= ST-1; i++ )
         {

/* CH4 au sol */
/* ---------- */

            if( ( strcmp(corps[i], "CH4") == 0 ) && ( j == 0 ) && ( Y[i][j] < *botCH4 ) )
            {
                fluxCH4 += (*botCH4 - Y[i][0]);
                Y[i][0] = *botCH4;
            }

         }
         
/* test evolution delta */
/* -------------------- */

         for( i = 0; i <= ST-1; i++ )
         {
            test = 1.0e-15; 
            if( ( Y[i][j] > test ) && ( ym1[i] > test ) )
            {
               conv = fabs( Y[i][j] - ym1[i] ) / ym1[i];

               if( conv > ts )
               {
/*
                  if( conv >= 0.1 )
                  {
                     out = fopen( outlog, "a" );
                     fprintf( out, "Lat no %d; declin:%e;", (*LAT)+1, (*DECLIN) );
                     fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
                     fclose( out );
                  }
*/
                  ts = conv;
               }
            }
         }

         if( ts < 0.1e0 )
         {
            for( i = 0; i <= ST-1; i++ )
                 if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
                 {
                  out = fopen( outlog, "a" );
                  fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
                           corps[i], ym1[i], Y[i][j], j );
                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
                  fclose( out );
                /*  exit(0); */
                  Y[i][j] = 1.e-20;
                 }
            for( i = 0; i <= NC-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
            time += delta;
            if(   ts < 1.00e-5 )                      delta *= 1.0e2;
            if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
            if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
            if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
            if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
         
            delta = min( deltamax, delta );
         }
         else
         {
            for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i];

            if(   ts > 0.8 )                    delta *= 1.e-6;
            if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
            if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
            if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
            if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
            if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
         }
         ts = 0.0e0;
      }               

/* +++++++++++++++++++ */
/*  end of time loop.  */
/* +++++++++++++++++++ */

      for( i = 0; i <= ST-1; i++ ) 
         if( ( strcmp(corps[i],"CH4") == 0 ) && ( j == 0 ) )
            fluxCH4 *= ( MASS[i]/(6.022e23*time) ); 

   }
   
/* **************** */        
/* end of main loop */
/* **************** */        

/* Plafond: maintient de la composition */
/* ------------------------------------ */
/* 
   for( i = 0; i <= ST-1; i++ ) 
      if( YTOP[i] != 0.0e0 )
      {
             fluxtop[i]  = (Y[i][NLEV-1] - YTOP[i]) * MASS[i]/(6.022e23*(*FIN));
             Y[i][NLEV-1]= YTOP[i];
      }
*/
/* Plafond: !! OU !!  flux vertical     */
/* ------------------------------------ */

   for( i = 0; i <= ST-1; i++ ) 
      if( FTOP[i] != 0.0e0 )
      {
             fluxtop[i]    = (- FTOP[i]/NB[NLEV-2]) * MASS[i]/6.022e23;
             Y[i][NLEV-2] += FTOP[i]/NB[NLEV-2]*(*FIN);
             Y[i][NLEV-2]  = max(Y[i][NLEV-2],0.0e0);
// on ajuste aussi le niveau dans la derniere couche...
// pour eviter les effets vers le haut
             Y[i][NLEV-1]  = Y[i][NLEV-2];
      }

/* Niveau de N2 */
/* ------------ */
      
   for( j = 0; j <= NLEV-1; j++ ) 
   {
      conv = 0.0e0;
      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
   }

   if( (*aerprod) == 1 )
   {
    frm2d( k_dep, 1, 5, 1 );
    frm2d(  faer, 1, 5, 1 );
    frm1d( productaer,  0 );
    frm1d( mmolaer,  0 );
    frm1d( csurn, 0 );
    frm1d( csurh, 0 );
   }

   frm1d(     ym1, 0 );
   frm1d(      fl, 0 );
   frm1d(      fp, 0 );
   fdm1d(      fd, 0 );
   frm1d(      mu, 0 );
   frm1d( fluxtop, 0 );
   frm2d(     jac, 0,   NC-1, 0 );
   fdm2d(    jacd, 0,   NC-1, 0 );
}