bottom_stress.c

波浪数值模拟

/*
 * Copyright (c) 2001-2005 Falk Feddersen
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */



/* bottom_stress.c   Falk Feddersen */       

#include <math.h>
#include "bottom_stress.h"

double c_d;


void  init_bottom_stress(double cd)  //bottom_stressinfo_t *S, int nx, int ny)
{
  c_d =  -1.0 * cd;   //0.003;   /* need the minus sign here */

}


void  deallocate_bottom_stress()
{

}


/*  TAUBX, TAUBY =  -(1/h) * cd * sqrt(V^2+U^2)*(U,V)  - note the normalization by depth and minus sign 

   TAUBX has dimensions  (NX-2,M), &  TAUBY = (NX-1,M)  
   U = (NX,M),   V(NX+1,M)
   HBX = (NX,M),   HBY(NX+1,M)
   ETABX = (NX,M),   ETABY(NX-1,M)

*/

void  rhs_uvmom_bottom_stress(field2D *TAUBX, field2D *TAUBY, const field2D *U, const field2D *V, const field2D *HBX, const field2D *HBY, 
                              const field2D *ETABX, const field2D *ETABY)
{
  const int M = TAUBX->M;
  const int MN = M-1;
  int NTX = TAUBX->N;
  //  int NTY = TAUBY->N;

  int i,j,jp;
  double Uavg, Vavg;
  double speedu, speedv;

  const double *ud; //, *udm;     //, *udp, *udmp;
  const double *vd; //, *vdp;     // *vddn;
  double Uij;  //, Uimj, Uijp, Uimjp;
  double Vij;  //, Vipj, Vijm, Vipjm;
  double depthx, depthy;

  const double *hbxd = &(HBX->data[M]);
  const double *hbyd = &(HBY->data[M]);      // because HBY has dimensions (NX+1,M)
  const double *etabxd = &(ETABX->data[M]);  // bring in one index in i     
  const double *etabyd = ETABY->data;        // this is ok

  double *taubxd = TAUBX->data;
  double *taubyd = TAUBY->data;


  ud = &(U->data[M]);
  //  udm = &(U->data[0]);

  vd = &(V->data[M]);
  //  vdp = &(V->data[2*M]);


  for (i=1;i<=NTX;i++) {
    //    j=0;
    //    jp=1;

    Uij = *ud;
    Vij = *vd;
    Uavg = 0.25*( *ud +  *(ud-M) + *(ud+1) + *(ud-M+1)) ; // Uij+ Uimj + Uijp + Uimjp);  // (Uij+ Uimj + Uijp + Uimjp);
    Vavg = 0.25*( *vd + *(vd+M) + *(vd+MN) + *(vd+MN) );  //(Vij+ Vipj +DR2(VV,i,jm)+DR2(VV,i+1,jm));  // (Vij+ Vipj +DR2(VV,i,jm)+DR2(VV,i+1,jm));
    speedu = sqrt(SQR(Vavg) + SQR(Uij));
    speedv = sqrt(SQR(Uavg) + SQR(Vij));
    depthx = ( *hbxd++ + *etabxd++);
    depthy = ( *hbyd++ + *etabyd++);
    *taubxd++  = c_d*speedu*Uij/ (depthx);           /* this is the pointer version */
    *taubyd++  = c_d*speedv*Vij/ (depthy);

    ud++;
    vd++;

    for (j=1;j<M-1;j++){

      Uij = *ud;
      Vij = *vd;
      Uavg = 0.25*( *ud +  *(ud-M) + *(ud+1) + *(ud-M+1)) ; // Uij+ Uimj + Uijp + Uimjp);  // (Uij+ Uimj + Uijp + Uimjp);
      Vavg = 0.25*( *vd + *(vd+M) + *(vd-1) + *(vd+M-1) );  //(Vij+ Vipj +DR2(VV,i,jm)+DR2(VV,i+1,jm));  // (Vij+ Vipj +DR2(VV,i,jm)+DR2(VV,i+1,jm));
      speedu = sqrt(SQR(Vavg) + SQR(Uij));
      speedv = sqrt(SQR(Uavg) + SQR(Vij));
      depthx = ( *hbxd++ + *etabxd++);
      depthy = ( *hbyd++ + *etabyd++);
      *taubxd++  = c_d*speedu*Uij/ (depthx);           /* this is the pointer version */
      *taubyd++  = c_d*speedv*Vij/ (depthy);

      ud++;
      vd++;

    }

    Uij = *ud;
    Vij = *vd;
    Uavg = 0.25*( *ud +  *(ud-M) + *(ud-MN) + *(ud-M-MN)) ; // Uij+ Uimj + Uijp + Uimjp);  // (Uij+ Uimj + Uijp + Uimjp);
    Vavg = 0.25*( *vd + *(vd+M) + *(vd-1) + *(vd+M-1) );  //(Vij+ Vipj +DR2(VV,i,jm)+DR2(VV,i+1,jm));  // (Vij+ Vipj +DR2(VV,i,jm)+DR2(VV,i+1,jm));
    speedu = sqrt(SQR(Vavg) + SQR(Uij));
    speedv = sqrt(SQR(Uavg) + SQR(Vij));
    depthx = ( *hbxd++ + *etabxd++);
    depthy = ( *hbyd++ + *etabyd++);
    *taubxd++  = c_d*speedu*Uij/ (depthx);           /* this is the pointer version */
    *taubyd++  = c_d*speedv*Vij/ (depthy);

    ud++;
    vd++;

  }

  i = NTX+1;    /* This is at the offshore boundary */

  for (j=0;j<M-1;j++){
    Vij = *vd;
    Uavg = 0.25*( *ud +  *(ud-M) + *(ud+1) + *(ud-M+1)) ; // Uij+ Uimj + Uijp + Uimjp);  // (Uij+ Uimj + Uijp + Uimjp);
    speedv = sqrt(SQR(Uavg) + SQR(Vij));
    depthy = ( *hbyd++ + *etabyd++);
    *taubyd++ = c_d*speedv*Vij/ (depthy);
    ud++;
    vd++;
  }
  j=M-1;
  jp=0;

  Vij = *vd;
  Uavg = 0.25*( *ud +  *(ud-M) + *(ud-MN) + *(ud-M-MN)) ; // Uij+ Uimj + Uijp + Uimjp);  // (Uij+ Uimj + Uijp + Uimjp);
  speedv = sqrt(SQR(Uavg) + SQR(Vij));
  depthy = ( *hbyd++ + *etabyd++);
  *taubyd++ = c_d*speedv*Vij/(depthy);



}





/* report on bottom stress stability parameters by doing a linearlized rayliegh friction analysis */

void  bottom_stress_stability_report( double Umax, double min_h0, double dt )
{

    // find the largest S_u, find smallest h  
 
  double Su;
 
  Su = fabs(c_d * Umax *dt/min_h0);     // < 0.96 (AB2)

  fprintf(stderr,"Su # = %g  : ",Su);
   if (Su>0.95)
      fprintf(stderr,"Warning: model may be Rayleigh Friction unstable\n");
   else
      fprintf(stderr,"Rayleigh friction OK\n");
  
    
}