ewSource.cpp 5.12 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187
#include <stdio.h>
#include <string.h>
#include <utilits.h>
#include <cOgrd.h>
#include "cOkadaEarthquake.h"
#include "easywave.h"

int Imin;
int Imax;
int Jmin;
int Jmax;

#define SRC_GRD 1
#define SRC_FLT 2

//====================================================
int ewSource()
{
  char dsaa_label[8];
  int i,j,ierr,srcType;
  double lon,lat,dz,absuzmax,absuzmin;
  FILE *fp;
  cOkadaEarthquake eq;
  cOgrd uZ;

  CNode& Node = *gNode;

  // check input file type: GRD or fault
  if( (fp = fopen( Par.fileSource, "rb" )) == NULL ) return Err.post( Err.msgOpenFile(Par.fileSource) );
  memset( dsaa_label, 0, 5 );
  ierr = fread( dsaa_label, 4, 1, fp );
  if( !strcmp( dsaa_label,"DSAA" ) || !strcmp( dsaa_label,"DSBB" ) )
    srcType = SRC_GRD;
  else
    srcType = SRC_FLT;
  fclose(fp);


  // load GRD file
  if( srcType == SRC_GRD) {
    ierr = uZ.readGRD( Par.fileSource ); if(ierr) return ierr;
  }

  // read fault(s) from file
  if( srcType == SRC_FLT) {
    int effSymSource = 0;
    long l;
    double dist,energy,factLat,effRad,effMax;
    
    ierr = eq.read( Par.fileSource ); if(ierr) return ierr;

    if( Par.adjustZtop ) {

		// check fault parameters
		Err.disable();
		ierr = eq.finalizeInput();
		while( ierr ) {
		  i = ierr/10;
		  ierr = ierr - 10*i;
		  if( ierr == FLT_ERR_STRIKE ) {
			Log.print( "No strike on input: Employing effective symmetric source model" );
			if( eq.nfault > 1 ) { Err.enable(); return Err.post("Symmetric source assumes only 1 fault"); }
			eq.fault[0].strike = 0.;
			effSymSource = 1;
		  }
		  else if( ierr == FLT_ERR_ZTOP ) {
			Log.print( "Automatic depth correction to fault top @ 10 km" );
			eq.fault[i].depth = eq.fault[i].width/2 * sindeg(eq.fault[i].dip) + 10.e3;
		  }
		  else {
			Err.enable();
			return ierr;
		  }
		  ierr = eq.finalizeInput();
		}
		Err.enable();

    } else {

		// check fault parameters
		Err.disable();
		ierr = eq.finalizeInput();
		if( ierr ) {
		  i = ierr/10;
		  ierr = ierr - 10*i;
		  if( ierr != FLT_ERR_STRIKE ) {
			Err.enable();
			ierr = eq.finalizeInput();
			return ierr;
		  }
		  Log.print( "No strike on input: Employing effective symmetric source model" );
		  Err.enable();
		  if( eq.nfault > 1 ) return Err.post("symmetric source assumes only 1 fault");
		  eq.fault[0].strike = 0.;
		  effSymSource = 1;
		  ierr = eq.finalizeInput(); if(ierr) return ierr;
		}
		Err.enable();

    }

    // calculate uplift on a rectangular grid
    // set grid resolution, grid dimensions will be set automatically
    uZ.dx = DLon; uZ.dy = DLat;
    ierr = eq.calculate( uZ ); if(ierr) return ierr;
    
    if( effSymSource ) {
      // integrate for tsunami energy
      energy = 0.;
      for( j=0; j<uZ.ny; j++ ) {
        factLat = Dx*cosdeg(uZ.getY(0,j))*Dy;
        for( i=0; i<uZ.nx; i++ )
          energy += pow(uZ(i,j),2.)*factLat;
      }
      energy *= (1000*9.81/2);
      effRad = eq.fault[0].length/sqrt(2*M_PI);
      effMax = 1./effRad / sqrt(M_PI/2) / sqrt(1000*9.81/2) * sqrt(energy);
      Log.print( "Effective source radius: %g km,  max height: %g m", effRad/1000, effMax );
  
      // transfer uplift onto tsunami grid and define deformed area for acceleration
      for( i=0; i<uZ.nx; i++ ) {
        for( j=0; j<uZ.ny; j++ ) {
          dist = GeoDistOnSphere( uZ.getX(i,j),uZ.getY(i,j), eq.fault[0].lon,eq.fault[0].lat ) * 1000;
          if( dist < effRad ) 
            uZ(i,j) = effMax*cos(M_PI/2*dist/effRad);
          else
            uZ(i,j) = 0.;
        }
      }
      
    } // effective source
  
  } // src_type == fault

  // remove noise in the source
  absuzmax = uZ.getMaxAbsVal();

  if( (Par.ssh0ThresholdRel + Par.ssh0ThresholdAbs) != 0 ) {

    absuzmin = RealMax;
    if( Par.ssh0ThresholdRel != 0 ) absuzmin = Par.ssh0ThresholdRel*absuzmax;
    if( Par.ssh0ThresholdAbs != 0 && Par.ssh0ThresholdAbs < absuzmin ) absuzmin = Par.ssh0ThresholdAbs;

    for( i=0; i<uZ.nx; i++ ) {
      for( j=0; j<uZ.ny; j++ ) {
        if( fabs(uZ(i,j)) < absuzmin ) uZ(i,j) = 0;
      }
    }

  }

  // calculated (if needed) arrival threshold (negative value means it is relative)
  if( Par.sshArrivalThreshold < 0 ) Par.sshArrivalThreshold = absuzmax * fabs(Par.sshArrivalThreshold);

  // transfer uplift onto tsunami grid and define deformed area for acceleration
  Imin = NLon; Imax = 1; Jmin = NLat; Jmax = 1;
  /* FIXME: change loops */
  for( i=1; i<=NLon; i++ ) {
    for( j=1; j<=NLat; j++ ) {

      lon = getLon(i);
      lat = getLat(j);

      if( Node(idx(j,i), iD) != 0. )
        dz = Node(idx(j,i), iH) = uZ.getVal( lon,lat );
      else
        dz = Node(idx(j,i), iH) = 0.;

      if( fabs(dz) > Par.sshClipThreshold ) {
        Imin = My_min( Imin, i );
        Imax = My_max( Imax, i );
        Jmin = My_min( Jmin, j );
        Jmax = My_max( Jmax, j );
      }

    }
  }

  if( Imin == NLon ) return Err.post( "Zero initial displacement" );

  Imin = My_max( Imin - 2, 2 );
  Imax = My_min( Imax + 2, NLon-1 );
  Jmin = My_max( Jmin - 2, 2 );
  Jmax = My_min( Jmax + 2, NLat-1 );

  return 0;
}