1 import sys; sys.path.append('../py_src')
   2 from pyOM_gui import pyOM_gui as pyOM
   3 #from pyOM_cdf import pyOM_cdf as pyOM
   4 from numpy import *
   5 
   6 class eady1(pyOM):
   7    """ Eady (1941) solution
   8    """
   9    def set_parameter(self):
  10      """set main parameter
  11      """
  12      M=self.fortran.main_module   
  13      
  14      (M.nx,M.ny,M.nz)    = (32,32,20)
  15      M.dt_tracer      = 1200.0 
  16      M.dt_mom         = 1200.0 
  17      
  18      M.congr_epsilon = 1e-12
  19      M.congr_max_iterations = 5000
  20      M.enable_streamfunction = 1
  21      
  22      M.enable_hydrostatic      = 1
  23      M.enable_cyclic_x         = 1
  24      
  25      M.enable_superbee_advection  = 1
  26      M.enable_explicit_vert_friction = 1
  27      M.enable_hor_friction = 1
  28      M.a_h = (20e3)**3*2e-11    
  29      M.kappam_0   = 1.e-4 
  30      M.kappah_0   = 1.e-4 
  31      
  32      M.enable_conserve_energy = 0
  33      M.coord_degree            = 0
  34      M.eq_of_state_type        = 1 
  35      M.enable_tempsalt_sources = 1
  36      return
  37 
  38 
  39    def set_grid(self):
  40        M=self.fortran.main_module   
  41        M.dxt[:]= 20e3
  42        M.dyt[:]= 20e3
  43        M.dzt[:]= 100.0
  44        return
  45 
  46    def set_topography(self):
  47        M=self.fortran.main_module   
  48        M.kbot[:]=0
  49        M.kbot[:,2:-2]=1
  50        return
  51 
  52    
  53    def set_coriolis(self):
  54      M=self.fortran.main_module   
  55      for j in range( M.yt.shape[0] ): M.coriolis_t[:,j] = 1e-4+0e-11*M.yt[j]
  56      return
  57 
  58    def set_forcing(self):
  59       M=self.fortran.main_module   
  60       
  61       # update density, etc of last time step
  62       M.temp[:,:,:,M.tau-1] = M.temp[:,:,:,M.tau-1] + self.t0
  63       self.fortran.calc_eq_of_state(M.tau)
  64       M.temp[:,:,:,M.tau-1] = M.temp[:,:,:,M.tau-1] - self.t0
  65       
  66       # advection of background temperature
  67       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) 
  68       M.temp_source[:] = (1.5+M.ab_eps)*self.dt0[...,M.tau-1] - ( 0.5+M.ab_eps)*self.dt0[...,M.taum1-1]
  69       return
  70      
  71 
  72    def set_initial_conditions(self):
  73      """ setup all initial conditions
  74      """
  75      M=self.fortran.main_module   
  76      U_0 = 0.5
  77      N_0 = 0.004
  78      f=M.coriolis_t[M.ny/2]
  79      h = (M.nz-2)*M.dzt[0]
  80      kx=1.6*f/(N_0*h)
  81      ky=pi/((M.ny-2)*M.dxt[0])
  82      d=f/N_0/(kx**2+ky**2)**0.5
  83 
  84      fxa=(exp(h/d)+exp(-h/d))/(exp(h/d)-exp(-h/d))
  85      c1= (1+0.25*(h/d)**2-h/d*fxa )*complex(1,0)
  86      c1=(sqrt(c1)*d/h+0.5)*U_0
  87      A=(U_0-c1)/U_0*h/d
  88      
  89      alpha = self.fortran.linear_eq_of_state.linear_eq_of_state_drhodt()
  90      grav = 9.81; rho0 = 1024.0
  91         
  92      # zonal velocity 
  93      for k in range(M.nz):
  94        M.u[:,:,k,M.tau-1]= (U_0/2+U_0*M.zt[k]/(M.nz*M.dzt[0]))*M.masku[:,:,k] 
  95      M.u[...,M.taum1-1] = M.u[...,M.tau-1]
  96      
  97      self.t0  = zeros( (M.i_blk+2*M.onx,M.j_blk+2*M.onx,M.nz) , 'd', order='F')
  98      self.dt0 = zeros( (M.i_blk+2*M.onx,M.j_blk+2*M.onx,M.nz,3) , 'd', order='F')
  99 
 100      # rho = alpha T ,  N^2 = b_z = - g/rho0 rho_z = - g/rho0 alpha T_z,  T = - N^2 z rho0/(g alpha)
 101      for k in range(M.nz):
 102         self.t0[:,:,k]=-N_0**2*M.zt[k]/grav/alpha*rho0*M.maskt[:,:,k]
 103      
 104      # fu = -p_y, p_z = -g rho,  f u_z = -g rho_y,  rho_y = - f u_z/g = alpha T_y
 105      for k in range(M.nz):
 106       for j in range(1,M.ny+1):
 107         jj = self.jf2py(j)
 108         uz = U_0/M.ht[:,jj]
 109         self.t0[:,jj+1,k]=(self.t0[:,jj,k]+M.dyt[jj]*uz*f/grav/alpha*rho0)*M.maskt[:,jj,k]
 110         
 111         
 112      # perturbation buoyancy
 113      for k in range(M.nz):
 114        for j in range(1,M.ny+1):
 115         jj = self.jf2py(j)
 116         phiz=A/d*sinh(M.zt[k]/d)+cosh(M.zt[k]/d)/d
 117         M.temp[:,jj,k,M.tau-1] =0.1*sin(kx*M.xt)*sin(ky*M.yt[jj])*abs(phiz)*M.maskt[:,jj,k]*rho0/grav/alpha
 118         #M.temp[:,jj,k,M.tau-1] =1e-5*random.randn(M.nx+2*M.onx)*M.maskt[:,jj,k]
 119      M.temp[...,M.taum1-1] = M.temp[...,M.tau-1]
 120      return
 121 
 122 
 123    def make_plot(self):
 124      """ make a plot using methods of self.figure
 125      """
 126      if hasattr(self,'figure'):
 127        M=self.fortran.main_module         # fortran module with model variables
 128        k=M.nz*3/4
 129        i=self.if2py(M.nx/2)
 130        j=self.jf2py(M.ny/2)
 131        x=M.xt[2:-2]/1e3
 132        y=M.yt[2:-2]/1e3
 133        z=M.zt[:]
 134 
 135        self.figure.clf()
 136        ax=self.figure.add_subplot(221)
 137        a=M.temp[i,2:-2,:,M.tau-1] +self.t0[i,2:-2,:]
 138        co=ax.contourf(y,z,a.transpose())
 139        ax.set_ylabel('y [km]')
 140        self.figure.colorbar(co)
 141        a=M.u[i,2:-2,:,M.tau-1]
 142        ax.contour(y,z,a.transpose(),10,colors='black')
 143        ax.set_title('total Temperature');  
 144        ax.axis('tight') 
 145          
 146               
 147        
 148        ax=self.figure.add_subplot(222)
 149        a= M.temp[2:-2,2:-2,k,M.tau-1]  
 150        co=ax.contourf(x,y,a.transpose())
 151        self.figure.colorbar(co)
 152        a=M.temp[2:-2,2:-2,k,M.tau-1] +self.t0[2:-2,2:-2,k]
 153        ax.contour(x,y,a.transpose(),10,colors='black')
 154        ax.axis('tight')
 155 
 156                        
 157        
 158        ax=self.figure.add_subplot(223)
 159        a=M.v[2:-2,j,:,M.tau-1] 
 160        co=ax.contourf(x,z,a.transpose())
 161        ax.set_title('meridional velocity'); 
 162        ax.set_xlabel('x [km]'); 
 163        self.figure.colorbar(co)
 164        ax.axis('tight')
 165        
 166        
 167        ax=self.figure.add_subplot(224)
 168        a=M.temp[2:-2,j,:,M.tau-1] 
 169        co=ax.contourf(x,z,a.transpose())
 170        ax.set_title('temperature perturbation'); 
 171        ax.set_xlabel('x [km]'); 
 172        self.figure.colorbar(co)
 173        ax.axis('tight')
 174        
 175        
 176      return
 177 
 178 if __name__ == "__main__": 
 179     model = eady1()
 180     M=model.fortran.main_module
 181     model.run( snapint = 86400.0)