tracer.c

波浪数值模拟

  }
  else {
    field2D_mult_const(DQDX, kappa);    /* if breaking is off then just use kappa */
    field2D_mult_const(DQDY, kappa);    /* DQDY *= kappa */
  }



  field2D_diffx(QDISS, DQDX, idx);         /* QDISS is now on eta points again (NX-1,M) */  
  field2D_diffy_v_to_eta(TMPQ, DQDY, idy);       /*  TMPQ  is now on Q points again  */  
  field2D_self_add(QDISS, TMPQ);                 /*  QDISS +=  d/dy ( hby * dQ/dy )   */  



}




void update_passive_tracerAB3(field2D *Q, field2D_time_derivative_var *DTVARS, const field2D_depth_var *Depth_Vars, double dt)
{

  const double *qt1d = DTVARS->QT1->data;
  const double *qt2d = DTVARS->QT2->data;
  const double *qt3d = DTVARS->QT3->data;

  const double *hd = REF2(Depth_Vars->H, 1);

  double *qd, *qd0;

  int i, j, M, N, N1;
  double ih, hdt;

  M = Q->M;
  N = Q->N;
  N1 = DTVARS->QT1->N;
  hdt = dt/12.0;

  qd = REF2(Q, 1);  


  /* loop over all the U and most V points */

  for (i=0;i<N1;i++) {
    for (j=0;j<M;j++) {
      ih = 1.0/(*hd++);
      *qd++ += hdt*(23.0* (*qt1d++) - 16.0* (*qt2d++) +5.0*(*qt3d++) ) * ih;
    }
  }


  /* next update the boundary conditions on passive tracer Q */
  /* No flux boundary conditions */

  qd0 = REF2(Q, 0);    /* for the shoreline boundary condition */
  for (j=0;j<M;j++) {
    qd0[j] = qd0[j+M];   // shore boundary condition
    *qd =  *(qd - M);   // far slip BC
    qd++;
  }


}


void update_passive_tracerAB2(field2D *Q, const field2D *QT1, const field2D *QT2,
                         const field2D_depth_var *Depth_Vars, const double dt)
{

  const double *qt1d = QT1->data;
  const double *qt2d = QT2->data;

  const double *hd = REF2(Depth_Vars->H, 1);

  double *qd, *qd0;

  int i, j, M, N, N1;
  double ih, hdt;

  M = Q->M;
  N = Q->N;
  N1 = QT1->N;
  hdt = 0.5*dt;

  /* loop over all the U and most V points */

  qd = REF2(Q, 1); 

  for (i=0;i<N1;i++) {
    for (j=0;j<M;j++) {
      ih = 1.0/(*hd++);
      *qd++ += hdt*(3.0* (*qt1d++) - *qt2d++) * ih;
    }
  }

  /* next update the boundary conditions on passive tracer Q */
  /* No flux boundary conditions */

  qd0 = Q->data;         // this is for the boundary conditions stuff

  for (j=0;j<M;j++) {
    qd0[j] = qd0[j+M];   // shore boundary condition
    *qd =  *(qd - M);   // far slip BC
    qd++;
  }


}



void update_passive_tracer_euler(field2D *Q, const field2D *QT, const field2D_depth_var *Depth_Vars, const double dt)
{

  const double *qtd = QT->data;
  const double *hd = REF2(Depth_Vars->H, 1);

  double *qd, *qd0;

  int i, j, M, N, N1;
  double ih;

  M = Q->M;
  N = Q->N;
  N1 = QT->N;

  qd = REF2(Q, 1);   


  /* loop over all the U and most V points */

  for (i=0;i<N1;i++) {
    for (j=0;j<M;j++) {
      ih = 1.0/(*hd++);
      *qd++ += dt * ( *qtd++ ) * ih;
    }
  }


  /* next update the boundary conditions on passive tracer Q */
  /* No flux boundary conditions */

  qd0 = Q->data;         // this is for the boundary conditions stuff
  for (j=0;j<M;j++) {
    qd0[j] = qd0[j+M];   // shore boundary condition
    *qd =  *(qd - M);   // far slip BC
    qd++;
  }

}



     

/* -------- OLD ROUTINES ------- */

/* Calculates   kappa * [ d/dx( h dq/dx)  + d/dy( h dq/dy) ]  and returns in QD */


/* void passive_newtonian1(field2D *QD, const field2D *QQ, const double *H, const double *HBAR, double kappa) */
/* { */

/*   int i,j,im,M,MN; */
/*   double tmp0, tmp1, tmp2, tmp3; */
/*   double *nvd; */
/*   double *qd, *qdup, *qddn; */
/*   double *qdd; */
/*   double him, hbim, hbim2; */
/*   double qq, qplus, qminus, qtmp0; */

/*   M = QD->M; */
/*   MN = M-1; */

/*   qdd = QD->data; */

/*   qd = &(QQ->data[M]); */
/*   qdup = &(QQ->data[2*M]); */
/*   qddn = &(QQ->data[M]); */

/*   /\* first do i=1 consider slip & no slip *\/ */

/*   i=1; */

/*   im = i-1; */
/*   hbim = HBAR[im]; */
/*   him = H[im]; */

/*   qminus = *(qd+MN); */
/*   qq = *qd++; */
/*   qplus = *qd++; */

/*   tmp1 = hbim * (*qdup++ - qq); */
/*   qtmp0 = 2.0*qq; */
/*   tmp2 = 0; */

/*   tmp3 = him*(qplus - qtmp0 + qminus); */
/*   *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */

/*   qminus = qq; */
/*   qq = qplus; */


/*   for (j=1;j<ny-1;j++) { */
/*     qplus = *qd++; */
/*     qtmp0 = 2.0*qq; */

/*     tmp1 = hbim * (*qdup++ - qq); */
/*     tmp2 = 0.0; */

/*     tmp3 = him*(qplus - qtmp0 + qminus); */
/*     *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */

/*     qminus = qq; */
/*     qq = qplus; */

/*   } */
/*   qplus = *(qd-M); */
/*   qtmp0 = 2.0*qq; */

/*   tmp1 = hbim * (*qdup++ - qq); */
/*   tmp2 = 0.0; */

/*   tmp3 = him*( qplus - qtmp0 + qminus); */
/*   *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */


/*   /\* away from BC do the interior *\/ */

/*   for (i=2;i<nx-1;i++) { */
/*     im = i-1; */
/*     hbim = HBAR[im]; */
/*     hbim2 = HBAR[im-1]; */
/*     him = H[im]; */


/*     qminus = *(qd+MN); */
/*     qq = *qd++; */
/*     qplus = *qd++; */

/*     tmp1 = hbim * (*qdup++ - qq); */
/*     qtmp0 = 2.0*qq; */
/*     tmp2 = hbim2 * (qq - *qddn++); */

/*     tmp3 = him*(qplus - qtmp0 + qminus); */
/*     *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */

/*     qminus = qq; */
/*     qq = qplus; */


/*     for (j=1;j<ny-1;j++) { */
/*       qplus = *qd++; */
/*       qtmp0 = 2.0*qq; */

/*       tmp1 = hbim * (*qdup++ - qq); */
/*       tmp2 = hbim2 * (qq - *qddn++); */

/*       tmp3 = him*(qplus - qtmp0 + qminus); */
/*       *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */

/*       qminus = qq; */
/*       qq = qplus; */

/*     } */
/*     qplus = *(qd-M); */
/*     qtmp0 = 2.0*qq; */

/*     tmp1 = hbim * (*qdup++ - qq); */
/*     tmp2 = hbim2 * (qq - *qddn++); */

/*     tmp3 = him*( qplus - qtmp0 + qminus); */
/*     *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */

/*   } */


/*   /\* last do i=nx-1  consider slip & no slip *\/ */



/*   i=nx-1; */
/*   im = i-1; */
/*   hbim2 = HBAR[im-1]; */
/*   him = H[im]; */

/*   qminus = *(qd+MN); */
/*   qq = *qd++; */
/*   qplus = *qd++; */
/*   qtmp0 = 2.0*qq; */

/*   tmp1 = 0.0; */
/*   tmp2 = hbim2 * ( qq - *qddn++); */

/*   tmp3 = him*(qplus - qtmp0 + qminus); */
/*   *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */

/*   qminus = qq; */
/*   qq = qplus; */

/*   for (j=1;j<ny-1;j++) { */
/*     qplus = *qd++; */
/*     qtmp0 = 2.0*qq; */

/*     tmp1 = 0.0; */
/*     tmp2 = hbim2 * ( qq - *qddn++); */

/*     tmp3 = him*(qplus - qtmp0 + qminus); */
/*     *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */

/*     qminus = qq; */
/*     qq = qplus; */

/*   } */

/*   qplus = *(qd-M); */
/*   qtmp0 = 2.0*qq; */

/*   tmp1 = 0.0; */
/*   tmp2 = hbim2 * ( qq - *qddn++); */


/*   tmp3 = him*( qplus - qtmp0 + qminus); */
/*   *qdd++ = kappa*(idx2*(tmp1-tmp2) + idy2*tmp3); */

/* } */

     


/* void laplacian_passive(field2D *LQ, field2D *QQ) */
/* { */

/*   int i,j,im,M,MN,N; */
/*   double qtmp0,qtmp1,qtmp2; */
/*   double *lqd, *qd, *qdup, *qddn; */

/*   double qminus, qq, qplus; */

/*   M = LQ->M; */
/*   N = QQ->N; */
/*   MN = M-1; */

/*   lqd = LQ->data; */

/*   qd = &(QQ->data[M]); */
/*   qdup = &(QQ->data[2*M]); */
/*   qddn = QQ->data; */

/*   for (i=1;i<N-1;i++) { */
/*     qminus = *(qd+MN); */
/*     qq = *qd++; */
/*     qplus = *qd++; */

/*     qtmp0 = qq * 2.0;  // here we do j=0 */
/*     qtmp1 =  *qdup++ - qtmp0 + *qddn++; */
/*     qtmp2 = qplus - qtmp0 + qminus; */
/*     *lqd++ = (idx2*qtmp1 + idy2*qtmp2); */

/*     qminus = qq; */
/*     qq = qplus; */

/*     for (j=1;j<ny-1;j++) { */
/*       qplus = *qd++; */
/*       qtmp0 =  qq * 2.0; */
/*       qtmp1 =  *qdup++ - qtmp0 + *qddn++; */
/*       qtmp2 =  qplus  - qtmp0 + qminus; */
/*       *lqd++ = (idx2*qtmp1 + idy2*qtmp2); */
/*       qminus = qq; */
/*       qq = qplus; */

/*     } */

/*     qplus = *(qd-M); */
/*     qtmp0 =  qq * 2.0; */
/*     qtmp1 =  *qdup++ - qtmp0 + *qddn++; */
/*     qtmp2 =  qplus  - qtmp0 + qminus; */
/*     *lqd++ = (idx2*qtmp1 + idy2*qtmp2); */
/*   } */

/* } */