1 import sys; sys.path.append('../py_src')
   2 
   3 
   4 from pyOM_gui import pyOM_gui as pyOM
   5 #from pyOM_cdf import pyOM_cdf as pyOM
   6 from numpy import *
   7 
   8 fac = 1.0
   9 mix = 2.5e-3
  10 
  11 
  12 class kelv1(pyOM):
  13    """ 
  14    """
  15    def set_parameter(self):
  16      """set main parameter
  17      """
  18      M=self.fortran.main_module   
  19 
  20      (M.nx,M.ny,M.nz)=(int(1.5*64*fac),  1, int(40*fac) )
  21      M.dt_mom     = 0.04/fac
  22      M.dt_tracer  = 0.04/fac
  23 
  24      M.enable_conserve_energy = 0
  25      M.coord_degree           = 0
  26      M.enable_cyclic_x        = 1
  27      M.enable_hydrostatic     = 0
  28      M.eq_of_state_type       = 1 
  29 
  30      M.congr_epsilon = 1e-12
  31      M.congr_max_iterations = 5000
  32      M.congr_epsilon_non_hydro=   1e-6
  33      M.congr_max_itts_non_hydro = 5000    
  34 
  35      M.enable_explicit_vert_friction = 1
  36      M.kappam_0 = mix/fac**2
  37      M.enable_hor_friction = 1
  38      M.a_h = mix/fac**2
  39 
  40      M.enable_tempsalt_sources = 1
  41      M.enable_momentum_sources = 1
  42      M.enable_superbee_advection = 1
  43      return
  44 
  45 
  46    def set_grid(self):
  47      M=self.fortran.main_module   
  48      M.dxt[:]=0.25/fac 
  49      M.dyt[:]=0.25/fac 
  50      M.dzt[:]=0.25/fac 
  51      return
  52 
  53    def t_star_fct(self,k):
  54      M=self.fortran.main_module
  55      return 9.85-6.5*tanh( (M.zt[k-1]-M.zt[M.nz/2-1] ) /M.zt[0]*100 )
  56 
  57    def u_star_fct(self,k):
  58      M=self.fortran.main_module   
  59      return 0.6+0.5*tanh( (M.zt[k-1]-M.zt[M.nz/2-1])/M.zt[0]*100)
  60 
  61          
  62    def set_initial_conditions(self):
  63      """ setup initial conditions
  64      """
  65      M=self.fortran.main_module   
  66 
  67      # target for restoring
  68      self.t_rest = zeros( M.u[:,:,:,0].shape, order = 'F' )
  69      self.t_star = zeros( M.u[:,:,:,0].shape, order = 'F' )
  70      self.u_star = zeros( M.u[:,:,:,0].shape, order = 'F' )
  71      
  72      for k in range(M.nz):
  73       self.t_star[:,:,k]=self.t_star_fct(k+1)
  74       self.u_star[:,:,k]=self.u_star_fct(k+1)
  75       for i in range(M.is_pe,M.ie_pe+1):
  76         if i < M.nx/8: self.t_rest[self.if2py(i),:,k] = 1./(15*M.dt_mom)
  77 
  78      # initial conditions
  79      from numpy.random import randn
  80      for k in range(M.nz):
  81        M.u[:,:,k,:]   = self.u_star_fct(k+1)
  82        for i in range(M.is_pe,M.ie_pe+1):
  83         ii = self.if2py(i)
  84         fxa=1e-3*M.zt[0]*sin(M.xt[ii]/(M.nx*M.dxt[1])*16*pi)*sin(M.zt[k]/( M.zt[0]-M.dzt[0]/2)*pi)
  85         M.temp[ii,:,k,:]= fxa+self.t_star_fct(k+1) 
  86      return
  87    
  88 
  89    def set_forcing(self):
  90      M=self.fortran.main_module   
  91      if M.enable_tempsalt_sources: M.temp_source[:]=self.t_rest*(self.t_star-M.temp[:,:,:,M.tau-1])#*M.maskt[:]
  92      if M.enable_momentum_sources: M.u_source[:]   =self.t_rest*(self.u_star-M.u[:,:,:,M.tau-1]   )#*M.masku[:]
  93      return
  94 
  95 
  96    def user_defined_signal(self):
  97        """ this routine must be called by all processors
  98        """
  99        M=self.fortran.main_module  
 100        a = zeros( (M.nx,M.ny), 'd', order = 'F')
 101        a[M.is_pe-1:M.ie_pe,0] = M.xt[2:-2]
 102        self.fortran.pe0_recv_2d(a)
 103        self.xt_gl = a[:,0].copy()
 104        
 105        self.temp_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 106        self.u_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 107        self.w_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 108        for k in range(M.nz):
 109          a[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] = where( M.maskt[2:-2,2:-2,k] >0,  M.temp[2:-2,2:-2,k,M.tau-1] , NaN) 
 110          self.fortran.pe0_recv_2d(a)
 111          self.temp_gl[:,:,k]=a.copy()
 112 
 113          a[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] = where( M.masku[2:-2,2:-2,k] >0,  M.u[2:-2,2:-2,k,M.tau-1] , NaN) 
 114          self.fortran.pe0_recv_2d(a)
 115          self.u_gl[:,:,k]=a.copy()
 116 
 117          a[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] = where( M.maskw[2:-2,2:-2,k] >0,  M.w[2:-2,2:-2,k,M.tau-1] , NaN) 
 118          self.fortran.pe0_recv_2d(a)
 119          self.w_gl[:,:,k]=a.copy()
 120 
 121        return
 122 
 123    def make_plot(self):
 124        M=self.fortran.main_module
 125        
 126        self.set_signal('user_defined') # following routine is called by all PEs
 127        self.user_defined_signal()
 128        
 129        self.figure.clf()
 130        ax=self.figure.add_subplot(211)
 131        
 132        co=ax.contourf(self.xt_gl,M.zt, self.temp_gl[:,0,:].transpose())
 133        self.figure.colorbar(co)
 134        ax.quiver(self.xt_gl[::2],M.zt[::2],self.u_gl[::2,0,::2].transpose(),self.w_gl[::2,0,::2].transpose() )
 135        ax.set_title('Temperature [deg C]')
 136        ax.set_xlabel('x [m]')
 137        ax.set_ylabel('z [m]')
 138        #ax.axis('tight')
 139        return
 140 
 141 
 142 if __name__ == "__main__":
 143       m=kelv1()
 144       dt=m.fortran.main_module.dt_tracer
 145       m.run( snapint = 2.0 ,runlen = 50.0)