1 #
   2 # solving the linear stability problem
   3 # for U(z), V(z), and B(z)
   4 #
   5 # version for primitive equation is pe
   6 # version for quasi-geostrophic approx. is qg
   7 #
   8 # constant grid spacing is assumed for both
   9 
  10 def pe(U,V,B0,dz,kx,ky,betax,betay,fh,f0,hx,hy,f_filter=0.5):
  11   # solution for primitive equations
  12   import numpy as np
  13   import sys
  14   cs = 150 # speed of sound, artificially reduced
  15   N=U.shape[0]
  16   # derivatives of U,V and B
  17   UZ=np.zeros(N,'d');VZ=np.zeros(N,'d');BZ=np.zeros(N,'d')
  18   for n in range(1,N-1):
  19    UZ[n]=(U[n+1]-U[n-1])/(2*dz)
  20    VZ[n]=(V[n+1]-V[n-1])/(2*dz)
  21    BZ[n]=(B0[n+1]-B0[n-1])/(2*dz)
  22   UZ[0]=UZ[1];UZ[-1]=UZ[-2];
  23   VZ[0]=VZ[1];VZ[-1]=VZ[-2]; 
  24   BZ[0]=BZ[1];BZ[-1]=BZ[-2];
  25   # allocate some variables
  26   I=complex(0,1)
  27   A  = np.zeros((5,5,N),'Complex64')
  28   B  = np.zeros((5,5,N),'Complex64')
  29   C  = np.zeros((5,5,N),'Complex64')
  30   AA = np.zeros((5*N,5*N),'Complex64')
  31   om_max= np.zeros((kx.shape[0],ky.shape[0]),'Complex64');
  32   omax  = complex(0,0)
  33   AAmax = np.zeros((5*N,5*N),'Complex64')
  34   kmax  = 0; lmax  = 0
  35   # enter main loop
  36   for j in range(ky.shape[0]): # loop over meridional wavelength
  37     sys.stdout.write('\b'*21+'calculating j=%3i/%3i'%(j,ky.shape[0]) )
  38     sys.stdout.flush()
  39     for i in range(kx.shape[0]): # loop over zonal wavelength
  40 
  41       Uk   = kx[i]*U+ky[j]*V   # k \cdot \v U
  42       UZk  = kx[i]*UZ+ky[j]*VZ # k \cdot \v U'
  43       UZhk =-kx[i]*VZ+ky[j]*UZ # k \cdot \rvec{\v U'}
  44 
  45       kh2 = ky[j]**2 + kx[i]**2 + 1e-18
  46       fhk = ( -betay*kx[i]+betax*ky[j] )/kh2
  47       fk  = (  betax*kx[i]+betay*ky[j] )/kh2
  48 
  49       for n in range(N): # loop over depth
  50         np1 = min(N-1,n+1)
  51         nm  = max(0,n-1)
  52 
  53         B[0,2,n]=  (ky[j]*fh+UZhk[n])/(2*kh2)
  54         B[1,2,n]= -(kx[i]*fh+UZk[n])/(2*kh2)
  55         B[3,2,n]= I*BZ[n]/2
  56         B[4,2,n]= -cs**2*I/dz
  57 
  58         C[2,0,n]=  ky[j]*fh/2
  59         C[2,1,n]= -kx[i]*fh/2
  60         C[2,3,n]= -I/2
  61         C[2,4,n]=  I/dz
  62 
  63         A[0,:,n]= [ Uk[n]+fhk,I*f0,(ky[j]*fh+UZhk[n])/(2*kh2),0,0 ]
  64         A[1,:,n]= [-I*f0 ,Uk[n], -(kx[i]*fh+UZk[n])/(2*kh2),0,I]
  65         A[2,:,n]= [ ky[j]*fh/2 ,-kx[i]*fh/2 ,(Uk[n]+Uk[np1])/2,-I/2,-I/dz ]
  66         A[3,:,n]= [ -f0*UZk[n],-f0*UZhk[n],I*BZ[n]/2,Uk[n], 0  ]
  67         A[4,:,n]= [ 0, -I*cs**2*kh2 ,  I*cs**2/dz, 0   , 0 ]
  68 
  69       # upper boundary condition
  70       A[2,:,-1]=0; B[2,:,-1]=0; C[2,:,-1]=0;
  71       A[:,2,-1]=0
  72       C[:,2,-1]=0
  73       C[:,2,-2]=0
  74 
  75       # lower boundary condition
  76       #B[:,3,1]=0; 
  77       #A[:,0,0]=A[:,0,0]-I*(-kx[i]*hy+ky[j]*hx)*B[:,2,0]
  78       #A[:,1,0]=A[:,1,0]+I*( kx[i]*hx+ky[j]*hy)*B[:,2,0]
  79 
  80       A[:,0,0]=A[:,0,0]-2*I*(-kx[i]*hy+ky[j]*hx)*B[:,2,0]
  81       #A[:,1,0]=A[:,1,0]+2*I*( kx[i]*hx+ky[j]*hy)*B[:,2,0]
  82       A[:,2,0]=A[:,2,0]-B[:,2,0]
  83 
  84       # build large matrix
  85       n1 = 0; n2 = 5
  86       AA[ n1:n2,(n1+5):(n2+5) ] = C[:,:,0]
  87       AA[ n1:n2,n1:n2         ] = A[:,:,0]
  88       for n in range(1,N-1):
  89         n1 = 5*n; n2 = 5*(n+1)
  90         AA[ n1:n2, n1  :n2  ]     = A[:,:,n]
  91         AA[ n1:n2, (n1-5):(n2-5)] = B[:,:,n]
  92         AA[ n1:n2, (n1+5):(n2+5)] = C[:,:,n]
  93       n1 = 5*(N-1); n2 = 5*N
  94       AA[ n1:n2, n1    :n2     ] = A[:,:,-1];
  95       AA[ n1:n2, (n1-5):(n2-5) ] = B[:,:,-1];
  96 
  97       # eigenvalues of matrix
  98       om = np.linalg.eigvals(AA)
  99       kh=np.sqrt( kx[i]**2 + ky[j]**2 )
 100       # search minimal imaginary eigenvalue
 101       if kh>0: 
 102         om = np.extract( np.abs( np.real( om/kh )) <f_filter*cs, om )
 103         n=np.argmin( np.imag(om) )
 104         om_max[i,j]=om[n]
 105         # look for global minimum
 106         if np.imag(om[n]) < np.imag(omax):
 107            omax=om[n]; kmax=kx[i]; lmax=ky[j]; AAmax[:,:] = AA
 108   sys.stdout.write('\n')
 109 
 110   #eigenvectors for global minimum
 111   om, phi=np.linalg.eig(AAmax)
 112   n=np.argmin( np.imag(om) )
 113   om=om[n]
 114   phi = phi[:,n]
 115 
 116   # scale 
 117   A=2*np.pi*np.abs(np.imag(om))/(kmax**2+lmax**2)
 118   s=np.max( np.abs(np.real(phi[::5]))+np.abs(np.imag(phi[::5])) )
 119   phi=A*phi/s
 120 
 121   #complete solution
 122   str= phi[0::5]
 123   pot= phi[1::5]
 124   u= complex(0,1)*lmax*str-complex(0,1)*kmax*pot
 125   v=-complex(0,1)*kmax*str-complex(0,1)*lmax*pot
 126   w= phi[2::5]
 127   b= phi[3::5]
 128   p= phi[4::5]
 129 
 130   return om_max,om,kmax,lmax,u,v,w,b,p
 131  
 132 
 133 def qg(U,V,B0,dz,kx,ky,beta,f0,hx,hy):
 134   import numpy as np
 135   import sys
 136   # quasi geostrophic vertical eigen value problem after
 137   # Smith (2007) The Geography of Linear Baroclinic 
 138   # Instability in Earth's Oceans, J. Mar. Res. 
 139   N=U.shape[0]
 140   # derivatives of U and B
 141   UZ=np.zeros(N,'d');VZ=np.zeros(N,'d');BZ=np.zeros(N,'d')
 142   for n in range(1,N-1):
 143    UZ[n]=(U[n+1]-U[n-1])/(2*dz)
 144    VZ[n]=(V[n+1]-V[n-1])/(2*dz)
 145    BZ[n]=(B0[n+1]-B0[n-1])/(2*dz)
 146   UZ[0]=UZ[1];UZ[-1]=UZ[-2]; BY=-f0*UZ
 147   VZ[0]=VZ[1];VZ[-1]=VZ[-2]; BX= f0*VZ
 148   BZ[0]=BZ[1];BZ[-1]=BZ[-2]
 149 
 150   # vertical differencing operator
 151   G      = np.zeros((N,N),'d')
 152   G[0,0] = -f0**2*( 1.0/(B0[1]-B0[0]))/dz
 153   G[0,1] = -f0**2*(-1.0/(B0[1]-B0[0]))/dz
 154   for n in range(1,N-1): 
 155     G[n,n-1] = -f0**2*(-1./(B0[n]-B0[n-1]))/dz
 156     G[n,n]   = -f0**2*( 1./(B0[n]-B0[n-1])+1./(B0[n+1]-B0[n]))/dz
 157     G[n,n+1] = -f0**2*(-1./(B0[n+1]-B0[n]))/dz
 158   G[-1,-1] = -f0**2*( 1./(B0[-1]-B0[-2]))/dz
 159   G[-1,-2] = -f0**2*(-1./(B0[-1]-B0[-2]))/dz
 160 
 161   # background PV gradient
 162   Qy = beta-np.dot(G,U) 
 163   Qx = np.dot(G,V)
 164 
 165   # topography
 166   Qx[-1]+=f0*hx/dz
 167   Qy[-1]+=f0*hy/dz
 168 
 169 
 170   om_max=np.zeros((kx.shape[0],ky.shape[0]),'Complex64')
 171   omax = complex(0,0); kmax=0.; lmax=0.
 172   Amax = np.zeros((N,N),'Complex64')
 173 
 174   # loop over wavenumbers
 175   for j in  range(ky.shape[0]):
 176     sys.stdout.write('\b'*22+' calculating j=%3i/%3i'%(j,ky.shape[0]) )
 177     sys.stdout.flush()
 178     for i in  range(kx.shape[0]):
 179       # construct matrixes and solve eigenvalue problem
 180       B = G-(kx[i]**2+ky[j]**2)*np.eye(N)
 181       A = kx[i]*np.diag(Qy)-ky[j]*np.diag(Qx) \
 182             +np.dot( kx[i]*np.diag(U)+ky[j]*np.diag(V) ,B)
 183       A = np.dot(np.linalg.inv(B),A) # this way is faster than gen. problem
 184       om= np.linalg.eigvals(A)
 185       n = np.argmin( np.imag(om) )
 186       om_max[i,j]=om[n]
 187       # look for global minimum
 188       if np.imag(om[n]) < np.imag(omax):
 189          omax=om[n]; kmax=kx[i]; lmax=ky[j]; Amax = A.copy()
 190   sys.stdout.write('\n')
 191 
 192   #eigenvector for global minimum
 193   om, phi = np.linalg.eig(Amax)
 194   n=np.argmin( np.imag(om) )
 195   om  = om[n]
 196   phi = phi[:,n]
 197 
 198   # scale 
 199   A=2*np.pi*np.abs(np.imag(om))/(kmax**2+lmax**2)
 200   s=np.max( np.abs(np.real(phi))+np.abs(np.imag(phi)) )
 201   phi=A*phi/s
 202 
 203   #complete solution
 204   u= complex(0,1.)*lmax*phi
 205   v=-complex(0,1.)*kmax*phi
 206   phiz=np.zeros( (N,), 'Complex64' )
 207   phiz[1:-1]= (phi[2:]-phi[:-2])/(2*dz)
 208   phiz[0]   = (phi[1]-phi[0])/dz
 209   phiz[-1]  = (phi[-1]-phi[-2])/dz
 210   b=f0*phiz
 211   p=f0*phi
 212   # N^2 w = - ( b_t + U b_x + V b_y - psi_y B_x + psi_x B_y)
 213   # w = -(i omega b - i k U b - i l V b + i l psi B_x - i k psi B_y)/BZ 
 214   w = -complex(0,1)*(om*b-kmax*U*b-lmax*V*b+lmax*phi*BX-kmax*phi*BY)/BZ 
 215 
 216   return om_max,om,kmax,lmax,u,v,w,b,p