1 import sys; sys.path.append('../py_src')
   2 
   3 from pyOM_gui import pyOM_gui as pyOM
   4 #from pyOM_ave import pyOM_ave as pyOM
   5 
   6 from numpy import *
   7 
   8 hresol = 0.25
   9 vresol = 1.0
  10 
  11 class acc1(pyOM):
  12    """  idealised Southern Ocean, similar to Viebahn and Eden (2010) Ocean modeling
  13    """
  14    def set_parameter(self):
  15      """set main parameter
  16      """
  17      M=self.fortran.main_module   
  18 
  19      (M.nx,M.ny,M.nz)    = ( int(128*hresol),  int(128*hresol), int(18*vresol) )
  20      M.dt_mom    = 1200./hresol
  21      M.dt_tracer = 1200.0/hresol*5
  22 
  23 
  24      M.enable_conserve_energy = 0
  25      M.coord_degree           = 0
  26      M.enable_cyclic_x        = 1
  27      M.enable_hydrostatic     = 1
  28      M.eq_of_state_type       = 1
  29      
  30      M.congr_epsilon = 1e-6
  31      M.congr_max_iterations = 5000
  32      M.enable_streamfunction = 1
  33      
  34      M.enable_superbee_advection = 1
  35 
  36      M.enable_implicit_vert_friction = 1
  37      I=self.fortran.isoneutral_module   
  38      I.enable_TEM_friction = 1
  39      I.k_gm_0 = 1000.0
  40 
  41      M.enable_hor_friction  = 1
  42      M.a_h = 5e4 
  43 
  44      M.enable_biharmonic_friction  = 0 # for eddy resolving version
  45      M.a_hbi  = 5e11/hresol**2
  46 
  47      M.enable_bottom_friction = 1
  48      M.r_bot = 1e-5*vresol
  49 
  50      M.kappah_0=1.e-4/vresol
  51      M.kappam_0=1.e-3/vresol
  52      #M.a_h = (1./hresol)**2*5e4     # 1/T = A_1 /dx^2 = A_2 /(f dx)^2
  53      return
  54 
  55 
  56    def set_grid(self):
  57      M=self.fortran.main_module
  58      M.dxt[:] = 20e3/hresol
  59      M.dyt[:] = 20e3/hresol
  60      M.dzt[:] = 50.0/vresol
  61      return
  62 
  63 
  64    def set_coriolis(self):
  65      M=self.fortran.main_module   
  66      phi0 =- 25.0 /180. *pi
  67      betaloc = 2*M.omega*cos(phi0)/M.radius
  68      for j in range( M.yt.shape[0] ): 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      L_y = 0.0
  76      if M.my_blk_j == M.n_pes_j: L_y = M.yu[-(1+M.onx)]
  77      self.fortran.global_max(L_y)
  78      L_x = 0.0
  79      if M.my_blk_i == M.n_pes_i: L_x = M.xu[-(1+M.onx)]
  80      self.fortran.global_max(L_x)
  81      if M.my_pe==0: print ' domain size is ',L_x,' m x ',L_y,' m'
  82      self.L_x = L_x; self.L_y = L_y    
  83      (X,Y)= meshgrid(M.xt,M.yt); X=X.transpose(); Y=Y.transpose()
  84      M.kbot[:]=1
  85      M.kbot[ logical_and( Y>L_y*0.5 , logical_or(  X>L_x*0.75 ,  X<L_x*0.25 )  ) ]=0
  86      return
  87 
  88    def set_initial_conditions(self):
  89      M=self.fortran.main_module
  90      alpha = self.fortran.linear_eq_of_state.linear_eq_of_state_drhodt()
  91      (XT,YT)= meshgrid(M.xt,M.yt); XT=XT.transpose(); YT=YT.transpose()
  92      (XU,YU)= meshgrid(M.xu,M.yu); XU=XU.transpose(); YU=YU.transpose()
  93 
  94      n=YT < self.L_y*0.5
  95      M.surface_taux[n] = .1e-3*sin(2*pi*YU[n] /self.L_y)
  96      self.t_rest = M.dzt[-1]/(30.*86400)
  97      db=-30e-3 *M.rho_0/M.grav/alpha
  98      self.t_star = YT*0 + db
  99      n=YT < self.L_y*0.5
 100      self.t_star[n] = db*YT[n]/(self.L_y/2.0)
 101      n=YT > self.L_y*0.75
 102      self.t_star[n] = db*(1-(YT[n]-self.L_y*0.75)/(self.L_y*0.25) )
 103      return
 104 
 105  
 106    def set_forcing(self):
 107      M=self.fortran.main_module   
 108      M.forc_temp_surface[:]=self.t_rest*(self.t_star-M.temp[:,:,-1,M.tau-1])
 109      return
 110 
 111    def set_diagnostics(self):
 112      M=self.fortran.main_module   
 113      self.register_average(name='temp',long_name='Temperature',         units = 'deg C' , grid = 'TTT', var = M.temp)
 114      self.register_average(name='salt',long_name='Salinity',            units = 'g/kg' ,  grid = 'TTT', var = M.salt)
 115      self.register_average(name='u',   long_name='Zonal velocity',      units = 'm/s' ,   grid = 'UTT', var = M.u)
 116      self.register_average(name='v',   long_name='Meridional velocity', units = 'm/s' ,   grid = 'TUT', var = M.v)
 117      self.register_average(name='w',   long_name='Vertical velocity',   units = 'm/s' ,   grid = 'TTU', var = M.w)
 118      self.register_average(name='taux',long_name='wind stress',         units = 'm^2/s' , grid = 'UT',  var = M.surface_taux)
 119      self.register_average(name='tauy',long_name='wind stress',         units = 'm^2/s' , grid = 'TU',  var = M.surface_tauy)
 120      self.register_average(name='psi' ,long_name='Streamfunction',      units = 'm^3/s' , grid = 'UU',  var = M.psi)
 121      return
 122 
 123    def user_defined_signal(self):
 124        """ this routine must be called by all processors
 125        """
 126        M=self.fortran.main_module  
 127        a = zeros( (M.nx,M.ny), 'd', order = 'F')
 128        a[M.is_pe-1:M.ie_pe,0] = M.xt[2:-2]
 129        self.fortran.pe0_recv_2d(a)
 130        self.xt_gl = a[:,0].copy()
 131        
 132        a[0,M.js_pe-1:M.je_pe] = M.yt[2:-2]
 133        self.fortran.pe0_recv_2d(a)
 134        self.yt_gl = a[0,:].copy()
 135        
 136        self.psi_gl = zeros( (M.nx,M.ny), 'd', order = 'F')
 137        self.psi_gl[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] =  where( M.maskt[2:-2,2:-2,-1] >0,  M.psi[2:-2,2:-2,M.tau-1] , NaN) 
 138 
 139        self.fortran.pe0_recv_2d(self.psi_gl)
 140        
 141        self.temp_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 142        for k in range(M.nz):
 143          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) 
 144          self.fortran.pe0_recv_2d(a)
 145          self.temp_gl[:,:,k]=a.copy()
 146 
 147        self.u_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 148        for k in range(M.nz):
 149          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) 
 150          self.fortran.pe0_recv_2d(a)
 151          self.u_gl[:,:,k]=a.copy()  
 152 
 153        return
 154    
 155    def make_plot(self):
 156        from scipy import stats
 157        M=self.fortran.main_module 
 158        
 159        self.set_signal('user_defined') # following routine is called by all PEs
 160        self.user_defined_signal()
 161        self.figure.clf()
 162        
 163        ax=self.figure.add_subplot(221)
 164        co=ax.contourf(self.xt_gl/1e3,self.yt_gl/1e3,self.temp_gl[:,:,-1].transpose())
 165        self.figure.colorbar(co)
 166        ax.set_title('temperature')
 167        
 168        ax=self.figure.add_subplot(222)
 169        co=ax.contourf(self.xt_gl/1e3,self.yt_gl/1e3,self.psi_gl.transpose()*1e-6)
 170        self.figure.colorbar(co)
 171        ax.set_title('Streamfunction [Sv]')
 172        
 173        ax=self.figure.add_subplot(223)
 174        co=ax.contourf(self.yt_gl/1e3,M.zt,stats.nanmean(self.u_gl,axis=0).transpose())
 175        self.figure.colorbar(co)
 176        ax.set_title('zonal velocity [m/s]')
 177        
 178        ax=self.figure.add_subplot(224)
 179        co=ax.contourf(self.yt_gl/1e3,M.zt,stats.nanmean(self.temp_gl,axis=0).transpose())
 180        self.figure.colorbar(co)
 181        ax.set_title('temperature [deg C]')
 182        
 183        return
 184 
 185 if __name__ == "__main__": acc1().run(snapint= 86400*40.0, runlen = 365*86400.*200)