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 
   6 from numpy import *
   7 from scipy.io.netcdf import netcdf_file as NF
   8 
   9 class flame_E7(pyOM):
  10    """ the 4/3 FLAME Atlantic model
  11    """
  12    def set_parameter(self):
  13      """set main parameter
  14      """
  15      M=self.fortran.main_module   
  16 
  17      (M.nx,M.ny,M.nz)   = (87,89,45)
  18      M.dt_mom    = 3600.0
  19      M.dt_tracer = 3600.0
  20 
  21      M.runlen = 86400*365.0
  22      
  23      M.coord_degree     = 1
  24 
  25      M.congr_epsilon = 1e-6
  26      M.congr_max_iterations = 20000
  27      M.enable_streamfunction = 1
  28 
  29      I=self.fortran.isoneutral_module   
  30      I.enable_neutral_diffusion = 1
  31      I.enable_skew_diffusion = 1
  32      I.k_iso_0 = 1000.0
  33      I.k_iso_steep = 200.0
  34      I.iso_dslope = 1./1000.0
  35      I.iso_slopec = 4./1000.0
  36 
  37      M.enable_hor_friction = 1
  38      M.a_h = 5e4 
  39      M.enable_hor_friction_cos_scaling = 1
  40      M.hor_friction_cosPower = 3
  41      M.enable_tempsalt_sources = 1
  42 
  43      M.enable_implicit_vert_friction = 1
  44      T=self.fortran.tke_module   
  45      T.enable_tke = 1
  46      T.c_k = 0.1
  47      T.c_eps = 0.7
  48      T.alpha_tke = 30.0
  49      T.mxl_min = 1e-8
  50      T.tke_mxl_choice = 2
  51 
  52      M.k_gm_0 = 1000.0
  53      #E=self.fortran.eke_module   
  54      #E.enable_eke = 1
  55 
  56      I=self.fortran.idemix_module   
  57      I.enable_idemix = 1
  58      I.enable_idemix_hor_diffusion = 1
  59 
  60      M.eq_of_state_type = 5
  61      return
  62 
  63 
  64    def set_grid(self):
  65        M=self.fortran.main_module   
  66        fid = NF('forcing.cdf','r')
  67        (i1,i2) = ( self.if2py(M.is_pe), self.if2py(M.ie_pe) )
  68        M.dxt[i1:i2+1] = fid.variables['dxtdeg'][M.is_pe-1:M.ie_pe]
  69        (j1,j2) = ( self.jf2py(M.js_pe), self.jf2py(M.je_pe) )
  70        M.dyt[j1:j2+1] = fid.variables['dytdeg'][M.js_pe-1:M.je_pe] 
  71        M.dzt[:] = fid.variables['dzt'][::-1] /100.0
  72        M.x_origin = fid.variables['xu'][0]
  73        M.y_origin = fid.variables['yu'][0]
  74        return
  75 
  76 
  77    def set_coriolis(self):
  78      M=self.fortran.main_module
  79      for j in range( M.yt.shape[0] ): M.coriolis_t[:,j] = 2*M.omega*sin(M.yt[j]/180.*pi)
  80      return
  81 
  82    def set_topography(self):
  83        """ setup topography
  84        """
  85        M=self.fortran.main_module   
  86        fid = NF('forcing.cdf','r')
  87        kmt = fid.variables['kmt'][:]
  88        M.kbot[:]=1
  89        for j in range(M.js_pe,M.je_pe+1):
  90          for i in range(M.is_pe,M.ie_pe+1):
  91            (ii,jj) = ( self.if2py(i), self.jf2py(j) )
  92            M.kbot[ii,jj]=min(M.nz,M.nz-kmt[j-1,i-1]+1)
  93            if  kmt[j-1,i-1] == 0: M.kbot[ii,jj] =0
  94        return
  95    
  96    def set_initial_conditions(self):
  97        """ setup initial conditions
  98        """
  99        M=self.fortran.main_module   
 100        fid = NF('forcing.cdf','r')
 101        
 102        # initial conditions
 103        t = fid.variables['temp_ic'][::-1,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe]
 104        s = fid.variables['salt_ic'][::-1,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe]
 105        for k in range(M.nz):
 106           M.temp[ self.if2py(M.is_pe):self.if2py(M.ie_pe)+1, self.jf2py(M.js_pe):self.jf2py(M.je_pe)+1, k,M.tau-1] = t[k,:,:].transpose()
 107           M.salt[ self.if2py(M.is_pe):self.if2py(M.ie_pe)+1, self.jf2py(M.js_pe):self.jf2py(M.je_pe)+1, k,M.tau-1] = s[k,:,:].transpose()*1000+35.0
 108        for d in (M.temp,M.salt):
 109            d[...,M.tau-1] = d[...,M.tau-1]*M.maskt
 110            d[...,M.taum1-1] = d[...,M.tau-1]
 111 
 112        # wind stress
 113        self.tx = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 114        self.ty = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 115        for k in range(12):
 116            self.tx[2:-2,2:-2,k] = fid.variables['taux'][k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()/10. /M.rho_0
 117            self.ty[2:-2,2:-2,k] = fid.variables['tauy'][k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()/10. /M.rho_0
 118        # zero out masks and do boundary exchange
 119        for d in (self.tx, self.ty) :
 120          d[ d <= -1e10 ] = 0.0
 121          for k in range(12):
 122            self.fortran.border_exchg_xy(M.is_pe-M.onx,M.ie_pe+M.onx,M.js_pe-M.onx,M.je_pe+M.onx,d[:,:,k]) 
 123            self.fortran.setcyclic_xy   (M.is_pe-M.onx,M.ie_pe+M.onx,M.js_pe-M.onx,M.je_pe+M.onx,d[:,:,k])
 124        
 125        # interpolate from B grid location to C grid
 126        for k in range(12):
 127            self.tx[:,1:,k] = (self.tx[:,:-1,k] + self.tx[:,1:,k])/(M.maskz[:,:-1,-1]+M.maskz[:,1:,-1]+1e-20)
 128            self.ty[1:,:,k] = (self.ty[:-1,:,k] + self.ty[1:,:,k])/(M.maskz[:-1,:,-1]+M.maskz[1:,:,-1]+1e-20)
 129            self.tx[:,:,k] = self.tx[:,:,k]*M.masku[:,:,-1]
 130            self.ty[:,:,k] = self.ty[:,:,k]*M.maskv[:,:,-1]
 131 
 132        # heat flux and salinity restoring
 133        self.sst_clim = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 134        self.sss_clim = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 135        self.sst_rest = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 136        self.sss_rest = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, 12 ), 'd' , order = 'F')
 137        for k in range(12):
 138            self.sst_clim[2:-2,2:-2,k] = fid.variables['sst_clim'][k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()
 139            self.sss_clim[2:-2,2:-2,k] = fid.variables['sss_clim'][k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()*1000+35
 140            self.sst_rest[2:-2,2:-2,k] = fid.variables['sst_rest'][k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()*41868.
 141            self.sss_rest[2:-2,2:-2,k] = fid.variables['sss_rest'][k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()/100.0
 142        for d in (self.sst_clim, self.sss_clim, self.sst_rest, self.sss_rest): d[ d  <= -1e10 ] = 0.0
 143 
 144        # sponge layers
 145        self.t_star    = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, M.nz, 12 ), 'd' , order = 'F')
 146        self.s_star    = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, M.nz, 12 ), 'd' , order = 'F')
 147        self.rest_tscl = zeros( ( M.i_blk+2*M.onx, M.j_blk+2*M.onx, M.nz ), 'd' , order = 'F')
 148        fid = NF('restoring_zone.cdf','r')
 149        (i1,j1) = ( self.if2py(M.is_pe  ), self.jf2py(M.js_pe) )
 150        (i2,j2) = ( self.if2py(M.ie_pe+1), self.jf2py(M.je_pe+1) )
 151        for k in range(M.nz):
 152          self.rest_tscl[i1:i2,j1:j2,k] = fid.variables['tscl'][0,k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()
 153        for n in range(12):
 154          for k in range(M.nz):
 155            self.t_star[i1:i2,j1:j2,k,n] = fid.variables['t_star'][n,k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()
 156            self.s_star[i1:i2,j1:j2,k,n] = fid.variables['s_star'][n,k,M.js_pe-1:M.je_pe,M.is_pe-1:M.ie_pe].transpose()
 157 
 158        I=self.fortran.idemix_module   
 159        if I.enable_idemix:
 160          f = reshape( fromfile('tidal_energy.bin', dtype='>f'), (M.nx,M.ny), order='F' )/M.rho_0
 161          I.forc_iw_bottom[2:-2,2:-2]  = f[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe]
 162          f = reshape( fromfile('wind_energy.bin', dtype='>f'), (M.nx,M.ny), order='F' )/M.rho_0*0.2
 163          I.forc_iw_surface[2:-2,2:-2] = f[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe]
 164        return
 165 
 166    def get_periodic_interval(self,currentTime,cycleLength,recSpacing,nbRec):
 167        """  interpolation routine taken from mitgcm
 168        """
 169        locTime = currentTime - recSpacing*0.5 + cycleLength*( 2 - round(currentTime/cycleLength) )
 170        tmpTime = locTime % cycleLength 
 171        tRec1 = 1 + int( tmpTime/recSpacing )
 172        tRec2 = 1 + tRec1 % int(nbRec) 
 173        wght2 = ( tmpTime - recSpacing*(tRec1 - 1) )/recSpacing
 174        wght1 = 1.0 - wght2
 175        return (wght1,wght2,tRec1,tRec2)
 176    
 177 
 178    def set_forcing(self):
 179      M=self.fortran.main_module   
 180      fxa = 365*86400.0
 181      (f1,f2,n1,n2) = self.get_periodic_interval((M.itt-1)*M.dt_tracer,fxa,fxa/12.,12)
 182      
 183      M.surface_taux[:]=(f1*self.tx[:,:,n1-1] + f2*self.tx[:,:,n2-1])
 184      M.surface_tauy[:]=(f1*self.ty[:,:,n1-1] + f2*self.ty[:,:,n2-1])
 185 
 186      T=self.fortran.tke_module   
 187      if T.enable_tke:
 188        T.forc_tke_surface[1:-1,1:-1] = sqrt( (0.5*(M.surface_taux[1:-1,1:-1]+M.surface_taux[:-2,1:-1]) )**2  \
 189                                             +(0.5*(M.surface_tauy[1:-1,1:-1]+M.surface_tauy[1:-1,:-2]) )**2 )**(3./2.) 
 190      cp_0 = 3991.86795711963
 191      M.forc_temp_surface[:]=(f1*self.sst_rest[:,:,n1-1]+f2*self.sst_rest[:,:,n2-1])* \
 192                             (f1*self.sst_clim[:,:,n1-1]+f2*self.sst_clim[:,:,n2-1]-M.temp[:,:,-1,M.tau-1])*M.maskt[:,:,-1]/cp_0/M.rho_0
 193      M.forc_salt_surface[:]=(f1*self.sss_rest[:,:,n1-1]+f2*self.sss_rest[:,:,n2-1])* \
 194                             (f1*self.sss_clim[:,:,n1-1]+f2*self.sss_clim[:,:,n2-1]-M.salt[:,:,-1,M.tau-1])*M.maskt[:,:,-1]
 195      # apply simple ice mask                     
 196      n=nonzero( logical_and( M.temp[:,:,-1,M.tau-1]*M.maskt[:,:,-1] <= -1.8 , M.forc_temp_surface <= 0.0 ) )                            
 197      M.forc_temp_surface[n]=0.0
 198      M.forc_salt_surface[n]=0.0
 199 
 200      if M.enable_tempsalt_sources:
 201         M.temp_source[:] = M.maskt*self.rest_tscl*(f1*self.t_star[:,:,:,n1-1]+f2*self.t_star[:,:,:,n2-1] - M.temp[:,:,:,M.tau-1] )
 202         M.salt_source[:] = M.maskt*self.rest_tscl*(f1*self.s_star[:,:,:,n1-1]+f2*self.s_star[:,:,:,n2-1] - M.salt[:,:,:,M.tau-1] )
 203      return
 204 
 205    def set_diagnostics(self):
 206      M=self.fortran.main_module   
 207      self.register_average(name='temp',long_name='Temperature',         units = 'deg C' , grid = 'TTT', var = M.temp)
 208      self.register_average(name='salt',long_name='Salinity',            units = 'g/kg' ,  grid = 'TTT', var = M.salt)
 209      self.register_average(name='u',   long_name='Zonal velocity',      units = 'm/s' ,   grid = 'UTT', var = M.u)
 210      self.register_average(name='v',   long_name='Meridional velocity', units = 'm/s' ,   grid = 'TUT', var = M.v)
 211      self.register_average(name='w',   long_name='Vertical velocity',   units = 'm/s' ,   grid = 'TTU', var = M.w)
 212      self.register_average(name='taux',long_name='wind stress',         units = 'm^2/s' , grid = 'UT',  var = M.surface_taux)
 213      self.register_average(name='tauy',long_name='wind stress',         units = 'm^2/s' , grid = 'TU',  var = M.surface_tauy)
 214      self.register_average(name='psi' ,long_name='Streamfunction',      units = 'm^3/s' , grid = 'UU',  var = M.psi)
 215      return
 216 
 217 
 218    def user_defined_signal(self):
 219        """ this routine must be called by all processors
 220        """
 221        M=self.fortran.main_module  
 222        a = zeros( (M.nx,M.ny), 'd', order = 'F')
 223        a[M.is_pe-1:M.ie_pe,0] = M.xt[2:-2]
 224        self.fortran.pe0_recv_2d(a)
 225        self.xt_gl = a[:,0].copy()
 226        
 227        a[0,M.js_pe-1:M.je_pe] = M.yt[2:-2]
 228        self.fortran.pe0_recv_2d(a)
 229        self.yt_gl = a[0,:].copy()
 230        
 231        self.psi_gl = zeros( (M.nx,M.ny), 'd', order = 'F')
 232        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) 
 233        self.fortran.pe0_recv_2d(self.psi_gl)
 234        
 235        self.temp_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 236        for k in range(M.nz):
 237          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) 
 238          self.fortran.pe0_recv_2d(a)
 239          self.temp_gl[:,:,k]=a.copy()
 240 
 241        self.salt_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 242        for k in range(M.nz):
 243          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.salt[2:-2,2:-2,k,M.tau-1] , NaN) 
 244          self.fortran.pe0_recv_2d(a)
 245          self.salt_gl[:,:,k]=a.copy()
 246 
 247        #self.kappa_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
 248        #for k in range(M.nz):
 249        #  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.kappah[2:-2,2:-2,k] , NaN) 
 250        #  self.fortran.pe0_recv_2d(a)
 251        #  self.kappa_gl[:,:,k]=a.copy()
 252 
 253        return
 254 
 255 
 256    def make_plot(self):
 257        M=self.fortran.main_module         # fortran module with model variables
 258        self.figure.clf()
 259        
 260        self.set_signal('user_defined') # following routine is called by all PEs
 261        self.user_defined_signal()
 262        
 263        ax=self.figure.add_subplot(221)
 264        co=ax.contourf(self.yt_gl,M.zt,self.temp_gl[M.nx/2-1,:,:].transpose())
 265        self.figure.colorbar(co)
 266        ax.set_title('Temperature [deg C]')
 267        ax.set_ylabel('z [m]')
 268        ax.axis('tight')
 269 
 270        ax=self.figure.add_subplot(223)
 271        co=ax.contourf(self.yt_gl,M.zt,self.salt_gl[M.nx/2-1,:,:].transpose())
 272        self.figure.colorbar(co)
 273        ax.set_title('Salinity [g/kg]')
 274        ax.set_ylabel('z [m]')
 275        ax.axis('tight')
 276 
 277        #ax=self.figure.add_subplot(223)
 278        #try:
 279        # co=ax.contourf(self.yt_gl,M.zw,log10(self.kappa_gl[M.nx/2-1,:,:].transpose()) )
 280        #except:
 281        # pass
 282        #self.figure.colorbar(co)
 283        #ax.set_title('Diffusivity')
 284        #ax.set_xlabel('Latitude [deg N]')
 285        #ax.set_ylabel('z [m]')
 286        #ax.axis('tight')
 287 
 288        ax=self.figure.add_subplot(222)
 289        co=ax.contourf(self.xt_gl,self.yt_gl,self.psi_gl.transpose()*1e-6)
 290        self.figure.colorbar(co)
 291        ax.set_title('Streamfunction [Sv]')
 292        #ax.set_xlabel('Longitude [deg E]')
 293        ax.set_ylabel('Latitude [deg N]')
 294        ax.axis('tight')
 295        
 296        ax=self.figure.add_subplot(224)
 297        co=ax.contourf(self.xt_gl,self.yt_gl,self.temp_gl[:,:,-1].transpose())
 298        self.figure.colorbar(co)
 299        ax.set_title('Sea Surface Temperature [deg C]')
 300        #ax.set_xlabel('Longitude [deg E]')
 301        ax.set_ylabel('Latitude [deg N]')
 302        ax.axis('tight')
 303        
 304        return
 305 
 306 if __name__ == "__main__": flame_E7().run(snapint= 86400.0)