forcing.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
 *
 */



/* forcing.c   Falk Feddersen */       


#include <stdlib.h>
#include <math.h>

#include "forcing.h"
#include "depth.h"
#include "bottom_stress.h"


field2D *FXX, *FYY;


/* Function to load the forcomg file.
   The flag is set to 0 if we are loading Fx and
   set to 1 if we are loading Fy      */


void load_full_forcing(char *fname,field2D *FF)
{
  int val;
  val = load_field2D_ascii(fname,FF);
  if (val != 1)
    funwaveC_error_message("Bad when loading the full forcing file!\n");

}



void load_line_forcing(char *fname,field2D *FF)
{

   FILE *fp;
   int i,j;
   gdouble a;

   if ((fp=fopen(fname,"r"))==NULL) {
       fprintf(stderr,"Error loading line forcing file: %s\n",fname);
       exit(-1);
   }

   for (i=0;i<FF->N;i++) {
      fscanf(fp,"%lf",&a);
      for (j=0;j<FF->M;j++)
          DR2(FF,i,j) = a;
   }
   fclose(fp);
}



void set_const_forc(gdouble f,field2D *FF)
{  

  field2D_set_to_value(FF,f);

}


/* This peturbs the forcing as described in ANH96 */

void peturb_forcing(field2D *FY, gdouble e, gdouble b,int N)
{

   int i,j,k,ny;
   gdouble *by;
   gdouble phase;

   ny = FY->M;

   by = (gdouble * ) g_malloc(sizeof(gdouble)*ny);
   
   if (by==NULL) 
     funwaveC_error_message("Bad memory allocation in peturb_forcing()\n");

   for (j=0;j<ny;j++) {
     by[j] = 0.0;
     for (k=1;k<=N;k++) {
         phase= 2.0*pi * (double ) (rand()/(RAND_MAX+1.0));
         by[j] += cos(2*pi*k*j/ny - phase);
     }
     by[j] *= e;
   }

   for (j=0;j<FY->M;j++)
     for (i=0;i<FY->N;i++)
        DR2(FY,i,j) *= (1+by[j]);

   g_free(by);
}


void balance_forcing()
{

  field2D *TMP_ETABX;
  field2D *TMP_ETABY;

  TMP_ETABX = allocate_field2D(nx,ny);
  TMP_ETABY = allocate_field2D(nx-1,ny);

  field2D_set_to_value(TMP_ETABX, 0.0);
  field2D_set_to_value(TMP_ETABY, 0.0);

  rhs_uvmom_bottom_stress(FXX, FYY, UU, VV, Depth_Vars->HBX, Depth_Vars->HBY, TMP_ETABX, TMP_ETABY);
  field2D_set_to_value(FXX, 0.0);

  deallocate_field2D(TMP_ETABX);
  deallocate_field2D(TMP_ETABY);

}


void  init_forcing(forcing_info_t  *F, int nx, int ny)
{

  double eee;

  FXX = allocate_field2D(nx-2, ny);
  FYY = allocate_field2D(nx-1, ny);


   /* load the forcing quantities */

  switch(F->dfx) {
   case 0:   set_const_forc(atof(F->fxfile->str),FXX);
             break;       // this is the zero quantity 
   case 1:   load_line_forcing(F->fxfile->str,FXX);
             break;
   case 2:   load_full_forcing(F->fxfile->str,FXX);
             break;
   default:  funwaveC_error_message("Error for x-forcing term in init file");
   }

   switch(F->dfy) {
   case 0:   set_const_forc(atof(F->fyfile->str),FYY);
             break;     // this is for zero forcing 
   case 1:   load_line_forcing(F->fyfile->str,FYY);
             break;
   case 2:   load_full_forcing(F->fyfile->str,FYY);
             break;
   case 3:   balance_forcing();   // this function sits in bottom_stress.c 
             break;
   case 4:   balance_forcing();  
   //             sscanf(sfy,"%*s  %lf %lf %d",&e,&b1,&N1);    
     //     save_field2D_ascii("Fyy0.dat",FYY);
             eee = 0.01;
             printf("*** Peturbing the forcing!!  e=%g\n",eee);
             peturb_forcing(FYY,eee,1,12);              // need to fix this!!!
	 //     save_field2D_ascii("Fyyp.dat",FYY);
             break;
   default:  funwaveC_error_message("Error for y-forcing term in init file");
   }


}

void  deallocate_forcing()
{

  deallocate_field2D(FXX);
  deallocate_field2D(FYY);

}



void  rhs_uvmom_forcing_add(field2D *UT, field2D *VT)
{

  int M = UT->M;
  int numFX = FXX->num;
  int numFY = FYY->num;
  int i;

  const double *fxxd = FXX->data;
  const double *fyyd = FYY->data;

  double *utd = &(UT->data[M]);
  double *vtd = &(VT->data[M]);

  for (i=0;i<numFX;i++) 
    *utd++ += *fxxd++;

  for (i=0;i<numFY;i++) 
    *vtd++ += *fyyd++;


}