1 import sys; sys.path.append('../py_src')
   2 
   3 from pyOM_gui import pyOM_gui as pyOM
   4 #from pyOM_cdf import pyOM_cdf as pyOM
   5 from numpy import *
   6 
   7 
   8 HRESOLVE = 0.5
   9 VRESOLVE = 0.5
  10 N_0     = 0.004
  11 M_0     = sqrt(1e-5*0.1/1024.0*9.801)
  12 spg_width = int(3*HRESOLVE)
  13 t_rest=1./(5.*86400)
  14 
  15 #DX  = 32094.1729769
  16 DX  = 30e3
  17 Lx  = DX*128
  18 H   = 1800.0
  19 
  20 class jets1(pyOM):
  21    """ a wide channel model with zonal jets
  22    """
  23    def set_parameter(self):
  24      """set main parameter
  25      """
  26      M=self.fortran.main_module   
  27      (M.nx,M.ny,M.nz)  = ( 128*HRESOLVE , 128*HRESOLVE , 18*VRESOLVE )
  28      M.dt_tracer  = 1800.0/HRESOLVE
  29      M.dt_mom     = 1800.0/HRESOLVE
  30 
  31      M.enable_conserve_energy = 0
  32      M.coord_degree           = 0
  33      M.enable_cyclic_x        = 1
  34      M.enable_hydrostatic     = 1
  35      M.eq_of_state_type       = 1
  36      
  37      M.congr_epsilon = 1e-12
  38      M.congr_max_iterations = 5000
  39      M.enable_streamfunction = 1
  40      
  41      M.kappah_0=1e-4/VRESOLVE**2
  42      M.kappam_0=1e-4/VRESOLVE**2
  43      
  44      M.enable_ray_friction      = 1
  45      M.r_ray = 1e-7 
  46      
  47      #M.enable_hor_friction = 1
  48      #M.a_h = 100/HRESOLVE**2  
  49      M.enable_biharmonic_friction  = 1
  50      M.a_hbi  = 5e11/HRESOLVE**2
  51      
  52      M.enable_superbee_advection = 1
  53      M.enable_tempsalt_sources = 1
  54      return
  55  
  56    def set_grid(self):
  57      M=self.fortran.main_module
  58      M.dxt[:]  = Lx/M.nx
  59      M.dyt[:]  = Lx/M.ny
  60      M.dzt[:]  = H/M.nz
  61      return
  62 
  63    def set_coriolis(self):
  64      M=self.fortran.main_module   
  65      phi0 = 10. /180. *pi
  66      betaloc = 2*M.omega*cos(phi0)/M.radius
  67      for j in range( M.yt.shape[0] ):
  68        M.coriolis_t[:,j] = 2*M.omega*sin(phi0)+betaloc*M.yt[j]
  69      return
  70 
  71    def set_topography(self):
  72      """ setup topography
  73      """
  74      M=self.fortran.main_module   
  75      M.kbot[:]=1
  76      return
  77 
  78    def set_initial_conditions(self):
  79      """ setup all initial conditions
  80      """
  81      M=self.fortran.main_module   
  82      alpha = self.fortran.linear_eq_of_state.linear_eq_of_state_drhodt() # use here general eq. of state !!!!!!!!!!
  83      self.T_star = zeros( M.temp[:,:,:,0].shape, order = 'F' )
  84      self.T_rest = zeros( M.temp[:,:,:,0].shape, order = 'F' )
  85 
  86      for j in range(M.yt.shape[0]):
  87        for k in range(M.nz):
  88          fxa = 0.5e-3*sin(M.xt*8.5/Lx*pi) *1024./9.81 /alpha  # drho = drho/dt  dt
  89          self.T_star[:,j,k]  = (  32+ (M_0**2*M.yt[j] - N_0**2*M.zt[k])*1024./9.81 /alpha )*M.maskt[:,j,k] 
  90          M.temp[:,j,k,M.tau-1]     = (fxa+ self.T_star[:,j,k])*M.maskt[:,j,k]
  91      M.temp[:,:,:,M.taum1-1] = M.temp[:,:,:,M.tau-1]
  92        
  93      for j in range(1,spg_width+1):
  94        if j<=M.je_pe and j>= M.js_pe:
  95          jj = self.jf2py(j)
  96          self.T_rest[:,jj,:]  =  t_rest/j*M.maskt[:,jj,:]
  97      for j in range(M.ny-1,M.ny-spg_width-1,-1):
  98        if j<=M.je_pe and j>= M.js_pe:
  99          jj = self.jf2py(j)
 100          self.T_rest[:,jj+1,:]  =  t_rest/(M.ny-j) *M.maskt[:,jj+1,:]
 101      return
 102   
 103 
 104    def set_forcing(self):
 105      M=self.fortran.main_module   
 106      if M.enable_tempsalt_sources: M.temp_source[:]=self.T_rest*(self.T_star-M.temp[:,:,:,M.tau-1])*M.maskt
 107      return
 108 
 109    def user_defined_signal(self):
 110        """ this routine must be called by all processors
 111        """
 112        M=self.fortran.main_module  
 113        a = zeros( (M.nx,M.ny), 'd', order = 'F')
 114        a[M.is_pe-1:M.ie_pe,0] = M.xt[2:-2]
 115        self.fortran.pe0_recv_2d(a)
 116        self.xt_gl = a[:,0].copy()
 117        
 118        a[0,M.js_pe-1:M.je_pe] = M.yt[2:-2]
 119        self.fortran.pe0_recv_2d(a)
 120        self.yt_gl = a[0,:].copy()
 121        
 122        self.psi_gl = zeros( (M.nx,M.ny), 'd', order = 'F')
 123        self.psi_gl[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] = where( M.maskz[2:-2,2:-2,-1] >0,  M.psi[2:-2,2:-2,M.tau-1] , NaN) 
 124        self.fortran.pe0_recv_2d(self.psi_gl)
 125        
 126        self.temp_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 127        for k in range(M.nz):
 128          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) 
 129          self.fortran.pe0_recv_2d(a)
 130          self.temp_gl[:,:,k]=a.copy()
 131 
 132        self.u_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 133        for k in range(M.nz):
 134          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) 
 135          self.fortran.pe0_recv_2d(a)
 136          self.u_gl[:,:,k]=a.copy()
 137         
 138        self.v_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 139        for k in range(M.nz):
 140          a[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] = where( M.maskv[2:-2,2:-2,k] >0,  M.v[2:-2,2:-2,k,M.tau-1] , NaN) 
 141          self.fortran.pe0_recv_2d(a)
 142          self.v_gl[:,:,k]=a.copy()
 143         
 144        return
 145    
 146 
 147    def make_plot(self):
 148      """ diagnose the model variables, could be replaced by other version
 149      """
 150      M=self.fortran.main_module         # fortran module with model variables
 151      self.set_signal('user_defined') # following routine is called by all PEs
 152      self.user_defined_signal()
 153      k=M.nz*3/4
 154      i=int(M.nx/2)
 155      self.figure.clf()
 156      
 157      ax=self.figure.add_subplot(221)
 158      co=ax.contourf(self.xt_gl/1e3,self.yt_gl/1e3,self.temp_gl[:,:,k].transpose())
 159      self.figure.colorbar(co)
 160      #ax.set_title('temperature')
 161      #ax.set_ylabel('y [km]')
 162      a=self.u_gl[::2,::2,k] 
 163      b=self.v_gl[::2,::2,k] 
 164      ax.quiver(self.xt_gl[::2]/1e3,self.yt_gl[::2]/1e3,a.transpose(),b.transpose() )
 165      
 166      ax=self.figure.add_subplot(222)
 167      co=ax.contourf(self.yt_gl/1e3,M.zt,self.temp_gl[i,:,:].transpose())
 168      self.figure.colorbar(co)
 169      #ax.set_title('temperature')
 170      
 171      ax=self.figure.add_subplot(223)
 172      co=ax.contourf(self.xt_gl/1e3,self.yt_gl/1e3,self.psi_gl.transpose()/1e6)
 173      self.figure.colorbar(co)
 174      #ax.set_title('[Sv]')
 175      #ax.set_ylabel('y [km]')
 176      
 177      ax=self.figure.add_subplot(224)
 178      co=ax.contourf(self.yt_gl/1e3,M.zt,self.u_gl[i,:,:].transpose())
 179      self.figure.colorbar(co)
 180      #ax.set_title('temperature')
 181      return
 182 
 183   
 184 if __name__ == "__main__":
 185     model=jets1()
 186     model.run(snapint=3*86400.0,runlen=365*86400.)