1 import sys; sys.path.append('../py_src')
   2 from pyOM_gui import pyOM_gui as pyOM
   3 from numpy import *
   4 
   5 fac=1
   6 N_0 = 2*pi/10.
   7 OM0 = 1./(1.5*10)
   8 
   9 class internal_wave1(pyOM):
  10    
  11    """ internal wave maker
  12    """
  13    
  14    def set_parameter(self):
  15      """set main parameter
  16      """
  17      M=self.fortran.main_module
  18      (M.nx,M.ny,M.nz) = (64*fac,1,64*fac)
  19      
  20      M.dt_mom   =0.025/fac
  21      M.dt_tracer=0.025/fac
  22      
  23      M.enable_conserve_energy = 0
  24      M.coord_degree           = 0
  25      M.enable_cyclic_x        = 1
  26      M.enable_hydrostatic     = 0
  27      M.eq_of_state_type       = 1
  28 
  29      M.congr_epsilon = 1e-12
  30      M.congr_max_iterations = 5000
  31      M.congr_epsilon_non_hydro=   1e-7
  32      M.congr_max_itts_non_hydro = 5000    
  33 
  34      M.enable_explicit_vert_friction = 1
  35      M.kappam_0 = 5e-3/fac**2   
  36      M.enable_hor_friction = 1
  37      M.a_h      = 5e-3/fac**2
  38      M.enable_superbee_advection = 1
  39 
  40      M.enable_tempsalt_sources = 1
  41      M.enable_momentum_sources = 1
  42      return
  43 
  44    def set_grid(self):
  45        M=self.fortran.main_module   
  46        M.dxt[:]= 0.25/fac
  47        M.dyt[:]= 0.25/fac
  48        M.dzt[:]= 0.25/fac
  49        return
  50 
  51    def set_topography(self):
  52        M=self.fortran.main_module   
  53        M.kbot[:]=0
  54        M.kbot[:,2:-2]=1
  55        return
  56 
  57 
  58    def set_initial_conditions(self):
  59      """ setup all initial conditions
  60      """
  61      M=self.fortran.main_module
  62      alpha = self.fortran.linear_eq_of_state.linear_eq_of_state_drhodt()
  63      grav = 9.81; rho0 = 1024.0
  64 
  65      self.t0  = zeros( (M.i_blk+2*M.onx,M.j_blk+2*M.onx,M.nz) , 'd', order='F')
  66      self.dt0 = zeros( (M.i_blk+2*M.onx,M.j_blk+2*M.onx,M.nz,3) , 'd', order='F')
  67 
  68      # prescribed background stratification
  69      for k in range(M.nz):
  70         self.t0[:,:,k]=-N_0**2*M.zt[k]/grav/alpha*rho0*M.maskt[:,:,k]
  71      return
  72 
  73    def set_forcing(self):
  74       M=self.fortran.main_module   
  75       # implement effect of background state 
  76       # update density, etc of last time step
  77       M.temp[:,:,:,M.tau-1] = M.temp[:,:,:,M.tau-1] + self.t0
  78       self.fortran.calc_eq_of_state(M.tau)
  79       M.temp[:,:,:,M.tau-1] = M.temp[:,:,:,M.tau-1] - self.t0
  80       
  81       # advection of background temperature
  82       self.fortran.advect_tracer(M.is_pe-M.onx,M.ie_pe+M.onx,M.js_pe-M.onx,M.je_pe+M.onx,self.t0,self.dt0[...,M.tau-1],M.nz) 
  83       M.temp_source[:] = (1.5+M.ab_eps)*self.dt0[...,M.tau-1] - ( 0.5+M.ab_eps)*self.dt0[...,M.taum1-1]
  84 
  85       i=M.nx/2
  86       if i >= M.is_pe and i<= M.ie_pe:
  87          ii = self.if2py(i); k=M.nz/2-1
  88          jj= self.jf2py(1)
  89          M.u_source[ii,jj,k]= M.masku[ii,jj,k]*1./(100*60.*M.dt_tracer)*sin(2*pi*OM0*M.itt*M.dt_tracer)
  90       return
  91 
  92   
  93    def make_plot(self):
  94        M=self.fortran.main_module   
  95        self.figure.clf()
  96        ax=self.figure.add_subplot(111)
  97        x=M.xt[2:-2];z=M.zt[:];
  98        jj= self.jf2py(1); ii = self.if2py(M.nx/2); k=M.nz/2-1
  99        a= M.temp[:,jj,:,M.tau-1]
 100        a[ii-2:ii+3,k-2:k+3]=0
 101 
 102        co=ax.contourf(x,z,a[2:-2,:].transpose() )
 103        self.figure.colorbar(co)
 104        a = a+self.t0[:,jj,:]
 105        ax.contour(x,z,a[2:-2,:].transpose(),10,colors='red')
 106 
 107        a= M.u[:,jj,:,M.tau-1]
 108        b= M.w[:,jj,:,M.tau-1]
 109        a[ii-3:ii+4,k-3:k+4]=0
 110        b[ii-3:ii+4,k-3:k+4]=0
 111        ax.quiver(x[::2],z[::2],a[2:-2:2,::2].transpose(),b[2:-2:2,::2].transpose() )
 112  
 113        ax.set_title('Temperature')
 114        ax.set_xlabel('x [m]')
 115        ax.set_ylabel('z [m]')
 116        ax.axis('tight')
 117        return
 118   
 119 if __name__ == "__main__": 
 120      internal_wave1().run(snapint = 1.0 ,runlen = 50.0)